大数阶乘

最新推荐文章于 2024-09-25 10:36:42 发布

一直活在明天

最新推荐文章于 2024-09-25 10:36:42 发布

阅读量2.4k

点赞数

文章标签：算法 fft 优化 c# delete n2

大数阶乘算法

2008-07-14 13:52

*************************************(1)****************************************************

      假如需要计算n+16的阶乘，n+16接近10000，已经求得n!（共有m个单元），（每个单元用一个long数表示，表示1-100000000）

第一种算法（传统算法）
计算(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个单元
该算法的复杂度：共需:m+(m+8)+(m+1+m+7)*7=16m+64次乘法，16m+64次求余，16m+64次除法

第二种算法：
1.将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
2.   n1   与   n2相乘，结果叫做p2，结果为2个单元，需要1次乘法。
        p2   与   n3相乘，结果叫做p3，需要2次乘法，1次加法(加法速度快，下面省略），2次除法，2次求余
        p3   与   n4相乘，结果叫做p4，需要3次乘法，3次除法，3次求余
        p4   与   n4相乘，结果叫做p5，需要4次乘法，4次除法，4次求余
        p5   与   n6相乘，结果叫做p6，需要5次乘法，5次除法，5次求余
        p6   与   n7相乘，结果叫做p7，需要6次乘法，6次除法，6次求余
        p7   与   n8相乘，结果叫做p8，结果为8个单元，需要7次乘法，7次除法，7次求余
        这一过程的复杂度为（1+2+3+...+7)共计28次乘法，28次求余，28次除法。
3.将   n!（m个单元)   与p8(8个单元)相乘，需要m*8次乘法，m次除法，m次求余，（注意不是m*8次求余，m*8次除法,原因请大家思考)
该算法的复杂度为   8m次乘法，m次求余，m次除法。
假定求余运算和除法运算和乘法的复杂度相同，则第第一种算法需要48m+192次乘法，
而第2种运算需要10m+84次乘法，当m很大时，第二种算法的速度是第一种算法的4.8倍。

      第二种算法表明，在计算阶乘时，通常的方法（先计算出n的阶乘，再用一位数乘以多位数的方法计算(n+1)的阶乘，再计算n+2的阶乘）不是最优的，更优化的算法是，计算出相邻的几个数的积（称之为部分积），用   部分积   乘以   部分积   的这种多位数   乘以   多位数的算法来计算。如果两个多位数均为n位，使用FFT算法可以将乘法运算的速度   由   n*n   提高   n*log2(n),当n很大时，FFT算法的速度大大快与常规算法。

#include   <math.h>
#include   <memory.h>
#include   <stdio.h>

#define   PI     3.1415926535897932384626433832795
#define   E 2.7182818284590452353602874713527
#define   TEN9     1000000000

typedef   unsigned   __int64   UINT64;
typedef   unsigned   long   DWORD;

int   calcResultLen(DWORD   n)
{
    return   log10(2*PI*(double)n)/2   +   n   *     log10((double)n/E)+2;
}
void   printfResult(DWORD   *buff,int   buffLen)
{
    DWORD   *p=buff;
    while   (   *p==0) p++;
    printf("/nn!=/n");
    for   (   int   i=0;i<buffLen;i++)
    {
    if   (i==0)
    printf("%9d ",p[i]);
    else
    printf("%09d ",p[i]);
    if   (i   %   8==7)
    printf("/n");
    }
    printf("/n");
}
void   calcFac(DWORD   *buff,   int   buffLen,DWORD   n)
{
    int     len=1;
    DWORD   carry,*p,*pEnd,*pBeg;
    unsigned   __int64   prod;

    memset(buff,0,   sizeof(DWORD)*buffLen);
    pBeg=buff+buffLen-1;
    *pBeg=1;
    for   (DWORD   i=2;i<=n;i++)
    {
    for   (p=pBeg,pEnd=p-len,carry=0;p>pEnd;p--)
    {
        prod=   (UINT64)(*p)   *   (UINT64)i   +   (UINT64)carry;
        *p=       prod   %   TEN9;
        carry=prod   /   TEN9;
    }
    if   (carry>0)
    {
        *(pBeg-len)=carry;
        len++;
    }
    }
}
int   main(int   argc,   char*   argv[])
{
    DWORD   n=1;
    printf("n=?");
    scanf("%d",&n);
    int   buffLen=(calcResultLen(n)+8)   /   9;
    DWORD   *buff=new   DWORD[buffLen];
    calcFac(buff,   buffLen,   n);
    printfResult(buff,buffLen);
    delete[]   buff;
    return   0;
}

**************************************(2)***************************************************
用数组的方法解决大数、巨数的阶乘结果越界的问题。
具体算法中有最朴实的乘法运算思想，请各位细细体味。
#include <stdio.h>
int main()
{
    int n; //阶乘大小
    printf("请输入n的大小:");
    scanf("%d",&n); //从键盘接收阶乘大小
    int a[200]; //确保保存最终运算结果的数组足够大
    int carry; //进位
    int digit = 1; //位数
    a[0] = 1; //将结果先初始化为1
    int temp; //阶乘的任一元素与临时结果的某位的乘积结果

    for(int i = 2; i <= n; ++i) //开始阶乘，阶乘元素从2开始依次“登场”
    {
        //按最基本的乘法运算思想来考虑，将临时结果的每位与阶乘元素相乘
        for(int j = 1, carry = 0; j <= digit; ++j)
        {
             //相应阶乘中的一项与当前所得临时结果的某位相乘（加上进位）
            temp = a[j-1] * i + carry;
            a[j-1] = temp % 10; //更新临时结果的位上信息
            carry = temp / 10; //看是否有进位
        }
        while(carry) //如果有进位
        {
            a[++digit-1] = carry % 10; //新加一位，添加信息。位数增1
            carry /= 10; //看还能不能进位
        }
    }
    printf("结果是:/n%d ! = ",n); //显示结果
    for(int i = digit; i >=1; --i)
    {
        printf("%d",a[i-1]);
    }
    return 0;
}

    阶乘，是求一组数列的乘积，其效率的高低，一、是取决于高精度乘法算法，二、是针对阶乘自身特点算法的优化。
    前者，是一种广为习知的技术，虽然不同算法效率可天壤之别，但终究可以通过学习获得，甚至可用鲁迅先生的“拿来主义”；
    后者，则有许多的发展空间，这里我将介绍一种由我完全独创的一种方法，欢迎大家讨论，如能起到“抛砖引玉”的效果，那就真正达到了目的。

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

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

    言归正转，下面以精确计算 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) = 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)，
即可转化成了可充分利用高效的平方算法。

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

阶乘的定义

阶乘是数学中的一个术语。对于一个非负整数n，n的阶乘指的是所有小于等于n的正整数的乘积，记为n!。例如，

符号n!是由Christian Kramp(1760 – 1826)于1808年引入的。

阶乘的严格定义为：

并且0！＝1，因为阶乘是针对所有的非负整数。
后者基于一个事实：0个数的乘积为1。这个是很有用的：
递归关系

适用于n = 0的情况；

这个定义使得组合学（combinatorics）中许多包含0的计算能够有效。

阶乘的概念相当简单、直接，但它的应用很广泛。在排列、组合、微积分（如泰勒级数）、概率论中都有它的身影。

但我这里最想说的是（与本文主题相关），在计算机科学的教学中，阶乘与斐波那契数列一道经常被选为递归算法的素材，因为阶乘满足下面的递归关系（如果n ≥ 1）：

言归正传

下面来考虑如何在程序中计算阶乘。根据阶乘的定义和它满足的递归关系，我们很容易得到这样的算法：

C# code

      
      
       
       
       
       public
       
        
       
       static
       
        
       
       long
       
        Calculate(
       
       int
       
        n)
{
    
       
       if
       
        (n 
       
       <
       
        
       
       0
       
       ) { 
       
       throw
       
        
       
       new
       
        ArgumentOutOfRangeException(
       
       "
       
       n必须为非负数。
       
       "
       
       ); }
    
       
       if
       
        (n 
       
       ==
       
        
       
       0
       
       ) { 
       
       return
       
        
       
       1
       
       ; }
    
       
       return
       
        n 
       
       *
       
        Calculate(n 
       
       -
       
        
       
       1
       
       );
}

随着n的增大，n!会迅速增大，其速度可能会超出你的想象。如果n不大，这种算法还可以，但对long类型来说，很快就会溢出。对计算器来说，大多数可以计算到69!，因为70! > 10^100。

上面这种累积相乘的算法的主要问题在于普通类型所容纳的数值太小，即使是double类型，它的最大值不过是1.79769313486232e308，即一个309位的数字。

我们来考虑另外一种方案，将乘积的每一位数字都存放在数组中，这样的话一个长度为10000的数组可以存放任何一个10000位以内的数字。

假设数组为uint[] array = new uint[10000]，因为1! = 1，所以首先置a[0] = 1，分别乘以2、3，得到3! = 6，此时仍只需要一个元素a[0]；然后乘以4得到24，我们把个位数4放在a[0]，十位数2放在a[1]，这样存放结果就需要两个元素；乘以5的时候，我们可以这样进行：用5与各元素由低到高逐一相乘，先计算个位数（a[0]）4 × 5，结果为20，这样将a[0]置为0，注意要将2进到十位数，然后计算原来的十位数（a[1]）2 × 5，结果为10加上刚才进的2 为12，这样十位数是2，而1则进到百位，这样就得到5! = 120；以此类推……

下面给出上述方案的一个实现：

C# code

      
      
       
       
       
       public
       
        
       
       static
       
        
       
       uint
       
       [] CalculateLargeNumber(
       
       int
       
        n)
{
    
       
       if
       
        (n 
       
       <
       
        
       
       0
       
       ) { 
       
       throw
       
        
       
       new
       
        ArgumentOutOfRangeException(
       
       "
       
       n必须为非负数。
       
       "
       
       ); }

    
       
       if
       
        (n 
       
       ==
       
        
       
       0
       
        
       
       ||
       
        n 
       
       ==
       
        
       
       1
       
       ) { 
       
       return
       
        
       
       new
       
        
       
       uint
       
       [] { 
       
       1
       
        }; }
    
       
       //
       
        数组的最大长度
       
       

       
           
       
       const
       
        
       
       int
       
        MaxLength 
       
       =
       
        
       
       100000
       
       ;
    
       
       uint
       
       [] array 
       
       =
       
        
       
       new
       
        
       
       uint
       
       [MaxLength];
    
       
       //
       
        1! = 1
       
       

       
           array[
       
       0
       
       ] 
       
       =
       
        
       
       1
       
       ;
    
       
       int
       
        i 
       
       =
       
        
       
       0
       
       ;
    
       
       int
       
        j 
       
       =
       
        
       
       0
       
       ;
    
       
       //
       
        valid为当前阶乘值的位数（如5! = 120，此时valid = 3）
       
       

       
           
       
       int
       
        valid 
       
       =
       
        
       
       1
       
       ;

    
       
       for
       
        (i 
       
       =
       
        
       
       2
       
       ; i 
       
       <=
       
        n; i
       
       ++
       
       )
    {
        
       
       long
       
        carry 
       
       =
       
        
       
       0
       
       ;
        
       
       for
       
        (j 
       
       =
       
        
       
       0
       
       ; j 
       
       <
       
        valid; j
       
       ++
       
       )
        {
            
       
       long
       
        multipleResult 
       
       =
       
        array[j] 
       
       *
       
        i 
       
       +
       
        carry;
            
       
       //
       
        计算当前位的数值
       
       

       
                   array[j] 
       
       =
       
        (
       
       uint
       
       )(multipleResult 
       
       %
       
        
       
       10
       
       );
            
       
       //
       
        计算进到高位的数值
       
       

       
                   carry 
       
       =
       
        multipleResult 
       
       /
       
        
       
       10
       
       ;
        }
        
       
       //
       
        为更高位赋值
       
       

       
               
       
       while
       
        (carry 
       
       !=
       
        
       
       0
       
       )
        {
            array[valid
       
       ++
       
       ] 
       
       =
       
        (
       
       uint
       
       )(carry 
       
       %
       
        
       
       10
       
       );
            carry 
       
       /=
       
        
       
       10
       
       ;
        }
    }
    
       
       //
       
        截取有效元素
       
       

       
           
       
       uint
       
       [] result 
       
       =
       
        
       
       new
       
        
       
       uint
       
       [valid];
    Array.Copy(array, result, valid);
    
       
       return
       
        result;
}