求倒数和的算法

Question

是否有任何快速算法来评估形式(a * n + b) / (c * n + d)的a, b, c, d固定和n范围从1到大约10^14左右?

显然，由于总和的大小，单独对每一项求和是行不通的。

编辑:对1 / (c * n + d)求和的算法就足够了，因为您可以将分数拆分并在O(1)时间内对每个分子求和。

Answer 1

你可以将求和化简为n + α的逆求和，这就得到了一个移位的调和数，对应于Digamma函数。该值渐近于ln(n) + Γ(欧拉常数)，并对从1到第1层(α)的缺失/额外初始项进行了校正。

对于Digamma函数的近似，请检查https://en.wikipedia.org/wiki/Digamma_function Computation_and_approximation

Answer 2

这里有一个近似函数，对于较大的x值特别准确，正如这里所解释的。

f = lambda x: 7/(6*x**15) - 691/(2730*x**13) + 5/(66*x**11) - 1/(30*x**9) + 1/(42*x**7) - 1/(30*x**5) + 1/(6*x**3) + 1/(2*x**2) + 1/x
# or 
f = lambda x: 0.97354060901812177/x**2 + 1/(x + 0.4849142940227510) # relative low accuracy and optimized for absolute error
eulers_mascheroni = 0.5772156649015328606065120900824024310421
print(f((eulers_mascheroni+(12*5)+7))) # obviously use an another language if you had speed bottlenecks

尽管你可以在大多数语言中使用库来为你做最好的实现

对于python，您可以使用scipy或mpmath(也用于浮点运算)

from mpmath import *
mp.dps = 50; mp.pretty = True
print(psi(1,eulers_mascheroni+(12*5)+7))

如果你要使用它，请查看那里的文档，他们对每个函数都有很好的解释。

from scipy import special
from numpy import arange
x = [*map(lambda c,n,d: c * n + d + eulers_mascheroni,
range(172),arange(24,42,.04),range(96))]
print(special.polygamma(1, x))