高精度快速阶乘算法

返回索引

关键字:高精度 阶乘 算法

    我在业余时间开发了一套《超大整数完全精度快速算法库》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),
即可转化成了可充分利用高效的平方算法


    最后,我再提供大家一个确保“小整数累乘不溢出”的技巧,这是我自创的,相信会对大家有借鉴作用:
  1. 应采用“倒序”,比方说,应按 999 * 997 * ... 的顺序;
  2. 可预先设定一个数组,记录如果当前起始数组的最大值不大于某个值,则连乘多少个数不会溢出;
  3. 可利用“几何平均值≤算术平均值”定理,可推导出“ k 个正整数连乘不大于某个定值,其充分条件是其和不大于某个临界值”,我们只需要事先计算出它们的对应关系并保存下来。

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

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



备注:
n!digitsA2-1A2-0A1-1A1-0A0-0BC
10,000!35,6600.004586 s0.004682 s0.006145 s0.006254 s0.007348 s2.168x10^-19 s0.095 s
100,000!456,5740.074358 s0.081094 s0.110150 s0.121145 s0.333286 s0.047 s0.188 s
200,000!973,3510.159463 s0.179471 s0.251844 s0.280748 s1.026817 s0.156 s0.234 s
400,000!2,067,1100.372189 s0.408471 s0.580490 s0.649423 s3.142477 s0.36 s0.578 s
800,000!4,375,0400.857702 s0.934181 s1.348855 s1.526333 s9.304887 s0.875 s1.391 s
1,000,000!5,565,7091.087945 s1.208318 s1.781194 s2.023373 s13.692116 s1.265 s0.891 s
10,000,000!65,657,06014.922132 s17.024912 s24.660277 s28.985707 s-----------20.141 s57.687 s
20,000,000!137,334,71530.775132 s35.551279 s51.123602 s61.061514 s-----------46.484 s--------
40,000,000!286,710,62567.570732 s78.600010 s112.555688 s136.195050 s-----------108.891 s--------

CPU: AMD Athlon 64 X2 Dual Core Processor 4800+, 2.512GHz(201MHz x12.5), L1 Cache 64KB, L2 Cache 512KB (参见:CPU-Z 1.42 报告)
OS : Windows XP SP2
RAM: 2GB DDR2 - 800MHz



  A -- HugeCalc V8.0.0.0
     A2-1 --> HugeCalc.ini: NumOfCores = 2; SSE2Support = 1; (测试双核、SSE2指令加速)
     A2-0 --> HugeCalc.ini: NumOfCores = 2; SSE2Support = 0; (测试双核、ALU 指令加速)
     A1-1 --> HugeCalc.ini: NumOfCores = 1; SSE2Support = 1; (测试单核、SSE2指令加速)
     A1-0 --> HugeCalc.ini: NumOfCores = 1; SSE2Support = 0; (测试单核、ALU 指令加速)
     A0-0 --> HugeCalc.ini: NumOfCores = 0; SSE2Support = 0; (测试无高级算法参与)
  B -- Mathematica V6.0.1.0
  C -- Maple V11.0 ( February 17 2007, Build ID 277223 )


  1. 作者:郭先强初稿:2004-06-14;本站发布日期:2004-06-15;
  2. 本文的初稿发表于著名的“CSDN - 技术社区 - 专题开发 数据结构与算法问题”版,原稿请参见本站的原版收藏贴(原贴已不幸被 csdn “私吞”);
  3. 相关下载:超大整数完全精度快速计算器/算法库
    超大整数高精度快速计算器。
  4. 版权所有,未经原作者授权,严禁转载!

返回索引