Problem 995: A Particular Pair of Polynomials

Project Euler Problem 995 Solution

EulerSolve provides an optimized solution for Project Euler Problem 995, A Particular Pair of Polynomials, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 995 looks like a question about divisibility of two sparse polynomials, but the useful object is not a polynomial expansion. For a prime $p$, $f_p(x)=1+x+\cdots+x^{p-1}$ is the cyclotomic polynomial $\Phi_p(x)$. Therefore divisibility by $f_p$ can be tested at a primitive $p$-th root of unity. This turns the problem into a question about how the divisor set of an integer $s$ is distributed among residue classes modulo $p$. The main reduction is that $f_p(x)\mid g_s(x)$ holds exactly when the positive divisors of $s$ form a perfect tiling of $\mathbb{F}_p^\times$ by residues. Consequently $p\nmid s$, $\tau(s)=p-1$, and each non-zero residue class modulo $p$ must be hit exactly once by a divisor of $s$. The implementation then searches for the numerically smallest such $s$, not by enumerating divisors of candidate integers, but by building a subgroup chain inside the cyclic group $\mathbb{F}_p^\times$. The final product $T(20000)$ is much too large to store. The program computes each $\log_{10} S(p)$, sums those logarithms over all primes $p\lt 20000$, and formats the final mantissa and exponent. Exact integer reconstruction is used only for small validation cases such as $S(5)=8$, $T(20)=1348422598656$, and the stated $T(100)$ checkpoint....

Detailed mathematical approach

Problem Summary

Problem 995 looks like a question about divisibility of two sparse polynomials, but the useful object is not a polynomial expansion. For a prime $p$, $f_p(x)=1+x+\cdots+x^{p-1}$ is the cyclotomic polynomial $\Phi_p(x)$. Therefore divisibility by $f_p$ can be tested at a primitive $p$-th root of unity. This turns the problem into a question about how the divisor set of an integer $s$ is distributed among residue classes modulo $p$.

The main reduction is that $f_p(x)\mid g_s(x)$ holds exactly when the positive divisors of $s$ form a perfect tiling of $\mathbb{F}_p^\times$ by residues. Consequently $p\nmid s$, $\tau(s)=p-1$, and each non-zero residue class modulo $p$ must be hit exactly once by a divisor of $s$. The implementation then searches for the numerically smallest such $s$, not by enumerating divisors of candidate integers, but by building a subgroup chain inside the cyclic group $\mathbb{F}_p^\times$.

The final product $T(20000)$ is much too large to store. The program computes each $\log_{10} S(p)$, sums those logarithms over all primes $p\lt 20000$, and formats the final mantissa and exponent. Exact integer reconstruction is used only for small validation cases such as $S(5)=8$, $T(20)=1348422598656$, and the stated $T(100)$ checkpoint.

Mathematical Approach and Notation

For a prime $p$ and a positive integer $n$, the problem defines

$$f_p(x)=\sum_{i=0}^{p-1}x^i,\qquad g_n(x)=1+\sum_{d\mid n}x^d.$$

The term $1$ in $g_n$ is important: it is an extra constant term whose exponent is $0$. The divisors $d\mid n$ are positive divisors and contribute exponents $d$. We seek the least positive integer $s$ such that $f_p(x)\mid g_s(x)$, denoted $S(p)$.

Throughout the analysis, $G=\mathbb{F}_p^\times$ denotes the multiplicative group of non-zero residues modulo $p$. This group is cyclic of order $p-1$. If $g$ is a primitive root modulo $p$, every non-zero residue has a unique representation $g^e$, $0\le e\lt p-1$; the exponent $e$ is the discrete logarithm used by the code.

Lemma 1: The Cyclotomic Root Condition

Let $\zeta\ne 1$ be a primitive $p$-th root of unity. Since $p$ is prime,

$$f_p(x)=1+x+\cdots+x^{p-1}=\Phi_p(x),$$

and $\Phi_p(x)$ is the minimal polynomial of $\zeta$ over $\mathbb{Q}$. Therefore $f_p(x)\mid g_s(x)$ if and only if $g_s(\zeta)=0$.

Group the exponents of $g_s$ by their residue classes modulo $p$. Let $c_r$ be the number of terms in $g_s$ whose exponent is congruent to $r\pmod p$. Then

$$g_s(\zeta)=\sum_{r=0}^{p-1}c_r\zeta^r.$$

The powers $1,\zeta,\ldots,\zeta^{p-1}$ satisfy one rational linear relation, namely $1+\zeta+\cdots+\zeta^{p-1}=0$, and all rational relations are multiples of it. Hence $g_s(\zeta)=0$ is equivalent to

$$c_0=c_1=\cdots=c_{p-1}.$$

Thus the polynomial divisibility problem has become a balancing condition on residue-class counts. The coefficients of the polynomials themselves never need to be expanded.

Lemma 2: The Divisors of $s$ Tile $\mathbb{F}_p^\times$

First suppose $p\mid s$, and write $s=p^a u$, with $a\ge 1$ and $p\nmid u$. The divisors not divisible by $p$ are exactly the divisors of $u$, so all non-zero residue classes together receive $\tau(u)$ contributions. If the common residue count is $C$, then

$$\tau(u)=(p-1)C.$$

The zero residue class receives the extra constant term and every divisor containing at least one factor $p$, so its count is $1+a\tau(u)$. Equality of all $c_r$ would require

$$C=1+a\tau(u)=1+a(p-1)C,$$

which is impossible for positive $C$. Hence no valid minimal or non-minimal solution can have $p\mid s$.

Now $p\nmid s$. Every divisor of $s$ is non-zero modulo $p$, and the only contribution to residue class $0$ is the separate constant term in $g_s$. Therefore $c_0=1$. Since all $c_r$ must be equal, every non-zero class must also occur exactly once. The divisibility condition is therefore equivalent to the bijection

$$\{d:d\mid s\}\longrightarrow \mathbb{F}_p^\times,\qquad d\longmapsto d\bmod p.$$

In particular, the number of divisors must be

$$\tau(s)=p-1.$$

Subgroup-Chain Model

Write the unknown integer in prime-power form

$$s=\prod_{j=1}^k q_j^{a_j},\qquad q_j\ne p,$$

and define the divisor-box lengths $\ell_j=a_j+1$. Since $\tau(s)=\prod_j(a_j+1)$, the lengths must satisfy

$$\prod_{j=1}^k \ell_j=p-1.$$

Choosing a factor $q_j^{a_j}$ adds the finite progression $1,q_j,q_j^2,\ldots,q_j^{\ell_j-1}$ to the divisor box. Modulo $p$, these progressions should multiply together without collisions and eventually cover the whole cyclic group $G$. The implementation uses the canonical subgroup-chain formulation of this condition: after some steps the already constructed residues are the unique subgroup $H\le G$ of size $h\mid p-1$.

Assume the current subgroup has size $h$. The quotient group $G/H$ has order

$$Q={p-1\over h}.$$

A new prime residue $q$ may be used with length $\ell\mid Q$ precisely when the coset $qH$ has order $\ell$ in $G/H$. Then the cosets

$$H,\;qH,\;q^2H,\;\ldots,\;q^{\ell-1}H$$

are distinct, and multiplying the old subgroup by $1,q,\ldots,q^{\ell-1}$ gives the unique subgroup of size $h\ell$. This is exactly the no-collision condition needed by the divisor tiling.

Using a primitive root $g$, write $q\equiv g^e\pmod p$. In the quotient of order $Q$, the order of $qH$ is

$$\operatorname{ord}_{G/H}(qH)={Q\over \gcd(e,Q)}.$$

Therefore the transition condition used in the solver is

$$\operatorname{ord}_{G/H}(qH)=\ell \quad\Longleftrightarrow\quad \gcd(e,Q)={Q\over \ell}.$$

Optimization as Dynamic Programming

The value of $s$ should be minimal as an integer. If a transition uses a prime $q$ with length $\ell$, then it contributes the factor $q^{\ell-1}$ to $s$. The code minimizes logarithms, so the transition weight is

$$w(\ell,q)=(\ell-1)\log_{10}q.$$

For a fixed state $h$ and a fixed length $\ell$, the future state depends only on $h\ell$, not on the path by which $h$ was reached. Because $G$ is cyclic, there is a unique subgroup of each size $h\mid p-1$. Hence the best prime for that transition is simply the smallest prime $q\ne p$ satisfying the discrete-log condition above.

The dynamic program stores $D(h)$, the smallest known value of $\log_{10}$ of the partial integer that constructs the subgroup of size $h$. It starts with

$$D(1)=0,$$

and for every divisor $\ell\gt 1$ of $Q=(p-1)/h$, it relaxes

$$D(h\ell)=\min\left(D(h\ell),\;D(h)+(\ell-1)\log_{10}q\right),$$

where $q$ is the smallest transition prime with $\gcd(\log_g q,Q)=Q/\ell$. The answer for one prime is $D(p-1)=\log_{10}S(p)$. Parent pointers record the chosen pairs $(\ell,q)$, which lets the program reconstruct exact small values for validation.

Correctness Argument

The root-condition lemma proves that polynomial divisibility is equivalent to equal residue counts. The divisor-tiling lemma then proves that any solution must have $p\nmid s$ and that its divisors must hit the elements of $G$ exactly once. Thus solving the original problem is equivalent to finding the smallest integer whose divisor box is a collision-free product decomposition of $G$.

Each DP transition is valid because the quotient-order condition makes the new powers $1,q,\ldots,q^{\ell-1}$ represent $\ell$ distinct cosets of the current subgroup. The old subgroup has $h$ elements, so the enlarged product has $h\ell$ distinct residues and is exactly the subgroup of that size. By induction, every DP path from $1$ to $p-1$ constructs a valid divisor tiling and hence a valid integer $s$.

Conversely, in the cyclic group setting relevant here, an exact divisor-box tiling can be read as a sequence of quotient complements: at each stage one chooses a divisor-box direction whose images are distinct cosets over the subgroup already constructed. Since the group has only one subgroup of each divisor size, the state $h$ loses no subgroup identity information. Therefore the DP considers the necessary canonical transitions, and taking the minimum over them yields $S(p)$.

Worked Examples and Checkpoints

For $p=2$, $f_2(x)=1+x$, and $g_1(x)=1+x$, so $S(2)=1$. This is the degenerate group case $G=\mathbb{F}_2^\times$, whose order is $1$; the DP starts and finishes at the same state.

For $p=5$, $G$ has order $4$. The smallest primitive residue is $2$. A single transition with $\ell=4$ gives

$$S(5)=2^{4-1}=8,$$

matching the problem statement. Its divisors $1,2,4,8$ reduce modulo $5$ to $1,2,4,3$, all non-zero residues exactly once.

For $p=7$, the optimal chain can be expressed as $(\ell,q)=(3,2)$ followed by $(2,3)$. The resulting integer is

$$S(7)=2^{3-1}3^{2-1}=12.$$

The divisors $1,2,4,3,6,12$ reduce modulo $7$ to $1,2,4,3,6,5$, again a complete tiling of $\mathbb{F}_7^\times$. The full implementation checks these local examples and also verifies the product checkpoints $T(20)=1348422598656$ and $T(100)=1.37451\mathrm{e}123$.

Implementation Details

The C++, Python, and Java versions follow the same pipeline. They sieve primes up to $2\,000\,000$, factor $p-1$, enumerate all divisors of $p-1$, find a primitive root modulo $p$, and build a discrete-log table mapping each non-zero residue to its exponent with respect to that root.

The function that searches a transition prime scans the precomputed prime list and rejects $q=p$. If no suitable prime appears inside the sieve range, the implementation continues by trial primality testing beyond the sieve. In practice the bound is generous for the $p\lt 20000$ instances, but the fallback keeps the logic complete rather than dependent on a guessed cutoff.

All large products are handled logarithmically. Exact multiplication is retained only for reconstructed small $S(p)$ values and for $T(20)$. This separation is important: the final $T(20000)$ has hundreds of thousands of decimal digits, while its scientific notation only needs the fractional part and integer part of the summed base-10 logarithm.

Complexity and Numerical Stability

For a single prime $p$, the DP state count is $\tau(p-1)$, the number of divisors of $p-1$. Each state tests divisor lengths $\ell\mid (p-1)/h$, and each transition scans primes until it finds a residue with the required quotient order. The memory use for the discrete-log table is $O(p)$, and it is discarded after that prime is solved.

The full target has only $2262$ primes below $20000$. The algorithm therefore spends its time on small divisor lattices of $p-1$, not on polynomial coefficients and not on divisors of a gigantic candidate integer. Long-double logarithms are sufficient because the output requires a rounded mantissa with five digits after the decimal point; the exact small checkpoints guard against structural errors in the reduction.

References

Problem page: Project Euler 995
Cyclotomic polynomial: Wikipedia - Cyclotomic polynomial
Finite field: Wikipedia - Finite field
Primitive root modulo $n$: Wikipedia - Primitive root modulo n
Dynamic programming: Wikipedia - Dynamic programming

Problem 995 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <sstream>
#include <string>
#include <utility>
#include <vector>

namespace {

constexpr int kLimit = 20000;
constexpr int kPrimeSearchLimit = 2000000;

std::vector<int> sieve_primes(int n) {
    std::vector<bool> is_prime(n + 1, true);
    if (n >= 0) {
        is_prime[0] = false;
    }
    if (n >= 1) {
        is_prime[1] = false;
    }
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (!is_prime[i]) {
            continue;
        }
        for (long long j = 1LL * i * i; j <= n; j += i) {
            is_prime[static_cast<std::size_t>(j)] = false;
        }
    }

    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[i]) {
            primes.push_back(i);
        }
    }
    return primes;
}

std::vector<std::pair<int, int>> factorize(int n, const std::vector<int>& primes) {
    std::vector<std::pair<int, int>> factors;
    int x = n;
    for (int q : primes) {
        if (1LL * q * q > x) {
            break;
        }
        if (x % q != 0) {
            continue;
        }
        int e = 0;
        while (x % q == 0) {
            x /= q;
            ++e;
        }
        factors.push_back({q, e});
    }
    if (x > 1) {
        factors.push_back({x, 1});
    }
    return factors;
}

std::vector<int> divisors_from_factors(const std::vector<std::pair<int, int>>& factors) {
    std::vector<int> divisors{1};
    for (auto [prime, exponent] : factors) {
        const std::vector<int> previous = divisors;
        int power = 1;
        for (int e = 1; e <= exponent; ++e) {
            power *= prime;
            for (int d : previous) {
                divisors.push_back(d * power);
            }
        }
    }
    std::sort(divisors.begin(), divisors.end());
    return divisors;
}

int primitive_root(int p, const std::vector<int>& small_primes) {
    const auto factors = factorize(p - 1, small_primes);
    for (int g = 2; g < p; ++g) {
        bool ok = true;
        for (auto [q, exponent] : factors) {
            (void)exponent;
            long long value = 1;
            long long base = g;
            int power = (p - 1) / q;
            while (power > 0) {
                if (power & 1) {
                    value = (value * base) % p;
                }
                base = (base * base) % p;
                power >>= 1;
            }
            if (value == 1) {
                ok = false;
                break;
            }
        }
        if (ok) {
            return g;
        }
    }
    return -1;
}

bool is_prime_by_trial(long long n, const std::vector<int>& primes) {
    if (n < 2) {
        return false;
    }
    for (int q : primes) {
        if (1LL * q * q > n) {
            return true;
        }
        if (n % q == 0) {
            return n == q;
        }
    }
    return true;
}

int smallest_transition_prime(int p, int quotient_order, int target_gcd,
                              const std::vector<int>& discrete_log,
                              const std::vector<int>& primes) {
    auto works = [&](long long q) {
        if (q == p) {
            return false;
        }
        const int residue = static_cast<int>(q % p);
        if (residue == 0) {
            return false;
        }
        const int e = discrete_log[residue];
        return std::gcd(e, quotient_order) == target_gcd;
    };

    for (int q : primes) {
        if (works(q)) {
            return q;
        }
    }

    long long q = primes.back() + 1LL;
    if ((q & 1LL) == 0) {
        ++q;
    }
    for (;; q += 2) {
        if (is_prime_by_trial(q, primes) && works(q)) {
            return static_cast<int>(q);
        }
    }
}

struct PrimeSolution {
    long double log10_value = 0.0L;
    std::vector<std::pair<int, int>> factors;
};

PrimeSolution solve_prime(int p, const std::vector<int>& primes) {
    if (p == 2) {
        return {};
    }

    const int n = p - 1;
    const int g = primitive_root(p, primes);
    assert(g > 0);

    std::vector<int> discrete_log(p, -1);
    int x = 1;
    for (int e = 0; e < n; ++e) {
        discrete_log[x] = e;
        x = static_cast<int>(1LL * x * g % p);
    }

    const auto divisors = divisors_from_factors(factorize(n, primes));
    const int d_count = static_cast<int>(divisors.size());
    const long double inf = std::numeric_limits<long double>::infinity();

    std::vector<long double> dp(d_count, inf);
    std::vector<int> parent(d_count, -1);
    std::vector<std::pair<int, int>> edge(d_count, {0, 0});

    auto index_of = [&](int value) {
        return static_cast<int>(std::lower_bound(divisors.begin(), divisors.end(), value) - divisors.begin());
    };

    dp[index_of(1)] = 0.0L;
    for (int i = 0; i < d_count; ++i) {
        const int subgroup_size = divisors[i];
        if (!std::isfinite(dp[i])) {
            continue;
        }

        const int quotient_order = n / subgroup_size;
        for (int length : divisors) {
            if (length <= 1 || quotient_order % length != 0) {
                continue;
            }
            const int target_gcd = quotient_order / length;
            const int q = smallest_transition_prime(p, quotient_order, target_gcd, discrete_log, primes);
            const int next_size = subgroup_size * length;
            const int j = index_of(next_size);
            const long double candidate = dp[i] + static_cast<long double>(length - 1) * std::log10(static_cast<long double>(q));
            if (candidate + 1e-24L < dp[j]) {
                dp[j] = candidate;
                parent[j] = i;
                edge[j] = {length, q};
            }
        }
    }

    const int finish = index_of(n);
    PrimeSolution solution;
    solution.log10_value = dp[finish];

    for (int at = finish; parent[at] != -1; at = parent[at]) {
        solution.factors.push_back(edge[at]);
    }
    std::reverse(solution.factors.begin(), solution.factors.end());
    return solution;
}

std::uint64_t exact_value(const PrimeSolution& solution) {
    std::uint64_t value = 1;
    for (auto [length, q] : solution.factors) {
        for (int i = 1; i < length; ++i) {
            value *= static_cast<std::uint64_t>(q);
        }
    }
    return value;
}

std::string scientific(long double log10_value) {
    long long exponent = static_cast<long long>(std::floor(log10_value));
    long double mantissa = std::pow(10.0L, log10_value - static_cast<long double>(exponent));
    mantissa = std::round(mantissa * 100000.0L) / 100000.0L;
    if (mantissa >= 10.0L) {
        mantissa /= 10.0L;
        ++exponent;
    }

    std::ostringstream out;
    out << std::fixed << std::setprecision(5) << static_cast<double>(mantissa) << 'e' << exponent;
    return out.str();
}

struct Solver {
    std::vector<int> primes = sieve_primes(kPrimeSearchLimit);

    PrimeSolution S(int p) const {
        return solve_prime(p, primes);
    }

    long double log_T(int limit) const {
        long double total = 0.0L;
        for (int p : primes) {
            if (p >= limit) {
                break;
            }
            total += S(p).log10_value;
        }
        return total;
    }
};

void validate() {
    Solver solver;

    assert(exact_value(solver.S(2)) == 1);
    assert(exact_value(solver.S(5)) == 8);

    std::uint64_t t20 = 1;
    for (int p : solver.primes) {
        if (p >= 20) {
            break;
        }
        t20 *= exact_value(solver.S(p));
    }
    assert(t20 == 1348422598656ULL);

    assert(scientific(solver.log_T(100)) == "1.37451e123");
    std::cerr << "Validation checkpoints passed.\n";
}

}  // namespace

int main() {
    validate();

    Solver solver;
    std::cout << scientific(solver.log_T(kLimit)) << '\n';
    return 0;
}

Python

import math
import sys
from bisect import bisect_left
from functools import reduce
from math import gcd


LIMIT = 20_000
PRIME_SEARCH_LIMIT = 2_000_000


def sieve_primes(n):
    is_prime = bytearray(b"\x01") * (n + 1)
    if n >= 0:
        is_prime[0] = 0
    if n >= 1:
        is_prime[1] = 0

    r = int(math.isqrt(n))
    for i in range(2, r + 1):
        if is_prime[i]:
            start = i * i
            is_prime[start : n + 1 : i] = b"\x00" * (((n - start) // i) + 1)

    return [i for i in range(2, n + 1) if is_prime[i]]


def factorize(n, primes):
    factors = []
    x = n
    for q in primes:
        if q * q > x:
            break
        if x % q:
            continue
        exponent = 0
        while x % q == 0:
            x //= q
            exponent += 1
        factors.append((q, exponent))
    if x > 1:
        factors.append((x, 1))
    return factors


def divisors_from_factors(factors):
    divisors = [1]
    for prime, exponent in factors:
        previous = divisors[:]
        power = 1
        for _ in range(exponent):
            power *= prime
            for d in previous:
                divisors.append(d * power)
    divisors.sort()
    return divisors


def primitive_root(p, small_primes):
    factors = factorize(p - 1, small_primes)
    for g in range(2, p):
        ok = True
        for q, _ in factors:
            if pow(g, (p - 1) // q, p) == 1:
                ok = False
                break
        if ok:
            return g
    return -1


def is_prime_by_trial(n, primes):
    if n < 2:
        return False
    for q in primes:
        if q * q > n:
            return True
        if n % q == 0:
            return n == q
    return True


def smallest_transition_prime(p, quotient_order, target_gcd, discrete_log, primes):
    def works(q):
        if q == p:
            return False
        residue = q % p
        if residue == 0:
            return False
        return gcd(discrete_log[residue], quotient_order) == target_gcd

    for q in primes:
        if works(q):
            return q

    q = primes[-1] + 1
    if q % 2 == 0:
        q += 1
    while True:
        if is_prime_by_trial(q, primes) and works(q):
            return q
        q += 2


class PrimeSolution:
    def __init__(self, log10_value=0.0, factors=None):
        self.log10_value = log10_value
        self.factors = factors or []


def solve_prime(p, primes):
    if p == 2:
        return PrimeSolution()

    n = p - 1
    root = primitive_root(p, primes)
    assert root > 0

    discrete_log = [-1] * p
    x = 1
    for exponent in range(n):
        discrete_log[x] = exponent
        x = x * root % p

    divisors = divisors_from_factors(factorize(n, primes))
    index = {value: i for i, value in enumerate(divisors)}
    inf = float("inf")

    dp = [inf] * len(divisors)
    parent = [-1] * len(divisors)
    edge = [(0, 0)] * len(divisors)
    transition_cache = {}

    dp[index[1]] = 0.0
    for i, subgroup_size in enumerate(divisors):
        if not math.isfinite(dp[i]):
            continue

        quotient_order = n // subgroup_size
        for length in divisors:
            if length <= 1 or quotient_order % length != 0:
                continue

            target_gcd = quotient_order // length
            key = (quotient_order, target_gcd)
            q = transition_cache.get(key)
            if q is None:
                q = smallest_transition_prime(
                    p, quotient_order, target_gcd, discrete_log, primes
                )
                transition_cache[key] = q

            next_size = subgroup_size * length
            j = index[next_size]
            candidate = dp[i] + (length - 1) * math.log10(q)
            if candidate + 1e-24 < dp[j]:
                dp[j] = candidate
                parent[j] = i
                edge[j] = (length, q)

    finish = index[n]
    factors = []
    at = finish
    while parent[at] != -1:
        factors.append(edge[at])
        at = parent[at]
    factors.reverse()
    return PrimeSolution(dp[finish], factors)


def exact_value(solution):
    value = 1
    for length, q in solution.factors:
        value *= q ** (length - 1)
    return value


def scientific(log10_value):
    exponent = math.floor(log10_value)
    mantissa = 10.0 ** (log10_value - exponent)
    mantissa = math.floor(mantissa * 100000.0 + 0.5) / 100000.0
    if mantissa >= 10.0:
        mantissa /= 10.0
        exponent += 1
    return f"{mantissa:.5f}e{exponent}"


class Solver:
    def __init__(self):
        self.primes = sieve_primes(PRIME_SEARCH_LIMIT)
        self._cache = {}

    def S(self, p):
        solution = self._cache.get(p)
        if solution is None:
            solution = solve_prime(p, self.primes)
            self._cache[p] = solution
        return solution

    def log_T(self, limit):
        values = []
        for p in self.primes:
            if p >= limit:
                break
            values.append(self.S(p).log10_value)
        return math.fsum(values)


def validate():
    solver = Solver()

    assert exact_value(solver.S(2)) == 1
    assert exact_value(solver.S(5)) == 8

    t20 = 1
    for p in solver.primes:
        if p >= 20:
            break
        t20 *= exact_value(solver.S(p))
    assert t20 == 1_348_422_598_656

    assert scientific(solver.log_T(100)) == "1.37451e123"
    print("Validation checkpoints passed.", file=sys.stderr)
    return solver


def main():
    solver = validate()
    print(scientific(solver.log_T(LIMIT)))


if __name__ == "__main__":
    main()

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;

public class Euler995 {
    private static final int LIMIT = 20_000;
    private static final int PRIME_SEARCH_LIMIT = 2_000_000;

    private static int[] sievePrimes(int n) {
        boolean[] isPrime = new boolean[n + 1];
        Arrays.fill(isPrime, true);
        if (n >= 0) {
            isPrime[0] = false;
        }
        if (n >= 1) {
            isPrime[1] = false;
        }

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

        int count = 0;
        for (int i = 2; i <= n; ++i) {
            if (isPrime[i]) {
                ++count;
            }
        }

        int[] primes = new int[count];
        int at = 0;
        for (int i = 2; i <= n; ++i) {
            if (isPrime[i]) {
                primes[at++] = i;
            }
        }
        return primes;
    }

    private static ArrayList<int[]> factorize(int n, int[] primes) {
        ArrayList<int[]> factors = new ArrayList<>();
        int x = n;
        for (int q : primes) {
            if ((long) q * q > x) {
                break;
            }
            if (x % q != 0) {
                continue;
            }
            int exponent = 0;
            while (x % q == 0) {
                x /= q;
                ++exponent;
            }
            factors.add(new int[] {q, exponent});
        }
        if (x > 1) {
            factors.add(new int[] {x, 1});
        }
        return factors;
    }

    private static int[] divisorsFromFactors(ArrayList<int[]> factors) {
        ArrayList<Integer> divisors = new ArrayList<>();
        divisors.add(1);
        for (int[] factor : factors) {
            int prime = factor[0];
            int exponent = factor[1];
            ArrayList<Integer> previous = new ArrayList<>(divisors);
            int power = 1;
            for (int e = 1; e <= exponent; ++e) {
                power *= prime;
                for (int d : previous) {
                    divisors.add(d * power);
                }
            }
        }
        int[] out = new int[divisors.size()];
        for (int i = 0; i < divisors.size(); ++i) {
            out[i] = divisors.get(i);
        }
        Arrays.sort(out);
        return out;
    }

    private static long modPow(long base, int exponent, int mod) {
        long value = 1;
        long b = base % mod;
        int e = exponent;
        while (e > 0) {
            if ((e & 1) != 0) {
                value = value * b % mod;
            }
            b = b * b % mod;
            e >>= 1;
        }
        return value;
    }

    private static int primitiveRoot(int p, int[] smallPrimes) {
        ArrayList<int[]> factors = factorize(p - 1, smallPrimes);
        for (int g = 2; g < p; ++g) {
            boolean ok = true;
            for (int[] factor : factors) {
                int q = factor[0];
                if (modPow(g, (p - 1) / q, p) == 1) {
                    ok = false;
                    break;
                }
            }
            if (ok) {
                return g;
            }
        }
        return -1;
    }

    private static int gcd(int a, int b) {
        int x = Math.abs(a);
        int y = Math.abs(b);
        while (y != 0) {
            int t = x % y;
            x = y;
            y = t;
        }
        return x;
    }

    private static boolean isPrimeByTrial(long n, int[] primes) {
        if (n < 2) {
            return false;
        }
        for (int q : primes) {
            if ((long) q * q > n) {
                return true;
            }
            if (n % q == 0) {
                return n == q;
            }
        }
        return true;
    }

    private static int smallestTransitionPrime(
            int p,
            int quotientOrder,
            int targetGcd,
            int[] discreteLog,
            int[] primes) {
        for (int q : primes) {
            if (works(p, quotientOrder, targetGcd, discreteLog, q)) {
                return q;
            }
        }

        long q = primes[primes.length - 1] + 1L;
        if ((q & 1L) == 0) {
            ++q;
        }
        while (true) {
            if (isPrimeByTrial(q, primes)
                    && works(p, quotientOrder, targetGcd, discreteLog, q)) {
                return (int) q;
            }
            q += 2;
        }
    }

    private static boolean works(
            int p,
            int quotientOrder,
            int targetGcd,
            int[] discreteLog,
            long q) {
        if (q == p) {
            return false;
        }
        int residue = (int) (q % p);
        if (residue == 0) {
            return false;
        }
        return gcd(discreteLog[residue], quotientOrder) == targetGcd;
    }

    private static final class PrimeSolution {
        private double log10Value;
        private final ArrayList<int[]> factors = new ArrayList<>();
    }

    private static PrimeSolution solvePrime(int p, int[] primes) {
        PrimeSolution solution = new PrimeSolution();
        if (p == 2) {
            return solution;
        }

        int n = p - 1;
        int root = primitiveRoot(p, primes);
        if (root <= 0) {
            throw new IllegalStateException("primitive root not found");
        }

        int[] discreteLog = new int[p];
        Arrays.fill(discreteLog, -1);
        int x = 1;
        for (int exponent = 0; exponent < n; ++exponent) {
            discreteLog[x] = exponent;
            x = (int) ((long) x * root % p);
        }

        int[] divisors = divisorsFromFactors(factorize(n, primes));
        HashMap<Integer, Integer> index = new HashMap<>();
        for (int i = 0; i < divisors.length; ++i) {
            index.put(divisors[i], i);
        }

        double[] dp = new double[divisors.length];
        Arrays.fill(dp, Double.POSITIVE_INFINITY);
        int[] parent = new int[divisors.length];
        Arrays.fill(parent, -1);
        int[][] edge = new int[divisors.length][2];
        HashMap<Long, Integer> transitionCache = new HashMap<>();

        dp[index.get(1)] = 0.0;
        for (int i = 0; i < divisors.length; ++i) {
            int subgroupSize = divisors[i];
            if (!Double.isFinite(dp[i])) {
                continue;
            }

            int quotientOrder = n / subgroupSize;
            for (int length : divisors) {
                if (length <= 1 || quotientOrder % length != 0) {
                    continue;
                }

                int targetGcd = quotientOrder / length;
                long key = (((long) quotientOrder) << 32) ^ (targetGcd & 0xffffffffL);
                Integer cached = transitionCache.get(key);
                int q;
                if (cached == null) {
                    q = smallestTransitionPrime(
                            p, quotientOrder, targetGcd, discreteLog, primes);
                    transitionCache.put(key, q);
                } else {
                    q = cached;
                }

                int nextSize = subgroupSize * length;
                int j = index.get(nextSize);
                double candidate = dp[i] + (length - 1) * Math.log10(q);
                if (candidate + 1e-24 < dp[j]) {
                    dp[j] = candidate;
                    parent[j] = i;
                    edge[j][0] = length;
                    edge[j][1] = q;
                }
            }
        }

        int finish = index.get(n);
        solution.log10Value = dp[finish];

        ArrayList<int[]> reversed = new ArrayList<>();
        for (int at = finish; parent[at] != -1; at = parent[at]) {
            reversed.add(new int[] {edge[at][0], edge[at][1]});
        }
        for (int i = reversed.size() - 1; i >= 0; --i) {
            solution.factors.add(reversed.get(i));
        }
        return solution;
    }

    private static long exactValue(PrimeSolution solution) {
        long value = 1L;
        for (int[] factor : solution.factors) {
            int length = factor[0];
            int q = factor[1];
            for (int i = 1; i < length; ++i) {
                value *= q;
            }
        }
        return value;
    }

    private static String scientific(double log10Value) {
        long exponent = (long) Math.floor(log10Value);
        double mantissa = Math.pow(10.0, log10Value - exponent);
        mantissa = Math.floor(mantissa * 100000.0 + 0.5) / 100000.0;
        if (mantissa >= 10.0) {
            mantissa /= 10.0;
            ++exponent;
        }
        return String.format(java.util.Locale.ROOT, "%.5fe%d", mantissa, exponent);
    }

    private static final class Solver {
        private final int[] primes = sievePrimes(PRIME_SEARCH_LIMIT);
        private final HashMap<Integer, PrimeSolution> cache = new HashMap<>();

        private PrimeSolution s(int p) {
            PrimeSolution solution = cache.get(p);
            if (solution == null) {
                solution = solvePrime(p, primes);
                cache.put(p, solution);
            }
            return solution;
        }

        private double logT(int limit) {
            double total = 0.0;
            for (int p : primes) {
                if (p >= limit) {
                    break;
                }
                total += s(p).log10Value;
            }
            return total;
        }
    }

    private static Solver validate() {
        Solver solver = new Solver();

        assert exactValue(solver.s(2)) == 1L;
        assert exactValue(solver.s(5)) == 8L;

        long t20 = 1L;
        for (int p : solver.primes) {
            if (p >= 20) {
                break;
            }
            t20 *= exactValue(solver.s(p));
        }
        assert t20 == 1_348_422_598_656L;

        assert scientific(solver.logT(100)).equals("1.37451e123");
        System.err.println("Validation checkpoints passed.");
        return solver;
    }

    public static void main(String[] args) {
        Solver solver = validate();
        System.out.println(scientific(solver.logT(LIMIT)));
    }
}