调和级数的高精度快速计算

mathe

可以转化为${Gamma'(n+1)}/{Gamma(n+1)}-{Gamma'(1)}/{Gamma(1)}$

wayne

官方资料Mathematica technology guide 证实了Mathematica在很多地方用到了AAS： (Automatic Algorithm Selection)

mathe

http://en.wikipedia.org/wiki/Digamma_function

wayne

竟然可以用解析式来表达 http://en.wikipedia.org/wiki/Polygamma 在Mathematica里面是函数PolyGamma PolyGamma[n + 1] - PolyGamma[1] $\psi (n+1)-\psi (1)$

云梦

当n很大时：Sum[1/k,{k,1,n}]~Log[n]+EulerGamma
n不是很大时：用以下公式：

云梦

以1024位有效数字为例：
n=894:  k=400:       BernoulliB[800]
n=1000,k=377:       BernoulliB[654]
n=2000,k=288:       BernoulliB[576]
n=10000:k=193:     BernoulliB[386]
n=100000:k=134:    BernoulliB[268]
n=1000000:k=104:  BernoullliB[208]
n=10000000:k=85:  BernoulliB[170]
n=10^10:  k=55
n=10^15:  k=35
n=10^20:  k=26
n=10^30:  k=17
n=10^50:  k=10
n=10^100: k=5
n=10^256: k=0
从上面计算可知：
（1） 当 n<894 时，用 ∑（1/k,{k,1,n}）
（2） 当 894<=n<10^256时：用Hn公式
（3） 当 n>10^256时：用Ln(n)+γ+1/2n
（4） 当 n>10^1024时：可直接用 Ln(n) +γ
γ=0.57721566490153286060651209008240......  

我用VB计算器计算：
n=10^1000 用时0.42秒
n=10^100         0.46秒
n=10^10          0.63秒
n=10^8            0.69秒
n=894              2.00秒
     
n=893              0.63秒
m=800             0.56秒

n=2399 用 (1)计算      1.59秒
n=2400 用(2)计算       1.59秒
真对上述分析：把 n设定在2400作为算法转换点是合理的。计算n=2399 和 n=2400 用（1）（2）式时间相同，不会造成算法转换后时间的差异。

云梦

当 n>2400 时，n可以为非整数。
而 n<2400 时,且为非整数时或者n<0时用PolyGamma[n+1]+EulerGamma替代计算。

云梦

vb调和级数（HarmonicNumber）计算程序：

Dim Ex As Chen
Dim N2 As Chen
Dim qx As Chen
Dim FX As Chen
Dim Hx As Chen
Dim Px As Chen
Dim i, n, Fz As Long
Ex = Ax


Ax = Ex
Bx.Bz = ""
Bx.St = "24"
Bx.zs = 3
hsfx = "-"
Call Addition
If Ax.Bz = "" And Ex.Bz = "" Then
    'n^-2
    Cx = Ex
    Ax = Ex
    Call Mult_
    Ax = Cx
    Cx.Bz = ""
    Cx.St = "1"
    Cx.zs = 0
    Call Multx_
    N2 = Cx
    'Ln(n)+EulerGamma+1/2n
    Ax = Ex
    Call Lnx
    Bx.Bz = ""
    Bx.St = Euler.Tag
    Bx.zs = -1
    hsfx = "+"
    Call Addition
    FX = Ax
    '1/2n
    Ax = Ex
    Bx = Ax
    hsfx = "+"
    Call Addition
    Cx.Bz = ""
    Cx.St = "1"
    Cx.zs = 0
    Call Multx_
    Ax = FX
    Bx = Cx
    hsfx = "+"
    Call Addition
    FX = Ax
   
    Hx.Bz = ""
    Hx.St = "0"
    Hx.zs = 0
   
    Px = Hx
    Px.St = "1"
    For i = 1 To 400
        'p/n^2->p   1/n^2k
        Ax = Px
        Cx = N2
        Call Mult_
        Px = Cx
        'B2k/2k
        Ax = pz2(2 * i)
        Call Mult_
        Ax.Bz = ""
        Ax.St = CStr(2 * i)
        Ax.zs = Len(Ax.St) - 1
        Call Multxs_
        qx = Cx
        
        Ax = Hx
        Bx = qx
        hsfx = "+"
        Call Addition
        Hx = Ax
        If qx.zs < -xChen Then Exit For
    Next
    Bx = Ax
    Ax = FX
    hsfx = "-"
    Call Addition
Else
    If Ex.zs > -1 Then
        '是否为整数
        Fz = 0
        Ax = Ex
        Bx = Ax
        Bx.St = Left(Bx.St, Bx.zs + 1)
        hsfx = "-"
        Call Addition
        If Ax.St > "0" Then Fz = 1
    Else
        Fz = 1
    End If
   
    If Fz = 0 And Ex.Bz = "" Then
        'n=整数时
        Hx.Bz = ""
        Hx.St = "0"
        Hx.zs = 0
        
        Px = Hx
        Px.St = "1"
        
        n = Val(Left(Ex.St, Ex.zs + 1))
        
        For i = 1 To n
            Cx = Px
            Ax.Bz = ""
            Ax.St = CStr(i)
            Ax.zs = Len(Ax.St) - 1
            Call Multxs_
            Ax = Hx
            Bx = Cx
            hsfx = "+"
            Call Addition
            Hx = Ax
        Next
    Else
        If Fz = 1 Then
            Ax = Ex
            Bx.Bz = ""
            Bx.St = "1"
            Bx.zs = 0
            hsfx = "+"
            Call Addition
            Call psi
            Bx.Bz = ""
            Bx.St = Euler.Tag
            Bx.zs = -1
            hsfx = "+"
            Call Addition
        Else
           Text12.Text = "ComplexInfinity"
        End If
    End If
End If

王守恩

调和级数前任意项求和