高精度快速阶乘算法

大家好，又见面了，我是你们的朋友全栈君。
我在业余时间开发了一套《超大整数完全精度快速算法库》HugeCalc，可快速计算超大整数的加、减、乘、除（商/余）、乘方、开方，也可快速计算大数的 Fibonacci 数列、（双）阶乘、排列、组合等，还可完成超大整数数组的最大公约数、最小公倍数等数论运算，现在，该套软件已被华军、天空、电脑之家、天天等下载站点收录。

自在网上公开以来，广受网友关注，经常有网友来联系，想交流一些算法心得。其中涉及最多的是关于“阶乘”算法，部分是在校大学生，也许是他们的毕业设计？:)

这里，我就把关于该算法模块的核心部分，也就是一些关键点，整理出来，以供大家参考。

    阶乘，是求一组数列的乘积，其效率的高低，一、是取决于高精度乘法算法，二、是针对阶乘自身特点算法的优化。

    前者，是一种广为习知的技术，虽然不同算法效率可天壤之别，但终究可以通过学习获得，甚至可用鲁迅先生的“拿来主义”；

    后者，则有许多的发展空间，这里我将介绍一种由我完全独创的一种方法，欢迎大家讨论，如能起到“抛砖引玉”的效果，那就真正达到了目的。

我在开发“阶乘”类算法时，始终遵循如下原则：

参与高精度乘法算法的两数，大小应尽可能地相近；
尽可能将乘法转化为乘方；
尽可能地优先调用平方；

言归正转，下面以精确计算 1000! 为例，阐述该算法：

    记 F1(n) = n*(n-1)*…*1；

       F2(n) = n*(n-2)*…*(2 or 1)；

       Pow2(r) = 2^r；

    有 F1(2k+1) = F2(2k+1) * F2(2k)

           = Pow2(k) * F2(2k+1) * F1(k),

       F1(2k) = Pow2(k) * F2(2k-1) * F1(k),

    及 Pow2(u) * Pow2(v) = Pow2(u+v),

∴ 1000! =
F1(1000)

         = Pow2(500)*F2(999)*
F1(500)

         = Pow2(750)*F2(999)*F2(499)*
F1(250)

         = Pow2(875)*F2(999)*F2(499)*F2(249)*
F1(125)

         = Pow2(937)*F2(999)*F2(499)*F2(249)*F2(125)*
F1(62)

         = Pow2(968)*F2(999)*F2(499)*F2(249)*F2(125)*F2(61)*
F1(31)

         = Pow2(983)*F2(999)*F2(499)*F2(249)*F2(125)*F2(61)*F2(31)*
F1(15)

         = …

如果你预存了某些小整数阶乘（比如这里的“
F1(15)”），则可提前终止分解，否则直至右边最后一项为“
F1(1)”为止；这样，我们
将阶乘转化为2的整数次幂与一些连续奇数的积（或再乘以一个小整数的阶乘）；

再定义：F2(e,f) = e*(e-2)*…*(f+2)，这里仍用“F2”，就当是“函数重载”好了:)，

则 F2(e) = F2(e,-1) = F2(e,f)*F2(f,-1) （e、f为奇数，0≤f≤e）

∴ F2(999) = F2(999,499)*F2(499,249)*F2(249,125)*F2(125,61)*F2(61,31)*F2(31),

   F2(499) = ____________F2(499,249)*F2(249,125)*F2(125,61)*F2(61,31)*F2(31),

   F2(249) = ________________________F2(249,125)*F2(125,61)*F2(61,31)*F2(31),

   F2(125) = ____________________________________F2(125,61)*F2(61,31)*F2(31),

   F2( 61) = _______________________________________________F2(61,31)*F2(31),

   F2( 31) = _________________________________________________________F2(31),

∴ F1(1000) =
F1(15) * Pow2(983) * F2(999,499) /

                 * [F2(499,249)^2] * [F2(249,125)^3] /

                 * [F2(61,31)^4] * [F2(31)^5]

    这样，我们又
将阶乘转化为了乘方运算。

    上式实际上是个形如 a * b^2 * c^3 * d^4 * e^5 的式子；我们再将指数转化为二进制，可得到如下公式：

        a * b^2 * c^3 * d^4 * e^5 = (a*c*e)*[(b*c)^2]*[(d*e)^4]

                                  =
(((e*d)^2)*(c*b))^2*(e*c*a)，

即可
转化成了可充分利用高效的平方算法。

最后，我再提供大家一个
确保“小整数累乘不溢出”的技巧，这是我自创的，相信会对大家有借鉴作用：

应采用“倒序”，比方说，应按 999 * 997 * … 的顺序；
可预先设定一个数组，记录如果当前起始数组的最大值不大于某个值，则连乘多少个数不会溢出；
可利用“几何平均值≤算术平均值”定理，可推导出“ k 个自然数连乘不大于某个定值，其充分条件是其和不大于某个临界值”，我们只需要事先计算出它们的对应关系并保存下来。

上述算法是
HugeCalc Ver1.2.0.1 的算法关键点，其效率已略高于liangbch(宝宝)的高级算法3.0版。

在最新的
HugeCalc Ver2.1.0.1 中，采用的是更彻底的分解成质因数的方法，技巧性倒反不如上面的强（但却需要有高效的素数筛法），但里面的很多思想仍在延续。

备注：

Factorial.exe Ver2.1.0.1	Celeron(tm) 466MHz 64MB RAM, Win98Sec	Pentium(R) 4 1.70GHz 256MB RAM, WinXP	Result
10,000!	0.069s	0.031s	2.8462… x 10^35,659
20,000!	0.236s	0.109s	1.8192… x 10^77,337
40,000!	0.795s	0.390s	2.0916… x 10^166,713
80,000!	2.661s	1.328s	3.0977… x 10^357,506
100,000!	4.177s	1.969s	2.8242… x 10^456,573
200,000!	13.663s	6.438s	1.4202… x 10^973,350
400,000!	43.818s	20.828s	2.5344… x 10^2,067,109
800,000!	139.337s	66.921s	5.6846… x 10^4,375,039