WV.29-大数阶乘算法9-求N!的高精度算法

求N!的高精度算法

?

本文是张一飞2001年写的论文，原文可从http://oibh.kuye.cn/download/thesis/thesis2001_zhangyifei.zip处下载

?求N!的高精度算法

本文中的算法主要针对Pascal语言

这篇文章的内容

你了解高精度吗？

你曾经使用过哪些数据结构？

你仔细思考过如何优化算法吗？

在这里，你将看到怎样成倍提速求N!的高精度算法

Pascal中的标准整数类型

Comp虽然属于实型，实际上是一个64位的整数

高精度算法的基本思想

Pascal中的标准整数类型最多只能处理在-263～263之间的整数。如果要支持更大的整数运算，就需要使用高精度

高精度算法的基本思想，就是将无法直接处理的大整数，分割成若干可以直接处理的小整数段，把对大整数的处理转化为对这些小整数段的处理

数据结构的选择

每个小整数段保留尽量多的位

使用Comp类型

采用二进制表示法

每个小整数段保留尽量多的位

一个例子：计算两个15位数的和

?方法一

?分为15个小整数段，每段都是1位数，需要15次1位数加法

?方法二

?分为5个小整数段，每段都是3位数，需要5次3位数加法

?方法三

?Comp类型可以直接处理15位的整数，故1次加法就可以了

?比较

?用Integer计算1位数的加法和3位数的加法是一样快的

?故方法二比方法一效率高

?虽然对Comp的操作要比Integer慢，但加法次数却大大减少

?实践证明，方法三比方法二更快

使用Comp类型

高精度运算中，每个小整数段可以用Comp类型表示

Comp有效数位为19～20位

求两个高精度数的和，每个整数段可以保留17位

求高精度数与不超过m位整数的积，每个整数段可以保留18–m位

求两个高精度数的积，每个整数段可以保留9位

如果每个小整数段保留k位十进制数，实际上可以认为其只保存了1位10^k进制数，简称为高进制数，称1位高进制数为单精度数

采用二进制表示法

采用二进制表示，运算过程中时空效率都会有所提高，但题目一般需要以十进制输出结果，所以还要一个很耗时的进制转换过程。因此这种方法竞赛中一般不采用，也不在本文讨论之列.

算法的优化

高精度乘法的复杂度分析

连乘的复杂度分析

设置缓存

分解质因数求阶乘

二分法求乘幂

分解质因数后的调整

高精度乘法的复杂度分析

计算n位高进制数与m位高进制数的积

?需要n*m次乘法

?积可能是n+m–1或n+m位高进制数

连乘的复杂度分析（1）

一个例子：计算5*6*7*8

?方法一：顺序连乘

?5*6=30，1*1=1次乘法

?30*7=210，2*1=2次乘法

?210*8=1680，3*1=3次乘法 ?共6次乘法

?方法二：非顺序连乘

?5*6=30，1*1=1次乘法

?7*8=56，1*1= 1次乘法

?30*56=1680，2*2=4次乘法???共6次乘法

连乘的复杂度分析（2）

若“n位数*m位数=n+m位数”，则n个单精度数，无论以何种顺序相乘，乘法次数一定为n(n-1)/2次

?证明：

?设F(n)表示乘法次数，则F(1)=0，满足题设

?设k<n时，F(k)=k(k-1)/2，现在计算F(n)

?设最后一次乘法计算为“k位数*(n-k)位数”，则

?F(n)=F(k)+F(n-k)+k (n-k)=n(n-1)/2（与k的选择无关）

设置缓存（1）

一个例子：计算9*8*3*2

?方法一：顺序连乘

?9*8=72，1*1=1次乘法

?72*3=216，2*1=2次乘法

?216*2=432，3*1=3次乘法

?方法二：非顺序连乘

?9*8=72，1*1=1次乘法

?3*2=6，1*1=1次乘法

?72*6=432，2*1=2次乘法

共6次乘法

?

特点：n位数*m位数可能是n+m-1位数

设置缓存（2）

考虑k+t个单精度数相乘a1*a2*…*ak*ak+1*…*ak+t

?设a1*a2*…*ak结果为m位高进制数（假设已经算出）

?ak+1*…*ak+t结果为1位高进制数

?若顺序相乘，需要t次“m位数*1位数”，共mt次乘法

?可以先计算ak+1*…*ak+t，再一起乘，只需要m+t次乘法

在设置了缓存的前提下，计算m个单精度数的积，如果结果为n位数，则乘法次数约为n(n–1)/2次，与m关系不大?

–设S=a1a2… am，S是n位高进制数
–可以把乘法的过程近似看做，先将这m个数分为n组，每组的积仍然是一个单精度数，最后计算后面这n个数的积。时间主要集中在求最后n个数的积上，这时基本上满足“n位数*m位数=n+m位数”，故乘法次数可近似的看做n(n-1)/2次?

设置缓存（3）

缓存的大小

?设所选标准数据类型最大可以直接处理t位十进制数

?设缓存为k位十进制数，每个小整数段保存t–k位十进制数

?设最后结果为n位十进制数，则乘法次数约为

?k/(n–k)?∑_(i=1..n/k)i=(n+k)n/(2k(t–k))，其中k远小于n

?要乘法次数最少，只需k (t–k)最大，这时k=t/2

?因此，缓存的大小与每个小整数段大小一样时，效率最高

?故在一般的连乘运算中，可以用Comp作为基本整数类型，每个小整数段为9位十进制数，缓存也是9位十进制数

分解质因数求阶乘

例：10!=2⁸*3⁴*5²*7

?n!分解质因数的复杂度远小于nlogn，可以忽略不计

?与普通算法相比，分解质因数后，虽然因子个数m变多了，但结果的位数n没有变，只要使用了缓存，乘法次数还是约为n(n-1)/2次

?因此，分解质因数不会变慢（这也可以通过实践来说明）

?分解质因数之后，出现了大量求乘幂的运算，我们可以优化求乘幂的算法。这样，分解质因数的好处就体现出来了

二分法求乘幂

二分法求乘幂，即：

?a²ⁿ⁺¹=a²ⁿ*a

?a²ⁿ=(aⁿ)²

?其中，a是单精度数

复杂度分析

?假定n位数与m位数的积是n+m位数

?设用二分法计算aⁿ需要F(n)次乘法

?F(2n)=F(n)+n²，F(1)=0

?设n=2^k，则有F(n)=F(2^k)=∑_(i=0..k–1)4ⁱ=(4^k–1)/3=(n²–1)/3

与连乘的比较

?用连乘需要n(n-1)/2次乘法，二分法需要(n²–1)/3

?连乘比二分法耗时仅多50%

?采用二分法，复杂度没有从n²降到nlogn

二分法求乘幂之优化平方算法

怎样优化

?(a+b)²=a²+2ab+b²

?例：12345²=123²*10000+45²+2*123*45*100

?把一个n位数分为一个t位数和一个n-t位数，再求平方

怎样分

?设求n位数的平方需要F(n)次乘法

?F(n)=F(t)+F(n-t)+t(n-t)，F(1)=1

?用数学归纳法，可证明F(n)恒等于n(n+1)/2

?所以，无论怎样分，效率都是一样

?将n位数分为一个1位数和n–1位数，这样处理比较方便

二分法求乘幂之复杂度分析

复杂度分析

?前面已经求出F(n)=n(n+1)/2，下面换一个角度来处理

?S²=(∑_(0≤i<n)a_i10ⁱ)²=∑_(0≤i<n)a_i²10²ⁱ+2∑_(0≤i<j<n)a_ia_j10^i+j

?一共做了n+C(n,2)=n(n+1)/2次乘法运算

?普通算法需要n²次乘法，比改进后的慢1倍

改进求乘幂的算法

?如果在用改进后的方法求平方，则用二分法求乘幂，需要(n+4)(n–1)/6次乘法，约是连乘算法n(n–1)/2的三分之一

分解质因数后的调整（1）

为什么要调整

?计算S=2¹¹3¹⁰，可以先算2¹¹，再算3¹⁰，最后求它们的积

?也可以根据S=2¹¹3¹⁰=6¹⁰*2，先算6¹⁰，再乘以?2即可

?两种算法的效率是不同的

分解质因数后的调整（2）

什么时候调整

?计算S=a^x+kb^x=(ab)^xa^k

?当k<xlog_ab时，采用(ab)^xa^k比较好，否则采用a^x+kb^x更快

?证明：

?可以先计算两种算法的乘法次数，再解不等式，就可以得到结论

?也可以换一个角度来分析。其实，两种算法主要差别在最后一步求积上。由于两种方法，积的位数都是一样的，所以两个因数的差越大，乘法次数就越小

?∴当a^xb^x–a^k>a^x+k–b^x时，选用(ab)^xa^k，反之，则采用a^x+kb^x。

?∴a^xb^x–a^k>a^x+k–b^x

?∴(b^x–a^k)(a^x+1)>0

?∴b^x>a^k

?这时k<xlog_ab

总结

内容小结

?用Comp作为每个小整数段的基本整数类型

?采用二进制优化算法

?高精度连乘时缓存和缓存的设置

?改进的求平方算法

?二分法求乘幂

?分解质因数法求阶乘以及分解质因数后的调整

应用

?高精度求乘幂（平方）

?高精度求连乘（阶乘）

?高精度求排列组合

结束语

求N!的高精度算法本身并不难，但我们仍然可以从多种角度对它进行优化。

其实，很多经典算法都有优化的余地。我们自己编写的一些程序也不例外。只要用心去优化，说不准你就想出更好的算法来了。

也许你认为本文中的优化毫无价值。确实是这样，竞赛中对高精度的要求很低，根本不需要优化。而我以高精度算法为例，不过想谈谈如何优化一个算法。我想说明的只有一点：算法是可以优化的。