Detailed mathematical approach

Problem Summary

For each positive integer \(n\), define

$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}.$

The task is to evaluate

$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$

A direct search would require checking huge numbers of residue tuples for every modulus. The successful approach is to turn \(A(n)\) into a multiplicative arithmetic function and then compute its prime-power factors explicitly.

Mathematical Approach

The key observation is that the eight-variable congruence can be rewritten as a second moment of a four-variable residue count. After that, the Chinese remainder theorem and quadratic Gauss sums do the heavy lifting.

Step 1: Rewrite the congruence as a residue-distribution square

For a fixed modulus \(n\), let

$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}.$

If two quadruples produce the same residue \(t\), then together they form one admissible octuple. Therefore

$A(n)=\sum_{t\bmod n} r_n(t)^2.$

So the original count is the second moment of the distribution of four-square residues modulo \(n\).

Step 2: Use the Chinese remainder theorem to prove multiplicativity

Suppose \(\gcd(m,n)=1\). A residue class modulo \(mn\) corresponds uniquely to a pair of residue classes modulo \(m\) and modulo \(n\). Under this identification, a quadruple modulo \(mn\) is just an independent pair of quadruples modulo \(m\) and modulo \(n\), so

$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$

Squaring and summing over all residue pairs gives

$A(mn)=A(m)\,A(n).$

Hence it is enough to determine \(A(p^e)\) for prime powers, and then multiply the local factors together.

Step 3: Convert the residue count to quadratic Gauss sums

Introduce the additive characters

$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$

By orthogonality of additive characters,

$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut).$

Now take the second moment in \(t\). Parseval's identity for the finite Fourier transform yields

$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8.$

This formula is the decisive reduction: instead of counting octuples directly, we only need the sizes of local quadratic Gauss sums.

Step 4: Evaluate the odd-prime-power factor

Let \(n=p^e\) with \(p\) odd. Write \(u=p^jv\) where \(0\le j<e\) and \(p\nmid v\). Then

$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$

So every frequency with valuation \(j\) contributes

$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}.$

There are \((p-1)p^{e-j-1}\) such frequencies, while \(u=0\) contributes \(p^{8e}\). Substituting into the character formula gives

$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}.$

The same formula can be written recursively as

$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}.$

This is exactly the update rule used for odd prime powers in the implementations.

Step 5: Handle powers of \(2\) separately

The modulus \(2^e\) needs its own branch because quadratic Gauss sums over powers of \(2\) behave differently. For odd frequency modulo \(2^m\) with \(m\ge 2\), one has

$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2},$

whereas the remaining modulus \(2\) case vanishes. Therefore only valuations \(j=0,1,\dots,e-2\) contribute when \(u=2^jv\) with \(v\) odd. The resulting local factor is

$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}.$

In the iterative form used by the code, this becomes

$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$

Step 6: Reconstruct the value for a general modulus

If

$n=\prod_{i=1}^{k} p_i^{e_i},$

then multiplicativity gives

$A(n)=\prod_{i=1}^{k} A(p_i^{e_i}).$

So the entire problem reduces to

$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$

Worked Example: \(n=10\) and the checkpoint up to \(10\)

Because \(10=2\cdot 5\), multiplicativity immediately gives

$A(10)=A(2)\,A(5).$

From the local formulas,

$A(2)=2^7=128,$

$A(5)=5^7+4\cdot 5^3=78\,625.$

Hence

$A(10)=128\cdot 78\,625=10\,064\,000.$

Applying the same formulas to all \(n\le 10\) gives

$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$

and therefore

$S(10)=18\,573\,381.$

This is the small checkpoint verified by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They first build a smallest-prime-factor table up to \(12\,345\,678\). That table makes it possible to factor every \(n\) in a fast incremental scan.

For each \(n\), the implementation strips equal prime factors to recover the prime-power decomposition. Every local factor \(A(p^e)\) is then evaluated modulo \(1\,001\,961\,001\): odd primes use the recurrence from Step 4, and the factor \(2^e\) uses the special recurrence from Step 5.

The local factors are multiplied to obtain \(A(n)\bmod 1\,001\,961\,001\), and a running sum accumulates

$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}.$

The C++ implementation also performs tiny exact sanity checks on small moduli before the full computation, including \(A(4)=18\,432\) and \(S(10)=18\,573\,381\).

Complexity Analysis

Let \(N=12\,345\,678\). Building the smallest-prime-factor table with a linear sieve costs \(O(N)\) time and \(O(N)\) memory. Factoring all numbers up to \(N\) from that table requires total work proportional to the total number of prime factors encountered, which is \(O(N\log\log N)\) on average. The whole method is therefore near-linear in practice and uses \(O(N)\) memory.

Footnotes and References

1. Problem page: Project Euler 875
2. Chinese remainder theorem: Wikipedia — Chinese remainder theorem
3. Gauss sum: Wikipedia — Gauss sum
4. Parseval's theorem: Wikipedia — Parseval's theorem
5. Four-square theorem: Wikipedia — Lagrange's four-square theorem