Detailed mathematical approach

Problem Summary

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.

Mathematical Approach

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.

Step 1: Pair \(i\) with \(n-i\)

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.

Step 2: Count the divisible residues prime by prime

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)}.$

Step 3: Split off the easy closed form

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.

Step 4: Rewrite \(\frac{n}{m_k(n)}\) as a totient-weighted divisor sum

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\).

Step 5: Compress all odd \(k\) with the same ceiling pattern

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.

Worked Example: \(n=72\) and \(k=3\)

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.

How the Code Works

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.

Complexity Analysis

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.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=850
2. Fractional part: Wikipedia — Fractional part
3. \(p\)-adic valuation: Wikipedia — p-adic valuation
4. Euler's totient function: Wikipedia — Euler's totient function
5. Ceiling function: Wikipedia — Floor and ceiling functions