Project Euler 851

发表于 分类于 阅读次数： 本文字数： 3.5k 阅读时长 ≈ 3 分钟

Project Euler 851

题目

SOP and POS

Let \(n\) be a positive integer and let \(E_n\) be the set of \(n\)-tuples of strictly positive integers.

For \(u = (u_1, \cdots, u_n)\) and \(v = (v_1, \cdots, v_n)\) two elements of \(E_n\), we define:

  • the Sum Of Products of \(u\) and \(v\), denoted by \(\langle u, v\rangle\), as the sum \(\displaystyle\sum_{i = 1}^n u_i v_i\);
  • the Product Of Sums of \(u\) and \(v\), denoted by \(u \star v\), as the product \(\displaystyle\prod_{i = 1}^n (u_i + v_i)\).

Let \(R_n(M)\) be the sum of \(u \star v\) over all ordered pairs \((u, v)\) in \(E_n\) such that \(\langle u, v\rangle = M\).

For example: \(R_1(10) = 36\), \(R_2(100) = 1873044\), \(R_2(100!) \equiv 446575636 \bmod 10^9 + 7\).

Find \(R_6(10000!)\). Give your answer modulo \(10^9+7\).

解决方案

令\(N=10^4,p=10^9+7\)。当 \(n=1\) 时，约束是 \(uv=M\)，权重是 \(u+v\)。于是

\[ R_1(M)=\sum_{\substack{u,v\ge 1\\uv=M}}(u+v) =\sum_{d\mid M}\left(d+\frac{M}{d}\right) =2\sum_{d\mid M}d =2\sigma_1(M). \]

对一般 \(n\)，令 \(m_i=u_i v_i\)。约束 \(\langle u,v\rangle=M\) 等价于\(m_1+m_2+\cdots+m_n=M, m_i\ge 1.\)同时权重可按坐标分离：\(\displaystyle u\star v=\prod_{i=1}^n (u_i+v_i).\)

固定一个分拆 \((m_1,\dots,m_n)\) 后，各坐标 \((u_i,v_i)\) 的选择彼此独立，且第 \(i\) 坐标的总贡献恰是 \(R_1(m_i)\)。因此 \[ R_n(M)=\sum_{\substack{m_1+\cdots+m_n=M\\m_i\ge 1}}\prod_{i=1}^n R_1(m_i). \] 这说明 \(R_n\) 是 \(R_1\) 的 \(n\) 次加性卷积。

令生成函数\(\displaystyle\mathcal{R}_n(q)=\sum_{M\ge 1} R_n(M)q^M,\)则卷积对应乘法：\(\displaystyle\mathcal{R}_n(q)=\bigl(\mathcal{R}_1(q)\bigr)^n.\)此外有\(\displaystyle\mathcal{R}_1(q)=\sum_{M\ge 1}2\sigma_1(M)q^M=2\sum_{m\ge 1}\sigma_1(m)q^m.\)

论文给出Eisenstein 级数：\(\displaystyle E_2(q)=1-24\sum_{n\ge 1}\sigma_1(n)q^n,\)因此\(\displaystyle\sum_{n\ge 1}\sigma_1(n)q^n=\frac{1-E_2(q)}{24},\mathcal{R}_1(q)=2\cdot\frac{1-E_2(q)}{24}=\frac{1-E_2(q)}{12}.\)于是\(\displaystyle\mathcal{R}_6(q)=\left(\frac{1-E_2(q)}{12}\right)^6,R_6(M)=[q^M]\left(\frac{1-E_2(q)}{12}\right)^6.\)

此外，论文介绍了使用Ramanujan \(\tau\) 函数\(\tau(n)\)。同时，论文在引言里给出：对满足条件的卷积和\(\displaystyle S_{e,f}(n)=\sum_{k=1}^{n-1}\sigma_e(k)\sigma_f(n-k),\) 可以写成若干 \(\sigma_j(n)\)、\(n\sigma_{e+f-1}(n)\) 以及 \(\tau(n)\) 的有理线性组合（对应其公式 1.4、1.5 的叙述）。

最关键的是论文公式 1.14：\(\displaystyle S_{5,5}(n)=\sum_{k=1}^{n-1}\sigma_5(k)\sigma_5(n-k)=\frac{65}{174132}\sigma_{11}(n)+\frac{1}{252}\sigma_5(n)-\frac{3}{691}\tau(n).\)因此，当推导中出现 \(S_{5,5}(n)\) 时，其化简结果必然包含 \(\tau(n)\) 分量，且该分量不能仅通过 \(\sigma\) 函数项表示。由此，\(R_6(n)\) 的闭式中含有 \(\tau(n)\) 是由卷积恒等式直接推出的结论，并非额外引入的假设。

把 \(\mathcal R_1^6\) 做 Cauchy 乘积，会产生多层卷积和；每一层都可用论文中已给出的卷积恒等式继续化简。在这个过程中，所有项最终都会落到有限类基函数上：$ n^a_{2b-1}(n),(b=1,,6,0a-b)\(和\)(n).$

于是可设待定系数模板 \[ R_6(n)=c_1(n)\sigma_1(n)+c_3(n)\sigma_3(n)+c_5(n)\sigma_5(n)+c_7(n)\sigma_7(n)+c_9(n)\sigma_9(n)+c_{11}\sigma_{11}(n)+c_\tau,\tau(n), \] 其中 \(c_1,\dots,c_9\) 是关于 \(n\) 的多项式（次数分别不超过 \(5,4,3,2,1\)），\(c_{11},c_\tau\) 为常数。

再用前若干个 \(n\) 的精确值做线性方程组求解，即可唯一确定所有有理系数。化简后得到：

\[ \begin{aligned} c_1(n)&=\frac{1}{20736}-\frac{5n}{3456}+\frac{5n^2}{432}-\frac{5n^3}{144}+\frac{n^4}{24}-\frac{n^5}{60},\\ c_3(n)&=\frac{25}{20736}-\frac{25n}{1728}+\frac{5n^2}{96}-\frac{5n^3}{72}+\frac{5n^4}{168},\\ c_5(n)&=\frac{35}{10368}-\frac{35n}{1728}+\frac{5n^2}{144}-\frac{5n^3}{288},\\ c_7(n)&=\frac{25}{10368}-\frac{25n}{3456}+\frac{25n^2}{5184},\\ c_9(n)&=\frac{11}{20736}-\frac{11n}{17280},\\ c_{11}&=\frac{455}{14328576},\\ c_\tau&=-\frac{1}{15671880}. \end{aligned} \]

令\(v_q(N!)\)表示\(N!\)的质因子\(q\)的数量，那么由于\(\sigma_k\) 是乘法函数，且对素幂\(q^e\)有\(\displaystyle\sigma_k(q^{e})=1+q^k+q^{2k}+\cdots+q^{ek}=\frac{q^{k(e+1)}-1}{q^k-1}.\) 于是

\[ \sigma_k(N!)=\prod_{q\le N}\frac{q^{k(v_p(N!)+1)}-1}{q^k-1}. \]

Ramanujan \(\tau\) 函数的页面说明，\(\tau\)同样是乘法函数。对素幂满足二阶递推：\(\tau(q^{e+1})=\tau(q)\tau(q^{e})-q^{11}\tau(q^{e-1}).\) 因此只要知道所有素数 \(q\le N\) 的 \(\tau(q)\) 就能用矩阵快速幂在 \(O(\log e_q)\) 内算出 \(\tau(q^{e_q})\)，再相乘得到 \(\tau(N!)\)。

由此合并结果最终可以计算得到\(R_6(N!)\)的值。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
from tools import Factorizer, int_sqrt

MOD = 1000000007
N = 10000

SIGMA_POWS = [1, 3, 5, 7, 9, 11]

COEFF_TERMS = [
    [(1, 20736, 0, 1), (5, 3456, 1, -1), (5, 432, 2, 1), (5, 144, 3, -1), (1, 24, 4, 1), (1, 60, 5, -1)],
    [(25, 20736, 0, 1), (25, 1728, 1, -1), (5, 96, 2, 1), (5, 72, 3, -1), (5, 168, 4, 1)],
    [(35, 10368, 0, 1), (35, 1728, 1, -1), (5, 144, 2, 1), (5, 288, 3, -1)],
    [(25, 10368, 0, 1), (25, 3456, 1, -1), (25, 5184, 2, 1)],
    [(11, 20736, 0, 1), (11, 17280, 1, -1)],
    [(455, 14328576, 0, 1)],
    [(1, 15671880, 0, -1)],
]

def inv_upto(n):
    inv = [0] * (n + 1)
    inv[1] = 1
    for i in range(2, n + 1):
        inv[i] = (MOD - (MOD // i) * inv[MOD % i] % MOD) % MOD
    return inv

def tau_upto(N):
    inv = inv_upto(N)
    tau = [0] * (N + 1)
    if N >= 1:
        tau[1] = 1
    for n in range(2, N + 1):
        bn = (int_sqrt(8 * n + 1) - 1) // 2
        s = 0
        sign = 1
        for m in range(1, bn + 1):
            t = m * (m + 1) // 2
            term = tau[n - t] * ((n - 1 - 9 * t) % MOD) % MOD
            term = term * (2 * m + 1) % MOD
            if sign:
                s += term
            else:
                s -= term
            s %= MOD
            sign ^= 1
        tau[n] = s * inv[n - 1] % MOD
    return tau

def factorial_prime_exponents(n, primes):
    exp = {}
    for p in primes:
        e = 0
        m = n
        while m:
            m //= p
            e += m
        exp[p] = e
    return exp

def mat_mul(a, b):
    a00, a01, a10, a11 = a
    b00, b01, b10, b11 = b
    return ((a00 * b00 + a01 * b10) % MOD,
            (a00 * b01 + a01 * b11) % MOD,
            (a10 * b00 + a11 * b10) % MOD,
            (a10 * b01 + a11 * b11) % MOD)

def mat_pow(mat, e):
    res = (1, 0, 0, 1)
    base = mat
    while e:
        if e & 1:
            res = mat_mul(res, base)
        base = mat_mul(base, base)
        e >>= 1
    return res

def tau_prime_power(p, e, tp):
    if e == 0:
        return 1
    if e == 1:
        return tp % MOD
    p11 = pow(p, 11, MOD)
    A = (tp % MOD, (-p11) % MOD, 1, 0)
    M = mat_pow(A, e - 1)
    return (M[0] * (tp % MOD) + M[1]) % MOD

def eval_terms(terms, pows, inv_den):
    s = 0
    for num, den, pw, sg in terms:
        v = num % MOD
        v = v * pows[pw] % MOD
        v = v * inv_den[den] % MOD
        if sg == 1:
            s += v
        else:
            s -= v
        s %= MOD
    return s

def solve(N):
    primes = Factorizer(N).primes
    exp = factorial_prime_exponents(N, primes)

    n_mod = 1
    for i in range(2, N + 1):
        n_mod = n_mod * i % MOD

    pows = [1] * 6
    for i in range(1, 6):
        pows[i] = pows[i - 1] * n_mod % MOD

    dens = set()
    for ts in COEFF_TERMS:
        for _, den, _, _ in ts:
            dens.add(den)
    inv_den = {d: pow(d, MOD - 2, MOD) for d in dens}

    sig = []
    for k in SIGMA_POWS:
        prod = 1
        for p in primes:
            e = exp[p]
            pk = pow(p, k, MOD)
            num = (pow(pk, e + 1, MOD) - 1) % MOD
            den = (pk - 1) % MOD
            prod = prod * num % MOD * pow(den, MOD - 2, MOD) % MOD
        sig.append(prod)

    tau_tbl = tau_upto(N)
    t = 1
    for p in primes:
        t = t * tau_prime_power(p, exp[p], tau_tbl[p]) % MOD

    coeffs = [eval_terms(ts, pows, inv_den) for ts in COEFF_TERMS]

    vals = sig + [t]
    ans = 0
    for c, v in zip(coeffs, vals):
        ans = (ans + c * v) % MOD
    return ans

print(solve(N))