For odd \(k\), define
$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$
where \(\{x\}\) denotes the fractional part of \(x\). The problem then asks for
$$S(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N} f_k(n),$$
and specifically for \(\lfloor S(33557799775533)\rfloor \bmod 977676779\). A direct double loop over all odd \(k\) and all \(n\) is hopeless at this scale, so the solution turns the fractional-part sum into a much more structured divisibility problem.
The C++, Python, and Java implementations all follow the same chain of reductions: first remove the fractional parts, then express the remaining term through prime exponents, and finally reorganize the whole computation by divisor profiles.
Because \(k\) is odd,
$$ (n-i)^k\equiv -\,i^k \pmod n. $$
So for \(1\le i\le n\), the two terms \(\left\{\frac{i^k}{n}\right\}\) and \(\left\{\frac{(n-i)^k}{n}\right\}\) add up to \(1\), except when \(n\mid i^k\), in which case both fractional parts are \(0\). Define
$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$
Then the pairing identity becomes
$$2f_k(n)=n-z_k(n).$$
This is the first crucial simplification: the problem is no longer about fractional parts themselves, only about how many residues produce exact divisibility.
Write
$$n=\prod_{p} p^{e_p}.$$
The condition \(n\mid i^k\) is equivalent to
$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad\text{for every prime }p\mid n,$$
where \(v_p\) is the \(p\)-adic valuation. Therefore \(i\) must be a multiple of
$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}.$$
Among the integers \(1,2,\dots,n\), exactly \(n/m_k(n)\) are such multiples. Hence
$$z_k(n)=\frac{n}{m_k(n)}$$
and therefore
$$2f_k(n)=n-\frac{n}{m_k(n)}.$$
Define
$$T(N)=2S(N).$$
Using the previous identity,
$$T(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N),$$
where
$$A(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$
for the actual odd input \(N\), and
$$H(N)=\sum_{\substack{k=1\\k\text{ odd}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}.$$
So the entire problem has been reduced to evaluating \(H(N)\) efficiently.
The implementation does not sum \(\frac{n}{m_k(n)}\) directly. Instead, for odd \(k\ge 3\), it introduces
$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}.$$
Then one has the identity
$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$
where \(\varphi\) is Euler's totient function.
Why is this true? It is enough to check one prime power. If \(n=p^E\), then the admissible exponents \(a\) in \(d=p^a\) are exactly those satisfying
$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E.$$
The largest such \(a\) is \(E-\left\lceil E/k\right\rceil\), so
$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$
Since both sides are multiplicative in \(n\), the formula follows for general \(n\). For \(k=1\), the situation is simpler: \(m_1(n)=n\), hence \(\frac{n}{m_1(n)}=1\).
Substituting the divisor identity gives
$$H(N)=N+\sum_{\substack{k=3\\k\text{ odd}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$
Now fix one divisor profile
$$d=\prod_{j=1}^{r} p_j^{a_j}.$$
Only profiles with
$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$
can contribute, because \(\rho_k(d)\ge \operatorname{rad}(d)\) for every odd \(k\ge 3\). This immediately implies that every prime in \(d\) is at most \(\sqrt N\), which is why the implementations only sieve primes up to \(\sqrt N\).
For a fixed \(d\), the values \(\left\lceil a_j/(k-1)\right\rceil\) do not change at every odd \(k\); they change only when \(k-1\) crosses one of finitely many exponent thresholds. Once \(k\) is larger than the biggest exponent \(a_j\), every ceiling becomes \(1\), so
$$\rho_k(d)=\operatorname{rad}(d).$$
That means the whole tail of large odd \(k\) can be summed in one block, while the small odd \(k\) values are handled individually. This is the compression that makes the problem tractable.
Take
$$72=2^3\cdot 3^2.$$
For \(k=3\),
$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6,$$
so
$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$
The divisor-profile formula gives the same result. Here
$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$
For the prime \(2\), the admissible exponents \(a\) satisfy \(a+\lceil a/2\rceil\le 3\), so \(a\in\{0,1,2\}\). For the prime \(3\), the condition is \(a+\lceil a/2\rceil\le 2\), so \(a\in\{0,1\}\). Therefore
$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$
which matches \(\frac{72}{6}\) exactly.
The implementations first generate all primes up to \(\sqrt N\) with a linear sieve. That is enough because every contributing divisor profile \(d\) must satisfy \(d\,\operatorname{rad}(d)\le N\), so no prime factor larger than \(\sqrt N\) can appear in a nontrivial profile.
They then enumerate divisor profiles \(d=\prod p^a\) recursively in increasing prime order. During that recursion they maintain four multiplicative pieces of state: the value of \(d\), the radical \(\operatorname{rad}(d)\), the totient \(\varphi(d)\), and the largest exponent currently present in the profile. The recursion stops immediately once \(d\,\operatorname{rad}(d)>N\), which is the main pruning rule.
For each profile, the code evaluates
$$\sum_{\substack{k=3\\k\text{ odd}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$
without iterating over every odd \(k\). Large odd exponents share the same denominator \(d\,\operatorname{rad}(d)\), so they are aggregated in one arithmetic block; only the finitely many smaller breakpoints are processed one by one. After multiplying by \(\varphi(d)\) and summing all profiles, the implementation obtains \(H(N)\), computes \(T(N)=A(N)-H(N)\) modulo \(2M\) with \(M=977676779\), and finally converts that residue into
$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M.$$
The C++ and Java implementations also split the top-level prime branches across threads, while the Python implementation uses the same arithmetic in a single-threaded recursive traversal.
The prime sieve up to \(\sqrt N\) costs \(O(\sqrt N)\) time and \(O(\sqrt N)\) memory. The harder part is the recursive enumeration of divisor profiles. If
$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$
then the recursion visits exactly the admissible profiles in \(\mathcal{D}(N)\), and for each profile it evaluates only the finitely many odd-\(k\) breakpoint intervals where some ceiling \(\left\lceil a/(k-1)\right\rceil\) changes. So the running time is roughly
$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right),$$
where \(B(d)\) is the number of distinct odd-\(k\) breakpoints induced by the exponents of \(d\). There is no simple closed elementary form for this quantity, but it is dramatically smaller than the original \(O(N^2)\) scan over all \((k,n)\). Memory usage stays at \(O(\sqrt N)\) for the prime table plus \(O(\log N)\) recursion depth.
Für ungerades \(k\) sei
$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$
wobei \(\{x\}\) den Nachkommateil von \(x\) bezeichnet. Gesucht ist dann
$$S(N)=\sum_{\substack{k=1\\k\text{ ungerade}}}^{N}\sum_{n=1}^{N} f_k(n),$$
genauer \(\lfloor S(33557799775533)\rfloor \bmod 977676779\). Ein direkter Doppelloop über alle ungeraden \(k\) und alle \(n\) ist völlig unbrauchbar, daher formt die Lösung die Summe in ein strukturiertes Teilbarkeitsproblem um.
Die C++-, Python- und Java-Implementierungen benutzen dieselbe Ideenfolge: erst die Bruchteile beseitigen, dann die verbleibenden Terme über Primexponenten beschreiben und schließlich alles nach Teilerprofilen gruppieren.
Da \(k\) ungerade ist, gilt
$$ (n-i)^k\equiv -\,i^k \pmod n. $$
Also addieren sich \(\left\{\frac{i^k}{n}\right\}\) und \(\left\{\frac{(n-i)^k}{n}\right\}\) zu \(1\), außer wenn \(n\mid i^k\); dann sind beide Bruchteile gleich \(0\). Definiere
$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$
Damit wird die Paarungsidentität zu
$$2f_k(n)=n-z_k(n).$$
Das ist die erste große Vereinfachung: Statt Bruchteilen müssen wir nur noch zählen, wie viele Restklassen eine exakte Teilbarkeit liefern.
Schreibe
$$n=\prod_{p} p^{e_p}.$$
Die Bedingung \(n\mid i^k\) ist äquivalent zu
$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad\text{für alle Primzahlen }p\mid n.$$
Also muss \(i\) ein Vielfaches von
$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}$$
sein. Unter \(1,2,\dots,n\) gibt es genau \(n/m_k(n)\) solche Vielfache. Folglich
$$z_k(n)=\frac{n}{m_k(n)}$$
und damit
$$2f_k(n)=n-\frac{n}{m_k(n)}.$$
Setze
$$T(N)=2S(N).$$
Mit der vorigen Identität folgt
$$T(N)=\sum_{\substack{k=1\\k\text{ ungerade}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N),$$
wobei
$$A(N)=\sum_{\substack{k=1\\k\text{ ungerade}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$
für den tatsächlichen ungeraden Eingabewert \(N\) gilt, und
$$H(N)=\sum_{\substack{k=1\\k\text{ ungerade}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}.$$
Damit bleibt als eigentliche Aufgabe nur noch die effiziente Berechnung von \(H(N)\).
Die Implementierung summiert \(\frac{n}{m_k(n)}\) nicht direkt. Für ungerades \(k\ge 3\) wird stattdessen
$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}$$
eingeführt. Dann gilt die Identität
$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$
wobei \(\varphi\) die Eulersche Phi-Funktion ist.
Warum stimmt das? Es reicht, eine Primzahlpotenz zu betrachten. Für \(n=p^E\) sind genau die Exponenten \(a\) in \(d=p^a\) zulässig, die
$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E$$
erfüllen. Der größte solche Exponent ist \(E-\left\lceil E/k\right\rceil\). Also
$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$
Da beide Seiten multiplikativ in \(n\) sind, folgt die Formel allgemein. Für \(k=1\) ist der Fall noch einfacher: \(m_1(n)=n\), also \(\frac{n}{m_1(n)}=1\).
Einsetzen der Teileridentität liefert
$$H(N)=N+\sum_{\substack{k=3\\k\text{ ungerade}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$
Fixiere nun ein Teilerprofil
$$d=\prod_{j=1}^{r} p_j^{a_j}.$$
Es können nur Profile mit
$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$
beitragen, weil stets \(\rho_k(d)\ge \operatorname{rad}(d)\) für ungerades \(k\ge 3\) gilt. Daraus folgt sofort: Jede Primzahl in \(d\) ist höchstens \(\sqrt N\). Genau deshalb reicht im Programm ein Sieb nur bis \(\sqrt N\).
Für festes \(d\) ändern sich die Werte \(\left\lceil a_j/(k-1)\right\rceil\) nicht bei jedem ungeraden \(k\), sondern nur an endlich vielen Schwellen. Sobald \(k\) größer ist als der größte Exponent \(a_j\), werden alle diese Ceiling-Ausdrücke gleich \(1\), also
$$\rho_k(d)=\operatorname{rad}(d).$$
Damit lässt sich der ganze Schwanz großer ungerader \(k\) in einem Block aufsummieren, während nur die wenigen kleinen ungeraden \(k\) einzeln behandelt werden müssen. Genau diese Kompression macht die Lösung schnell.
Nimm
$$72=2^3\cdot 3^2.$$
Für \(k=3\) ist
$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6,$$
also
$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$
Die Teilerprofil-Formel liefert dasselbe Resultat. Hier gilt
$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$
Für die Primzahl \(2\) ist \(a+\lceil a/2\rceil\le 3\), also \(a\in\{0,1,2\}\). Für die Primzahl \(3\) ist \(a+\lceil a/2\rceil\le 2\), also \(a\in\{0,1\}\). Deshalb
$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$
genau gleich \(\frac{72}{6}\).
Die Implementierungen erzeugen zunächst mit einem linearen Sieb alle Primzahlen bis \(\sqrt N\). Das genügt, weil jedes beitragende Teilerprofil \(d\) die Schranke \(d\,\operatorname{rad}(d)\le N\) erfüllen muss und somit kein Primfaktor größer als \(\sqrt N\) in einem nichttrivialen Profil auftreten kann.
Danach werden die Teilerprofile \(d=\prod p^a\) rekursiv in aufsteigender Primzahlreihenfolge erzeugt. Während dieser Rekursion werden vier multiplikative Zustandsgrößen mitgeführt: der Wert von \(d\), der Radikalteil \(\operatorname{rad}(d)\), der Totient \(\varphi(d)\) und der größte bisher auftretende Exponent. Sobald \(d\,\operatorname{rad}(d)>N\) gilt, wird der Ast sofort abgeschnitten.
Für jedes Profil berechnet der Code
$$\sum_{\substack{k=3\\k\text{ ungerade}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$
ohne über alle ungeraden \(k\) zu iterieren. Große ungerade Exponenten teilen denselben Nenner \(d\,\operatorname{rad}(d)\) und werden daher in einem einzigen Block zusammengefasst; nur die endlich vielen kleineren Schwellenwerte werden einzeln ausgewertet. Nach Multiplikation mit \(\varphi(d)\) und Summation über alle Profile erhält die Implementierung \(H(N)\), berechnet \(T(N)=A(N)-H(N)\) modulo \(2M\) mit \(M=977676779\) und wandelt dieses Ergebnis anschließend in
$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M$$
um. Die C++- und Java-Implementierungen parallelisieren zusätzlich die obersten Primzahlzweige; die Python-Implementierung verwendet dieselbe Arithmetik in einer einzelnen rekursiven Traversierung.
Das Primzahlsieb bis \(\sqrt N\) benötigt \(O(\sqrt N)\) Zeit und \(O(\sqrt N)\) Speicher. Der schwierigere Teil ist die rekursive Erzeugung der Teilerprofile. Definiert man
$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$
dann besucht die Rekursion genau die zulässigen Profile in \(\mathcal{D}(N)\), und für jedes Profil werden nur die endlich vielen Odd-\(k\)-Schwellen ausgewertet, an denen sich ein Ceiling-Term \(\left\lceil a/(k-1)\right\rceil\) ändert. Die Laufzeit ist daher ungefähr
$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right),$$
wobei \(B(d)\) die Anzahl der durch die Exponenten von \(d\) induzierten Schwellen bezeichnet. Dafür gibt es keine einfache geschlossene Elementarform, aber diese Größe ist drastisch kleiner als der ursprüngliche \(O(N^2)\)-Ansatz über alle Paare \((k,n)\). Der Speicherbedarf bleibt bei \(O(\sqrt N)\) für die Primtabelle plus \(O(\log N)\) Rekursionstiefe.
Tek \(k\) için
$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\}$$
tanımlanır; burada \(\{x\}\), \(x\)'in kesir kısmıdır. Problem daha sonra
$$S(N)=\sum_{\substack{k=1\\k\text{ tek}}}^{N}\sum_{n=1}^{N} f_k(n)$$
toplamını ele alır ve özellikle \(\lfloor S(33557799775533)\rfloor \bmod 977676779\) değerini ister. Bu ölçekte tüm tek \(k\) ve tüm \(n\) üzerinde doğrudan çift döngü kurmak imkansızdır; bu yüzden çözüm, kesirli kısımları düzenli bir bölünebilme problemine çevirir.
C++, Python ve Java implementasyonları aynı indirgeme zincirini kullanır: önce kesir kısımlarını ortadan kaldırır, sonra kalan ifadeyi asal üslerle yazar ve en sonunda tüm toplamı bölen profilleri üzerinden yeniden organize eder.
\(k\) tek olduğu için
$$ (n-i)^k\equiv -\,i^k \pmod n. $$
Dolayısıyla \(\left\{\frac{i^k}{n}\right\}\) ile \(\left\{\frac{(n-i)^k}{n}\right\}\) toplamı normalde \(1\) eder; yalnızca \(n\mid i^k\) olduğunda her iki terim de \(0\) olur. Şunu tanımlayalım:
$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$
O zaman eşleştirme özdeşliği
$$2f_k(n)=n-z_k(n)$$
şeklini alır. Bu ilk büyük sadeleştirmedir: artık kesir kısımlarını değil, tam bölünebilirlik üreten kalıntıların sayısını hesaplıyoruz.
\(n\)'yi
$$n=\prod_{p} p^{e_p}$$
şeklinde yazalım. \(n\mid i^k\) koşulu şu koşulla denktir:
$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad\text{her }p\mid n\text{ için}.$$
Yani \(i\),
$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}$$
sayısının katı olmalıdır. \(1,2,\dots,n\) arasında tam olarak \(n/m_k(n)\) tane böyle sayı vardır. Demek ki
$$z_k(n)=\frac{n}{m_k(n)}$$
ve dolayısıyla
$$2f_k(n)=n-\frac{n}{m_k(n)}.$$
Şunu tanımlayalım:
$$T(N)=2S(N).$$
Önceki eşitlikten
$$T(N)=\sum_{\substack{k=1\\k\text{ tek}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N)$$
elde edilir; burada
$$A(N)=\sum_{\substack{k=1\\k\text{ tek}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$
gerçek tek girdi \(N\) için kapalı formdur ve
$$H(N)=\sum_{\substack{k=1\\k\text{ tek}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}$$
zor kısımdır. Böylece bütün problem \(H(N)\)'yi hızlı hesaplama işine indirgenmiş olur.
Implementasyon \(\frac{n}{m_k(n)}\)'yi doğrudan toplamıyor. Bunun yerine tek \(k\ge 3\) için
$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}$$
tanımını yapıyor. Sonra şu özdeşlik kullanılıyor:
$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$
burada \(\varphi\), Euler phi fonksiyonudur.
Bu neden doğru? Tek bir asal kuvvette kontrol etmek yeterlidir. \(n=p^E\) ise, \(d=p^a\) içindeki uygun üsler tam olarak
$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E$$
koşulunu sağlayanlardır. En büyük uygun \(a\), \(E-\left\lceil E/k\right\rceil\) olur. Böylece
$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$
Her iki taraf da çarpımsal olduğu için genel \(n\) için de formül geçerlidir. \(k=1\) durumu daha da basittir: \(m_1(n)=n\), yani \(\frac{n}{m_1(n)}=1\).
Bölen özdeşliğini yerine koyunca
$$H(N)=N+\sum_{\substack{k=3\\k\text{ tek}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$
elde edilir. Şimdi bir bölen profili sabitleyelim:
$$d=\prod_{j=1}^{r} p_j^{a_j}.$$
Yalnızca
$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$
koşulunu sağlayan profiller katkı verebilir; çünkü her tek \(k\ge 3\) için \(\rho_k(d)\ge \operatorname{rad}(d)\) olur. Bu da şu anlama gelir: \(d\)'de görünen her asal en fazla \(\sqrt N\) olabilir. Implementasyonun yalnızca \(\sqrt N\)'ye kadar asal elemesi yapmasının nedeni tam olarak budur.
Sabit bir \(d\) için \(\left\lceil a_j/(k-1)\right\rceil\) değerleri her tek \(k\)'de değişmez; sadece sonlu sayıda eşikte değişirler. \(k\), en büyük üs \(a_j\)'den büyük olduğunda tüm tavan ifadeleri \(1\)'e düşer ve
$$\rho_k(d)=\operatorname{rad}(d)$$
olur. Böylece büyük tek \(k\)'lerin tamamı tek bir blok halinde toplanabilir; yalnızca küçük tek \(k\)'ler tek tek işlenir. Problemi uygulanabilir kılan sıkıştırma tam olarak budur.
Şunu alalım:
$$72=2^3\cdot 3^2.$$
\(k=3\) için
$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6$$
olur. O halde
$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$
Bölen profili formülü de aynı sonucu verir. Burada
$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$
\(2\) asalı için uygun üsler \(a+\lceil a/2\rceil\le 3\) koşulundan \(a\in\{0,1,2\}\) olur. \(3\) asalı için ise \(a+\lceil a/2\rceil\le 2\) olduğundan \(a\in\{0,1\}\). Bu yüzden
$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$
ve bu da tam olarak \(\frac{72}{6}\)'ya eşittir.
Implementasyonlar önce lineer bir elekle \(\sqrt N\)'ye kadar tüm asalları üretir. Bu yeterlidir; çünkü katkı veren her bölen profili \(d\,\operatorname{rad}(d)\le N\) koşulunu sağlamak zorundadır ve bu nedenle \(\sqrt N\)'den büyük bir asal, önemsiz olmayan bir profilde yer alamaz.
Daha sonra bölen profilleri \(d=\prod p^a\), asallar artan sırada olacak şekilde özyinelemeli olarak gezilir. Bu gezi sırasında dört çarpımsal durum tutulur: \(d\)'nin değeri, \(\operatorname{rad}(d)\), \(\varphi(d)\) ve profildeki en büyük üs. \(d\,\operatorname{rad}(d)>N\) olduğu anda dal anında budanır.
Her profil için kod
$$\sum_{\substack{k=3\\k\text{ tek}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$
ifadesini bütün tek \(k\)'leri dolaşmadan hesaplar. Büyük tek üsler aynı paydayı \(d\,\operatorname{rad}(d)\) paylaştığı için tek bir aritmetik blokta toplanır; yalnızca sonlu sayıdaki küçük eşik değerler tek tek işlenir. Tüm profiller \(\varphi(d)\) ile ağırlıklandırılıp toplandıktan sonra implementasyon \(H(N)\)'yi elde eder, \(M=977676779\) olmak üzere \(T(N)=A(N)-H(N)\) değerini \(2M\) modunda hesaplar ve son olarak bunu
$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M$$
şeklinde sonuca çevirir. C++ ve Java sürümleri ayrıca en üst düzey asal dallarını iş parçacıkları arasında paylaştırır; Python sürümü aynı aritmetiği tek iş parçacığında uygular.
\(\sqrt N\)'ye kadar asal eleği \(O(\sqrt N)\) zaman ve \(O(\sqrt N)\) bellek kullanır. Daha zor kısım, bölen profillerinin özyinelemeli olarak gezilmesidir. Şunu tanımlarsak
$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$
özyineleme tam olarak \(\mathcal{D}(N)\) içindeki uygun profilleri ziyaret eder; her profil için de yalnızca \(\left\lceil a/(k-1)\right\rceil\) değerlerinden biri değiştiğinde ortaya çıkan sonlu sayıdaki tek-\(k\) eşikleri incelenir. Bu yüzden çalışma süresi yaklaşık olarak
$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right)$$
şeklindedir; burada \(B(d)\), \(d\)'nin üslerinden doğan farklı tek-\(k\) eşik sayısıdır. Bu büyüklüğün basit bir kapalı formu yoktur, fakat yine de tüm \((k,n)\) çiftleri üzerinde \(O(N^2)\) dolaşmaktan dramatik biçimde küçüktür. Bellek kullanımı asal tablosu için \(O(\sqrt N)\), özyineleme derinliği için ise \(O(\log N)\) seviyesinde kalır.
Para \(k\) impar se define
$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$
donde \(\{x\}\) denota la parte fraccionaria de \(x\). El problema pide entonces
$$S(N)=\sum_{\substack{k=1\\k\text{ impar}}}^{N}\sum_{n=1}^{N} f_k(n),$$
y en particular \(\lfloor S(33557799775533)\rfloor \bmod 977676779\). Un doble bucle directo sobre todos los \(k\) impares y todos los \(n\) es imposible a esta escala, así que la solución transforma la suma de partes fraccionarias en un problema de divisibilidad mucho más rígido.
Las implementaciones en C++, Python y Java siguen la misma secuencia: eliminar las partes fraccionarias, describir el término restante mediante exponentes primos y después reorganizar todo por perfiles de divisores.
Como \(k\) es impar,
$$ (n-i)^k\equiv -\,i^k \pmod n. $$
Por tanto, \(\left\{\frac{i^k}{n}\right\}\) y \(\left\{\frac{(n-i)^k}{n}\right\}\) suman \(1\), salvo cuando \(n\mid i^k\), caso en el que ambas partes fraccionarias son \(0\). Definimos
$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$
Entonces la identidad de emparejamiento es
$$2f_k(n)=n-z_k(n).$$
Esta es la primera simplificación importante: el problema deja de tratar sobre partes fraccionarias y pasa a tratar sobre cuántos residuos producen divisibilidad exacta.
Escribimos
$$n=\prod_{p} p^{e_p}.$$
La condición \(n\mid i^k\) equivale a
$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad\text{para todo primo }p\mid n.$$
Así, \(i\) debe ser múltiplo de
$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}.$$
Entre \(1,2,\dots,n\), hay exactamente \(n/m_k(n)\) múltiplos de ese número. Por consiguiente,
$$z_k(n)=\frac{n}{m_k(n)}$$
y, por tanto,
$$2f_k(n)=n-\frac{n}{m_k(n)}.$$
Definimos
$$T(N)=2S(N).$$
Usando la identidad anterior,
$$T(N)=\sum_{\substack{k=1\\k\text{ impar}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N),$$
donde
$$A(N)=\sum_{\substack{k=1\\k\text{ impar}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$
para la entrada real, que es impar, y
$$H(N)=\sum_{\substack{k=1\\k\text{ impar}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}.$$
Así, toda la dificultad se concentra en calcular \(H(N)\) con rapidez.
La implementación no suma \(\frac{n}{m_k(n)}\) de forma directa. Para \(k\) impar y \(k\ge 3\), introduce
$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}.$$
Entonces vale la identidad
$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$
donde \(\varphi\) es la función phi de Euler.
¿Por qué? Basta revisar una sola potencia prima. Si \(n=p^E\), los exponentes \(a\) admisibles en \(d=p^a\) son exactamente los que satisfacen
$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E.$$
El mayor \(a\) posible es \(E-\left\lceil E/k\right\rceil\). Por ello,
$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$
Como ambos lados son multiplicativos en \(n\), la fórmula vale en general. Para \(k=1\), el caso es todavía más sencillo: \(m_1(n)=n\), así que \(\frac{n}{m_1(n)}=1\).
Sustituyendo la identidad de divisores se obtiene
$$H(N)=N+\sum_{\substack{k=3\\k\text{ impar}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$
Fijemos ahora un perfil de divisor
$$d=\prod_{j=1}^{r} p_j^{a_j}.$$
Solo pueden contribuir los perfiles que cumplan
$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$
porque siempre se tiene \(\rho_k(d)\ge \operatorname{rad}(d)\) para todo \(k\) impar con \(k\ge 3\). Eso implica inmediatamente que todo primo que aparezca en \(d\) es a lo sumo \(\sqrt N\), y por eso la criba de primos solo llega hasta \(\sqrt N\).
Para un \(d\) fijo, las cantidades \(\left\lceil a_j/(k-1)\right\rceil\) no cambian en cada \(k\) impar, sino solo cuando \(k-1\) cruza ciertos umbrales finitos. En cuanto \(k\) supera el mayor exponente \(a_j\), todos esos techos pasan a ser \(1\), de modo que
$$\rho_k(d)=\operatorname{rad}(d).$$
Eso permite sumar toda la cola de \(k\) impares grandes en un solo bloque, mientras que los pocos \(k\) impares pequeños se tratan individualmente. Esa compresión es la clave de la eficiencia.
Tomemos
$$72=2^3\cdot 3^2.$$
Para \(k=3\),
$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6,$$
luego
$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$
La fórmula de perfiles de divisores produce exactamente lo mismo. Aquí
$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$
Para el primo \(2\), la condición \(a+\lceil a/2\rceil\le 3\) da \(a\in\{0,1,2\}\). Para el primo \(3\), la condición \(a+\lceil a/2\rceil\le 2\) da \(a\in\{0,1\}\). Por tanto,
$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$
que coincide con \(\frac{72}{6}\).
Las implementaciones generan primero todos los primos hasta \(\sqrt N\) con una criba lineal. Eso basta porque todo perfil \(d\) que realmente contribuya debe satisfacer \(d\,\operatorname{rad}(d)\le N\), de modo que ningún primo mayor que \(\sqrt N\) puede aparecer en un perfil no trivial.
Después enumeran recursivamente los perfiles \(d=\prod p^a\) en orden creciente de primos. Durante esa recursión se mantienen cuatro cantidades multiplicativas: el valor de \(d\), su radical \(\operatorname{rad}(d)\), su totiente \(\varphi(d)\) y el exponente máximo presente en el perfil. La poda principal es inmediata: si \(d\,\operatorname{rad}(d)>N\), esa rama se abandona.
Para cada perfil, el código evalúa
$$\sum_{\substack{k=3\\k\text{ impar}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$
sin recorrer todos los \(k\) impares. Los \(k\) grandes comparten el mismo denominador \(d\,\operatorname{rad}(d)\) y se acumulan en un único bloque aritmético; solo los pocos puntos de ruptura pequeños se procesan individualmente. Tras multiplicar por \(\varphi(d)\) y sumar todos los perfiles, la implementación obtiene \(H(N)\), calcula \(T(N)=A(N)-H(N)\) módulo \(2M\) con \(M=977676779\) y finalmente transforma ese residuo en
$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M.$$
Las versiones en C++ y Java además reparten entre hilos las ramas superiores definidas por el primer primo; la versión en Python usa la misma aritmética en un recorrido recursivo de un solo hilo.
La criba de primos hasta \(\sqrt N\) cuesta \(O(\sqrt N)\) tiempo y \(O(\sqrt N)\) memoria. La parte más delicada es la enumeración recursiva de perfiles de divisores. Si definimos
$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$
la recursión visita exactamente los perfiles admisibles de \(\mathcal{D}(N)\), y para cada uno solo evalúa los intervalos finitos de \(k\) impares donde cambia algún techo \(\left\lceil a/(k-1)\right\rceil\). Por eso el tiempo total es aproximadamente
$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right),$$
donde \(B(d)\) es el número de puntos de ruptura impares inducidos por los exponentes de \(d\). No existe una forma cerrada elemental sencilla para esta cantidad, pero sigue siendo muchísimo menor que el recorrido original \(O(N^2)\) sobre todos los pares \((k,n)\). La memoria se mantiene en \(O(\sqrt N)\) para la tabla de primos y \(O(\log N)\) para la profundidad de la recursión.
对奇数 \(k\)，定义
$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$
其中 \(\{x\}\) 表示 \(x\) 的小数部分。题目接着要求
$$S(N)=\sum_{\substack{k=1\\k\text{ 为奇数}}}^{N}\sum_{n=1}^{N} f_k(n),$$
并且特别关心 \(\lfloor S(33557799775533)\rfloor \bmod 977676779\)。如果直接对所有奇数 \(k\) 与所有 \(n\) 做双重枚举，规模完全不可接受，所以解法必须把这个分数和改写成一个结构更强的整除计数问题。
C++、Python 和 Java 的实现采用的是同一条思路链：先把小数部分消掉，再把剩余项写成素因子指数形式，最后按照“约数轮廓”重新组织求和。
因为 \(k\) 是奇数，所以
$$ (n-i)^k\equiv -\,i^k \pmod n. $$
因此 \(\left\{\frac{i^k}{n}\right\}\) 与 \(\left\{\frac{(n-i)^k}{n}\right\}\) 通常会凑成 \(1\)；唯一的例外是 \(n\mid i^k\)，这时两个小数部分都等于 \(0\)。定义
$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$
于是就得到配对恒等式
$$2f_k(n)=n-z_k(n).$$
这是第一步关键化简。问题不再是“如何直接处理小数部分”，而是“有多少个剩余类会让 \(i^k\) 被 \(n\) 整除”。
把 \(n\) 分解为
$$n=\prod_{p} p^{e_p}.$$
条件 \(n\mid i^k\) 等价于
$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad\text{对每个 }p\mid n.$$
也就是说，\(i\) 必须是
$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}$$
的倍数。在 \(1,2,\dots,n\) 里，这样的数恰好有 \(n/m_k(n)\) 个，因此
$$z_k(n)=\frac{n}{m_k(n)}$$
从而
$$2f_k(n)=n-\frac{n}{m_k(n)}.$$
定义
$$T(N)=2S(N).$$
由上式可得
$$T(N)=\sum_{\substack{k=1\\k\text{ 为奇数}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N),$$
其中
$$A(N)=\sum_{\substack{k=1\\k\text{ 为奇数}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$
对题目的实际奇数输入 \(N\) 是一个直接闭式，而
$$H(N)=\sum_{\substack{k=1\\k\text{ 为奇数}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}$$
才是真正困难的部分。于是整道题就被压缩成“怎样快速求出 \(H(N)\)”这一件事。
实现并不是直接对 \(\frac{n}{m_k(n)}\) 求和。对奇数 \(k\ge 3\)，它先定义
$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}.$$
然后使用恒等式
$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$
其中 \(\varphi\) 是欧拉函数。
这个恒等式为什么成立？只要先看一个素数幂即可。若 \(n=p^E\)，那么 \(d=p^a\) 中允许出现的指数 \(a\) 正好满足
$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E.$$
其中最大的可行 \(a\) 是 \(E-\left\lceil E/k\right\rceil\)。于是
$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$
左右两边对不同素数都是乘法性的，所以一般的 \(n\) 也成立。至于 \(k=1\)，情况更简单：\(m_1(n)=n\)，因此 \(\frac{n}{m_1(n)}=1\)。
把约数恒等式代回去以后得到
$$H(N)=N+\sum_{\substack{k=3\\k\text{ 为奇数}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$
现在固定一个约数轮廓
$$d=\prod_{j=1}^{r} p_j^{a_j}.$$
只有满足
$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$
的轮廓才可能有贡献，因为对所有奇数 \(k\ge 3\) 都有 \(\rho_k(d)\ge \operatorname{rad}(d)\)。这立刻说明：凡是在 \(d\) 中出现的素数都不可能超过 \(\sqrt N\)。这就是实现只需要筛到 \(\sqrt N\) 的根本原因。
对固定的 \(d\) 而言，\(\left\lceil a_j/(k-1)\right\rceil\) 不会在每一个奇数 \(k\) 都改变，它只会在有限个阈值上改变。一旦 \(k\) 大于轮廓中的最大指数 \(a_j\)，所有这些天花板值都会退化成 \(1\)，于是
$$\rho_k(d)=\operatorname{rad}(d).$$
这意味着所有足够大的奇数 \(k\) 都可以放进同一个尾段里一次累加，而只有前面少量较小的奇数 \(k\) 需要单独处理。整道题之所以能跑得动，靠的就是这一步压缩。
取
$$72=2^3\cdot 3^2.$$
当 \(k=3\) 时，
$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6,$$
所以
$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$
约数轮廓公式会给出完全相同的结果。此时
$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$
对素数 \(2\)，条件 \(a+\lceil a/2\rceil\le 3\) 给出 \(a\in\{0,1,2\}\)。对素数 \(3\)，条件 \(a+\lceil a/2\rceil\le 2\) 给出 \(a\in\{0,1\}\)。因此
$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$
恰好等于 \(\frac{72}{6}\)。这说明实现中的约数重排与原始定义是完全一致的。
实现的第一步是用线性筛生成所有不超过 \(\sqrt N\) 的素数。这已经足够，因为任何真正有贡献的约数轮廓 \(d\) 都必须满足 \(d\,\operatorname{rad}(d)\le N\)，于是其中不会出现大于 \(\sqrt N\) 的素数。
接下来，程序按素数递增顺序递归枚举所有可能的轮廓 \(d=\prod p^a\)。递归过程中持续维护四个乘法性状态：\(d\) 本身、\(\operatorname{rad}(d)\)、\(\varphi(d)\) 以及当前轮廓里的最大指数。只要发现 \(d\,\operatorname{rad}(d)>N\)，这一整条分支就立刻剪掉。
对每个轮廓，程序并不会逐个遍历所有奇数 \(k\)，而是计算
$$\sum_{\substack{k=3\\k\text{ 为奇数}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$
当 \(k\) 足够大时，分母统一变成 \(d\,\operatorname{rad}(d)\)，这一大段可以一次性求和；只有前面有限个真正发生变化的断点需要逐个处理。所有轮廓乘上 \(\varphi(d)\) 后累加，就得到 \(H(N)\)。然后实现按模 \(2M\)（其中 \(M=977676779\)）求出 \(T(N)=A(N)-H(N)\)，最后再把这个结果转成
$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M.$$
C++ 和 Java 实现还会把最外层以首个素数划分的分支分配到多个线程上；Python 实现则在单线程里执行同样的递归与模运算。
筛到 \(\sqrt N\) 的素数需要 \(O(\sqrt N)\) 时间和 \(O(\sqrt N)\) 空间。更难的部分是约数轮廓的递归枚举。若记
$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$
那么递归恰好访问 \(\mathcal{D}(N)\) 中所有可行轮廓；对每个轮廓，只需要处理有限个让 \(\left\lceil a/(k-1)\right\rceil\) 发生变化的奇数 \(k\) 断点。因此总时间大致可写成
$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right),$$
其中 \(B(d)\) 表示该轮廓由指数诱导出的不同奇数断点个数。这个量没有简单的初等闭式，但和原始的 \(O(N^2)\) 双重枚举相比已经小了很多个数量级。空间方面主要是 \(O(\sqrt N)\) 的素数表，以及 \(O(\log N)\) 的递归深度。
Для нечётного \(k\) определим
$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$
где \(\{x\}\) означает дробную часть числа \(x\). Далее требуется вычислить
$$S(N)=\sum_{\substack{k=1\\k\text{ нечётно}}}^{N}\sum_{n=1}^{N} f_k(n),$$
а именно \(\lfloor S(33557799775533)\rfloor \bmod 977676779\). На таком масштабе прямой двойной перебор по всем нечётным \(k\) и всем \(n\) невозможен, поэтому решение превращает сумму дробных частей в более жёсткую задачу о делимости.
Во всех трёх реализациях, на C++, Python и Java, используется одна и та же схема: сначала убрать дробные части, затем выразить оставшийся вклад через показатели простых множителей и, наконец, перегруппировать вычисление по профилям делителей.
Поскольку \(k\) нечётно, имеем
$$ (n-i)^k\equiv -\,i^k \pmod n. $$
Следовательно, \(\left\{\frac{i^k}{n}\right\}\) и \(\left\{\frac{(n-i)^k}{n}\right\}\) в сумме дают \(1\), кроме случая \(n\mid i^k\), когда обе дробные части равны \(0\). Обозначим
$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$
Тогда получается тождество
$$2f_k(n)=n-z_k(n).$$
Это первый ключевой шаг: вместо самих дробных частей нужно считать лишь количество остатков, дающих точную делимость.
Разложим
$$n=\prod_{p} p^{e_p}.$$
Условие \(n\mid i^k\) равносильно системе
$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad\text{для каждого простого }p\mid n.$$
Значит, \(i\) должен быть кратен числу
$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}.$$
Среди чисел \(1,2,\dots,n\) таких кратных ровно \(n/m_k(n)\). Поэтому
$$z_k(n)=\frac{n}{m_k(n)}$$
и, следовательно,
$$2f_k(n)=n-\frac{n}{m_k(n)}.$$
Положим
$$T(N)=2S(N).$$
Тогда из предыдущего равенства следует
$$T(N)=\sum_{\substack{k=1\\k\text{ нечётно}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N),$$
где
$$A(N)=\sum_{\substack{k=1\\k\text{ нечётно}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$
для реального нечётного входного значения \(N\), а
$$H(N)=\sum_{\substack{k=1\\k\text{ нечётно}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}.$$
Итак, всё сводится к быстрому вычислению \(H(N)\).
Реализация не суммирует \(\frac{n}{m_k(n)}\) напрямую. Для нечётного \(k\ge 3\) вводится функция
$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}.$$
После этого используется тождество
$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$
где \(\varphi\) обозначает функцию Эйлера.
Почему это верно? Достаточно проверить один простой множитель. Если \(n=p^E\), то допустимые показатели \(a\) в \(d=p^a\) — это в точности те, для которых
$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E.$$
Максимальное такое \(a\) равно \(E-\left\lceil E/k\right\rceil\). Поэтому
$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$
Поскольку обе стороны мультипликативны по \(n\), формула верна и в общем случае. Для \(k=1\) всё ещё проще: \(m_1(n)=n\), значит \(\frac{n}{m_1(n)}=1\).
Подстановка тождества по делителям даёт
$$H(N)=N+\sum_{\substack{k=3\\k\text{ нечётно}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$
Зафиксируем один профиль делителя
$$d=\prod_{j=1}^{r} p_j^{a_j}.$$
Вклад могут давать только профили, для которых
$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$
потому что всегда \(\rho_k(d)\ge \operatorname{rad}(d)\) при нечётном \(k\ge 3\). Отсюда сразу следует: любой простой множитель, появляющийся в \(d\), не превосходит \(\sqrt N\). Именно поэтому в коде достаточно просеять простые числа только до \(\sqrt N\).
Для фиксированного \(d\) значения \(\left\lceil a_j/(k-1)\right\rceil\) меняются не на каждом нечётном \(k\), а только на конечном наборе порогов. Как только \(k\) становится больше наибольшего показателя \(a_j\), все такие Ceiling-выражения становятся равными \(1\), и тогда
$$\rho_k(d)=\operatorname{rad}(d).$$
Значит, весь хвост больших нечётных \(k\) можно просуммировать одним блоком, а отдельно обрабатывать нужно только небольшое число малых нечётных \(k\). Именно эта компрессия и делает задачу вычислимой.
Возьмём
$$72=2^3\cdot 3^2.$$
Для \(k=3\)
$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6,$$
поэтому
$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$
Формула по профилям делителей даёт тот же результат. Здесь
$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$
Для простого \(2\) условие \(a+\lceil a/2\rceil\le 3\) даёт \(a\in\{0,1,2\}\). Для простого \(3\) условие \(a+\lceil a/2\rceil\le 2\) даёт \(a\in\{0,1\}\). Следовательно,
$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$
что в точности совпадает с \(\frac{72}{6}\).
Сначала реализации строят линейным решетом все простые числа до \(\sqrt N\). Этого достаточно, потому что любой профиль \(d\), который реально вносит вклад, обязан удовлетворять ограничению \(d\,\operatorname{rad}(d)\le N\), а значит, в нетривиальном профиле не может появиться простой множитель больше \(\sqrt N\).
Далее рекурсивно перечисляются все профили вида \(d=\prod p^a\) в порядке возрастания простых чисел. Во время рекурсии поддерживаются четыре мультипликативные характеристики состояния: значение \(d\), радикал \(\operatorname{rad}(d)\), тотиент \(\varphi(d)\) и максимальный показатель, уже встреченный в профиле. Как только выполняется \(d\,\operatorname{rad}(d)>N\), ветвь немедленно отсекается.
Для каждого профиля код вычисляет
$$\sum_{\substack{k=3\\k\text{ нечётно}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$
без прохода по всем нечётным \(k\). Большие нечётные показатели имеют один и тот же знаменатель \(d\,\operatorname{rad}(d)\), поэтому их вклад суммируется одним арифметическим блоком; отдельно рассматриваются лишь немногие малые точки излома. После умножения на \(\varphi(d)\) и суммирования по всем профилям получается \(H(N)\), затем вычисляется \(T(N)=A(N)-H(N)\) по модулю \(2M\), где \(M=977676779\), и в конце результат переводится в
$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M.$$
В версиях на C++ и Java верхние ветви, задаваемые первым простым множителем, дополнительно распараллеливаются; версия на Python выполняет ту же арифметику в одном рекурсивном проходе.
Решето до \(\sqrt N\) требует \(O(\sqrt N)\) времени и \(O(\sqrt N)\) памяти. Более сложная часть — рекурсивный перебор профилей делителей. Если обозначить
$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$
то рекурсия посещает ровно допустимые профили из \(\mathcal{D}(N)\), а для каждого профиля рассматривает только конечное число интервалов нечётных \(k\), на которых меняется хотя бы одно выражение \(\left\lceil a/(k-1)\right\rceil\). Поэтому время работы можно оценить как
$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right),$$
где \(B(d)\) — число различных нечётных точек излома, порождённых показателями профиля \(d\). У этой величины нет простой элементарной замкнутой формы, но она на порядки меньше исходного перебора \(O(N^2)\) по всем парам \((k,n)\). Память остаётся на уровне \(O(\sqrt N)\) для таблицы простых и \(O(\log N)\) для глубины рекурсии.
لكل \(k\) فردي نعرّف
$$f_k(n)=\sum_{i=1}^{n}\left\{\frac{i^k}{n}\right\},$$
حيث ترمز \(\{x\}\) إلى الجزء الكسري من \(x\). ثم تطلب المسألة
$$S(N)=\sum_{\substack{k=1\\k\ \mathrm{odd}}}^{N}\sum_{n=1}^{N} f_k(n),$$
وبالتحديد القيمة \(\lfloor S(33557799775533)\rfloor \bmod 977676779\). من المستحيل عمليًا إجراء حلقة مزدوجة مباشرة على جميع قيم \(k\) الفردية وجميع قيم \(n\) عند هذا الحجم، لذلك يحوّل الحل مجموع الأجزاء الكسرية إلى مسألة قسمة أكثر انتظامًا.
تستخدم تطبيقات C++ وPython وJava السلسلة نفسها من الاختزالات: إزالة الأجزاء الكسرية أولًا، ثم كتابة ما تبقى بلغة أسس العوامل الأولية، ثم إعادة تنظيم المجموع كله على أساس ملفات القواسم.
بما أن \(k\) فردي، فإن
$$ (n-i)^k\equiv -\,i^k \pmod n. $$
ومن ثم فإن \(\left\{\frac{i^k}{n}\right\}\) و\(\left\{\frac{(n-i)^k}{n}\right\}\) يجتمعان عادة إلى \(1\)، إلا إذا كان \(n\mid i^k\)، وعندها يكون الجزآن الكسريان كلاهما \(0\). لنعرّف
$$z_k(n)=\#\left\{1\le i\le n:\ n\mid i^k\right\}.$$
عندئذ نحصل على الهوية
$$2f_k(n)=n-z_k(n).$$
وهذه أول نقلة أساسية: لم نعد بحاجة إلى التعامل المباشر مع الأجزاء الكسرية، بل يكفي أن نعرف كم عدد البواقي التي تحقق القسمة التامة.
لنكتب
$$n=\prod_{p} p^{e_p}.$$
إن الشرط \(n\mid i^k\) يكافئ
$$v_p(i)\ge \left\lceil\frac{e_p}{k}\right\rceil\qquad(\forall\, p\mid n).$$
إذًا يجب أن يكون \(i\) من مضاعفات
$$m_k(n)=\prod_{p\mid n} p^{\left\lceil e_p/k\right\rceil}.$$
ومن بين الأعداد \(1,2,\dots,n\) يوجد بالضبط \(n/m_k(n)\) عددًا من هذا النوع. لذلك
$$z_k(n)=\frac{n}{m_k(n)}$$
ومن ثم
$$2f_k(n)=n-\frac{n}{m_k(n)}.$$
لنعرّف
$$T(N)=2S(N).$$
باستخدام العلاقة السابقة نحصل على
$$T(N)=\sum_{\substack{k=1\\k\ \mathrm{odd}}}^{N}\sum_{n=1}^{N}\left(n-\frac{n}{m_k(n)}\right)=A(N)-H(N),$$
حيث
$$A(N)=\sum_{\substack{k=1\\k\ \mathrm{odd}}}^{N}\sum_{n=1}^{N} n=\frac{N+1}{2}\cdot\frac{N(N+1)}{2}$$
لأن قيمة \(N\) الفعلية في المسألة فردية، بينما
$$H(N)=\sum_{\substack{k=1\\k\ \mathrm{odd}}}^{N}\sum_{n=1}^{N}\frac{n}{m_k(n)}.$$
وهكذا تنحصر الصعوبة الفعلية في حساب \(H(N)\) بكفاءة.
لا تقوم الشيفرة بجمع \(\frac{n}{m_k(n)}\) مباشرة. بل عند \(k\) الفردي مع \(k\ge 3\) تعرّف
$$\rho_k(d)=\prod_{p^a\parallel d} p^{\left\lceil a/(k-1)\right\rceil}.$$
ثم تستخدم الهوية
$$\frac{n}{m_k(n)}=\sum_{\substack{d\ge 1\\ d\,\rho_k(d)\mid n}} \varphi(d),$$
حيث \(\varphi\) هي دالة أويلر التاتية.
لماذا هذه الهوية صحيحة؟ يكفي فحص قوة أولية واحدة. إذا كان \(n=p^E\)، فإن الأسس المسموح بها \(a\) في \(d=p^a\) هي بالضبط تلك التي تحقق
$$a+\left\lceil\frac{a}{k-1}\right\rceil\le E.$$
وأكبر \(a\) ممكن هو \(E-\left\lceil E/k\right\rceil\). لذلك
$$\sum_{\substack{a\ge 0\\ a+\lceil a/(k-1)\rceil\le E}} \varphi(p^a)=\sum_{a=0}^{E-\lceil E/k\rceil}\varphi(p^a)=p^{E-\lceil E/k\rceil}=\frac{p^E}{p^{\lceil E/k\rceil}}.$$
وبما أن الطرفين ضربيان بالنسبة إلى \(n\)، فإن الصيغة تصح في الحالة العامة أيضًا. أما عند \(k=1\) فالوضع أبسط بكثير: \(m_1(n)=n\)، ومن ثم \(\frac{n}{m_1(n)}=1\).
بالتعويض نحصل على
$$H(N)=N+\sum_{\substack{k=3\\k\ \mathrm{odd}}}^{N}\ \sum_{d\ge 1}\varphi(d)\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor.$$
ثبّت الآن ملف قاسم واحدًا
$$d=\prod_{j=1}^{r} p_j^{a_j}.$$
لا يمكن أن يساهم إلا ما يحقق
$$d\,\operatorname{rad}(d)\le N,\qquad \operatorname{rad}(d)=\prod_{p\mid d}p,$$
لأن \(\rho_k(d)\ge \operatorname{rad}(d)\) دائمًا لكل \(k\) فردي مع \(k\ge 3\). وهذا يفرض فورًا أن كل أولي يظهر في \(d\) لا بد أن يكون \(\le \sqrt N\)، ولهذا تكتفي الشيفرة بغربلة الأعداد الأولية حتى \(\sqrt N\).
بالنسبة إلى \(d\) ثابت، فإن القيم \(\left\lceil a_j/(k-1)\right\rceil\) لا تتغير عند كل \(k\) فردي، بل فقط عندما يعبر \(k-1\) عددًا محدودًا من العتبات. وما إن تصبح \(k\) أكبر من أكبر أس \(a_j\) في الملف، حتى تنهار جميع تلك السقوف إلى \(1\)، أي يصبح
$$\rho_k(d)=\operatorname{rad}(d).$$
وبهذا يمكن جمع الذيل الكامل للقيم الفردية الكبيرة لـ \(k\) في كتلة واحدة، بينما تُعالج فقط القيم الفردية الصغيرة القليلة على نحو منفصل. هذا هو الضغط الذي يجعل الخوارزمية قابلة للتنفيذ.
لنأخذ
$$72=2^3\cdot 3^2.$$
عند \(k=3\) نحصل على
$$m_3(72)=2^{\lceil 3/3\rceil}3^{\lceil 2/3\rceil}=2\cdot 3=6,$$
ومن ثم
$$z_3(72)=\frac{72}{6}=12,\qquad 2f_3(72)=72-12=60,\qquad f_3(72)=30.$$
وصيغة ملفات القواسم تعطي النتيجة نفسها تمامًا. هنا
$$\rho_3(d)=\prod_{p^a\parallel d} p^{\lceil a/2\rceil}.$$
بالنسبة إلى الأولي \(2\)، فإن الشرط \(a+\lceil a/2\rceil\le 3\) يعطي \(a\in\{0,1,2\}\). وبالنسبة إلى الأولي \(3\)، فإن الشرط \(a+\lceil a/2\rceil\le 2\) يعطي \(a\in\{0,1\}\). لذلك
$$\sum_{\substack{d\,\rho_3(d)\mid 72}}\varphi(d)=\left(1+\varphi(2)+\varphi(4)\right)\left(1+\varphi(3)\right)=(1+1+2)(1+2)=12,$$
وهو بالضبط \(\frac{72}{6}\).
تبدأ التطبيقات بتوليد جميع الأعداد الأولية حتى \(\sqrt N\) باستخدام غربال خطي. وهذا يكفي لأن أي ملف قاسم \(d\) يساهم فعليًا لا بد أن يحقق \(d\,\operatorname{rad}(d)\le N\)، وبالتالي لا يمكن أن يحتوي ملف غير تافه على أولي أكبر من \(\sqrt N\).
بعد ذلك تُعدّد ملفات القواسم \(d=\prod p^a\) عوديًا بترتيب تصاعدي في الأعداد الأولية. وخلال العودية تُحفظ أربع كميات ضربية: قيمة \(d\)، وجذرها الخالي من التكرار \(\operatorname{rad}(d)\)، وقيمة \(\varphi(d)\)، وأكبر أس موجود في الملف الحالي. وكلما تحقق الشرط \(d\,\operatorname{rad}(d)>N\) تُقصّ تلك الشعبة مباشرة.
لكل ملف، تحسب الشيفرة
$$\sum_{\substack{k=3\\k\ \mathrm{odd}}}^{N}\left\lfloor\frac{N}{d\,\rho_k(d)}\right\rfloor$$
من دون المرور على كل \(k\) فردي. فالقيم الفردية الكبيرة تشترك في المقام نفسه \(d\,\operatorname{rad}(d)\)، ولذلك تُجمع في كتلة حسابية واحدة، بينما لا يُعالَج فرديًا إلا عدد محدود من نقاط الانكسار الصغيرة. وبعد ضرب كل ملف في \(\varphi(d)\) وجمع جميع الملفات، تحصل التنفيذات على \(H(N)\)، ثم تحسب \(T(N)=A(N)-H(N)\) بترديد \(2M\) حيث \(M=977676779\)، وأخيرًا تُحوِّل هذا المقدار إلى
$$\left\lfloor\frac{T(N)}{2}\right\rfloor \bmod M=\lfloor S(N)\rfloor \bmod M.$$
كما أن تطبيقي C++ وJava يوزعان الفروع العليا المحددة بواسطة الأولي الأول على عدة خيوط، بينما يستخدم تطبيق Python الحساب نفسه في عبور عودي أحادي الخيط.
توليد الأعداد الأولية حتى \(\sqrt N\) يكلف \(O(\sqrt N)\) زمنًا و\(O(\sqrt N)\) ذاكرة. أما الجزء الأصعب فهو تعداد ملفات القواسم عوديًا. إذا عرّفنا
$$\mathcal{D}(N)=\left\{d\ge 1:\ d\,\operatorname{rad}(d)\le N\right\},$$
فإن العودية تزور بالضبط الملفات المسموح بها داخل \(\mathcal{D}(N)\)، ولكل ملف لا تتعامل إلا مع العدد المحدود من فترات \(k\) الفردية التي يتغير عندها أحد الحدود \(\left\lceil a/(k-1)\right\rceil\). لذلك يمكن توصيف زمن التنفيذ تقريبًا بـ
$$O\!\left(\sqrt N+\sum_{d\in\mathcal{D}(N)}(1+B(d))\right),$$
حيث تمثل \(B(d)\) عدد نقاط الانكسار الفردية المختلفة التي تفرضها أسس الملف \(d\). لا توجد صيغة ابتدائية مغلقة بسيطة لهذه الكمية، لكنها ما تزال أصغر بكثير من المرور الأصلي \(O(N^2)\) على جميع الأزواج \((k,n)\). أما الذاكرة فتظل عند \(O(\sqrt N)\) لجدول الأعداد الأولية، مع عمق عودية مقداره \(O(\log N)\).