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

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

Equivalently, after swapping in the other direction,

$$ \nu_p(N!)= \sum_{j\ge 0} p^j \sum_{n=j+1}^{q} T_n. $$

This is the exact formula implemented by the code.

4) Why only the first $e$ layers matter modulo $p^e$. We only need the answer modulo $p^e$. If $j\ge e$, then $p^j$ is already divisible by $p^e$, so those terms vanish. Hence

$$ \nu_p(N!)\equiv \sum_{j=0}^{e-1} p^j S_j \pmod{p^e}, \qquad S_j=\sum_{n=j+1}^{q} T_n. $$

So the entire huge problem reduces to computing only $e$ suffix sums of the digit sequence.

5) Why the code stores only a short prefix. Let

$$T_{\mathrm{all}}=\sum_{n=0}^{q}T_n.$$

Then each suffix can be written as

$$ S_j=T_{\mathrm{all}}-\sum_{n=0}^{j}T_n. $$

Therefore the solver only needs:

$$\text{(a) the total sum of all digits, and}\qquad \text{(b) prefix sums up to index }e-1.$$

This is why it stores

$$\texttt{first\_len}=\min(e,q+1)$$

prefix values, not the entire sequence.

6) Sequence generation. The digits are generated from

$$s_{n+1}=s_n^2 \bmod 50515093,\qquad s_0=290797,$$

and then

$$T_n=s_n\bmod p.$$

The algorithm streams the generator exactly once, accumulating the total digit sum and the first few prefix sums.

Worked Example

Take a toy base-$5$ number with digits

$$T_0=2,\qquad T_1=4,\qquad T_2=1.$$

Then

$$N=2+4\cdot 5+1\cdot 25=47.$$

Legendre gives

$$\nu_5(47!)=\left\lfloor\frac{47}{5}\right\rfloor+\left\lfloor\frac{47}{25}\right\rfloor=9+1=10.$$

Now use the suffix formula:

$$S_0=T_1+T_2=5,\qquad S_1=T_2=1.$$

$$\nu_5(47!)=5^0S_0+5^1S_1=5+5=10,$$

exactly as expected.

Checks and Complexity

The source file includes the checkpoint

$$\mathrm{NF}(3,10000)\bmod 3^{20}=624955285,$$

and also compares the fast method against a brute-force big-integer computation for a smaller case.

Time complexity is

$$O(q+e),$$

because generating the digits costs $O(q)$ and the final valuation sum costs $O(e)$. Memory usage is

$$O(e),$$

since only a short prefix array and a few accumulators are stored.

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());
    }
}