Problem 851: SOP and POS

Project Euler Problem 851 Solution

EulerSolve provides an optimized solution for Project Euler Problem 851, SOP and POS, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $E_n=(\mathbb{Z}_{>0})^n$. For two positive-integer vectors $u=(u_1,\dots,u_n)$ and $v=(v_1,\dots,v_n)$, define $$\langle u,v\rangle=\sum_{i=1}^{n}u_iv_i,\qquad u\star v=\prod_{i=1}^{n}(u_i+v_i).$$ The quantity $R_n(M)$ is the sum of $u\star v$ over all ordered pairs $(u,v)\in E_n\times E_n$ satisfying $\langle u,v\rangle=M$. The target is $$R_6(10000!) \pmod{10^9+7}.$$ A direct enumeration of six-dimensional vector pairs is hopeless, so the solution converts the problem into coefficient extraction from a generating function and then reduces that generating function with quasimodular-form identities. Mathematical Approach The computation rests on two ideas: first, each coordinate contributes independently through a one-variable divisor sum; second, the sixth power of the resulting generating series can be rewritten in a weight-12 quasimodular basis whose coefficients are easy to evaluate at $10000!$. Step 1: Collapse One Coordinate to a Divisor Sum For $n=1$, the condition $\langle u,v\rangle=M$ is simply $uv=M$. Every divisor $d\mid M$ gives the ordered pair $(d,M/d)$, whose contribution to $u\star v$ is $$d+\frac{M}{d}.$$ Therefore $$R_1(M)=\sum_{d\mid M}\left(d+\frac{M}{d}\right)=2\sum_{d\mid M}d=2\sigma_1(M),$$ where $\sigma_1(M)$ is the usual sum-of-divisors function....

Detailed mathematical approach

Problem Summary

Let $E_n=(\mathbb{Z}_{>0})^n$. For two positive-integer vectors $u=(u_1,\dots,u_n)$ and $v=(v_1,\dots,v_n)$, define

$$\langle u,v\rangle=\sum_{i=1}^{n}u_iv_i,\qquad u\star v=\prod_{i=1}^{n}(u_i+v_i).$$

The quantity $R_n(M)$ is the sum of $u\star v$ over all ordered pairs $(u,v)\in E_n\times E_n$ satisfying $\langle u,v\rangle=M$. The target is

$$R_6(10000!) \pmod{10^9+7}.$$

A direct enumeration of six-dimensional vector pairs is hopeless, so the solution converts the problem into coefficient extraction from a generating function and then reduces that generating function with quasimodular-form identities.

Mathematical Approach

The computation rests on two ideas: first, each coordinate contributes independently through a one-variable divisor sum; second, the sixth power of the resulting generating series can be rewritten in a weight-12 quasimodular basis whose coefficients are easy to evaluate at $10000!$.

Step 1: Collapse One Coordinate to a Divisor Sum

For $n=1$, the condition $\langle u,v\rangle=M$ is simply $uv=M$. Every divisor $d\mid M$ gives the ordered pair $(d,M/d)$, whose contribution to $u\star v$ is

$$d+\frac{M}{d}.$$

Therefore

$$R_1(M)=\sum_{d\mid M}\left(d+\frac{M}{d}\right)=2\sum_{d\mid M}d=2\sigma_1(M),$$

where $\sigma_1(M)$ is the usual sum-of-divisors function.

Step 2: Turn Six Coordinates into an Additive Convolution

For a general pair of vectors, set

$$m_i=u_iv_i \qquad (1\le i\le n).$$

Then

$$m_1+\cdots+m_n=M,$$

and for a fixed decomposition $M=m_1+\cdots+m_n$, the contributions from the coordinates multiply. Hence

$$R_n(M)=\sum_{\substack{m_1+\cdots+m_n=M\\m_i\ge 1}}R_1(m_1)\cdots R_1(m_n).$$

If we define the generating series

$$G(q)=\sum_{m\ge 1}R_1(m)q^m,$$

then this is exactly the $n$-fold additive convolution identity

$$R_n(M)=[q^M]\,G(q)^n.$$

Step 3: Replace the Generating Series by Eisenstein Series

The normalized quasimodular Eisenstein series $E_2$ satisfies

$$E_2(q)=1-24\sum_{m\ge 1}\sigma_1(m)q^m.$$

Since $R_1(m)=2\sigma_1(m)$, we obtain

$$G(q)=2\sum_{m\ge 1}\sigma_1(m)q^m=\frac{1-E_2(q)}{12}.$$

For the required dimension $n=6$, this gives

$$R_6(M)=[q^M]\left(\frac{1-E_2(q)}{12}\right)^6.$$

Because $M=10000!>0$, the constant term contributes nothing, so the task becomes the extraction of the positive-$q$ coefficients of powers of $E_2$.

Step 4: Reduce Powers of $E_2$ in Weight 12

Write

$$D=q\frac{d}{dq}.$$

For any series $f(q)=\sum_{m\ge 0}a_mq^m$, we have

$$[q^M]D^r f=M^r a_M.$$

The implementation uses the standard quasimodular reductions

$$E_2^2=E_4+12DE_2,$$

$$E_2^3=E_6+9DE_4+72D^2E_2,$$

$$E_2^4=E_8+8DE_6+\frac{216}{5}D^2E_4+288D^3E_2,$$

$$E_2^5=E_{10}+\frac{15}{2}DE_8+\frac{240}{7}D^2E_6+144D^3E_4+864D^4E_2,$$

$$E_2^6=E_{12}-\frac{4608}{24185}\Delta+\frac{36}{5}DE_{10}+30D^2E_8+\frac{720}{7}D^3E_6+\frac{2592}{7}D^4E_4+\frac{10368}{5}D^5E_2.$$

Here $\Delta(q)=\sum_{m\ge 1}\tau(m)q^m$ is the discriminant cusp form. The positive-$q$ coefficients of the Eisenstein series are

$$[q^M]E_2=-24\sigma_1(M),\qquad [q^M]E_4=240\sigma_3(M),\qquad [q^M]E_6=-504\sigma_5(M),$$

$$[q^M]E_8=480\sigma_7(M),\qquad [q^M]E_{10}=-264\sigma_9(M),\qquad [q^M]E_{12}=\frac{65520}{691}\sigma_{11}(M),$$

and

$$[q^M]\Delta=\tau(M).$$

So every coefficient needed for $(1-E_2)^6$ becomes a linear combination of $M^r\sigma_{2k-1}(M)$ and $\tau(M)$.

Step 5: Evaluate the Arithmetic Data at $M=10000!$

Let $N=10000$. For every prime $p\le N$, Legendre's formula gives the exponent of $p$ in $N!$:

$$v_p(N!)=\sum_{t\ge 1}\left\lfloor\frac{N}{p^t}\right\rfloor.$$

Thus

$$10000!=\prod_{p\le 10000}p^{v_p(10000!)}.$$

For every odd $r\in\{1,3,5,7,9,11\}$, multiplicativity gives

$$\sigma_r(10000!)=\prod_{p\le 10000}\left(1+p^r+\cdots+p^{r\,v_p(10000!)}\right)=\prod_{p\le 10000}\frac{p^{r(v_p(10000!)+1)}-1}{p^r-1}.$$

The Ramanujan tau function is multiplicative as well, so

$$\tau(10000!)=\prod_{p^e\parallel 10000!}\tau(p^e),$$

with the prime-power recurrence

$$\tau(p^e)=\tau(p)\tau(p^{e-1})-p^{11}\tau(p^{e-2}).$$

To obtain the prime values $\tau(p)$, the implementation first builds the sequence $\tau(n)$ up to $10000$ from

$$\tau(1)=1,\qquad (n-1)\tau(n)=-24\sum_{k=1}^{n-1}\sigma_1(k)\tau(n-k)\qquad (n\ge 2).$$

Worked Example: $R_1(10)$ and $R_2(3)$

For $M=10$, the divisor pairs are $(1,10)$, $(2,5)$, $(5,2)$, and $(10,1)$. Their contributions are $11$, $7$, $7$, and $11$, so

$$R_1(10)=11+7+7+11=36=2(1+2+5+10).$$

Now use convolution for $n=2$ and $M=3$. The only decompositions into positive parts are $3=1+2$ and $3=2+1$. Since $R_1(1)=2$ and $R_1(2)=2(1+2)=6$, we get

$$R_2(3)=R_1(1)R_1(2)+R_1(2)R_1(1)=2\cdot 6+6\cdot 2=24.$$

This small case shows exactly why generating functions turn the problem into a power of one basic series.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they sieve the primes up to $10000$ and compute every $v_p(10000!)$. From those exponents they evaluate $\sigma_1,\sigma_3,\sigma_5,\sigma_7,\sigma_9,\sigma_{11}$ multiplicatively modulo $10^9+7$, and they also compute $10000! \bmod (10^9+7)$ so that the factors $M,M^2,\dots,M^5$ needed after differentiation are available.

Next the implementation builds $\sigma_1(n)$ for all $1\le n\le 10000$, uses the recurrence above to obtain $\tau(n)$ up to $10000$, and then extends from prime values to prime powers through the second-order recurrence for $\tau(p^e)$. A $2\times 2$ matrix power is used to evaluate that recurrence efficiently for each prime power.

Finally the program substitutes all divisor sums, the $\tau$ term, and the powers of $10000!$ into the weight-12 reduction of $(1-E_2)^6$. The last step is division by $12^6=2985984$, implemented modulo $10^9+7$ by multiplying with the modular inverse.

Complexity Analysis

Let $N=10000$. The prime sieve and the factorial valuations cost $O(N\log\log N)$ time. Building the table of $\sigma_1(n)$ by divisor accumulation costs $O(N\log N)$. The recurrence for $\tau(n)$ is quadratic, $O(N^2)$, and this is the dominant step in the actual implementations. The remaining prime-product evaluations and prime-power recurrences are lower-order, roughly $O(\pi(N)\log N)$. Memory usage is $O(N)$.

Footnotes and References

Problem page: https://projecteuler.net/problem=851
Divisor function: Wikipedia - Divisor function
Eisenstein series: Wikipedia - Eisenstein series
Quasimodular form: Wikipedia - Quasimodular form
Ramanujan tau function: Wikipedia - Ramanujan tau function
Legendre's formula: Wikipedia - Legendre's formula

Problem 851 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <map>
#include <mutex>
#include <thread>
#include <utility>
#include <vector>

using namespace std;

namespace {

constexpr long long kMod = 1'000'000'007LL;

long long mod_add(long long a, long long b) {
    a += b;
    if (a >= kMod) a -= kMod;
    return a;
}

long long mod_sub(long long a, long long b) {
    a -= b;
    if (a < 0) a += kMod;
    return a;
}

long long mod_mul(long long a, long long b) {
    return static_cast<long long>((__int128)a * b % kMod);
}

long long mod_pow(long long a, long long e) {
    long long r = 1 % kMod;
    a %= kMod;
    if (a < 0) a += kMod;
    while (e > 0) {
        if (e & 1) r = mod_mul(r, a);
        a = mod_mul(a, a);
        e >>= 1;
    }
    return r;
}

long long mod_inv(long long a) {
    return mod_pow(a, kMod - 2);
}

vector<int> sieve_primes(int n) {
    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; i * 1LL * i <= n; ++i) {
        if (is_prime[i]) {
            for (long long j = 1LL * i * i; j <= n; j += i) {
                is_prime[static_cast<size_t>(j)] = false;
            }
        }
    }
    vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[i]) primes.push_back(i);
    }
    return primes;
}

int vp_factorial(int n, int p) {
    int e = 0;
    while (n) {
        n /= p;
        e += n;
    }
    return e;
}

long long factorial_mod(int n) {
    long long r = 1;
    for (int i = 2; i <= n; ++i) {
        r = mod_mul(r, i);
    }
    return r;
}

long long geom_sum(long long base, long long e) {
    base %= kMod;
    if (base < 0) base += kMod;
    if (e < 0) return 0;
    if (base == 1) return (e + 1) % kMod;
    long long num = mod_sub(mod_pow(base, e + 1), 1);
    long long den = mod_sub(base, 1);
    return mod_mul(num, mod_inv(den));
}

long long sigma_k_factorial(const vector<pair<int, int>>& pe, int k) {
    long long res = 1;
    for (auto [p, e] : pe) {
        long long base = mod_pow(p, k);
        long long s = geom_sum(base, e);
        res = mod_mul(res, s);
    }
    return res;
}

struct Mat2 {
    long long a00, a01, a10, a11;
};

Mat2 mat_mul(const Mat2& A, const Mat2& B) {
    Mat2 C;
    C.a00 = (mod_mul(A.a00, B.a00) + mod_mul(A.a01, B.a10)) % kMod;
    C.a01 = (mod_mul(A.a00, B.a01) + mod_mul(A.a01, B.a11)) % kMod;
    C.a10 = (mod_mul(A.a10, B.a00) + mod_mul(A.a11, B.a10)) % kMod;
    C.a11 = (mod_mul(A.a10, B.a01) + mod_mul(A.a11, B.a11)) % kMod;
    return C;
}

Mat2 mat_pow(Mat2 base, long long e) {
    Mat2 r{1, 0, 0, 1};
    while (e > 0) {
        if (e & 1) r = mat_mul(r, base);
        base = mat_mul(base, base);
        e >>= 1;
    }
    return r;
}

long long tau_prime_power(long long tau_p, int p, int e) {
    if (e == 0) return 1;
    if (e == 1) return tau_p;
    long long p11 = mod_pow(p, 11);
    Mat2 M{tau_p, mod_sub(0, p11), 1, 0};
    Mat2 P = mat_pow(M, e - 1);
    return (mod_mul(P.a00, tau_p) + mod_mul(P.a01, 1)) % kMod;
}

vector<long long> compute_sigma1_upto(int N) {
    vector<long long> sig(N + 1, 0);
    for (int d = 1; d <= N; ++d) {
        for (int m = d; m <= N; m += d) {
            sig[m] += d;
        }
    }
    for (int i = 0; i <= N; ++i) sig[i] %= kMod;
    return sig;
}

vector<long long> compute_tau_upto(int N, const vector<long long>& sigma1) {
    vector<long long> tau(N + 1, 0);
    tau[0] = 0;
    tau[1] = 1;
    const long long neg24 = mod_sub(0, 24);
    for (int n = 2; n <= N; ++n) {
        long long s = 0;
        for (int k = 1; k < n; ++k) {
            s = (s + sigma1[k] * tau[n - k]) % kMod;
        }
        long long inv = mod_inv(n - 1);
        tau[n] = mod_mul(mod_mul(neg24, s), inv);
    }
    return tau;
}

long long compute_tau_factorial_parallel(const vector<pair<int, int>>& pe,
                                         const vector<long long>& tau_small) {
    unsigned hw = thread::hardware_concurrency();
    unsigned T = hw ? hw : 4;
    if (T > 8) T = 8;

    vector<long long> partial(T, 1);
    vector<thread> threads;
    threads.reserve(T);
    int m = static_cast<int>(pe.size());

    auto worker = [&](unsigned tid) {
        int L = static_cast<int>((long long)m * tid / T);
        int R = static_cast<int>((long long)m * (tid + 1) / T);
        long long prod = 1;
        for (int idx = L; idx < R; ++idx) {
            int p = pe[idx].first;
            int e = pe[idx].second;
            long long tau_p = tau_small[p];
            long long t = tau_prime_power(tau_p, p, e);
            prod = mod_mul(prod, t);
        }
        partial[tid] = prod;
    };

    for (unsigned t = 0; t < T; ++t) threads.emplace_back(worker, t);
    for (auto& th : threads) th.join();

    long long total = 1;
    for (unsigned t = 0; t < T; ++t) total = mod_mul(total, partial[t]);
    return total;
}

long long conv_Rn_small(int n, int M) {
    vector<long long> sig1 = compute_sigma1_upto(M);
    vector<long long> f(M + 1, 0);
    for (int m = 1; m <= M; ++m) {
        f[m] = (2 * sig1[m]) % kMod;
    }
    vector<long long> dp(M + 1, 0), ndp(M + 1, 0);
    dp[0] = 1;
    for (int iter = 0; iter < n; ++iter) {
        fill(ndp.begin(), ndp.end(), 0);
        for (int s = 0; s <= M; ++s) {
            if (!dp[s]) continue;
            for (int m = 1; s + m <= M; ++m) {
                if (!f[m]) continue;
                ndp[s + m] = (ndp[s + m] + dp[s] * f[m]) % kMod;
            }
        }
        dp.swap(ndp);
    }
    return dp[M];
}

}  // namespace

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    const int N = 10000;

    vector<int> primes = sieve_primes(N);
    vector<pair<int, int>> pe;
    pe.reserve(primes.size());
    for (int p : primes) {
        pe.push_back({p, vp_factorial(N, p)});
    }

    long long n_mod = factorial_mod(N);
    long long n_pow[6];
    n_pow[0] = 1;
    for (int i = 1; i <= 5; ++i) n_pow[i] = mod_mul(n_pow[i - 1], n_mod);

    map<int, long long> sig;
    vector<int> ks = {1, 3, 5, 7, 9, 11};
    mutex sig_mtx;
    vector<thread> sig_threads;
    for (int k : ks) {
        sig_threads.emplace_back([&, k]() {
            long long v = sigma_k_factorial(pe, k);
            lock_guard<mutex> lk(sig_mtx);
            sig[k] = v;
        });
    }
    for (auto& th : sig_threads) th.join();

    vector<long long> sigma1_small = compute_sigma1_upto(N);
    vector<long long> tau_small = compute_tau_upto(N, sigma1_small);

    auto norm = [](long long x) {
        x %= kMod;
        if (x < 0) x += kMod;
        return x;
    };

    if (norm(tau_small[1]) != 1 || norm(tau_small[2]) != norm(-24) ||
        norm(tau_small[3]) != 252 || norm(tau_small[5]) != 4830) {
        cerr << "[Validation failed] tau initial values mismatch\n";
        return 1;
    }

    long long tau_fact = compute_tau_factorial_parallel(pe, tau_small);

    long long c2 = mod_mul(norm(-24), sig[1]);
    long long c4 = mod_mul(240, sig[3]);
    long long c6 = mod_mul(norm(-504), sig[5]);
    long long c8 = mod_mul(480, sig[7]);
    long long c10 = mod_mul(norm(-264), sig[9]);
    long long factor12 = mod_mul(65520, mod_inv(691));
    long long c12 = mod_mul(factor12, sig[11]);

    long long inv2 = mod_inv(2);
    long long inv5 = mod_inv(5);
    long long inv7 = mod_inv(7);
    long long inv24185 = mod_inv(24185);

    long long a1 = c2;

    long long a2 = mod_add(c4, mod_mul(12, mod_mul(n_pow[1], c2)));

    long long a3 = c6;
    a3 = mod_add(a3, mod_mul(9, mod_mul(n_pow[1], c4)));
    a3 = mod_add(a3, mod_mul(72, mod_mul(n_pow[2], c2)));

    long long a4 = c8;
    a4 = mod_add(a4, mod_mul(8, mod_mul(n_pow[1], c6)));
    a4 = mod_add(a4, mod_mul(mod_mul(216, inv5), mod_mul(n_pow[2], c4)));
    a4 = mod_add(a4, mod_mul(288, mod_mul(n_pow[3], c2)));

    long long a5 = c10;
    a5 = mod_add(a5, mod_mul(mod_mul(15, inv2), mod_mul(n_pow[1], c8)));
    a5 = mod_add(a5, mod_mul(mod_mul(240, inv7), mod_mul(n_pow[2], c6)));
    a5 = mod_add(a5, mod_mul(144, mod_mul(n_pow[3], c4)));
    a5 = mod_add(a5, mod_mul(864, mod_mul(n_pow[4], c2)));

    long long a6 = c12;
    a6 = mod_add(a6, mod_mul(mod_mul(norm(-4608), inv24185), tau_fact));
    a6 = mod_add(a6, mod_mul(mod_mul(36, inv5), mod_mul(n_pow[1], c10)));
    a6 = mod_add(a6, mod_mul(30, mod_mul(n_pow[2], c8)));
    a6 = mod_add(a6, mod_mul(mod_mul(720, inv7), mod_mul(n_pow[3], c6)));
    a6 = mod_add(a6, mod_mul(mod_mul(2592, inv7), mod_mul(n_pow[4], c4)));
    a6 = mod_add(a6, mod_mul(mod_mul(10368, inv5), mod_mul(n_pow[5], c2)));

    long long S = 0;
    S = mod_add(S, mod_mul(norm(-6), a1));
    S = mod_add(S, mod_mul(15, a2));
    S = mod_add(S, mod_mul(norm(-20), a3));
    S = mod_add(S, mod_mul(15, a4));
    S = mod_add(S, mod_mul(norm(-6), a5));
    S = mod_add(S, a6);

    long long inv2985984 = mod_inv(2985984);
    long long ans = mod_mul(S, inv2985984);

    {
        long long r1_10 = (2 * (1 + 2 + 5 + 10)) % kMod;
        if (r1_10 != 36) {
            cerr << "[Validation failed] R1(10) expected 36\n";
            return 1;
        }
        long long r2_100 = conv_Rn_small(2, 100);
        if (r2_100 != 1873044) {
            cerr << "[Validation failed] R2(100) expected 1873044, got " << r2_100
                 << "\n";
            return 1;
        }
        long long r6_20_brut = conv_Rn_small(6, 20);

        auto sigma_small = [&](int n, int k) -> long long {
            long long res = 0;
            for (int d = 1; d * 1LL * d <= n; ++d) {
                if (n % d != 0) continue;
                res = (res + mod_pow(d, k)) % kMod;
                if (d * 1LL * d != n) {
                    res = (res + mod_pow(n / d, k)) % kMod;
                }
            }
            return res;
        };
        auto coeffE = [&](int n, int w) -> long long {
            if (w == 2) return mod_mul(norm(-24), sigma_small(n, 1));
            if (w == 4) return mod_mul(240, sigma_small(n, 3));
            if (w == 6) return mod_mul(norm(-504), sigma_small(n, 5));
            if (w == 8) return mod_mul(480, sigma_small(n, 7));
            if (w == 10) return mod_mul(norm(-264), sigma_small(n, 9));
            if (w == 12) return mod_mul(factor12, sigma_small(n, 11));
            return 0LL;
        };

        long long tau20 = tau_small[20];
        long long n20 = 20;
        long long n20p[6];
        n20p[0] = 1;
        for (int i = 1; i <= 5; ++i) n20p[i] = mod_mul(n20p[i - 1], n20);

        long long c2_20 = coeffE(20, 2);
        long long c4_20 = coeffE(20, 4);
        long long c6_20 = coeffE(20, 6);
        long long c8_20 = coeffE(20, 8);
        long long c10_20 = coeffE(20, 10);
        long long c12_20 = coeffE(20, 12);

        long long a1_20 = c2_20;
        long long a2_20 = mod_add(c4_20, mod_mul(12, mod_mul(n20p[1], c2_20)));

        long long a3_20 = c6_20;
        a3_20 = mod_add(a3_20, mod_mul(9, mod_mul(n20p[1], c4_20)));
        a3_20 = mod_add(a3_20, mod_mul(72, mod_mul(n20p[2], c2_20)));

        long long a4_20 = c8_20;
        a4_20 = mod_add(a4_20, mod_mul(8, mod_mul(n20p[1], c6_20)));
        a4_20 = mod_add(a4_20, mod_mul(mod_mul(216, inv5), mod_mul(n20p[2], c4_20)));
        a4_20 = mod_add(a4_20, mod_mul(288, mod_mul(n20p[3], c2_20)));

        long long a5_20 = c10_20;
        a5_20 =
            mod_add(a5_20, mod_mul(mod_mul(15, inv2), mod_mul(n20p[1], c8_20)));
        a5_20 =
            mod_add(a5_20, mod_mul(mod_mul(240, inv7), mod_mul(n20p[2], c6_20)));
        a5_20 = mod_add(a5_20, mod_mul(144, mod_mul(n20p[3], c4_20)));
        a5_20 = mod_add(a5_20, mod_mul(864, mod_mul(n20p[4], c2_20)));

        long long a6_20 = c12_20;
        a6_20 = mod_add(a6_20, mod_mul(mod_mul(norm(-4608), inv24185), tau20));
        a6_20 =
            mod_add(a6_20, mod_mul(mod_mul(36, inv5), mod_mul(n20p[1], c10_20)));
        a6_20 = mod_add(a6_20, mod_mul(30, mod_mul(n20p[2], c8_20)));
        a6_20 =
            mod_add(a6_20, mod_mul(mod_mul(720, inv7), mod_mul(n20p[3], c6_20)));
        a6_20 =
            mod_add(a6_20, mod_mul(mod_mul(2592, inv7), mod_mul(n20p[4], c4_20)));
        a6_20 =
            mod_add(a6_20, mod_mul(mod_mul(10368, inv5), mod_mul(n20p[5], c2_20)));

        long long S20 = 0;
        S20 = mod_add(S20, mod_mul(norm(-6), a1_20));
        S20 = mod_add(S20, mod_mul(15, a2_20));
        S20 = mod_add(S20, mod_mul(norm(-20), a3_20));
        S20 = mod_add(S20, mod_mul(15, a4_20));
        S20 = mod_add(S20, mod_mul(norm(-6), a5_20));
        S20 = mod_add(S20, a6_20);

        long long r6_20_form = mod_mul(S20, inv2985984);
        if (r6_20_form != r6_20_brut) {
            cerr << "[Validation failed] R6(20) mismatch: brut=" << r6_20_brut
                 << " form=" << r6_20_form << "\n";
            return 1;
        }
    }

    cout << ans << "\n";
    return 0;
}

Python

import sys

kMod = 1000000007

def mod_add(a, b):
    a += b
    if a >= kMod: a -= kMod
    return a

def mod_sub(a, b):
    a -= b
    if a < 0: a += kMod
    return a

def mod_mul(a, b):
    return (a * b) % kMod

def mod_pow(a, e):
    r = 1 % kMod
    a %= kMod
    if a < 0: a += kMod
    while e > 0:
        if e & 1: r = mod_mul(r, a)
        a = mod_mul(a, a)
        e >>= 1
    return r

def mod_inv(a):
    return mod_pow(a, kMod - 2)

def sieve_primes(n):
    is_prime = [True] * (n + 1)
    if n >= 0: is_prime[0] = False
    if n >= 1: is_prime[1] = False
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False
    primes = [i for i in range(2, n + 1) if is_prime[i]]
    return primes

def vp_factorial(n, p):
    e = 0
    while n:
        n //= p
        e += n
    return e

def factorial_mod(n):
    r = 1
    for i in range(2, n + 1):
        r = mod_mul(r, i)
    return r

def geom_sum(base, e):
    base %= kMod
    if base < 0: base += kMod
    if e < 0: return 0
    if base == 1: return (e + 1) % kMod
    num = mod_sub(mod_pow(base, e + 1), 1)
    den = mod_sub(base, 1)
    return mod_mul(num, mod_inv(den))

def sigma_k_factorial(pe, k):
    res = 1
    for p, e in pe:
        base = mod_pow(p, k)
        s = geom_sum(base, e)
        res = mod_mul(res, s)
    return res

def mat_mul(A, B):
    C = [0] * 4
    C[0] = (mod_mul(A[0], B[0]) + mod_mul(A[1], B[2])) % kMod
    C[1] = (mod_mul(A[0], B[1]) + mod_mul(A[1], B[3])) % kMod
    C[2] = (mod_mul(A[2], B[0]) + mod_mul(A[3], B[2])) % kMod
    C[3] = (mod_mul(A[2], B[1]) + mod_mul(A[3], B[3])) % kMod
    return C

def mat_pow(base, e):
    r = [1, 0, 0, 1]
    while e > 0:
        if e & 1: r = mat_mul(r, base)
        base = mat_mul(base, base)
        e >>= 1
    return r

def tau_prime_power(tau_p, p, e):
    if e == 0: return 1
    if e == 1: return tau_p
    p11 = mod_pow(p, 11)
    M = [tau_p, mod_sub(0, p11), 1, 0]
    P = mat_pow(M, e - 1)
    return (mod_mul(P[0], tau_p) + mod_mul(P[1], 1)) % kMod

def compute_sigma1_upto(N):
    sig = [0] * (N + 1)
    for d in range(1, N + 1):
        for m in range(d, N + 1, d):
            sig[m] += d
    for i in range(N + 1):
        sig[i] %= kMod
    return sig

def compute_tau_upto(N, sigma1):
    tau = [0] * (N + 1)
    tau[1] = 1
    neg24 = mod_sub(0, 24)
    for n in range(2, N + 1):
        s = 0
        for k in range(1, n):
            s = (s + sigma1[k] * tau[n - k]) % kMod
        inv = mod_inv(n - 1)
        tau[n] = mod_mul(mod_mul(neg24, s), inv)
    return tau

def compute_tau_factorial_serial(pe, tau_small):
    prod = 1
    for p, e in pe:
        tau_p = tau_small[p]
        t = tau_prime_power(tau_p, p, e)
        prod = mod_mul(prod, t)
    return prod

def solve():
    N = 10000
    primes = sieve_primes(N)
    pe = [(p, vp_factorial(N, p)) for p in primes]
    
    n_mod = factorial_mod(N)
    n_pow = [1] * 6
    for i in range(1, 6):
        n_pow[i] = mod_mul(n_pow[i - 1], n_mod)
        
    sig = {}
    ks = [1, 3, 5, 7, 9, 11]
    for k in ks:
        sig[k] = sigma_k_factorial(pe, k)
        
    sigma1_small = compute_sigma1_upto(N)
    tau_small = compute_tau_upto(N, sigma1_small)
    
    def norm(x):
        x %= kMod
        if x < 0: x += kMod
        return x
        
    tau_fact = compute_tau_factorial_serial(pe, tau_small)
    
    c2 = mod_mul(norm(-24), sig[1])
    c4 = mod_mul(240, sig[3])
    c6 = mod_mul(norm(-504), sig[5])
    c8 = mod_mul(480, sig[7])
    c10 = mod_mul(norm(-264), sig[9])
    factor12 = mod_mul(65520, mod_inv(691))
    c12 = mod_mul(factor12, sig[11])
    
    inv2 = mod_inv(2)
    inv5 = mod_inv(5)
    inv7 = mod_inv(7)
    inv24185 = mod_inv(24185)
    
    a1 = c2
    a2 = mod_add(c4, mod_mul(12, mod_mul(n_pow[1], c2)))
    
    a3 = c6
    a3 = mod_add(a3, mod_mul(9, mod_mul(n_pow[1], c4)))
    a3 = mod_add(a3, mod_mul(72, mod_mul(n_pow[2], c2)))
    
    a4 = c8
    a4 = mod_add(a4, mod_mul(8, mod_mul(n_pow[1], c6)))
    a4 = mod_add(a4, mod_mul(mod_mul(216, inv5), mod_mul(n_pow[2], c4)))
    a4 = mod_add(a4, mod_mul(288, mod_mul(n_pow[3], c2)))
    
    a5 = c10
    a5 = mod_add(a5, mod_mul(mod_mul(15, inv2), mod_mul(n_pow[1], c8)))
    a5 = mod_add(a5, mod_mul(mod_mul(240, inv7), mod_mul(n_pow[2], c6)))
    a5 = mod_add(a5, mod_mul(144, mod_mul(n_pow[3], c4)))
    a5 = mod_add(a5, mod_mul(864, mod_mul(n_pow[4], c2)))
    
    a6 = c12
    a6 = mod_add(a6, mod_mul(mod_mul(norm(-4608), inv24185), tau_fact))
    a6 = mod_add(a6, mod_mul(mod_mul(36, inv5), mod_mul(n_pow[1], c10)))
    a6 = mod_add(a6, mod_mul(30, mod_mul(n_pow[2], c8)))
    a6 = mod_add(a6, mod_mul(mod_mul(720, inv7), mod_mul(n_pow[3], c6)))
    a6 = mod_add(a6, mod_mul(mod_mul(2592, inv7), mod_mul(n_pow[4], c4)))
    a6 = mod_add(a6, mod_mul(mod_mul(10368, inv5), mod_mul(n_pow[5], c2)))
    
    S = 0
    S = mod_add(S, mod_mul(norm(-6), a1))
    S = mod_add(S, mod_mul(15, a2))
    S = mod_add(S, mod_mul(norm(-20), a3))
    S = mod_add(S, mod_mul(15, a4))
    S = mod_add(S, mod_mul(norm(-6), a5))
    S = mod_add(S, a6)
    
    inv2985984 = mod_inv(2985984)
    ans = mod_mul(S, inv2985984)
    
    return str(ans)

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

Java

import java.util.ArrayList;

public class Euler851 {

    static final long kMod = 1000000007L;

    static long modAdd(long a, long b) {
        a += b;
        if (a >= kMod)
            a -= kMod;
        return a;
    }

    static long modSub(long a, long b) {
        a -= b;
        if (a < 0)
            a += kMod;
        return a;
    }

    static long modMul(long a, long b) {
        long aMod = (a % kMod + kMod) % kMod;
        long bMod = (b % kMod + kMod) % kMod;
        return (aMod * bMod) % kMod;
    }

    static long modPow(long a, long e) {
        long r = 1 % kMod;
        a %= kMod;
        if (a < 0)
            a += kMod;
        while (e > 0) {
            if ((e & 1) == 1)
                r = modMul(r, a);
            a = modMul(a, a);
            e >>= 1;
        }
        return r;
    }

    static long modInv(long a) {
        return modPow(a, kMod - 2);
    }

    static ArrayList<Integer> sievePrimes(int n) {
        boolean[] isPrime = new boolean[n + 1];
        java.util.Arrays.fill(isPrime, true);
        if (n >= 0)
            isPrime[0] = false;
        if (n >= 1)
            isPrime[1] = false;
        for (int i = 2; i * i <= n; ++i) {
            if (isPrime[i]) {
                for (int j = i * i; j <= n; j += i) {
                    isPrime[j] = false;
                }
            }
        }
        ArrayList<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; ++i) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static int vpFactorial(int n, int p) {
        int e = 0;
        while (n > 0) {
            n /= p;
            e += n;
        }
        return e;
    }

    static long factorialMod(int n) {
        long r = 1;
        for (int i = 2; i <= n; ++i) {
            r = modMul(r, i);
        }
        return r;
    }

    static long geomSum(long base, long e) {
        base %= kMod;
        if (base < 0)
            base += kMod;
        if (e < 0)
            return 0;
        if (base == 1)
            return (e + 1) % kMod;
        long num = modSub(modPow(base, e + 1), 1);
        long den = modSub(base, 1);
        return modMul(num, modInv(den));
    }

    static class Pair {
        int p, e;

        Pair(int p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static long sigmaKFactorial(ArrayList<Pair> pe, int k) {
        long res = 1;
        for (Pair pair : pe) {
            long base = modPow(pair.p, k);
            long s = geomSum(base, pair.e);
            res = modMul(res, s);
        }
        return res;
    }

    static class Mat2 {
        long a00, a01, a10, a11;

        Mat2() {
        }

        Mat2(long a00, long a01, long a10, long a11) {
            this.a00 = a00;
            this.a01 = a01;
            this.a10 = a10;
            this.a11 = a11;
        }
    }

    static Mat2 matMul(Mat2 A, Mat2 B) {
        Mat2 C = new Mat2();
        C.a00 = (modMul(A.a00, B.a00) + modMul(A.a01, B.a10)) % kMod;
        C.a01 = (modMul(A.a00, B.a01) + modMul(A.a01, B.a11)) % kMod;
        C.a10 = (modMul(A.a10, B.a00) + modMul(A.a11, B.a10)) % kMod;
        C.a11 = (modMul(A.a10, B.a01) + modMul(A.a11, B.a11)) % kMod;
        return C;
    }

    static Mat2 matPow(Mat2 base, long e) {
        Mat2 r = new Mat2(1, 0, 0, 1);
        while (e > 0) {
            if ((e & 1) == 1)
                r = matMul(r, base);
            base = matMul(base, base);
            e >>= 1;
        }
        return r;
    }

    static long tauPrimePower(long tauP, int p, int e) {
        if (e == 0)
            return 1;
        if (e == 1)
            return tauP;
        long p11 = modPow(p, 11);
        Mat2 M = new Mat2(tauP, modSub(0, p11), 1, 0);
        Mat2 P = matPow(M, e - 1);
        return (modMul(P.a00, tauP) + modMul(P.a01, 1)) % kMod;
    }

    static long[] computeSigma1Upto(int N) {
        long[] sig = new long[N + 1];
        for (int d = 1; d <= N; ++d) {
            for (int m = d; m <= N; m += d) {
                sig[m] += d;
            }
        }
        for (int i = 0; i <= N; ++i)
            sig[i] %= kMod;
        return sig;
    }

    static long[] computeTauUpto(int N, long[] sigma1) {
        long[] tau = new long[N + 1];
        tau[0] = 0;
        tau[1] = 1;
        long neg24 = modSub(0, 24);
        for (int n = 2; n <= N; ++n) {
            long s = 0;
            for (int k = 1; k < n; ++k) {
                s = (s + sigma1[k] * tau[n - k]) % kMod;
            }
            long inv = modInv(n - 1);
            tau[n] = modMul(modMul(neg24, s), inv);
        }
        return tau;
    }

    static long computeTauFactorialSerial(ArrayList<Pair> pe, long[] tauSmall) {
        long prod = 1;
        for (Pair pair : pe) {
            int p = pair.p;
            int e = pair.e;
            long tauP = tauSmall[p];
            long t = tauPrimePower(tauP, p, e);
            prod = modMul(prod, t);
        }
        return prod;
    }

    static long norm(long x) {
        x %= kMod;
        if (x < 0)
            x += kMod;
        return x;
    }

    public static String solve() {
        int N = 10000;
        ArrayList<Integer> primes = sievePrimes(N);
        ArrayList<Pair> pe = new ArrayList<>();
        for (int p : primes) {
            pe.add(new Pair(p, vpFactorial(N, p)));
        }

        long nMod = factorialMod(N);
        long[] nPow = new long[6];
        nPow[0] = 1;
        for (int i = 1; i <= 5; ++i)
            nPow[i] = modMul(nPow[i - 1], nMod);

        long[] sig = new long[12];
        int[] ks = { 1, 3, 5, 7, 9, 11 };
        for (int k : ks) {
            sig[k] = sigmaKFactorial(pe, k);
        }

        long[] sigma1Small = computeSigma1Upto(N);
        long[] tauSmall = computeTauUpto(N, sigma1Small);

        long tauFact = computeTauFactorialSerial(pe, tauSmall);

        long c2 = modMul(norm(-24), sig[1]);
        long c4 = modMul(240, sig[3]);
        long c6 = modMul(norm(-504), sig[5]);
        long c8 = modMul(480, sig[7]);
        long c10 = modMul(norm(-264), sig[9]);
        long factor12 = modMul(65520, modInv(691));
        long c12 = modMul(factor12, sig[11]);

        long inv2 = modInv(2);
        long inv5 = modInv(5);
        long inv7 = modInv(7);
        long inv24185 = modInv(24185);

        long a1 = c2;

        long a2 = modAdd(c4, modMul(12, modMul(nPow[1], c2)));

        long a3 = c6;
        a3 = modAdd(a3, modMul(9, modMul(nPow[1], c4)));
        a3 = modAdd(a3, modMul(72, modMul(nPow[2], c2)));

        long a4 = c8;
        a4 = modAdd(a4, modMul(8, modMul(nPow[1], c6)));
        a4 = modAdd(a4, modMul(modMul(216, inv5), modMul(nPow[2], c4)));
        a4 = modAdd(a4, modMul(288, modMul(nPow[3], c2)));

        long a5 = c10;
        a5 = modAdd(a5, modMul(modMul(15, inv2), modMul(nPow[1], c8)));
        a5 = modAdd(a5, modMul(modMul(240, inv7), modMul(nPow[2], c6)));
        a5 = modAdd(a5, modMul(144, modMul(nPow[3], c4)));
        a5 = modAdd(a5, modMul(864, modMul(nPow[4], c2)));

        long a6 = c12;
        a6 = modAdd(a6, modMul(modMul(norm(-4608), inv24185), tauFact));
        a6 = modAdd(a6, modMul(modMul(36, inv5), modMul(nPow[1], c10)));
        a6 = modAdd(a6, modMul(30, modMul(nPow[2], c8)));
        a6 = modAdd(a6, modMul(modMul(720, inv7), modMul(nPow[3], c6)));
        a6 = modAdd(a6, modMul(modMul(2592, inv7), modMul(nPow[4], c4)));
        a6 = modAdd(a6, modMul(modMul(10368, inv5), modMul(nPow[5], c2)));

        long S = 0;
        S = modAdd(S, modMul(norm(-6), a1));
        S = modAdd(S, modMul(15, a2));
        S = modAdd(S, modMul(norm(-20), a3));
        S = modAdd(S, modMul(15, a4));
        S = modAdd(S, modMul(norm(-6), a5));
        S = modAdd(S, a6);

        long inv2985984 = modInv(2985984);
        long ans = modMul(S, inv2985984);

        return Long.toString(ans);
    }

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