计算阶乘的 18 位非零尾数

阶乘不用解释吧，阶乘的 18 位非零尾数需要简单解释一下。阶乘的十进制结果中末尾都会有很多个数字 0 (5以下的除外)，去掉这些连续的 0 之后最后 18 位数字即是这里所说的 18 位非零尾数。之所以称这为“非零”，是因为最后一个数字肯定不是 0。举两个例子，10! = 3,628,800，其 18 位非零尾数为 36,288；24! = 620,448,401,733,239,439,360,000，其 18 位非零尾数为 044,840,173,323,943,936

分析这个问题,可以发现主要在于计算
1*2*3*4*6*7*8*9*11*12*... (mod 5^18)
也就是去掉所有5的倍数后的乘积.

分析方法如下:
由于对于充分大的N,N!中因子2的数目要远远大于因子5的数目,我们可以认为对于充分大的N,2的因子数目
减去5的因子数目不小于18了,也就是,去掉N!末尾所有的0后面的数我们认为还是被2^18整除,
实际上,对于28!,
2的因子数目为:
[28/2]+[28/4]+[28/8]+[28/16]=25
5的因子数目为:
[28/5]+[28/25]=6
所以28!的2的因子数目同5的因子差值已经不小于18了,对于更大的N,这个差值也不会变小.
所以我们只需要对小于28的整数特殊处理(直接相乘就可以了),其余的,都认为N!去掉末尾的0后,
还是被2^18整除,所以我们只要计算出这个数对于5^18的余数,就可以通过中国剩余定理计算出关于10^18的余数了.

对于N!,我们可以将5的倍数和不是5的倍数的数分开处理:
所有的5的倍数构成5^[N/5]*[N/5]!
如果我们已经得到所有不是5的倍数的乘积关于5^18的余数,记为g(N)
所求的结果记为f(N),那么就有f(N)=g(N)*f([N/5]) (mod 5^18)
所以关键是求g(N),也就是1到N之间所有不是5的倍数的数的乘积关于5^18的余数
也就是最前面列出的.
后面只讨论如何计算:
g(N)=1*2*3*4*6*7*8*9*11*12*...*N (mod 5^18)

还可以扩展这个问题到一般情况,记f(N,L)是N!去掉末尾的0后最后L位，而g(N,L)是1,2,...,N中去掉5的倍数后所有数字乘积的最后L位。
对于T位数N,计算N!去掉末尾的0的最后L位的算法f(N,L)的时间复杂度可以达到:
O(L^3 log(L)^2 + L^3 Log(L) *T+T*log(T))
使用的空间复杂度为O(L^3+L^2 log(L)+T)
计算过程中,我们总假设一个长度为L的数占用的空间为O(L),两个长度为L的数相乘花费的时间为O(L*log(L)),而相加的时间复杂度为O(L).

方法如下:
i)计算杨辉三角形 关于5^L的余数,并保存,记为C(U,V) (mod 5^L) 其中1<=U<=L, 0<=V<=U.
空间复杂度为O(L^3) (L^2个长度为L的整数),需要花费O(L^2)次加法,每次加法需要O(L)的时间,所以时间复杂度也是O(L^3)
ii)记F(k,x)=(x*5^k+1)(x*5^k+2)(x*5^k+3)(x*5^k+4)(x*5^k+6)...(x*5^k+5^k-1) (mod 5^L)
其中连乘中所有5的倍数全部被去掉
那么F(1,x)=24+50*5*x+35*5^2*x^2+10*5^3*x^3+5^4*x^4 (mod 5^L)
假设F(n,x)=a0+a1*5*x+a2*5^2*x^2+....+a(L-1)*5^(L-1)*x^(L-1) (mod 5^L)
那么F(n+1,x)=F(n,5*x)*F(n,5*x+1)*F(n,5*x+2)*F(n,5*x+3)*F(n,5*x+4) (mod 5^L)
首先将每个F(n,5*x+t)展开,重新整理成关于(5x)的多项式,这个需要L^2次乘法运算,使用到i)中的杨辉三角形.花费时间复杂度为O(L^3*log(L))
然后将4个多项式俩俩相乘,(5x)次数大于L的都可以抛弃,每次两个多项式相乘最多花费O(L^2)次乘法(如果使用快速乘法可以只用O(L log(L))次乘法),所以花费的时间复杂度还是O(L^3*log(L))
这样,通过O(L^3 log(L)^2),我们可以计算出F(k,x),其中1<=k<=L
而保存所有的F(k,x),需要花费的空间是O(L^2 log(L)) (L *log(L)个长度为L的数)
而对于F(k,x), k>L,我们知道F(k,x)=1,不需要再计算了.
iii)对于f(N,L),f(N,L) (mod 5^L)=g(N,L)*f([N/5],L) (mod 5^L)
而g(N,L) (mod 5^L)我们可以将它划分成若干个长度不一的类似ii)中的连乘式,其中长度不超过L的有L个，长度大于L的由于对应F(k,x)总是1,可以统一处理。对于每个长度不超过L的连乘式,我们代入公式ii)后,
需要L次乘法,
所以计算g(N,L)(mod 5^L)共花费时间为O(L^3 log(L))
v) 在计算完g(N,L) (mod 5^L)以后,我们还需要计算f([N/5],L),同样的,如果我们把N看成T位5进制数(事先转化一个数为5进制只需要Tlog(T)的时间),那么N/5是T-1位5进制数,所以通过同样算法,可以再花费O(L^3 log(L))的时间递归到T-2位数,...,这样总共经过T步后就会得出最终结果.
这个递归过程中总共花费时间最多为 O(T* L^3 log(L))
而保存原式输入数据X (以及中间数据X/5, X/25,...等等)需要一个长度为T的空间,需要O(T)的空间
所以总共需要时间复杂度
O(L^3 log(L)^2 + L^3 Log(L) T + T log(T))
使用的空间复杂度为O(L^3+L^2 log(L)+T)

上面过程只是算出f(N,L) (mod 5^L),而对于充分大的N(超过10L/3肯定够了,而对于这么小的N,直接计算乘积就可以了),f(N,L) (mod 2^L)总是0,所以通过中国剩余定理就可以计算出f(N,L) (mod 10^L)
对于现在的计算机基本上在L<=100而且T<=1000 (也就是N<=10^1000)时不会有问题.
下面程序使用了GMP库,
但是由于下面的程序中，对于输入的T位整数N,没有通过事先转化为5进制的方法，所以实际上花费的时间复杂度为 O(L^3 log(L)^2 +L^3 log(L)T +T^2)
代码以及预编译过的可执行程序可以在下面链接下载到:

#include
#include
#include

#define TRIANGLE(x,y)  triangle[(x)*L+(y)]
#define F(x,y)         fun[(x)*L+(y)]
void PolyMul(mpz_t *src1, mpz_t *src2, mpz_t *dst, int L, mpz_t mod)
{
    int i,j;mpz_t tmp,tmpp;
    mpz_init(tmp);mpz_init(tmpp);
    for(i=0;i mpz_set_ui(dst[i],0);
    for(i=0;i for(j=0;j     mpz_mul(tmp,src1[i],src2[j]);
     mpz_add(tmpp,tmp,dst[i+j]);
     mpz_mod(dst[i+j],tmpp,mod);
 }
    }
    mpz_clear(tmp);mpz_clear(tmpp);
}
void Ft4(mpz_t *fun,//Use Polynomial F(s,*)
  mpz_t t4,  //F(s,*)[t4] (mod five_L)
  int s,
  mpz_t five_L,
  int L,
  mpz_t result)
{
    int i;mpz_t tmp1,tmp2;mpz_init(tmp1);mpz_init(tmp2);
    mpz_mul(tmp1,F(s,L-1),t4);
    mpz_mod(result,tmp1,five_L);
    for(i=L-2;i>=1;i--){
 mpz_add(tmp1,result,F(s,i));
 mpz_mul(tmp2,tmp1,t4);
 mpz_mod(result,tmp2,five_L);
    }
    mpz_add(tmp1,result,F(s,0));
    mpz_mod(result,tmp1,five_L);
    mpz_clear(tmp1);mpz_clear(tmp2);
}

void Get_Mult_No_5(mpz_t N,
  int L,
  mpz_t *fun,
  mpz_t five_L,
  mpz_t *five_i,
  mpz_t result)
{
    mpz_t t1,t2,t3,t4;
    int i;
    mpz_init(t1);mpz_init(t2);mpz_init(t3);mpz_init(t4);
    mpz_fdiv_qr(t1,t2,N,five_L);
    i=mpz_tstbit(t1,0);//If odd, multiply result by -1, otherwise 1
    if(i){//odd
 mpz_sub_ui(result,five_L,1);
    }else{
 mpz_set_ui(result,1);
    }
    mpz_set_ui(t4,0);
    for(i=L-1;i>=0;i--){
 int r,s;
 mpz_fdiv_qr(t1,t3,t2,five_i[i]);
 mpz_set(t2,t3);
 //Assert t1 is from 0 to 4
 if(mpz_cmp_ui(t1,5)>=0){
     fprintf(stderr,"Data error i=%d\n",i);
     exit(-1);
 }
 r=mpz_get_ui(t1);
 mpz_mul_ui(t4,t4,5);
 if(i>0){
   for(s=0;s      Ft4(fun,t4,i,five_L,L,t1);
      mpz_mul(t3,result,t1);
      mpz_mod(result,t3,five_L);
      mpz_add_ui(t4,t4,1);
   }
 }else{
   for(s=1;s<=r;s++){
      mpz_add_ui(t4,t4,1);
      mpz_mul(t3,result,t4);
      mpz_mod(result,t3,five_L);
   }
 }
    }
    mpz_clear(t4);
    mpz_clear(t3);
    mpz_clear(t2);
    mpz_clear(t1);
}

void Get_m5L(mpz_t N, int L, mpz_t m5L, mpz_t fiveL,mpz_t count_5)
{
    int i,j;
    mpz_t *triangle = (mpz_t *)malloc(sizeof(mpz_t)*L*L);
    mpz_t *fun = (mpz_t *)malloc(sizeof(mpz_t)*L*L);
    mpz_t *tmp1 = (mpz_t *)malloc(sizeof(mpz_t)*L);
    mpz_t *tmp2 = (mpz_t *)malloc(sizeof(mpz_t)*L);
    mpz_t *tmp3 = (mpz_t *)malloc(sizeof(mpz_t)*L);
    mpz_t *tmp4 = (mpz_t *)malloc(sizeof(mpz_t)*L);
    mpz_t *tmp5 = (mpz_t *)malloc(sizeof(mpz_t)*L);
    mpz_t tmp,tmpp,tmppp; mpz_init(tmp);mpz_init(tmpp);mpz_init(tmppp);

    for(i=0;i    for(i=0;i    for(i=0;i mpz_init(tmp1[i]);
 mpz_init(tmp2[i]);
 mpz_init(tmp3[i]);
 mpz_init(tmp4[i]);
 mpz_init(tmp5[i]);
    }
  
   //Initialize Pascal Triangle
    mpz_set_ui(TRIANGLE(0,0),1);
    mpz_set_ui(TRIANGLE(1,0),1);
    mpz_set_ui(TRIANGLE(1,1),1);
    for(i=2;i mpz_set_ui(TRIANGLE(i,0),1);
 mpz_set_ui(TRIANGLE(i,i),1);
 for(j=1;j     mpz_add(TRIANGLE(i,j),TRIANGLE(i-1,j-1),TRIANGLE(i-1,j));
 }
    }
    /*
    for(i=1;i        int j;
        for(j=0;j<=i;j++){
       mpz_out_str(stdout,10,TRIANGLE(i,j));
     printf(",");
 }
 printf("\n\n");
    }*/

    for(i=0;i    mpz_set_si(F(0,1),1);
    //Initialize F(1,x)
    mpz_set_ui(F(1,0),24);
    if(L>1)
       mpz_set_ui(F(1,1),50*5);
    if(L>2)
       mpz_set_ui(F(1,2),35*5*5);
    if(L>3)
       mpz_set_ui(F(1,3),10*5*5*5);
    if(L>4)
       mpz_set_ui(F(1,4),1*5*5*5*5);
    for(i=0;i<=4;i++)mpz_mod(F(1,i),F(1,i),fiveL);
    for(i=5;i

    for(i=2;i int j,k;

// fprintf(stderr, "i=%d\n",i);
 
 //set tmp1[i] to be 5^i
 mpz_set_ui(tmp1[0],1);
 for(j=1;j     mpz_mul_ui(tmp1[j],tmp1[j-1],5);
 }
 for(j=0;j     mpz_mul(tmp,F(i-1,j),tmp1[j]);
     mpz_mod(tmp4[j],tmp,fiveL);
 }
/* fprintf(stderr,"F(%d,*)5^*:",i-1);
 for(j=0;j     mpz_out_str(stderr,10,tmp4[j]);fprintf(stderr,",");
 }
 fprintf(stderr,"\n");
*/
 for(j=0;j for(j=0;j     for(k=0;k<=j;k++){
   mpz_mul(tmp,F(i-1,j),TRIANGLE(j,k));//F(i-1,j)*C(j,k)
/*   if(k==0){
       fprintf(stderr,"F(%d,%d) is ",i-1,j);
       mpz_out_str(stderr,10,F(i-1,j));
       fprintf(stderr,"\nTRIANGLE(%d,%d) is ",j,k);
       mpz_out_str(stderr,10,TRIANGLE(j,k));
       fprintf(stderr,"\ntmp is ");
       mpz_out_str(stderr,10,tmp);
       fprintf(stderr,"\n5^0 is ");
       mpz_out_str(stderr,10,tmp1[k]);
       fprintf(stderr,"\n");
   }
*/   mpz_mod(tmpp,tmp,fiveL);
   mpz_mul(tmp,tmpp,tmp1[k]);//*5^k
   mpz_add(tmpp,tmp3[k],tmp);
   mpz_mod(tmp3[k],tmpp,fiveL);
/*   if(k==0){
       fprintf(stderr,"tmp is ");
       mpz_out_str(stderr,10,tmp);
       fprintf(stderr,"\ntmpp is ");
       mpz_out_str(stderr,10,tmpp);
       fprintf(stderr,"\ntmp3[0] is ");
       mpz_out_str(stderr,10,tmp3[k]);
       fprintf(stderr,"\n");
   }
*/     }
 }

/* fprintf(stderr,"F(%d,*)(5x+1)^*:",i-1);
 for(j=0;j     mpz_out_str(stderr,10,tmp3[j]);fprintf(stderr,",");
 }
 fprintf(stderr,"\n");
*/ //Next get Multiplication of polynomial tmp4, tmp3 into tmp5
 PolyMul(tmp3,tmp4,tmp5,L,fiveL); //tmp5=tmp3*tmp4;
        mpz_set_ui(tmp2[0],1); //set tmp2 to 2^j
        for(j=1;j            mpz_mul_ui(tmp2[j],tmp2[j-1],2);
        }
  
 for(j=0;j for(j=0;j     for(k=0;k<=j;k++){
  mpz_mul(tmp,F(i-1,j),TRIANGLE(j,k));//F(i-1,j)*C(j,k)
  mpz_mod(tmpp,tmp,fiveL);
  mpz_mul(tmp,tmpp,tmp1[k]);//*5^k
  mpz_mod(tmpp,tmp,fiveL);
  mpz_mul(tmp,tmpp,tmp2[j-k]);//*2^(j-k)
  mpz_add(tmpp,tmp3[k],tmp);
  mpz_mod(tmp3[k],tmpp,fiveL);
     }
 }
/* fprintf(stderr,"F(%d,*)(5x+2)^*:",i-1);
 for(j=0;j     mpz_out_str(stderr,10,tmp3[j]);fprintf(stderr,",");
 }
 fprintf(stderr,"\n");
*/
 PolyMul(tmp3,tmp5,tmp4,L,fiveL);//tmp4=tmp3*tmp5

        mpz_set_ui(tmp2[0],1); //set tmp2 to 3^j
        for(j=1;j            mpz_mul_ui(tmp2[j],tmp2[j-1],3);
        }

 for(j=0;j for(j=0;j     for(k=0;k<=j;k++){
  mpz_mul(tmp,F(i-1,j),TRIANGLE(j,k));//F(i-1,j)*C(j,k)
  mpz_mod(tmpp,tmp,fiveL);
  mpz_mul(tmp,tmpp,tmp1[k]);//*5^k
  mpz_mod(tmpp,tmp,fiveL);
  mpz_mul(tmp,tmpp,tmp2[j-k]);//*3^(j-k)
  mpz_add(tmpp,tmp3[k],tmp);
  mpz_mod(tmp3[k],tmpp,fiveL);
     }
 }
/* fprintf(stderr,"F(%d,*)(5x+3)^*:",i-1);
 for(j=0;j     mpz_out_str(stderr,10,tmp3[j]);fprintf(stderr,",");
 }
 fprintf(stderr,"\n");
*/ PolyMul(tmp3,tmp4, tmp5, L, fiveL); //tmp5=tmp3*tmp4;
 
        mpz_set_ui(tmp2[0],1); //set tmp2 to 4^j
        for(j=1;j            mpz_mul_ui(tmp2[j],tmp2[j-1],4);
        }

 for(j=0;j for(j=0;j     for(k=0;k<=j;k++){
  mpz_mul(tmp,F(i-1,j),TRIANGLE(j,k));//F(i-1,j)*C(j,k)
  mpz_mod(tmpp,tmp,fiveL);
  mpz_mul(tmp,tmpp,tmp1[k]);//*5^k
  mpz_mod(tmpp,tmp,fiveL);
  mpz_mul(tmp,tmpp,tmp2[j-k]);//*4^(j-k)
  mpz_add(tmpp,tmp3[k],tmp);
  mpz_mod(tmp3[k],tmpp,fiveL);
     }
 }
/* fprintf(stderr,"F(%d,*)(5x+1)^*:",i-1);
 for(j=0;j     mpz_out_str(stderr,10,tmp3[j]);fprintf(stderr,",");
 }
 fprintf(stderr,"\n");
*/
 PolyMul(tmp3,tmp5, tmp4, L, fiveL); //tmp4=tmp3*tmp5;

 //Now tmp4 has result F(i,X)
 for(j=0;j    }

    //Now we could free PASCAL Triangle
    for(i=1;i    free(triangle);
    for(i=0;i mpz_clear(tmp2[i]);
 mpz_clear(tmp3[i]);
 mpz_clear(tmp4[i]);
 mpz_clear(tmp5[i]);
    }
    free(tmp2);free(tmp3);free(tmp4);free(tmp5);
    //tmp1 which holds 5^i reserved for later use

    //Check for F(L-1,i), F(L-1,0) should be -1 and F(L-1, i) should be 0 for all i>0
    mpz_add_ui(tmp,F(L-1,0),1);
    if(mpz_cmp(tmp,fiveL)!=0){
 for(i=0;i     int j;
     for(j=0;j  mpz_out_str(stdout,10,F(i,j));
  printf(",");
     }
     printf("\n\n");
 }
 fprintf(stderr,"Cannot pass verification for F(%d,0)=1\n",L-1);
 exit(-1);
    }
    for(i=1;i if(mpz_sgn(F(L-1,i))!=0){
     fprintf(stderr,"Cannot pass verification for F(%d,%d)=0\n",L-1,i);
     exit(-1);
 }
    }
 
    mpz_set_ui(m5L,1); //Initialize to 1;
    mpz_set_ui(count_5,0);
    while(mpz_sgn(N)>0){
 Get_Mult_No_5(N,L,fun,fiveL,tmp1,tmp);
 mpz_mul(tmpp,m5L,tmp);
 mpz_mod(m5L,tmpp,fiveL);
 mpz_fdiv_q_ui(N,N,5);
 mpz_add(count_5,count_5,N);
    }
    for(i=0;i    free(tmp1);
    for(i=1;i    free(fun);
    mpz_clear(tmp);mpz_clear(tmpp);mpz_clear(tmppp);
}

void short_calc(int N, int L){
    mpz_t result;
    mpz_t tenL;
    int i,count5;
    mpz_init(result);
    mpz_init(tenL);
    mpz_set_ui(tenL,10);
    mpz_pow_ui(tenL,tenL,L);
    mpz_set_ui(result,1);
    count5=0;
    i=N;
    while(i>0){i/=5;count5+=i;}
    for(i=1;i<=N;i++){
 int j=i;
 while(j%5==0)j/=5;
 while(count5>0&&j%2==0){j/=2;count5--;}
 mpz_mul_ui(result,result,j);
 mpz_mod(result,result,tenL);
    }
    mpz_out_str(stdout,10,result);
    printf("\n");
}
void calc(mpz_t N, int L){
     int smallN=L*4;
     mpz_t m5L;
     mpz_t two,twomL;
     mpz_t fiveL;
     mpz_t result;
     mpz_t count_5;
     if(mpz_cmp_ui(N,smallN)<0){
  short_calc(mpz_get_ui(N),L);
  return;
     }
     mpz_init(count_5);
     mpz_init(m5L);
     mpz_init(fiveL);
     mpz_set_ui(m5L,5);
     mpz_pow_ui(fiveL, m5L, L);//fiveL=5^L
     Get_m5L(N,L,m5L,fiveL,count_5);  //m5L = final_result % (5^L)
     mpz_add_ui(count_5,count_5,L);
     mpz_init(two);
     mpz_init(twomL);
     mpz_set_ui(two,2);
     mpz_add_ui(twomL,fiveL,1);
     mpz_fdiv_q_ui(twomL,twomL,2);
     mpz_powm(twomL,twomL,count_5,fiveL);// (twomL = 2^(-L) mod (5^L))
     mpz_init(result);
     mpz_mul(result,twomL,m5L);   
     mpz_mod(twomL,result, fiveL); //twomL = (m5L*2^(-L)) %(5^L)
     mpz_pow_ui(m5L,two,L); //m5L= 2^L
     mpz_mul(result,twomL, m5L); //2^L * (m5L*2^(-L))%(5^L), the final result
     mpz_out_str(stdout,10,result);
     printf("\n");
     mpz_clear(result);
     mpz_clear(twomL);
     mpz_clear(two);
     mpz_clear(fiveL);
     mpz_clear(m5L);
}

int
main(void)
{
    mpz_t N;
    int L;
    size_t len;
    mpz_init(N);
    printf("Please input large number N, so that we could process for N!\n");
    len=mpz_inp_str(N,stdin,10);
    if(len>1024){
     printf("Current only process number whose length no more than 1000\n");
     return -1;
    }
    printf("Please input number of digits to get (no more than 100):");
    scanf("%d",&L);
    if(L>100){
    printf("Input Out of range\n");
    return -1;
    }
    if(L<3)L=3;
    printf("Calcuate last %d non-zero digits of ",L);
    mpz_out_str(stdout,10,N);
    printf("\n");
    calc(N,L);
    mpz_clear(N);
    return 0;
} 

 

转自:http://blog.csdn.net/mathe/archive/2006/08/28/1132404.aspx