二分快速方法

----翻译,添改自 Mathematical Constants and computation 的

Binary splitting method

许多公式, 常数的计算, 用经典方法都要 O(n²) 的时间去计算 n 位有效数字 ( 例如 arctan的级数,e的级数见{无穷级数} ). 使用{FFT乘法}和以下介绍的二分法, 可以令计算在 O(n log(n)³) (或 O(n log(n)²) 计算 e) 的时间里完成.

为什么二分法可以加快计算速度？

看一个简单的例子: 计算 A的8次方. 即 X=AAAAAAAA

用最易理解的方法,即反复计算 X=XA 七次(X初值是A):

设 A 是一个2位十进制.即 3>lg(A)>=1,其中lg(z)=log(z)/log(10)

1. X=A 这一步消耗的时间可忽略

2. X=XA 需要 2*2 的时间, X 现在有 4 位数, 总消耗时间 4

3. X=XA 需要 4*2 的时间, X 现在有 6 位数, 总消耗时间 12

4. X=XA 需要 6*2 的时间, X 现在有 8 位数, 总消耗时间 24

5. X=XA 需要 8*2 的时间, X 现在有 10 位数, 总消耗时间 40

6. X=XA 需要 10*2 的时间, X 现在有 12 位数, 总消耗时间 60

7. X=XA 需要 12*2 的时间, X 现在有 14 位数, 总消耗时间 84

8. X=XA 需要 14*2 的时间, X 现在有 16 位数, 总消耗时间 112

再用二分法试一下:

假设 A 还是2位.

1. X=AAAA*AAAA=X1*X1 均分成两部分.

2. X1=AA*AA =X2*X2 其中一部分再均分.

3. X2=A*A 可以直接计算了, 耗时 2*2 , X2 现有 4 位数, 总耗时 4

4. X1=X2*X2 耗时 4*4 , X1 现有 8 位数, 总耗时 20

5. X =X1*X1 耗时 8*8 , X 现有 16 位数, 总耗时 84

计算时间之比: 112:84=1.33…, 或者说后者是前者的 84/112=75%.

更一般的, 设 A 有 L 位数(或者 L=lgA), 求Aⁿ，注意到用一般方法每步需要的时间是等差数列,公差为 L²,有n-1项, 所以耗时 n(n-1)L²/2.而二分法每步耗时 L², 4L², 16L² ... 是等比数列.所以耗时

后者是前者耗时的

→2/3. 约为 66.7%.

设 A 有 L 位数, 用O(nlogn)的乘法, 求Aⁿn(n-1)L²/2

可见不论是n增长还是L增长，最终 T=O(nlogn). 非一般法的 T=O(n²) 可比。

1.1 阶乘.

假如你要计算很多位数的 n!. 一个基本的方法是计算下列值.

x₁ = 1, x₂ = 2x₁, x₃ = 3x₂, ¼, x_n = nx_n-1.

值 x_k 有 O(klog(k)) 位数(斯特林公式: ), 因此计算 x_k+1 = kx_k 有 O(klog(k)). 最后它的总时间复杂度是:

O(2log(2)+¼+n log(n)) = O(n²log(n)).

改用二分法(binary splitting) 递归地对半分割 n 个数的乘积.

设

p(a,b) º (a+1)(a+2) ¼(b-1) b =

明显有递推关系:

p(a,b) = p

æ
ç
è

a+b

ö
÷
ø

æ
ç
è

a+b

ö
÷
ø

这样处理后, 当中涉及到的乘法乘数和被乘数的位数都相差不太大, 可以利用快速乘法提高效率. 当然最好是能做到两乘数的大小几乎一样, 因为这样通常可以减少计算步骤, 又充分利用了快速乘法. 最后n!=p(0,n).

1.2 推广

上例给了我们一个很好的启示. 事实上, 许多级数都可以用类似的办法来处理, 但是要注意必须使用 FFT乘法才有较好效果.

2 二分法在指数级数

我们要计算 e 的前 n 位小数, 用级数:

e = å
k ³ 0 1
k!

设有:

Q(a,b) = (a+1)(a+2)¼b

P(a,b) = b(b-1)¼(a+2) + b(b-1)¼(a+3) +¼+ (b-1)b + b + 1

整数 P(a,b) & Q(a,b) 满足:

P(a,b)
Q(a,b) =
1
a+1 +
1
(a+1)(a+2) + ¼+
1
(a+1)(a+2)¼b

特别地, 1+P(0,K)/Q(0,K) 是e的级数前 K 项和.

二分地计算 P(a,b) & Q(a,b) , m 是a,b的均值:

m = é
ê
ë a+b
2 ù
ú
û (整数部分).

有递推关系:

P(a,b) = P(a,m)Q(m,b)+P(m,b), Q(a,b) = Q(a,m)Q(m,b).

也可以用其它方法计算m, 原则是计算简单和使得P(a,m), Q(m,b)或Q(a,m), Q(m,b)的大小相近.

最后的除法计算方法见 {高精度求倒数} {}

{程序}

3 二分法用在 arctan 级数.

我们再来看 arctan 级数的计算

arctan æ
ç
è 1
q ö
÷
ø = 1
q - 1
3q³ + 1
5q⁵ +¼ =
å
k ³ 0
(-1)^k
(2k+1)q^2k+1

比指数级数稍有复杂, 设:

R(a,b) = (2a+3)(2a+5)¼(2b+1),

Q(a,b) = (2a+3)(2a+5)¼(2b+1) q^2(b-a),

P(a,b) = (-1)^a+1

R(a,b)

2a+3

q^2(b-a-1) + (-1)^a+2

R(a,b)

2a+5

q^2(b-a-2) +¼

+ (-1)^b-1

R(a,b)

2b-1

q² + (-1)^b

R(a,b)

2b+1

有:

P(a,b)

Q(a,b)

b
å
k=a+1

(-1)^k

(2k+1)q^2(k-a)

令 a = 0, b = K, 表达式

æ
ç
è

1 +

P(0,K)

Q(0,K)

ö
÷
ø

等于计算前 n 项 arctan 级数.

递归处理(m 是 a+b 均值):

P(a,b) = P(a,m)Q(m,b)+R(a,m)P(m,b), Q(a,b) = Q(a,m)Q(m,b), R(a,b) = R(a,m)R(m,b).

4 自己动手来二分

初看上面的二分计算式会感到很奇怪: 它们是怎样得到的?

暂且不直接入主题, 而去验算上面的式子, 也许会有一点收获:

对于阶乘, 我们要验证:

P(a,b) = P(a,m)P(m,b) 和P(0,n)=n!

右边 P(a,m) * P(m,b)=(a+1)(a+2)...(m-1)m * (m+1)(m+2)...(b-1)b = P(a,b)

当 a=0, b=n 时 P(a,b)=1*2*3*...*n=n!

第二例的验证就有点复杂了. 但你若亲自验证了它, 就表明你已正式踏入了二分法世界的大门.

P(a,m)Q(m,b)+P(m,b)

=[(a+2)(a+3)...m + (a+3)(a+4)...m + ... + 1 ] * (m+1)(m+2)...b + [(m+2)(m+3)...b + (m+3)(m+4)...b + ... + 1 ]

=(a+2)(a+3)...m * (m+1)(m+2)...b + (a+3)(a+4)...m * (m+1)(m+2)...b + ... + 1

=P(a,b)

剩下的步骤就很容易了.

至于第三例 atan(1/p) 的验算是同样的方法. 但你若想象不出中间过程, 最好还是纸上写一遍, 加深理解和记忆.

回想后两例, 发现 P(a,m)Q(m,b)+P(m,b) 和 P(a,b) = P(a,m)Q(m,b)+R(a,m)P(m,b) 实际上只是提出了公因子而得到的, 这就正构造二分计算式的关键.

下面以构造函数 e^(1/x)=sum_(0->∞) 1/x^()

　------未完待续