Problem 717: Summation of a Modular Formula

Project Euler Problem 717 Solution

EulerSolve provides an optimized solution for Project Euler Problem 717, Summation of a Modular Formula, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each odd prime $p$, the computation associates an integer $g(p)$ obtained from a modular expression involving powers of $2$, and then sums these values: $$G(n)=\sum_{\substack{3\le p \lt n\\ p\text{ prime}}} g(p).$$ The difficulty is that the exponent $2^p$ is already enormous, and the quantity for each prime is not taken merely modulo $p$ but depends on how a congruence lifts from modulo $p$ to modulo $p^2$. The implemented method therefore combines exponent reduction, modular inversion, a controlled lift to $p^2$, and a prime sieve. Mathematical Approach Fix an odd prime $p$. The implementation computes $g(p)$ by first finding a residue modulo $p$, then lifting it to a compatible residue modulo $p^2$, and finally extracting the multiple of $p$ that measures the difference between the two levels. Step 1: Compress the Tower Exponent We begin with the residue $$b_p \equiv 2^{2^p}\pmod{p}.$$ Computing $2^{2^p}$ directly is hopeless, but Fermat's little theorem gives $2^{p-1}\equiv 1 \pmod p$, so exponents can be reduced modulo $p-1$. Therefore $$2^{2^p}\equiv 2^{\,2^p\bmod (p-1)}\pmod p.$$ If we define $$a_p = 2^p \bmod (p-1),$$ then the needed residue becomes $$b_p \equiv 2^{a_p}\pmod p.$$ This reduces a gigantic exponent tower to two ordinary modular exponentiations....

Detailed mathematical approach

Problem Summary

For each odd prime $p$, the computation associates an integer $g(p)$ obtained from a modular expression involving powers of $2$, and then sums these values:

$$G(n)=\sum_{\substack{3\le p \lt n\\ p\text{ prime}}} g(p).$$

The difficulty is that the exponent $2^p$ is already enormous, and the quantity for each prime is not taken merely modulo $p$ but depends on how a congruence lifts from modulo $p$ to modulo $p^2$. The implemented method therefore combines exponent reduction, modular inversion, a controlled lift to $p^2$, and a prime sieve.

Mathematical Approach

Fix an odd prime $p$. The implementation computes $g(p)$ by first finding a residue modulo $p$, then lifting it to a compatible residue modulo $p^2$, and finally extracting the multiple of $p$ that measures the difference between the two levels.

Step 1: Compress the Tower Exponent

We begin with the residue

$$b_p \equiv 2^{2^p}\pmod{p}.$$

Computing $2^{2^p}$ directly is hopeless, but Fermat's little theorem gives $2^{p-1}\equiv 1 \pmod p$, so exponents can be reduced modulo $p-1$. Therefore

$$2^{2^p}\equiv 2^{\,2^p\bmod (p-1)}\pmod p.$$

If we define

$$a_p = 2^p \bmod (p-1),$$

then the needed residue becomes

$$b_p \equiv 2^{a_p}\pmod p.$$

This reduces a gigantic exponent tower to two ordinary modular exponentiations.

Step 2: Solve the Halving Congruence Modulo $p$

Because $p$ is odd, the number $2$ is invertible modulo $p$. Its inverse is

$$2^{-1}\equiv \frac{p+1}{2}\pmod p,$$

since $2\cdot \frac{p+1}{2}=p+1\equiv 1\pmod p$.

Now choose the unique residue $h_p\in\{0,1,\dots,p-1\}$ satisfying

$$2h_p\equiv b_p\pmod p.$$

Equivalently,

$$h_p \equiv b_p\cdot \frac{p+1}{2}\pmod p.$$

This is the modular analogue of dividing the residue $b_p$ by $2$.

Step 3: Lift from Modulo $p$ to Modulo $p^2$

Next form the residue

$$w_p \equiv h_p\,2^p \pmod{p^2},\qquad 0\le w_p \lt p^2.$$

The key observation is that $w_p$ has the same value as $b_p$ modulo $p$. Indeed, by Fermat's little theorem,

$$2^p \equiv 2 \pmod p,$$

$$w_p \equiv h_p\,2 \equiv b_p \pmod p.$$

Therefore the difference $w_p-b_p$ is divisible by $p$. If $w_p\ge b_p$, we can divide that difference directly by $p$. If $w_p<b_p$, we add one full modulus $p^2$ first, which preserves the congruence class modulo $p^2$ while making the representative nonnegative.

Thus $g(p)$ is defined by

$$g(p)=\frac{\delta_p}{p},$$

where $\delta_p$ is the unique integer with

$$\delta_p \equiv w_p-b_p \pmod{p^2},\qquad 0\le \delta_p \lt p^2.$$

Because $w_p\equiv b_p\pmod p$, the numerator is automatically a multiple of $p$.

Step 4: Interpret the Quotient

The residue $b_p$ captures the value of the large power modulo $p$. The residue $w_p$ is a compatible lift modulo $p^2$. Since the two agree modulo $p$, their difference must be of the form $p\cdot q$. The number $g(p)$ is exactly this quotient:

$$w_p=b_p+p\,g(p)\quad\text{inside } \mathbb{Z}/p^2\mathbb{Z}.$$

So $g(p)$ measures the first correction term that appears when passing from arithmetic modulo $p$ to arithmetic modulo $p^2$.

Step 5: Sum Over Odd Primes

Once $g(p)$ is available for one prime, the target quantity is just

$$G(n)=\sum_{\substack{3\le p \lt n\\ p\text{ prime}}} g(p).$$

The remaining algorithmic task is therefore to enumerate every odd prime below $n$ and apply the same modular procedure to each one.

Worked Example: $p=31$

This prime is large enough to show every ingredient while still fitting on paper.

First reduce the tower exponent:

$$a_{31}=2^{31}\bmod 30=8.$$

Hence

$$b_{31}\equiv 2^8\equiv 256\equiv 8 \pmod{31}.$$

Now solve the halving congruence:

$$2h_{31}\equiv 8\pmod{31},$$

so the least nonnegative solution is

$$h_{31}=4.$$

Next compute the lifted power modulo $31^2=961$:

$$2^{31}\equiv 374\pmod{961}.$$

Therefore

$$w_{31}\equiv 4\cdot 374=1496\equiv 535\pmod{961}.$$

Now subtract the modulo-$31$ residue:

$$535-8=527=31\cdot 17.$$

$$g(31)=17.$$

This is one of the checkpoints used by the implementations. Another small aggregate checkpoint is

$$G(100)=474.$$

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They first generate a primality table up to $n$ with the classical sieve of Eratosthenes, then iterate only over odd primes. For each such prime, they perform three modular exponentiations: one modulo $p-1$ to compress the exponent, one modulo $p$ to recover the residue $2^{2^p}\bmod p$, and one modulo $p^2$ to build the lifted quantity.

After that, the implementation multiplies by the modular inverse of $2$ modulo $p$ in order to solve $2x\equiv b_p\pmod p$. The lifted residue is compared with the modulo-$p$ residue, and if the subtraction would be negative, one copy of $p^2$ is added before dividing by $p$. That guarantees an integer quotient in the canonical range. The resulting $g(p)$ is then accumulated into the final sum.

The C++ and Java versions include dedicated modular multiplication routines to keep intermediate products safe under 64-bit arithmetic, while the Python version relies on built-in arbitrary-precision integers and modular exponentiation.

Complexity Analysis

The sieve up to $n$ costs $O(n\log\log n)$ time and $O(n)$ memory. For each prime $p<n$, the algorithm performs a constant number of fast modular exponentiations, each taking $O(\log p)$ modular multiplications. Summed over all primes, the arithmetic work is

$$O\!\left(\sum_{p \lt n,\ p\text{ prime}}\log p\right)=O(\pi(n)\log n),$$

which is smaller than the sieve term for large $n$. Thus the full method runs in $O(n\log\log n+\pi(n)\log n)$ time and uses $O(n)$ memory.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=717
Fermat's little theorem: Wikipedia - Fermat's little theorem
Modular arithmetic: Wikipedia - Modular arithmetic
Modular exponentiation: Wikipedia - Modular exponentiation
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
Fermat quotient: Wikipedia - Fermat quotient

Problem 717 source code

C++

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

namespace {

using i128 = __int128_t;
using u64 = std::uint64_t;

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>(static_cast<i128>(a) * static_cast<i128>(b) % static_cast<i128>(mod));
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    if (mod == 1) {
        return 0;
    }
    u64 result = 1 % mod;
    base %= mod;
    while (exp > 0) {
        if (exp & 1ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1ULL;
    }
    return result;
}

u64 g_value(const u64 p) {
    const u64 e = pow_mod(2, p, p - 1);
    const u64 r = pow_mod(2, e, p);

    const u64 t = (r * ((p + 1) / 2)) % p;

    const u64 mod = p * p;
    const u64 two_p_mod = pow_mod(2, p, mod);
    u64 u = mul_mod(t, two_p_mod, mod);
    if (u < r) {
        u += mod;
    }
    u -= r;

    assert(u % p == 0);
    return u / p;
}

u64 G(const int n) {
    std::vector<bool> is_prime(static_cast<std::size_t>(n), true);
    if (n > 0) {
        is_prime[0] = false;
    }
    if (n > 1) {
        is_prime[1] = false;
    }

    for (int i = 2; static_cast<long long>(i) * i < n; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j < n; j += i) {
            is_prime[static_cast<std::size_t>(j)] = false;
        }
    }

    u64 sum = 0;
    for (int p = 3; p < n; p += 2) {
        if (is_prime[static_cast<std::size_t>(p)]) {
            sum += g_value(static_cast<u64>(p));
        }
    }
    return sum;
}

}  // namespace

int main() {
    assert(g_value(31) == 17);
    assert(G(100) == 474);
    assert(G(10'000) == 2'819'236);

    std::cout << G(10'000'000) << '\n';
    return 0;
}

Python

def g_value(p):
    e = pow(2, p, p - 1)
    r = pow(2, e, p)
    t = (r * ((p + 1) // 2)) % p
    
    mod = p * p
    two_p_mod = pow(2, p, mod)
    u = (t * two_p_mod) % mod
    
    if u < r:
        u += mod
    u -= r
    return u // p

def solve():
    n = 10000000
    is_prime = bytearray([1] * n)
    if n > 0: is_prime[0] = 0
    if n > 1: is_prime[1] = 0
    
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            for j in range(i * i, n, i):
                is_prime[j] = 0
                
    total_sum = 0
    for p in range(3, n, 2):
        if is_prime[p]:
            total_sum += g_value(p)
            
    return str(total_sum)

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

Java

import java.math.BigInteger;

public class Euler717 {
    static long mulMod(long a, long b, long mod) {
        BigInteger bigA = BigInteger.valueOf(a);
        BigInteger bigB = BigInteger.valueOf(b);
        BigInteger bigMod = BigInteger.valueOf(mod);
        return bigA.multiply(bigB).remainder(bigMod).longValue();
    }

    static long powMod(long base, long exp, long mod) {
        if (mod == 1)
            return 0;
        long result = 1 % mod;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) != 0) {
                result = mulMod(result, base, mod);
            }
            base = mulMod(base, base, mod);
            exp >>= 1;
        }
        return result;
    }

    static long gValue(long p) {
        long e = powMod(2, p, p - 1);
        long r = powMod(2, e, p);
        long t = mulMod(r, (p + 1) / 2, p);

        long mod = p * p;
        long twoPMod = powMod(2, p, mod);
        long u = mulMod(t, twoPMod, mod);

        if (u < r) {
            u += mod;
        }
        u -= r;

        return u / p;
    }

    public static String solve() {
        int n = 10000000;
        byte[] isPrime = new byte[n];
        java.util.Arrays.fill(isPrime, (byte) 1);
        if (n > 0)
            isPrime[0] = 0;
        if (n > 1)
            isPrime[1] = 0;

        for (int i = 2; (long) i * i < n; ++i) {
            if (isPrime[i] != 0) {
                for (int j = i * i; j < n; j += i) {
                    isPrime[j] = 0;
                }
            }
        }

        long totalSum = 0;
        for (int p = 3; p < n; p += 2) {
            if (isPrime[p] != 0) {
                totalSum += gValue(p);
            }
        }

        return Long.toString(totalSum);
    }

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