证实了Mathematica在很多地方用到了AAS：
(Automatic Algorithm Selection) http://en.wikipedia.org/wiki/Digamma_function 竟然可以用解析式来表达
http://en.wikipedia.org/wiki/Polygamma
在Mathematica里面是函数PolyGamma
PolyGamma - PolyGamma
\psi (n+1)-\psi (1) 当n很大时：Sum~Log+EulerGamma
n不是很大时：用以下公式： 以1024位有效数字为例：
n=894:k=400: BernoulliB
n=1000,k=377: BernoulliB
n=2000,k=288: BernoulliB
n=10000:k=193: BernoulliB
n=100000:k=134: BernoulliB
n=1000000:k=104:BernoullliB
n=10000000:k=85:BernoulliB
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）式时间相同，不会造成算法转换后时间的差异。 本帖最后由 云梦 于 2013-8-15 16:18 编辑
当 n>2400 时，n可以为非整数。
而 n<2400 时,且为非整数时或者n<0时用PolyGamma+EulerGamma替代计算。
本帖最后由 云梦 于 2013-8-15 16:19 编辑
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 调和级数前任意项求和
页:
1
[2]