Problem 875: Quadruple Congruence

View on Project Euler

Project Euler Problem 875 Solution

EulerSolve provides an optimized solution for Project Euler Problem 875, Quadruple Congruence, with C++, Python, Java, and a step-by-step mathematical explanation.

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

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

Problem 875 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

using u64 = std::uint64_t;
using u128 = unsigned __int128;

static constexpr int LIMIT = 12'345'678;
static constexpr u64 MOD = 1'001'961'001ULL;

static u64 q_prime_power_mod(int p, int e) {
    if (e == 0) return 1;

    if (p == 2) {
        u64 q = 128 % MOD;
        u64 add = 2'048 % MOD;
        for (int k = 2; k <= e; ++k) {
            q = (q * 128ULL + add) % MOD;
            add = (add * 16ULL) % MOD;
        }
        return q;
    }

    u64 pm = static_cast<u64>(p) % MOD;
    u64 p2 = (pm * pm) % MOD;
    u64 p3 = (p2 * pm) % MOD;
    u64 p4 = (p2 * p2) % MOD;
    u64 p7 = (p4 * p3) % MOD;
    u64 coeff = (pm + MOD - 1) % MOD;

    u64 q = 1;
    u64 term = p3;
    for (int k = 1; k <= e; ++k) {
        q = (q * p7 + coeff * term) % MOD;
        term = (term * p4) % MOD;
    }
    return q;
}

static u128 q_prime_power_exact(int p, int e) {
    if (e == 0) return 1;

    if (p == 2) {
        u128 q = 128;
        u128 add = 2048;
        for (int k = 2; k <= e; ++k) {
            q = q * 128 + add;
            add *= 16;
        }
        return q;
    }

    u128 p7 = 1;
    for (int i = 0; i < 7; ++i) p7 *= static_cast<u128>(p);
    u128 p4 = 1;
    for (int i = 0; i < 4; ++i) p4 *= static_cast<u128>(p);
    u128 term = static_cast<u128>(p) * static_cast<u128>(p) * static_cast<u128>(p);

    u128 q = 1;
    for (int k = 1; k <= e; ++k) {
        q = q * p7 + static_cast<u128>(p - 1) * term;
        term *= p4;
    }
    return q;
}

static u128 q_exact_from_factorization(int n) {
    int x = n;
    u128 q = 1;
    for (int p = 2; 1LL * p * p <= x; ++p) {
        if (x % p != 0) continue;
        int e = 0;
        while (x % p == 0) {
            x /= p;
            ++e;
        }
        q *= q_prime_power_exact(p, e);
    }
    if (x > 1) q *= q_prime_power_exact(x, 1);
    return q;
}

static u64 q_bruteforce_small(int n) {
    std::vector<int> sq(n);
    for (int i = 0; i < n; ++i) sq[i] = static_cast<int>((1LL * i * i) % n);

    std::vector<u64> cnt(n, 0);
    for (int a = 0; a < n; ++a) {
        for (int b = 0; b < n; ++b) {
            for (int c = 0; c < n; ++c) {
                for (int d = 0; d < n; ++d) {
                    int s = sq[a] + sq[b] + sq[c] + sq[d];
                    s %= n;
                    ++cnt[s];
                }
            }
        }
    }

    u64 q = 0;
    for (u64 v : cnt) q += v * v;
    return q;
}

static u64 compute_Q_mod(int n) {
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 10);

    for (int i = 2; i <= n; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            long long v = 1LL * p * i;
            if (v > n) break;
            spf[static_cast<int>(v)] = p;
            if (p == spf[i]) break;
        }
    }

    u64 sum = 1;
    for (int i = 2; i <= n; ++i) {
        int x = i;
        u64 qn = 1;
        while (x > 1) {
            int p = spf[x];
            int e = 0;
            do {
                x /= p;
                ++e;
            } while (x > 1 && spf[x] == p);

            qn = (qn * q_prime_power_mod(p, e)) % MOD;
        }

        sum += qn;
        if (sum >= MOD) sum -= MOD;
    }

    return sum;
}

int main() {
    assert(q_bruteforce_small(4) == 18'432ULL);
    for (int n = 1; n <= 8; ++n) {
        assert(q_bruteforce_small(n) == static_cast<u64>(q_exact_from_factorization(n)));
    }

    u128 q10 = 0;
    for (int i = 1; i <= 10; ++i) q10 += q_exact_from_factorization(i);
    assert(static_cast<u64>(q10) == 18'573'381ULL);
    assert(compute_Q_mod(10) == 18'573'381ULL);

    std::cout << compute_Q_mod(LIMIT) << '\n';
    return 0;
}

Python

def solve():
    LIMIT = 12345678; MOD = 1001961001

    def q_prime_power(p, e):
        if e == 0: return 1
        if p == 2:
            q = 128 % MOD; add = 2048 % MOD
            for _ in range(2, e+1):
                q = (q * 128 + add) % MOD; add = add * 16 % MOD
            return q
        pm = p % MOD; p2 = pm*pm%MOD; p3 = p2*pm%MOD; p4 = p2*p2%MOD; p7 = p4*p3%MOD
        coeff = (pm + MOD - 1) % MOD
        q = 1; term = p3
        for _ in range(1, e+1):
            q = (q * p7 + coeff * term) % MOD; term = term * p4 % MOD
        return q

    spf = list(range(LIMIT+1))
    primes = []
    for i in range(2, LIMIT+1):
        if spf[i] == i: primes.append(i)
        for p in primes:
            if i*p > LIMIT: break
            spf[i*p] = p
            if p == spf[i]: break

    s = 1
    for i in range(2, LIMIT+1):
        x = i; qn = 1
        while x > 1:
            p = spf[x]; e = 0
            while x > 1 and spf[x] == p: x //= p; e += 1
            qn = qn * q_prime_power(p, e) % MOD
        s = (s + qn) % MOD

    return str(s)

if __name__ == '__main__':
    print(solve())

Java

public class Euler875 {

    static final int LIMIT = 12345678;
    static final long MOD = 1001961001L;

    static long qPrimePowerMod(int p, int e) {
        if (e == 0)
            return 1;

        if (p == 2) {
            long q = 128 % MOD;
            long add = 2048 % MOD;
            for (int k = 2; k <= e; ++k) {
                q = (q * 128L + add) % MOD;
                add = (add * 16L) % MOD;
            }
            return q;
        }

        long pm = (p % MOD + MOD) % MOD;
        long p2 = (pm * pm) % MOD;
        long p3 = (p2 * pm) % MOD;
        long p4 = (p2 * p2) % MOD;
        long p7 = (p4 * p3) % MOD;
        long coeff = (pm + MOD - 1) % MOD;

        long q = 1;
        long term = p3;
        for (int k = 1; k <= e; ++k) {
            q = (q * p7 + coeff * term) % MOD;
            term = (term * p4) % MOD;
        }
        return q;
    }

    static long computeQMod(int n) {
        int[] spf = new int[n + 1];
        int[] primes = new int[n / 10 + 1000];
        int numPrimes = 0;

        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                if (numPrimes == primes.length) {
                    int[] newPrimes = new int[primes.length * 2];
                    System.arraycopy(primes, 0, newPrimes, 0, primes.length);
                    primes = newPrimes;
                }
                primes[numPrimes++] = i;
            }
            for (int j = 0; j < numPrimes; ++j) {
                int p = primes[j];
                long v = (long) p * i;
                if (v > n)
                    break;
                spf[(int) v] = p;
                if (p == spf[i])
                    break;
            }
        }

        long sum = 1;
        for (int i = 2; i <= n; ++i) {
            int x = i;
            long qn = 1;
            while (x > 1) {
                int p = spf[x];
                int e = 0;
                while (x > 1 && spf[x] == p) {
                    x /= p;
                    ++e;
                }
                qn = (qn * qPrimePowerMod(p, e)) % MOD;
            }
            sum += qn;
            if (sum >= MOD)
                sum -= MOD;
        }

        return sum;
    }

    public static String solve() {
        return Long.toString(computeQMod(LIMIT));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}