大数阶乘算法

大家好，又见面了，我是你们的朋友全栈君。

一:精度要求较低的阶乘算法

如果只是要求算法的速度,而对精度要求比较低的话可以直接使用,斯特林公式计算n!

斯特林公式如下:

n!=sqrt(2*PI*n)*(n/e)^n*(1+1/12/n+1/288/n²–139/51840/n³-571/2488320/n⁴+…)

或

ln(n!)=0.5*ln(2*PI)+(n+0.5)*ln(n)-n+(1/12/n -1/360/n³+ 1/1260/n⁵)-…,

这里PI为圆周率,而最后一顼为雅格布·伯努力数是无穷的级数，这里我们取前5项即可得到接近16位有效数字的近似值，而精度的提高可由雅格布·伯努力数取的项数增加而得到。

另也可对斯特林公式进行修行，下面即为蔡永裕先生的《谈Stirling公式的改良》一文中得出的斯特林公式的修正版：

n! = n^n*sqrt(2*n*PI)/exp(n*(1-1/12sqrt(n)+0.4))

该公式可以得到较高的精度，较计算也较方便。

二：高精度阶乘算法

算法1:硬乘

最容易想到的自然是硬乘。模拟人工计算乘法的方法，一位位的乘。以AB*C为例，其中A，B，C各占一个基数位，以N进制为例。

第一步：

Y　=　B*C

看Y是否大于等于N，如果是的话进位。设Cn为进位量，即Cn =Y/N，而本位Y=Y%N

第二步：

X = A*C+Cn

看X是否大于等于N，如果是的话进位。而本位X=X%N

在这个过程中可以用数组来存放大数的每一个位，为了提高效率可以使用long型来存放大数的每一位，同时改10进制为100000，当然也可以更大一些但最好基数仍为了10的幕数，这样便于输出，当然取的基数更大些可以更好的发挥long型的效能，但是注意不要让X*X>=LONG_MAX,即基数最大只可为sqrt(LONG_MAX)。

第三步：

然后通过

      n! = n*(n-1)*(n-2)*(n-3)……3*2

得出结果，该方法计算10000！大概需3-4秒(CY1.7G ,512内存)

用一这种方法，效率并不好，假如已经求得n!（共有m个单元）,计算(n+1)!需要 m次乘法，m次加法(加法速度较快,可以不予考虑,下同)，m次求余（求本位)，m次除法(求进位），结果为m+1的单元计算(n+2)!需要m+1次乘法，m+1次求余，m+1次除法,结果为m+1个单元计算(n+3)!需要m+1次乘法，m+1次求余m+1次除法，结果为m+2个单元。 
      计算(n+4)!   需要   m+2次乘法，  m+2次求余   ，m+2次除法  ，结果为m+2个单元。  
      计算(n+5)!   需要   m+2次乘法，  m+2次求余   ，m+2次除法  ，结果为m+3个单元。  
     计算(n+6)!  …  
     计算(n+7)!  …  
     计算(n+8)!  …  
     计算(n+9)!  …  
     计算(n+10)!  …  
     计算(n+11)!  …  
     计算(n+12)!  …  
     计算(n+13)!  …  
      计算(n+14)!   需要  m+7次乘法，m+7次求余   ，m+7次除法  ，结果为m+7个单元  
     计算(n+15)!  需要   m+7次乘法，m+7次求余  ，m+7次除法   ，结果为m+8个单元 
     计算(n+16)!  需要   m+8次乘法，m+8次求余  ，m+8次除法   ，结果为m+8个单元 

由此可估计该算法的复杂度约为：T(n) =O(n!)，是一个NP难问题。   

算法2：分治法  

1  总体思想

假如已经求得n!(共有m个单元)，求(n+16)!

 第一步:

将n+1   与n+2   相乘,将n+3  与n+4   相乘,将n+5 与n+6…n+15与n+16,得到8个数，叫做n1,n2,n3,n4,n5,n6,n7,n8。

第二步：

n1与n2相乘,结果叫p2,结果为2个单元,需要1次乘法,1次除法,1次求余。

n3与n4相乘,结果叫p4,结果为2个单元,需要1次乘法,1次除法,1次求余。

n5与n6相乘,结果叫p6,结果为2个单元,需要1次乘法,1次除法,1次求余。

n7与n8相乘,结果叫p8,结果为2个单元,需要1次乘法,1次除法,1次求余。

第三步：

p2与p4相乘,结果叫q4,结果为4个单元,需要4次乘法,4次除法,4次求余,2次加法。

p6与p8相乘,结果叫q8,结果为4个单元,需要4次乘法,4次除法,4次求余,2次法。

第四步：

p6与p8相乘,结果叫f8,结果为8个单元,需要8次乘法,8次除法,8次求余,4次加法。这一过程的复杂度为20次乘法，20次除法，20次求余。

假定求余运算和除法运算和乘法的复杂度相同，则可知其符合分治法所需时间的计算公式,故可得:   T(n) = log(n^2)

因数学水平及时间有限不能给出算法1和算法2的精确 算法复杂度，只能给出算法1和算法2之间大致的复杂度。     

第二种算法表明，在计算阶乘时，通常的方法（先计算出n的阶乘，再用一位数乘以多位数的方法计算(n+1)的阶乘，再计算n+2的阶乘）不是最优的，更优化的算法是，计算出相邻的几个数的积（称之为部分积），用部分积乘以部分积的这种多位数乘以 多位数的算法来计算。并且根据平衡子问题的思想：在用分治法设计算法时，最好使子问题的规模大致相同。即在计算过程中为提高效率可在两相乘时，两个乘法的长度最为接近的优先进行。为此可以根据乘数的长度构造小根堆，将最短乘数和次短乘法相乘，并将结果插入小根堆中，继续运算，直到得到结果。

2  两个数相乘

设X和Y都是n（n是偶数）位的L进制的整数(当X，Y位数不等时，可在较小的数的左边补零来使其位数和较大的数相等。如X=123，Y=4325,则令X=0123）。在第一种算法中，两个大数相乘采用的是硬乘。效率较低，如果将每两个一位数的乘法或加法看作一步运算的话，那么这种方法要作O(n^2)步运算才能求出乘积XY。

这里我们用二分法来计算大数乘法。将n位的L进制的X和Y都分为2段，每段的长为n/2(如n为奇数可在左边加零)，如下图如示:

由此，X=A*L^(n/2)+B         Y=C*L^(n/2)+D

所以　X*Y = (A*L^(n/2)+B)(C*L^(n/2)+D)   = AC*L^n+(AD+BC)*L^(n/2)+BD

而AD+BC = (A-B)(D-C)+AC +BD

所以　X*Y= AC*L^n +((A-B)(D-C)+AC +BD）*L^(n/2)+BD

此时，此式仅需做3次n/2位整数的乘法(AC,BD和(A-B)(D-C)),6次加、减法和2次移位。由此可得：

Ｔ(n) = O(1)                 n= 1时

Ｔ(n) = 3T(n/2)+O(n)    n>1时

根据分治法效率计算公式可得：Ｔ（n）= O(n^log3) =O(n^1.59)

同理,若将n三等分,则

(Axx+Bx+C)*(Dxx+Ex+F)=ADxx+(AE+BD)xxx+(AF+BE+CD)xx+(BF+CE)x+CF

=ADxxxx+[(A-B)(E-D)+AD+BE]xxx+[BE+(A-C)(F-D)+CF+AD]xx+[(B-C)(F-E)+CF+BE]+CF

可知此式仅需做6次n/3位整数的乘法，

根据分治法效率计算公式可得：Ｔ（n）= O(n^log6/3) =O(n^2)

取自：http://blog.sina.com.cn/s/blog_40a3dcc1010008nf.html

相关实现：

#include<iostream>
#define MAX 1000
using namespace std;

int main(void)
{
    int n;
    while(scanf("%d",&n)==1&&n>=0)
    {
        int i,j;
        int a[MAX];      //存数运算结果
        int p,h;           //p存储当前结果的位数,h为进位
        a[1]=1;
        p=1;  
        for(i=2;i<=n;i++)   //循环与2,3,4.....n相乘
        {
            for(j=1,h=0;j<=p;j++)    //让a[]的每位与i相乘
            {
                a[j]=a[j]*i+h;
                h=a[j]/10;
                a[j]=a[j]%10;
            }
            while(h>0)         //如果h不为0
            {
                a[j]=h%10;
                h=h/10;
                j++;
            }
            p=j-1;            //将当前的位数赋给p
        }
        for(i=p;i>=2;i--)
        {
            printf("%d",a[i]);
        }
        printf("%d\n",a[i]);
    }
    return 0;
}

取自：http://www.cnblogs.com/dolphin0520/archive/2011/07/16/2108006.html

代码解读:

首先，定义两个整型的数组：
int fac[1000];//暂且先设定是1000位，我称之为“结果数组”
int add[1000];//我称之为“进位数组”

现在具体说明两个数组的作用：
1.fac[1000]
比如说，一个数5的阶乘是120，那么我就用这个数组存储它：
fac[0]=0
fac[1]=2
fac[2]=1
现在明白了数组fac的作用了吧。用这样的数组我们可以放阶乘后结果是1000位的数。

2.在介绍add[1000]之前，我介绍一下算法的思想，就以6！为例：
 从上面我们知道了5！是怎样存储的。
 就在5！的基础上来计算6！，演示如下：
fac[0]=fac[0]*6=0
fac[1]=fac[1]*6=12
fac[2]=fac[2]*6=6

3.现在就用到了我们的：“进位数组”add[1000].
 先得说明一下：add就是在第2步中用算出的结果中，第i位向第i+1位的进位数值。还是接上例：
add[0]=0;
add[1]=1;  // 计算过程：就是 (fac[1]+add[0])  %  10＝1
add[2]=0;  // 计算过程：就是 (fac[2]+add[1]) % 10=0
…….
…….

add[i+1] =( fac[i+1] + add ) % 10

现在就知道了我们的数组add[]是干什么的了吧，并且明白了add是如何求得的了吧。

4.知道了add[1000]的值，现在就在 1. 和 3 . 两步的基础上来计算最终的结果。（第2步仅作为我们理解的步骤，我们下面的计算已经包含了它）

fac[0] = ( fac[0]*6 )  mod 10 =0 
fac[1] = ( fac[1]*6 + add[0] ) mod 10 ＝2
fac[2] = ( fac[2]*6 + add[1] ) mod 10=7

到这里我们已经计算完了6！。然后就可以将数组fac[1000]中的数，以字符的形式按位输出到屏幕上了。

5.还有一点需要说明，就是我们需要一个变量来记录fac[1000]中实际用到了几位，比如5！用了前3位；
我们在这里定义为 top 
为了计算top，我们有用到了add[1000],还是以上面的为例：
 5！时，top=3,在此基础上我们来看6！时top=?
由于  add[2]=0      所以  top=3（没有变）
也就是说，如果最高位有进位时，我们的top=top+1,否则，top值不变。

6.总结一下，可以发现，我们把阶乘转化为 两个10以内的数的乘法，还有两个10以内的数的加法了。 因此，不管再大的数，基本上都能算出了，只要你的数组够大就行了。

取自：http://blog.csdn.net/jaylongli/article/details/4865035