Project Euler 833

发表于 更新于 分类于 阅读次数： 本文字数： 4.1k 阅读时长 ≈ 4 分钟

Project Euler 833

题目

Square Triangle Products

Triangle numbers \(T_k\) are integers of the form \(\dfrac{k(k+1)} 2\).

A few triangle numbers happen to be perfect squares like \(T_1=1\) and \(T_8=36\), but more can be found when considering the product of two triangle numbers. For example, \(T_2 \cdot T_{24}=3 \cdot 300=30^2\).

Let \(S(n)\) be the sum of \(c\) for all integers triples \((a, b, c)\) with \(0<c \le n\), \(c^2=T_a \cdot T_b\) and \(0<a<b\). For example, \(S(100)= \sqrt{T_1 T_8}+\sqrt{T_2 T_{24}}+\sqrt{T_1 T_{49}}+\sqrt{T_3 T_{48}}=6+30+35+84=155\).

You are given \(S(10^5)=1479802\) and \(S(10^9)=241614948794\).

Find \(S(10^{35})\). Give your answer modulo \(136101521\).

解决方案

1. 把平方条件化为“公共非平方核”

由定义\(T_aT_b=c^2\Longleftrightarrow\dfrac{a(a+1)}2\cdot\dfrac{b(b+1)}2=c^2\Longleftrightarrow a(a+1)b(b+1)=(2c)^2.\)令\(U=a(a+1), V=b(b+1),\)则条件等价为\(UV\)为完全平方数。

一个基本事实是：若 \(UV\) 为完全平方，则存在唯一的表示\(U=km^2, V=kn^2, \gcd(m,n)=1,\)其中 \(k=\gcd(U,V)\)。

证明思路很直接：把公共部分 \(k=\gcd(U,V)\) 提出来后，\(\dfrac Uk\) 与 \(\dfrac Vk\) 互素，而互素乘积若为平方则两者均为平方。

因此对于 \(T_aT_b\) 为平方的情形，可以理解成：存在某个公共核 \(k\)，使得\(T_a=km^2, T_b=kn^2, \gcd(m,n)=1.\)

从\(T_a=km^2\)出发：\(\dfrac{a(a+1)}2=km^2\Longleftrightarrow a(a+1)=2km^2.\)配方\((2a+1)^2=4a(a+1)+1=8km^2+1,\)从而得到佩尔型方程\((2a+1)^2-8km^2=1.\)

记\(X=2a+1, Y=m,\)则\(X^2-8kY^2=1.\)因此：对于固定的 \(k\)，所有满足 \(T_a=km^2\) 的 \(a\)，对应佩尔方程 \(X^2-8kY^2=1\) 的正整数解。

设佩尔方程\(X^2-8kY^2=1\)，存在最小正解（基本解）\((X_1,Y_1)\)，则所有正解满足\(X_i+Y_i\sqrt{8k}=(X_1+Y_1\sqrt{8k})^i.\)

由乘法展开可得递推

\[ \begin{aligned} X_{i+1}&=X_1X_i+8kY_1Y_i,\\ Y_{i+1}&=Y_1X_i+X_1Y_i. \end{aligned} \]

从第二式可见对所有 \(i\) 都有 \(Y_1\mid Y_i\)，这是一个显而易见的归纳：\(Y_1\mid Y_2\)，\(Y_{i+1}\)是\(Y_1\)和\(Y_i\)的线性组合。因此若 \(Y_1>1\)，则任意两项 \(Y_i,Y_j\) 都不可能互素。但我们真正需要的是\(T_a=km^2, T_b=kn^2, \gcd(m,n)=1,\)其中 \(m,n\) 就是对应解中的 \(Y_i,Y_j\)。

为了存在互素的 \(Y_i,Y_j\)，必须有\(Y_1=1.\)于是基本解为 \((X_1,1)\)，并且\(X_1^2-8k=1\Longleftrightarrow k=\dfrac{X_1^2-1}{8}.\)

注意 \(X_1\) 必为奇数，令\(X_1=2x+1(x\ge1),\)则\(k=\dfrac{(2x+1)^2-1}{8}=\dfrac{x(x+1)}2=T_x.\)

这说明：能产生有效三元组 \((a,b,c)\) 的公共核 \(k\)，必然是某个三角数 \(T_x\)。

令\(u=2x+1,\)则\(u^2-1=4x(x+1)=8T_x.\)考虑代数数\(u+\sqrt{u^2-1}.\)其幂展开为\((u+\sqrt{u^2-1})^t=A_t(u)+B_t(u)\sqrt{u^2-1},\)其中 \(A_t,B_t\) 是整系数多项式，并满足递推

\[ \begin{aligned} A_{t+1}&=uA_t+(u^2-1)B_t,\\ B_{t+1}&=A_t+uB_t, \end{aligned} \]

其中\(A_0=1,B_0=0.\)

对比标准的 Chebyshev 第二类多项式 \(U_n\)：\(U_0(u)=1, U_1(u)=2u, U_{n+1}(u)=2uU_n(u)-U_{n-1}(u),\)可以得到

\[ B_t(u)=U_{t-1}(u) (t\ge1). \]

因此当我们取两个指数 \(i<j\) 时，若它们对应同一条序列，则\(T_{a_i}=k\cdot B_i(u)^2, T_{a_j}=k\cdot B_j(u)^2.\)于是\(c=\sqrt{T_{a_i}T_{a_j}}=k\cdot B_i(u)B_j(u).\)

代入 \(k=T_x=\dfrac{u^2-1}{8}\)，得\(c=\dfrac{u^2-1}{8}\cdot B_i(u)B_j(u)=\dfrac{u^2-1}{8}\cdot U_{i-1}(u)U_{j-1}(u).\)因此只要枚举互素指数对 \((i,j)\)，并对不同的 \(x\) 求和即可。

对固定的互素对 \((i,j)\)，定义\(F_{i,j}(x)=(u^2-1),U_{i-1}(u),U_{j-1}(u), u=2x+1.\)那么\(c=\dfrac{F_{i,j}(x)}8.\)因此目标变成

\[ S(n)=\sum_{\substack{1\le i<j\gcd(i,j)=1}} \sum_{\substack{x\ge 1F_{i,j}(x)\le 8n}} \frac{F_{i,j}(x)}8. \]

对每一对 \((i,j)\)，由于 \(U_k(u)\) 随 \(x\) 增长极快，存在最大 \(x_{\max}\) 使得\(F_{i,j}(x)\le 8n \iff x\le x_{\max}.\)于是该对的贡献为\(\displaystyle{\frac18\sum_{x=1}^{x_{\max}}F_{i,j}(x).}\)

注意 \(F_{i,j}(x)\) 本质上是关于 \(x\) 的多项式（通过 \(u=2x+1\) 代入），而 \(U_{i-1}(u)\)、\(U_{j-1}(u)\) 均为多项式，因此 \(F_{i,j}\) 也是多项式；其次数大约为 \(i+j\) 量级。

如果已知某个次数为 \(d\) 的多项式 \(P(x)\) 在\(x=0,1,2,\dots,d\)这些点的取值，那么可以用拉格朗日插值在任意整数点求值。进一步，若定义前缀和\(\displaystyle{Q(m)=\sum_{x=1}^m P(x),}\)则 \(Q(m)\) 是次数 \(d+1\) 的多项式，只需得到 \(Q(0),Q(1),\dots,Q(d+1)\) 就能插值求任意 \(Q(m)\)。

因此每对 \((i,j)\) 的求和可以变成：

  1. 计算 \(F_{i,j}(t)\) 在 \(t=1,2,\dots,d+1\) 的值；
  2. 得到其前缀和数组 \(y[t]=\sum_{x=1}^t F_{i,j}(x)\)；
  3. 用插值直接求 \(y[x_{\max}]\)。

为了找到\(x_{\max}=\max\{x: F_{i,j}(x)\le 8n\},\)可以先用增长级别估计一个初始猜测，再二分修正。由于 \(U_k(u)\) 随 \(u\) 近似指数增长，\(F_{i,j}(x)\) 对 \(x\) 也增长极快，因此二分次数很少。

由于这里需要计算\(\displaystyle{Q(m)=\sum_{x=1}^m P(x)}\)，而我们掌握的是连续点 \(0,1,\dots,k\) 的值。

那么令 \(k=\deg(Q)\)，已知 \(Q(0),Q(1),\dots,Q(k)\)。则

\[ Q(n)=\sum_{i=0}^k Q(i)\cdot \frac{\prod_{t\ne i}(n-t)}{\prod_{t\ne i}(i-t)}. \]

其中分母可写成阶乘形式并预处理逆元：\(\displaystyle{\prod_{t\ne i}(i-t)=(-1)^{k-i} i!(k-i)!.}\)为了避免分母出现除\(0\)的情况，使用前缀和后缀进行结果的拼接。

最终，我们首先预处理出来多项式表\(U_{i-1}(u)\)。我们可以预先对小范围的 \(x\) 建表：

  • 因为对次数为 \(d\) 的多项式前缀和插值，需要 \(x=0,\dots,d+1\)；
  • 本题中 \(d\approx i+j\)，最大检查到 \(i,j\le L\)，其中\(L=60\)即可，所以只需预处理到 \(x\le 2L+2\)。
  • 然后建表递推求\(U_i(u)\)。

为了避免重复计数，只枚举\(1\le i<j\le L, \gcd(i,j)=1.\) 对每对\((i,j)\)：求 \(x_{\max}\)，并取 \(K=(i+j)+1\)，构造前缀和数组 \(y[0,\dots,K]\)；最终插值得到 \(\sum_{x=1}^{x_{\max}}F_{i,j}(x)\)。在最后乘上\(\dfrac{1}{8}\)得到结果。

代码

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
from math import gcd, log2
from tools import Combinatorics

N = 10**35
MOD = 136101521

def lagrange_consecutive(y, n, mod, invfact):
    k = len(y) - 1
    if n <= k:
        return y[n] % mod
    n_mod = n % mod
    pre = [1] * (k + 2)
    suf = [1] * (k + 2)
    for i in range(k + 1):
        pre[i + 1] = pre[i] * (n_mod - i) % mod
    for i in range(k, -1, -1):
        suf[i] = suf[i + 1] * (n_mod - i) % mod
    ans = 0
    for i in range(k + 1):
        num = pre[i] * suf[i + 1] % mod
        den = invfact[i] * invfact[k - i] % mod
        if (k - i) & 1:
            den = (-den) % mod
        ans = (ans + y[i] * num % mod * den) % mod
    return ans


def chebU_int(n, u):
    if n == 0:
        return 1
    if n == 1:
        return 2 * u
    a, b = 1, 2 * u
    for _ in range(2, n + 1):
        a, b = b, 2 * u * b - a
    return b


def eval_F_int(i, j, x):
    u = 2 * x + 1
    Ui = chebU_int(i - 1, u)
    Uj = chebU_int(j - 1, u)
    return (u * u - 1) * Ui * Uj


def find_xmax(i, j, n):
    bound = 8 * n
    if eval_F_int(i, j, 1) > bound:
        return 0
    d = i + j
    e = 2 * d - 2
    lg = log2(bound)
    approx_exp = (lg - e) / d
    guess = int(2 ** max(approx_exp, 0))
    l = max(0, guess - 200)
    r = max(1, guess + 200)
    while eval_F_int(i, j, r) <= bound:
        r *= 2
    while l < r:
        mid = (l + r + 1) // 2
        if eval_F_int(i, j, mid) <= bound:
            l = mid
        else:
            r = mid - 1
    return l


def solve(n):
    L = 60
    maxK = 2 * L + 5
    co = Combinatorics(maxK, MOD)
    fact, invfact = co.fact, co.inv_fact
    inv8 = pow(8, MOD - 2, MOD)

    max_need_x = 2 * L + 2
    Utab = [[0] * L for _ in range(max_need_x + 1)]
    for x in range(max_need_x + 1):
        u = 2 * x + 1
        Utab[x][0] = 1
        if L > 1:
            Utab[x][1] = (2 * u) % MOD
        two_u = (2 * u) % MOD
        for t in range(2, L):
            Utab[x][t] = (two_u * Utab[x][t - 1] - Utab[x][t - 2]) % MOD

    ans = 0
    for i in range(1, L + 1):
        for j in range(i + 1, L + 1):
            if gcd(i, j) != 1:
                continue
            x_max = find_xmax(i, j, n)
            if x_max == 0:
                continue
            deg = i + j
            K = deg + 1
            y = [0] * (K + 1)
            acc = 0
            for t in range(1, K + 1):
                u = 2 * t + 1
                Ui = Utab[t][i - 1]
                Uj = Utab[t][j - 1]
                f = ((u * u - 1) % MOD) * Ui % MOD * Uj % MOD
                acc = (acc + f) % MOD
                y[t] = acc
            s = lagrange_consecutive(y, x_max, MOD, invfact)
            ans = (ans + s) % MOD

    ans = ans * inv8 % MOD
    return ans

print(solve(N))