求n∈[1,10^14]时∑1/(c*n+d)的快速算法

核心思路 #

利用**Digamma函数（ψ函数）**的渐近展开式，将超大规模求和转化为常数项计算+大参数渐近近似，时间复杂度为O(1)，完全适配N≈10^14的场景。

公式转化 #

首先将求和式转化为Digamma函数的形式：

S(N) = ∑ₙ=1^N 1/(c*n + d) = (1/c) * [ψ(N + 1 + d/c) - ψ(1 + d/c)]

其中：

• ψ(z)是Digamma函数，为调和数H_{z-1}的解析延拓，满足递推关系ψ(z+1) = ψ(z) + 1/z，且ψ(1) = -γ（γ为欧拉-Mascheroni常数，≈0.57721566490153286）。

分步实现 #

1. 预处理常数项ψ(1 + d/c)

   • 若d/c为整数k（即d=kc）：ψ(1+k) = H_k - γ，其中H_k是第k个调和数，固定参数下只需计算一次。
   • 若d/c为既约分数p/q：利用Digamma函数的有理分式公式计算：

     ψ(1 + p/q) = -γ + (π/2)cot(πp/q) + ∑_{k=1}^{q-1} cos(2πkp/q) * ln(2sin(πk/q))
     

   • 通用场景：也可通过递推公式ψ(z) = ψ(z+n) - ∑_{k=0}^{n-1} 1/(z+k)，将z转化为接近1的数后用已知近似值计算。

2. 计算大参数项ψ(N + 1 + d/c)
   当N≥10^14时，z=N+1+d/c是极大值，使用Digamma函数的渐近展开式，取前3-5项即可达到双精度浮点数的精度要求：

   ψ(z) ≈ ln(z) - 1/(2z) - 1/(12z²) - 1/(120z⁴) - 1/(252z⁶)
   

   为避免大数字计算的精度损失，可将z拆分为N*(1 + (1+d/c)/N)，则：

   ln(z) = ln(N) + ln(1 + (1+d/c)/N) ≈ ln(N) + (1+d/c)/N - (1+d/c)²/(2N²) + (1+d/c)³/(3N³)
   

   代入渐近式后计算即可。

3. 合并结果
   将两步的结果代入公式，计算(ψ(大参数项) - ψ(常数项))/c，得到S(N)的高精度近似值。

优化技巧 #

• 先约分：计算g = gcd(c, d)，令c'=c/g，d'=d/g，则S(N) = (1/g)*∑ₙ=1^N 1/(c'*n + d')，减少后续计算的参数规模。
• 精度控制：若只需粗略估计，取渐近展开前两项（ln(z)-1/(2z)）即可；若需更高精度，可增加展开项数。

示例验证 #

以c=2，d=1，N=10^14为例：

• 常数项ψ(1+1/2)=ψ(3/2)=2 - γ - 2ln2≈0.036489974
• 大参数项z=10^14 + 1.5，ψ(z)≈ln(10^14+1.5)-1/(2*(10^14+1.5))≈32.2361913819
• S(N)=(32.2361913819 - 0.036489974)/2≈16.0998507039，与实际求和的精确值误差小于1e-14。

内容的提问来源于stack exchange，提问作者user19565054