Problem 288: An Enormous Factorial

Project Euler Problem 288 Solution

EulerSolve provides an optimized solution for Project Euler Problem 288, An Enormous Factorial, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A pseudo-random sequence produces digits $T_n\in\{0,1,\dots,p-1\}$. These digits define a huge base-$p$ integer $$N=\sum_{n=0}^{q}T_n p^n.$$ The required quantity is the $p$-adic valuation $$\nu_p(N!)$$ taken modulo $p^e$. For the actual problem, $p=61$, $q=10^7$, $e=10$, but the explanation below is written for general $p,q,e$. Mathematical Approach 1) Start from Legendre’s formula. For every positive integer $N$, $$\nu_p(N!)=\sum_{k\ge 1}\left\lfloor\frac{N}{p^k}\right\rfloor.$$ Since $N$ has only $q+1$ base-$p$ digits, the sum actually stops at $k=q$. 2) Expand each floor using the base-$p$ digits. Because $$N=\sum_{n=0}^{q}T_n p^n,$$ dividing by $p^k$ and taking the floor simply drops the lowest $k$ digits. Therefore $$\left\lfloor\frac{N}{p^k}\right\rfloor=\sum_{n=k}^{q}T_n p^{\,n-k}.$$ This already shows why $T_0$ never contributes to $\nu_p(N!)$: every term in Legendre’s sum starts with division by at least one power of $p$. 3) Swap the order of summation. Substitute the digit expansion into Legendre: $$ \nu_p(N!)= \sum_{k=1}^{q}\sum_{n=k}^{q}T_n p^{\,n-k}. $$ Now fix a digit position $n$. It contributes once for every $k=1,2,\dots,n$. Writing $j=n-k$, we obtain $$ \nu_p(N!)= \sum_{n=1}^{q} T_n \sum_{j=0}^{n-1} p^j....

Detailed mathematical approach

Problem Summary

For a prime $p$, a pseudo-random generator produces the numbers $S_0=290797$ and $S_{n+1}=S_n^2\bmod 50515093$, and from them the numbers $T_n=S_n\bmod p$. For an integer $q\ge0$ these define

$$N(p,q)=\sum_{n=0}^{q}T_n\,p^n.$$

Let $\operatorname{Nfac}(p,q)$ be the factorial $N(p,q)!$, and let $\operatorname{NF}(p,q)$ be the number of factors $p$ in it, that is, the exponent of the highest power of $p$ that divides $N(p,q)!$. The statement gives $\operatorname{NF}(3,10000)\bmod 3^{20}$${}=624955285$ and asks for $\operatorname{NF}(61,10^7)\bmod 61^{10}$, the remainder of $\operatorname{NF}(61,10^7)$ on division by $61^{10}$. ^[1] Already $N(61,10^7)$ has about $1.79\cdot10^7$ decimal digits, so its factorial can never be written down; its factors $p$ must be counted without it.

The whole solution in one sentence: the numbers $T_n$ are the digits of $N(p,q)$ in base $p$, Legendre's formula turns the number of factors $p$ in its factorial into a simple expression in these digits, and modulo $p^{10}$ this expression needs only the sum of all $10^7+1$ digits and the first ten of them, which one pass over the generator provides.

Mathematical Approach

We first count the factors of a prime in an ordinary factorial, then rewrite this count in terms of base-$p$ digits, and finally keep only what matters modulo $p^e$. Throughout, $e$ is the exponent of the modulus: $e=20$ in the checkpoint and $e=10$ in the task.

1. Factors of a prime in a factorial

Take $10!=3628800=2^8\cdot3^4\cdot5^2\cdot7$. Its factors $3$ come from the multiples of $3$ among $1,2,\dots,10$: each of $3$, $6$ and $9$ gives one factor, and $9=3^2$ gives a second one, so $10!$ has $3+1=4$ factors $3$. In general, for a prime $p$ and a positive integer $m$, let $\nu_p(m)$ be the exponent of the highest power of $p$ that divides $m$, the $p$-adic valuation of $m$. Since $p$ is prime, the valuation of a product is the sum of the valuations of its factors, ^[2] so

$$\nu_p(n!)=\sum_{m=1}^{n}\nu_p(m)=\sum_{k\ge1}\left\lfloor\frac{n}{p^k}\right\rfloor.$$

The second equality is Legendre's formula. ^[3] To prove it, count the pairs $(m,k)$ with $1\le m\le n$, $k\ge1$ and $p^k\mid m$ in two ways. For a fixed $m$ there are $\nu_p(m)$ such $k$, namely $k=1,\dots,\nu_p(m)$; for a fixed $k$ there are $\lfloor n/p^k\rfloor$ such $m$, the multiples $p^k,2p^k,\dots$ up to $n$. The sum is finite, because its terms are $0$ once $p^k\gt n$. For $10!$ it reads $\lfloor10/3\rfloor+\lfloor10/9\rfloor=3+1=4$, in agreement with the OEIS table of the factors $3$ in $n!$. ^[4] The primality of $p$ matters: $4^4=2^8$ divides $10!$, but $\lfloor10/4\rfloor+\lfloor10/16\rfloor=2$, because a factor $4$ can also be assembled from two factors $2$ of different numbers. In this notation $\operatorname{NF}(p,q)=\nu_p(N!)$ with $N=N(p,q)$.

2. The generator writes N in base p

In base $p$, every nonnegative integer is written with digits from $0$ to $p-1$, the digit at position $n$ counting multiples of $p^n$, just as decimal digits count units, tens and hundreds. Since $T_n=S_n\bmod p$ is a remainder on division by $p$, it always lies between $0$ and $p-1$. So the sum that defines $N(p,q)$ is exactly the base-$p$ expansion of $N(p,q)$, with the digit $T_n$ at position $n$; a top digit $T_q=0$ changes nothing.

Here is the example used throughout. The generator starts with $S_0=290797$, $S_1=629527$, $S_2=13339144$, $S_3=15552512$ and $S_4=17939732$, whose remainders on division by $3$ are $T_0,\dots,T_4=1,1,1,2,2$. Hence $N(3,4)=1+3+9+54+162$${}=229$, which is $22111_3$ in base $3$, with the top digit written first.

The recurrence for $S_n$ is the squaring map of the Blum Blum Shub generator, ^[5] and its modulus $50515093=5807\cdot8699$ is, as there, a product of two primes that both leave the remainder $3$ on division by $4$. Nothing about randomness is needed below, only $0\le T_n\le p-1$. From $S_1$ on, the sequence repeats with period $6308948$; the programs do not need this.

3. Dividing by a power of p deletes digits

Write $N=N(p,q)$ and split its expansion at a position $k$ with $1\le k\le q$:

$$N=p^k\sum_{n=k}^{q}T_n\,p^{n-k}+\sum_{n=0}^{k-1}T_n\,p^n.$$

Every digit is at most $p-1$, so the second sum is at most $(p-1)(1+p+\dots+p^{k-1})$${}=p^k-1$. Dividing by $p^k$ therefore gives the integer $\sum_{n=k}^{q}T_n\,p^{n-k}$ plus a fraction that is at least $0$ and less than $1$, and rounding down removes the fraction: ^[6]

$$\left\lfloor\frac{N}{p^k}\right\rfloor=\sum_{n=k}^{q}T_n\,p^{n-k}.$$

Dividing by $p^k$ and rounding down simply deletes the last $k$ digits, just as dividing by $1000$ deletes the last three decimal digits. The bound $T_n\le p-1$ is essential here; larger "digits" could carry into the quotient. For $N(3,4)=22111_3$ the quotients for $k=1,2,3,4$ are $76=2211_3$, $25=221_3$, $8=22_3$ and $2=2_3$, and Legendre's formula gives

$$\operatorname{NF}(3,4)=76+25+8+2=111.$$

4. What each digit contributes

Since $N\lt p^{q+1}$, the quotients vanish for $k\gt q$, and inserting those of Section 3 into Legendre's formula gives a double sum:

$$\operatorname{NF}(p,q)=\sum_{k=1}^{q}\ \sum_{n=k}^{q}T_n\,p^{n-k}.$$

Fix a position $n\ge1$. Its digit $T_n$ appears once for each $k=1,\dots,n$, with the weights $p^{n-1},\dots,p,1$, so it contributes $T_n$ times the geometric sum

$$1+p+\dots+p^{n-1}=\frac{p^n-1}{p-1},$$

while the lowest digit $T_0$ contributes nothing, since every quotient deletes it. In the example, $(3^n-1)/2$ equals $1$, $4$, $13$ and $40$ for $n=1,2,3,4$, so the five digits contribute $0$, $1$, $4$, $26$ and $80$, again $111$ in total. Adding all contributions and using $\sum_{n=0}^{q}T_n\,p^n=N$ gives the second form of Legendre's formula: ^[3]

$$\begin{aligned}\operatorname{NF}(p,q)&=\frac{N-s}{p-1},\\ s&=T_0+T_1+\dots+T_q,\end{aligned}$$

where $s$ is the sum of the base-$p$ digits of $N$. For the example $s=7$ and $(229-7)/2=111$. So $\operatorname{NF}(p,q)$ is a little less than $N/(p-1)$; already $\operatorname{NF}(3,10000)$ has $4771$ decimal digits, which is why only a remainder is asked for.

5. Grouping by powers of p

The double sum of Section 4 can also be grouped by the power $p^j$ with $j=n-k$. For a fixed $j\ge0$, the pairs with $1\le k\le n\le q$ and $n-k=j$ are given by the positions $n=j+1,\dots,q$, each once, so

$$\begin{aligned}\operatorname{NF}(p,q)&=\sum_{j=0}^{q-1}p^j\,R_j,\\ R_j&=\sum_{n=j+1}^{q}T_n.\end{aligned}$$

Thus $R_j$ is the sum of all digits above position $j$; the programs call it the suffix sum. For the digits $1,1,1,2,2$ of the example these sums are $R_0=6$, $R_1=5$, $R_2=4$ and $R_3=2$, and $6+3\cdot5+9\cdot4+27\cdot2=111$ once more. The numbers $R_j$ are not digits: $R_0=6$ exceeds $2$, so $R_3,R_2,R_1,R_0$ do not form the base-$3$ expansion of $111$, which is $11010_3$. The identity is exact all the same; it just leaves the carries undone.

6. Only e terms survive modulo p^e

For $j\ge e$ the power $p^j$ is a multiple of $p^e$, so all those terms leave the remainder $0$:

$$\operatorname{NF}(p,q)\equiv\sum_{j=0}^{e-1}p^j\,R_j\pmod{p^e}.$$

Here $x\equiv y\pmod m$ means that $x$ and $y$ leave the same remainder on division by $m$. However many digits $N$ has, only $e$ suffix sums are needed, and for $j\le q$ each of them is the total digit sum minus a short prefix:

$$R_j=s-(T_0+T_1+\dots+T_j).$$

For $j\ge q$ the sum $R_j$ is empty and equal to $0$; generator values after $T_q$ are not digits of $N$ and must not enter. So it suffices to know $s$ and the first $\min(e,q+1)$ prefix sums. For $p=3$ and $e=2$, the example gives $R_0+3R_1=6+15=21$, which leaves the remainder $3$ on division by $9$, and indeed $111=12\cdot9+3$.

7. A second route through the digit sum

The form of Section 4 needs the same data. The numbers $p-1$ and $p^e$ have no common factor, so $p-1$ has an inverse modulo $p^e$, a number $u$ with $(p-1)\,u\equiv1\pmod{p^e}$. ^[7] Multiplying $(p-1)\operatorname{NF}(p,q)=N-s$ by $u$ gives

$$\operatorname{NF}(p,q)\equiv(N-s)\,u\pmod{p^e},$$

and modulo $p^e$ the number $N$ reduces to its digits below position $e$, that is, to $\sum_{n=0}^{r}T_n\,p^n$ with $r=\min(q,e-1)$. In the example with $e=2$: $229\equiv1+1\cdot3=4\pmod9$, the inverse of $2$ modulo $9$ is $5$, since $2\cdot5=10$, and $(4-7)\cdot5=-15$ leaves the remainder $3$ on division by $9$, as before. Both routes need exactly the digit sum and the first $\min(e,q+1)$ digits. The programs follow Section 6, which needs no inverse; the second route served as an independent check of them.

8. The size of the real task

For $p=61$ and $q=10^7$, the number $N$ has $10^7+1$ digits in base $61$, about $5.9\cdot10^7$ bits. Storing it would take about 7.4 MB, but Legendre's formula evaluated directly would divide it by $61$ ten million times, roughly $5\cdot10^{12}$ operations on 64-bit words. Section 6 needs instead one pass that adds up $10^7+1$ digits of at most $60$ each, so every $R_j$ is at most $60\,(10^7+1)$, about $6\cdot10^8$, and fits easily into 32 bits. The modulus $61^{10}=713342911662882601$, about $7.1\cdot10^{17}$, fits into a signed 64-bit integer, whose largest value is $2^{63}-1\approx9.2\cdot10^{18}$. A product $p^jR_j$, however, does not: in the task the term with $j=9$ exceeds $10^{24}$, far above $2^{64}\approx1.8\cdot10^{19}$.

How the Code Works

The C++ program ^[8] keeps its parameters in an Options structure with the defaults $p=61$, $q=10^7$ and exponent $e=10$. parse_arguments, with the helper parse_int_after_prefix, accepts --p=, --q=, --exponent= and --skip-checkpoints, and ends the program with exit code $1$ on anything else. It requires $p\ge2$, $q\ge0$ and $1\le e\le24$, but it does not check that $p$ is prime, and it has no overflow checks: an oversized numeral overflows the int parser, $q+1$ must fit into an int, and $p^e$ into 64 bits. pow_u64 computes the modulus $p^e$ by repeated multiplication in unsigned 64-bit arithmetic, which is exact for the defaults, since $61^{10}\lt2^{64}$; for $p=61$ and --exponent=11 it would silently wrap around, and the results can be wrong.

solve(p, q, exponent) carries out Section 6 in two loops. It sets first_len to $\min(e,q+1)$, the number of prefix sums it needs, and creates the vector prefix with prefix[0] $=0$. The first loop runs over $n=0,\dots,q$: it takes the digit t = s % p, adds it to total, which ends as the digit sum $s$ of Section 4, stores prefix[n+1] $=T_0+\dots+T_n$ while $n$ is below first_len, and advances the generator by squaring s modulo $50515093$, in 128-bit arithmetic, although $S_n^2\lt2.6\cdot10^{15}$ fits into 64 bits. The digits themselves are never stored. The second loop runs over $j=0,\dots,e-1$: it takes pref = prefix[j+1], or total once $j$ passes the last digit, which happens only when $q+1\lt e$, so that suffix = total - pref is $R_j$. It then adds $p^jR_j$ to answer modulo $p^e$, keeping the running power p_pow $=p^j\bmod p^e$. Products and sums are formed in unsigned __int128, a GCC and Clang extension, ^[9] because by Section 8 they can exceed 64 bits. A third branch of this loop, which would use prefix[first_len], can never be reached.

brute_small builds $N(p,q)$ exactly as a Boost.Multiprecision cpp_int, an arbitrary-precision integer, ^[10] and evaluates Legendre's formula literally: it divides a copy of $N$ by $p$ again and again and adds up the quotients. These are the numbers $\lfloor N/p^k\rfloor$, because $\lfloor\lfloor x/p^k\rfloor/p\rfloor=\lfloor x/p^{k+1}\rfloor$. ^[6] run_checkpoints requires solve(3, 10000, 20) to give the value $624955285$ of the statement and solve(5, 200, 8) to agree with brute_small(5, 200, 8); if a check fails, the program reports it and exits with code $2$. main runs these checks first, unless --skip-checkpoints is given, and then prints the result of solve.

The Python solve is the same function with the default arguments $p=61$, $q=10^7$ and exponent $10$. Its integers have unlimited precision, ^[11] so nothing can overflow, and it returns the result as a string. The Java version hard-codes the parameters and uses long for the generator, the digit sum and the prefix sums, which is safe because $S_n^2\lt2^{63}$ and $s\le6.1\cdot10^8$. For the multiply-add modulo $p^e$ and for the update of the power it switches to BigInteger, ^[11] whose mod never returns a negative value. The same loop also computes a 64-bit product named term, which could overflow, but its value is never used. Neither port runs the checkpoints.

Complexity and Verification

The first loop makes $q+1$ steps of a few machine operations each, and the second loop makes $e$ steps with two multiplications and three reductions each, so the time is $O(q+e)$ arithmetic operations. For the task all intermediate values stay below $2^{128}$, so in all three languages every operation acts on numbers of bounded size. The memory is $O(e)$, at most $e+1$ prefix sums and a few counters, independent of $q$. For the task this means $10^7+1$ generator steps, and on a recent computer each of the three programs finishes in well under a second.

The C++ program checks itself on the value from the statement and against brute_small for $p=5$, $q=200$ and $e=8$, where $N(5,200)$ has $200$ base-$5$ digits, because $T_{200}=0$. The following checks were run for this article and are not part of the published solutions:

A brute force that builds $N(p,q)$ as an exact integer and evaluates Legendre's formula directly agrees with the formula of Section 6 and with the route of Section 7 for every prime $p\in\{2,3,5,7,11,13,61,97\}$, every $q\le59$ and $q\in\{100,257,500,1000\}$, and every $e\le12$, in all $6144$ cases.
The same brute force reproduces the checkpoint $624955285$; there $N(3,10000)$ has $4771$ decimal digits.
The compiled C++ program, called with command-line parameters, agrees with the brute force in $4144$ runs: the same eight primes, every $q\le24$ and $q\in\{50,100,300\}$, and every $e\le24$ with $p^e\lt2^{64}$, including all cases with $q+1\lt e$.
For the task, the three published programs and an independent program for the route of Section 7 print the same number.

Footnotes and References

The references below supply the background facts used above. The digit form of the quotients, the suffix sums $R_j$ and their reduction modulo $p^e$ are derived explicitly in this article.

Project Euler 288 — An Enormous Factorial. The official statement defines N(p, q) as the sum of Tₙ·pⁿ for n = 0, …, q, with S₀ = 290797, Sₙ₊₁ = Sₙ² mod 50515093 and Tₙ = Sₙ mod p, lets NF(p, q) be the number of factors p in Nfac(p, q), the factorial of N(p, q), gives NF(3, 10000) mod 3²⁰ = 624955285 and asks for NF(61, 10⁷) mod 61¹⁰.
Wikipedia — p-adic valuation. For a prime p, the p-adic valuation νₚ(n) of a nonzero integer n is the exponent of the highest power of p that divides n, and νₚ(rs) = νₚ(r) + νₚ(s).
Wikipedia — Legendre's formula. For a prime p, νₚ(n!) is the sum of ⌊n/pⁱ⌋ over i ≥ 1, proved by counting one factor p for each multiple of p, an additional one for each multiple of p², and so on. The article also proves the alternative form νₚ(n!) = (n − sₚ(n))/(p − 1), where sₚ(n) is the sum of the base-p digits of n, by writing each quotient in terms of the digits, as in Sections 3 and 4.
OEIS A054861 — Greatest k such that 3ᵏ divides n!. The entry gives a(n) = ⌊n/3⌋ + ⌊n/9⌋ + ⌊n/27⌋ + … and a(n) = (n − A053735(n))/2, where A053735(n) is the sum of the base-3 digits of n; its terms include a(10) = 4.
Wikipedia — Blum Blum Shub. The generator computes xₙ₊₁ = xₙ² mod M, where M is the product of two large primes that are both congruent to 3 modulo 4, and outputs the parity bit or a few least significant bits of each value.
Wikipedia — Floor and ceiling functions. For a positive y, x mod y = x − y⌊x/y⌋ is the remainder of x on division by y and lies in [0, y); for a positive integer n, ⌊⌊x/m⌋/n⌋ = ⌊x/(mn)⌋. Section 3 uses the first fact, brute_small the second.
Wikipedia — Modular multiplicative inverse. An inverse of a modulo m exists if and only if gcd(a, m) = 1, it is unique modulo m, and the extended Euclidean algorithm finds it. Section 7 uses the inverse of p − 1 modulo pᵉ.
C++ — Euler288.cpp. The linked immutable C++ revision contains parse_arguments, solve, brute_small and run_checkpoints with the two checks described above. It is the implementation record for this article.
GCC manual — 128-bit Integers. GCC supports the types __int128 and unsigned __int128 on targets with an integer mode wide enough to hold 128 bits; solve uses unsigned __int128 for the products of Section 8.
Boost.Multiprecision — cpp_int. The type cpp_int, declared in boost/multiprecision/cpp_int.hpp, is an arbitrary-precision signed integer; brute_small keeps N(p, q) and its quotients in it.
Java documentation — BigInteger describes immutable arbitrary-precision integers whose mod always returns a value between 0 and the modulus minus 1. Python — numeric types states that integers have unlimited precision.

Problem 288 source code

C++

#include <boost/multiprecision/cpp_int.hpp>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using boost::multiprecision::cpp_int;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    int p = 61;
    int q = 10000000;
    int exponent = 10;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--p=", options.p) ||
            parse_int_after_prefix(arg, "--q=", options.q) ||
            parse_int_after_prefix(arg, "--exponent=", options.exponent)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.p >= 2 && options.q >= 0 && options.exponent >= 1 && options.exponent <= 24;
}

u64 pow_u64(const u64 base, const int exp) {
    u64 value = 1;
    for (int i = 0; i < exp; ++i) {
        value *= base;
    }
    return value;
}

u64 solve(const int p, const int q, const int exponent) {
    const u64 mod = pow_u64(static_cast<u64>(p), exponent);
    const int first_len = std::min(exponent, q + 1);
    std::vector<u64> prefix(static_cast<std::size_t>(first_len + 1), 0ULL);

    u64 s = 290797ULL;
    u64 total = 0ULL;
    for (int n = 0; n <= q; ++n) {
        const u64 t = s % static_cast<u64>(p);
        total += t;
        if (n < first_len) {
            prefix[static_cast<std::size_t>(n + 1)] = prefix[static_cast<std::size_t>(n)] + t;
        }
        s = static_cast<u64>((static_cast<u128>(s) * s) % 50515093ULL);
    }

    u64 answer = 0ULL;
    u64 p_pow = 1ULL;
    for (int j = 0; j < exponent; ++j) {
        u64 pref = 0ULL;
        if (j + 1 <= first_len) {
            pref = prefix[static_cast<std::size_t>(j + 1)];
        } else if (first_len == q + 1) {
            pref = total;
        } else {
            pref = prefix[static_cast<std::size_t>(first_len)];
        }
        const u64 suffix = total - pref;
        answer = static_cast<u64>((static_cast<u128>(answer) +
                                   static_cast<u128>(p_pow) * (suffix % mod)) %
                                  mod);
        p_pow = static_cast<u64>((static_cast<u128>(p_pow) * static_cast<u64>(p)) % mod);
    }
    return answer;
}

u64 brute_small(const int p, const int q, const int exponent) {
    const u64 mod = pow_u64(static_cast<u64>(p), exponent);
    u64 s = 290797ULL;
    cpp_int n_value = 0;
    cpp_int p_pow = 1;
    for (int i = 0; i <= q; ++i) {
        const int t = static_cast<int>(s % static_cast<u64>(p));
        n_value += cpp_int(t) * p_pow;
        p_pow *= p;
        s = static_cast<u64>((static_cast<u128>(s) * s) % 50515093ULL);
    }

    cpp_int valuation = 0;
    cpp_int current = n_value;
    while (current > 0) {
        current /= p;
        valuation += current;
    }
    return static_cast<u64>(valuation % mod);
}

bool run_checkpoints() {
    if (solve(3, 10000, 20) != 624955285ULL) {
        std::cerr << "Checkpoint failed for stated sample NF(3,10000) mod 3^20" << '\n';
        return false;
    }
    if (solve(5, 200, 8) != brute_small(5, 200, 8)) {
        std::cerr << "Checkpoint failed for brute cross-check at p=5 q=200" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.p, options.q, options.exponent) << '\n';
    return 0;
}

Python

def solve(p=61, q=10000000, exponent=10):
    mod = p ** exponent
    first_len = min(exponent, q + 1)
    prefix = [0] * (first_len + 1)

    s = 290797
    total = 0
    for n in range(q + 1):
        t = s % p
        total += t
        if n < first_len:
            prefix[n + 1] = prefix[n] + t
        s = (s * s) % 50515093

    answer = 0
    p_pow = 1
    for j in range(exponent):
        pref = 0
        if j + 1 <= first_len:
            pref = prefix[j + 1]
        elif first_len == q + 1:
            pref = total
        else:
            pref = prefix[first_len]
            
        suffix = total - pref
        answer = (answer + p_pow * (suffix % mod)) % mod
        p_pow = (p_pow * p) % mod

    return str(answer)

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

Java

public class Euler288 {
    static long powU64(long base, int exp) {
        long value = 1;
        for (int i = 0; i < exp; ++i) {
            value *= base;
        }
        return value;
    }

    public static String solve() {
        int p = 61;
        int q = 10000000;
        int exponent = 10;

        long mod = powU64(p, exponent);
        int firstLen = Math.min(exponent, q + 1);
        long[] prefix = new long[firstLen + 1];

        long s = 290797L;
        long total = 0L;
        for (int n = 0; n <= q; ++n) {
            long t = s % p;
            total += t;
            if (n < firstLen) {
                prefix[n + 1] = prefix[n] + t;
            }
            s = (s * s) % 50515093L;
        }

        long answer = 0L;
        long pPow = 1L;
        for (int j = 0; j < exponent; ++j) {
            long pref = 0L;
            if (j + 1 <= firstLen) {
                pref = prefix[j + 1];
            } else if (firstLen == q + 1) {
                pref = total;
            } else {
                pref = prefix[firstLen];
            }

            long suffix = total - pref;

            // answer = (answer + pPow * (suffix % mod)) % mod
            long term = (pPow % mod) * (suffix % mod); // pPow < mod, suffix%mod < mod. Could overflow long if mod is
                                                       // large.
            // But mod is 61^10 = 7.19 * 10^17.
            // We need u128 multiplication equivalent or BigInteger.
            // Using BigInteger to be safe.
            java.math.BigInteger bigTerm = java.math.BigInteger.valueOf(pPow)
                    .multiply(java.math.BigInteger.valueOf(suffix % mod));
            java.math.BigInteger bigAnswer = java.math.BigInteger.valueOf(answer)
                    .add(bigTerm).mod(java.math.BigInteger.valueOf(mod));
            answer = bigAnswer.longValue();

            pPow = java.math.BigInteger.valueOf(pPow)
                    .multiply(java.math.BigInteger.valueOf(p))
                    .mod(java.math.BigInteger.valueOf(mod)).longValue();
        }

        return String.valueOf(answer);
    }

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