擂台：计算阶乘的 18  位非零尾数

求多少大的数的阶乘的 18 位非零尾数？不指明这一点，怎么行？

至少百位以上的吧。硬算的话，10亿以内的数都能在可接受的时间内硬算出来。

good!不过我想知道为什么一定要18位呢?只有18位才有好算法吗?

这主要是为了内部数据表示比较方便罢了。

    一个 32 位整数，最多可以表示到 10 的 9 次方，两个 32 位整数就可以表示到 10 的 18 次方。当然，两个 32 位整数最多可以表示到 10 的 19 次方，但是拿一个 64 位整数做一个基本运算单元在 32 位 CPU 上进行运算是不合算且不方便的。所以，就只计算到了 18 位。

    另外，尾数的个数可以再适当的增加五六位，不过增加之后，对于很大的一个数字 N，计算量的增加也是比较大的。这也是原因之一。

    BTW：最新版本中，对于 N 等于 10 的 10000 次方减 1, 在我机器上 0.3 秒多一点点就可以算出 N! 的 18 位非零尾数。

不是很难的,主要在于计算
1*2*3*4*6*7*8*9*11*12*... (mod 5^18)
也就是去掉所有5的倍数后的乘积.

方法如下:
由于对于充分大的N,N!中因子2的数目要远远大于因子5的数目,我们可以认为对于充分大,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]) (去掉了所有的因子5)
所以关键是求g(N),也就是1到N之间所有不是5的倍数的数的乘积关于5^18的余数
也就是最前面列出的.
现在讨论如何计算:
1*2*3*4*6*7*8*9*11*12*...*N (mod 5^18)
可以知道的是,
1*2*3*4*6*7*8*9*11*12*...*(5^17-1)=-1 (mod 5^18)
同样
(5^17*k+1)*(5^17*k+2)*...*(5^17*k+5^17-1)=-1 (mod 5^18) 所有5的倍数不包括
所以对于N>5^17,我们可以将前面的若干段去掉,全部用-1替换,最后变成一个长度不超过
5^17的连乘
(k*5^17+1)*(k*5^17+2)*...*(k*5^17+s) (mod 5^18)  其中s<5^17,5的倍数全部被去掉
然后剩下的数据我们再根据5^16的长度来分段,其中完整的段可以写成
(k*5^16+1)*(k*5^16+2)*....*(k*5^16+5^16-1) (mod 5^18)
这个可以看成关于k的函数,是个周期为5的函数,可以事先算好,保存起来.
同样剩下还有一段不完整的,我们再根据5^15长度来分段,完整段关于模5^18是周期为25的函数,
同样可以事先计算. 
然后再划分余下不完整的段...
当最后长度为5^10时,是个周期为5^7的函数,也可以事先计算,不过一个70多K长的数据,还可以忍受的
最后长度不超过5^10的部分,直接计算它的值就可以了.

再给个运行时候直接计算
1*2*3*4*6*7*8*9*11*12*... (mod 5^18)
中任何一个长度为5^10段乘积的方法
这样就不需要预先计算了.
首先直接计算
1*2*3*4*6*7*8*9*11*...*(5^10-1)(mod 5^18)
然后对于
(k*5^10+1)(k*5^10+2)....(k*5^10+5^10-1) (mod 5^18) (去掉所有5的倍数)
=1*2*3*...*(5^10-1)+1*2*3*...*(5^10-1)*(1/1+1/2+...+1/(5^10-1))*k*5^10
(mod 5^18)
其中
1*2*3*...*(5^10-1)*(1/1+1/2+...+1/(5^10-1))项中
由于1*2*3*...*(5^10-1)已经知道,所以我们可以通过计算1,2,...,5^10-1关于5^18的数论倒数的方法
来得到,这个不难

现在我们已经有公式
(k*5^10+1)(k*5^10+2)....(k*5^10+5^10-1) (mod 5^18)=(a+bk*5^10) (mod 5^18)
我们可以将k=5s,5s+1,5s+2,5s+3,5s+4带入相乘,得到
(s*5^11+1)(s*5^11+2)....(s*5^11+5^11-1) (mod 5^18)
=(a+5bs*5^10)(a+(5bs+b)*5^10)(a+(5bs+2b)*5^10)(a+(5bs+3b)*5^10)(a+(5bs+4b)*5^10) (mod 5^18)
=(a^5+a^4*5^10*(5bs+5bs+b+5bs+2b+5bs+3b+5bs+4b)) (mod 5^18)
=(a^5+a^4*5^11*(5bs+2b)) (mod 5^18)
用同样的方法可以推导出连续5^12项,5^13项,...的值.

从上面方法我们可以看出,由于只要求计算18位,我们可以从连续5^9项出发,这样开始连乘
可以只需要计算1/5数据的连乘,可以快5倍的速度.

另外还有一种可以考虑的方法是先只计算5^6项的乘积
(k*5^6+1)(k*5^6+2)...(k*5^6+5^6-1)(mod 5^18)
=1*2*...*(5^6-1) * (1+ (1/1+1/2+...+1/(5^6-1)) *5^6*k+ (1/(1*2)+1/(1*3)+...+1/((5^6-2)(5^6-1))*5^12*k^2) (mod 5^18)
这里主要的问题在于计算
=1*2*...*(5^6-1) *(1/(1*2)+1/(1*3)+...+1/((5^6-2)(5^6-1)) (mod 5^18)
同样,对于1/i我们可以用i关于5^18的数论倒数替换
然后
(1/(1*2)+1/(1*3)+...+1/((5^6-2)(5^6-1))=
1/2 *{ (1/1+1/2+...+1/(5^6-1))^2 - 1/1^2-1/2^2-...-1/(5^6-1)^2)}
于是我们只需要计算
1*2*...*(5^6-1)
1/1+1/2+...+1/(5^6-1)
1/1^2+1/2^2+...+1/(5^6-1)^2
关于5^18的数论倒数,这样,结果就可以变成:
(k*5^6+1)(k*5^6+2)...(k*5^6+5^6-1)(mod 5^18)
=a+b*k*5^6+c*k^2*5^12 (mod 5^18)
然后用类似上面的方法,我们可以得出从5^6到5^7的系数递推方法,
同样方法可以用于5^7到5^8,...
通过这种方法,计算所需的时间会非常非常小了,因为无论我们需要计算的连乘
还是连续和,都不超过5^6项,对现在的计算机,应该是只需要毫秒级的时间

呵呵,这个新算法够2000分了吧!!!

想了一下扩展问题:
对于N位数X,计算X!去掉末尾的0的最后L位的算法的时间复杂度可以达到:
O(L^3 log(L) log(N) + L^2 N Log(L) Log(N))
使用的空间复杂度为O(L^3+L^2 log(N)+ N)

计算过程中,我们总假设一个长度为L的数占用的空间为O(L),两个长度为L的数相乘花费的时间为O(L*log(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) log(N)),我们可以计算出F(k,x),其中1<=k<=ceil(log5(N))
   而保存所有的F(k,x),需要花费的空间是O(L^2 log(N)) (L *log(N)个数)
iii)对于X! (mod 5^L),我们首先只取X!中所有不是5的倍数的数,它们的乘积看成Y,那么X!=Y*floor(X/5)! 
     而Y (mod 5^L)我们可以将它划分成若干个长度不一的类似ii)中的连乘式,这样的式子有log5(N)个,对于每个连乘式,我们代入公式ii),需要L次乘法,
     所以计算Y(mod 5^L)共花费时间为O(log(N) L^2 log(L))
v)在计算完Y (mod 5^L)以后,我们还需要计算floor(X/5)! (mod 5^L),同样的,如果我们把X看成N位5进制数,那么X/5是N-1位5进制数,所以通过同样算法,可以再花费O(log(N-1) L^2 log(L))的时间递归到L-2位数,...,这样总共经过N步后就会得出最终结果.
   这个递归过程中总共花费时间最多为 O(N log(N) L^2 log(L))
而保存原式输入数据X (以及中间数据X/5, X/25,...等等)需要一个长度为N的空间,需要O(N)的空间
所以总共需要
时间复杂度: O(L^3 log(L) log(N) + L^2 N Log(L) Log(N))
空间复杂度: O(L^3+L^2 log(N)+ N)

上面过程只是算出X! (mod 5^L),而对于充分大的X(超过10L/3肯定够了,而对于这么小的X,直接计算乘积就可以了),X! (mod 2^L)总是0,所以通过中国剩余定理就可以计算出X! (mod 10^L)
对于现在的计算机基本上在L<=100而且N<=1000时不会有问题.

忘了一件事,时间空间复杂度可以达到:
对于N位数X,计算X!去掉末尾的0的最后L位的算法的时间复杂度可以达到:
O(L^3 log(L)^2 + L^2 N Log(L) Log(N))
使用的空间复杂度为O(L^3+L^2 log(L)+ N)
因为计算F(k,x)只要算1<=k<=L-1,对于k>=L,容易证明 F(k,x)=-1

给个程序,计算100位以内基本上很快就出来
使用了GMP来做大整数运算,在Linux上面编译通过

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

#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<L;i++)
mpz_set_ui(dst[i],0);
    for(i=0;i<L;i++){
for(j=0;j<L-i;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<r;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<L;i++)for(j=0;j<=i;j++)mpz_init(TRIANGLE(i,j));
    for(i=0;i<L;i++)for(j=0;j<L;j++)mpz_init(F(i,j));
    for(i=0;i<L;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<L;i++){
mpz_set_ui(TRIANGLE(i,0),1);
mpz_set_ui(TRIANGLE(i,i),1);
for(j=1;j<i;j++){
    mpz_add(TRIANGLE(i,j),TRIANGLE(i-1,j-1),TRIANGLE(i-1,j));
}
    }
    /*
    for(i=1;i<L;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<L;i++)mpz_set_si(F(0,i),0);
    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<L;i++)mpz_set_ui(F(1,i),0);

    for(i=2;i<L;i++){//Set F(i,x) now
int j,k;

//set tmp1[i] to be 5^i
mpz_set_ui(tmp1[0],1);
for(j=1;j<L;j++){
    mpz_mul_ui(tmp1[j],tmp1[j-1],5);
}
for(j=0;j<L;j++){ //tmp4 is F(i-1,j)*5^j mod 5^L
    mpz_mul(tmp,F(i-1,j),tmp1[j]);
    mpz_mod(tmp4[j],tmp,fiveL);
}
for(j=0;j<L;j++) mpz_set_ui(tmp3[j],0);
for(j=0;j<L;j++){//Adding F(i-1,j)*(5x+1)^j mod 5^L into tmp3
    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_add(tmpp,tmp3[k],tmp);
 mpz_mod(tmp3[k],tmpp,fiveL);
    }
}

//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<L;j++){
            mpz_mul_ui(tmp2[j],tmp2[j-1],2);
        }

for(j=0;j<L;j++) mpz_set_ui(tmp3[j],0);
for(j=0;j<L;j++){ //Adding F(i-1,j)*(5x+2)^j mod 5^L into tmp3
    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);
    }
}

PolyMul(tmp3,tmp5,tmp4,L,fiveL);//tmp4=tmp3*tmp5

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

for(j=0;j<L;j++) mpz_set_ui(tmp3[j],0);
for(j=0;j<L;j++){ //Adding F(i-1,j)*(5x+3)^j mod 5^L into tmp3
    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);
    }
}
PolyMul(tmp3,tmp4, tmp5, L, fiveL); //tmp5=tmp3*tmp4;

        mpz_set_ui(tmp2[0],1); //set tmp2 to 4^j
        for(j=1;j<L;j++){
            mpz_mul_ui(tmp2[j],tmp2[j-1],4);
        }

for(j=0;j<L;j++) mpz_set_ui(tmp3[j],0);
for(j=0;j<L;j++){ //Adding F(i-1,j)*(5x+4)^j mod 5^L into tmp3
    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);
    }
}

PolyMul(tmp3,tmp5, tmp4, L, fiveL); //tmp4=tmp3*tmp5;

//Now tmp4 has result F(i,X)
for(j=0;j<L;j++)mpz_set(F(i,j),tmp4[j]);
    }

    //Now we could free PASCAL Triangle
    for(i=1;i<L;i++)for(j=0;j<=i;j++)mpz_clear(TRIANGLE(i,j));
    free(triangle);
    for(i=0;i<L;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<L;i++){
    int j;
    for(j=0;j<L;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<L;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<L;i++)mpz_clear(tmp1[i]);
    free(tmp1);
    for(i=0;i<L;i++)for(j=0;j<L;j++)mpz_clear(F(i,j));
    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;
}

呵呵 这个paper挺不错。

强帖，虽然很难。

该问题的解决必须充分运用数论知识，mathe给出了非常理想的算法并给出了源代码；
但可惜的是是基于Linux平台下的，下面给出改编成Windows平台下的源代码。
该程序在 VC6.0 下 Win98 or WinXP 下均编译通过。

由于调用了最新的 HugeCalc V5.x，所以请按如下步骤操作：
1、从 http://maths.diy.myrice.com/software.htm#02 下载 HugeCalc 的完全版；
2、解压缩 HugeCalc.rar;
3、在 \HugeCalc_Dll_Import\ 的同级目录里新建目录，如：\FactorialTail\
4、在新目录新建一个cpp文件，将如下代码复制进去编译即可。


//#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>

#include "../HugeCalc_Dll_Import/HugeCalc.h"    // 公共接口
#include "../HugeCalc_Dll_Import/HugeInt.h"        // 10进制系统
#include "../HugeCalc_Dll_Import/HugeIntX.h"    // 16进制系统

#ifndef _UNICODE
    #pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalc.lib" )
#else
    #pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalcU.lib" )
#endif

//    Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
//        Win32 Debug:    Debug Multithreaded DLL
//        Win32 Release:    Multithreaded DLL


#define TRIANGLE(x,y)  triangle[(x)*L+(y)]
#define F(x,y)         fun[(x)*L+(y)]

void PolyMul(const CHugeInt * const src1, const CHugeInt * const src2, CHugeInt * const dst, int L, const CHugeInt &mod)
{
    int i,j;
    CHugeInt tmp;
    for(i=0;i<L;++i) dst[i]=0;

    for(i=0;i<L;++i)
    for(j=0;j<L-i;++j)
        ( dst[i+j] += tmp.Mul( src1[i],src2[j] )) %= mod;
}

void Ft4(const CHugeInt * const fun,//Use Polynomial F(s,*)
     const CHugeInt &t4,  //F(s,*)[t4] (mod five_L)
     int s,
     const CHugeInt &five_L,
     int L,
     CHugeInt &result)
{
    int i;
    CHugeInt tmp1;
    tmp1.Mul( F(s, L-1), t4 );
    result = tmp1 % five_L;
    for(i=L-2;i>=1;i--){
        (tmp1 = result ) +=F(s, i);
        result.Mul( tmp1, t4 ) %= five_L;
    }
    ( result.Swap( tmp1 ) += F(s, 0)) %= five_L;
}

int GetVal( const CHugeInt &HugeInt )
{
    SINT64 s64Num = -1;
    HugeInt.CanConvertSINT64( s64Num );
    return (int)s64Num;
}

void Get_Mult_No_5(const CHugeInt &N,
        int L,
        const CHugeInt *const fun,
        const CHugeInt &five_L,
        const CHugeInt *const five_i,
        CHugeInt &result)
{
    CHugeInt t1,t2,t3,t4;
    int i;
    t1.Div(N,five_L,&t2);
    if(t1.IsOdd()){//odd
        --(result = five_L);
    }else{
        result = 1;
    }

    int r,s;
    for(i=L-1;i>=0;--i){
        t1.Div(t2,five_i[i],&t3);
        t2=t3;
        //Assert t1 is from 0 to 4
        if(Compare(t1,5)>=0){
            fprintf(stderr,"Data error i=%d\n",i);
            exit(-1);
        }
        r=GetVal( t1 );
        t4 *= 5;
        if(i>0){
            for(s=0;s<r;++s){
                Ft4(fun,t4,i,five_L,L,t1);
                t3.Mul(result,t1);
                result=t3 % five_L;
                ++t4;
            }
        }else{
            for(s=1;s<=r;++s){
                ++t4;
                t3.Mul(result,t4);
                result = t3 % five_L;
            }
        }
    }
}

void Get_m5L( CHugeInt &N, int L, CHugeInt &m5L, const CHugeInt &fiveL, CHugeInt &count_5)
{
    int i,j;
    int L2 = L * L;
    CHugeInt *triangle = new CHugeInt[L2];
    CHugeInt *fun = new CHugeInt[L2];
    CHugeInt *tmp1 = new CHugeInt[L];
    CHugeInt *tmp2 = new CHugeInt[L];
    CHugeInt *tmp3 = new CHugeInt[L];
    CHugeInt *tmp4 = new CHugeInt[L];
    CHugeInt *tmp5 = new CHugeInt[L];
    CHugeInt tmp,tmpp,tmppp;

    for(i=0;i<L;++i)for(j=0;j<=i;++j) TRIANGLE(i,j) = 0;
    for(i=0;i<L;++i)for(j=0;j<L;++j) F(i,j)=0;
    for(i=0;i<L;++i){
        tmp1[i]=0;
        tmp2[i]=0;
        tmp3[i]=0;
        tmp4[i]=0;
        tmp5[i]=0;
    }

    //Initialize Pascal Triangle
    TRIANGLE(0,0) = 1;
    TRIANGLE(1,0) = 1;
    TRIANGLE(1,1) = 1;
    for(i=2;i<L;++i){
        TRIANGLE(i,0) = 1;
        TRIANGLE(i,i) = 1;
        for(j=1;j<i;++j)
            (TRIANGLE(i,j)=TRIANGLE(i-1,j))+=TRIANGLE(i-1,j-1);
    }

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

    for(i=2;i<L;++i){//Set F(i,x) now
        int j,k;

        //set tmp1[i] to be 5^i
        tmp1[0]=1;
        for(j=1;j<L;++j) ( tmp1[j] = tmp1[j-1] ) *= 5;
        for(j=0;j<L;++j){ //tmp4 is F(i-1,j)*5^j mod 5^L
            tmp.Mul( F(i-1,j),tmp1[j] );
            tmp4[j] = tmp % fiveL;
        }
        for(j=0;j<L;++j) tmp3[j]=0;
        for(j=0;j<L;++j){//Adding F(i-1,j)*(5x+1)^j mod 5^L into tmp3
            for(k=0;k<=j;++k){
                tmp.Mul( F(i-1,j),TRIANGLE(j,k) );//F(i-1,j)*C(j,k)
                tmpp = tmp % fiveL;
                tmp.Mul( tmpp, tmp1[k]);//*5^k
                ( tmpp = tmp3[k] ) += tmp;
                tmp3[k] = tmpp % fiveL;
            }
        }
        //Next get Multiplication of polynomial tmp4, tmp3 into tmp5
        PolyMul(tmp3,tmp4,tmp5,L,fiveL); //tmp5=tmp3*tmp4;

        tmp2[0]=1; //set tmp2 to 2^j
        for(j=1;j<L;++j) ( tmp2[j] = tmp2[j-1] ) *= 2;

        for(j=0;j<L;++j) tmp3[j]=0;
        for(j=0;j<L;++j){ //Adding F(i-1,j)*(5x+2)^j mod 5^L into tmp3
            for(k=0;k<=j;++k){
                tmp.Mul( F(i-1,j),TRIANGLE(j,k));//F(i-1,j)*C(j,k)
                tmpp = tmp % fiveL;
                tmp.Mul( tmpp,tmp1[k]);//*5^k
                tmpp = tmp % fiveL;
                tmp.Mul( tmpp,tmp2[j-k]);//*2^(j-k)
                (tmpp=tmp3[k])+=tmp;
                tmp3[k]= tmpp % fiveL;
            }
        }

        PolyMul(tmp3,tmp5,tmp4,L,fiveL);//tmp4=tmp3*tmp5

        tmp2[0]=1; //set tmp2 to 3^j
        for(j=1;j<L;++j) ( tmp2[j] = tmp2[j-1] ) *= 3;

        for(j=0;j<L;++j) tmp3[j]=0;
        for(j=0;j<L;++j){ //Adding F(i-1,j)*(5x+3)^j mod 5^L into tmp3
            for(k=0;k<=j;++k){
                tmp.Mul(F(i-1,j),TRIANGLE(j,k));//F(i-1,j)*C(j,k)
                tmpp = tmp % fiveL;
                tmp.Mul(tmpp,tmp1[k]);//*5^k
                tmpp = tmp % fiveL;
                tmp.Mul(tmpp,tmp2[j-k]);//*3^(j-k)
                (tmpp=tmp3[k])+=tmp;
                tmp3[k] = tmpp % fiveL;
            }
        }
        PolyMul( tmp3, tmp4, tmp5, L, fiveL); //tmp5=tmp3*tmp4;

        tmp2[0] = 1; //set tmp2 to 4^j
        for(j=1;j<L;++j) (tmp2[j]=tmp2[j-1]) *= 4;

        for(j=0;j<L;++j) tmp3[j]=0;
        for(j=0;j<L;++j){ //Adding F(i-1,j)*(5x+4)^j mod 5^L into tmp3
            for(k=0;k<=j;++k){
                tmp.Mul( F(i-1,j),TRIANGLE(j,k));//F(i-1,j)*C(j,k)
                tmpp = tmp % fiveL;
                tmp.Mul( tmpp,tmp1[k]);//*5^k
                tmpp = tmp % fiveL;
                tmp.Mul( tmpp,tmp2[j-k]);//*4^(j-k)
                ( tmpp=tmp3[k] ) += tmp;
                tmp3[k] = tmpp % fiveL;
            }
        }

        PolyMul(tmp3,tmp5, tmp4, L, fiveL); //tmp4=tmp3*tmp5;

        //Now tmp4 has result F(i,X)
        for(j=0;j<L;++j) F(i,j) = tmp4[j];
    }

    //Now we could free PASCAL Triangle
    delete []triangle;
    delete []tmp2;
    delete []tmp3;
    delete []tmp4;
    delete []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
    ++(tmp = F(L-1,0));
    if( tmp != fiveL ){
        for(i=0;i<L;++i){
            int j;
            for(j=0;j<L;++j){
                printf( F(i,j).ConvertToStr(FS_NORMAL) );
                printf(",");
            }
            printf("\n\n");
        }
        fprintf(stderr,"Cannot pass verification for F(%d,0)=1\n",L-1);
        exit(-1);
    }
    for(i=1;i<L;++i){
        if( F(L-1,i).GetSign() !=0 ){
            fprintf(stderr,"Cannot pass verification for F(%d,%d)=0\n",L-1,i);
            exit(-1);
        }
    }

    m5L=1; //Initialize to 1;
    count_5 = 0;
    while(N > 0){
        Get_Mult_No_5(N,L,fun,fiveL,tmp1,tmp);
        tmpp.Mul(m5L,tmp);
        m5L = tmpp % fiveL;
        N/=5;
        count_5+=N;
    }
    delete []tmp1;
    delete []fun;
}

void short_calc(int N, int L){
    CHugeInt result;
    CHugeInt tenL;
    int i,count5;
    tenL=10;
    tenL.Pow(L);
    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;}
        result *= j;
        result %= tenL;
    }
    printf( result.ConvertToStr(FS_NORMAL) );
    printf("\n");
}
void calc( CHugeInt &N, int L){
    int smallN=L*4;
    CHugeInt m5L;
    CHugeInt two,twomL;
    CHugeInt fiveL;
    CHugeInt result;
    CHugeInt count_5;
    if(N<smallN){
        short_calc(GetVal(N),L);
        return;
    }

    m5L=5;
    fiveL.Pow(m5L, L);//fiveL=5^L
    Get_m5L(N,L,m5L,fiveL,count_5);  //m5L = final_result % (5^L)
    count_5+=L;
    two=2;
    ++(twomL = fiveL);
    twomL/=2;

    // mpz_powm(twomL,twomL,count_5,fiveL);// (twomL = 2^(-L) mod (5^L))
    if ( fiveL<0x10000000 )
    {
        int mod = GetVal(fiveL);
        twomL = HugeCalc::PowMod( (UINT32)(twomL % mod), (UINT32)(count_5 % (mod*4/5)), mod );
    }
    else
    {
        CHugeInt hugeBase(twomL % fiveL);
        CHugeIntX r(count_5 % (fiveL*4/5));    // Euler
        for( UINT32 k = r.GetBits() - 1; 0!=k; )
        {
            twomL.Pow( 2 ) %= fiveL;
            if ( r.IsBitOne(--k) )
                ( twomL *= hugeBase ) %= fiveL;
        }
    }

    result.Mul(twomL,m5L);
    twomL=result % fiveL; //twomL = (m5L*2^(-L)) %(5^L)
    m5L.Pow(two,L); //m5L= 2^L
    result.Mul(twomL, m5L); //2^L * (m5L*2^(-L))%(5^L), the final result
    printf( result.ConvertToStr(FS_NORMAL) );
    printf("\n");
}

int
main(void)
{
    if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
    {
        printf( "\r\n您未通过 HugeCalc 的认证许可！\r\n\r\n" \
                "解决方案可选下列方案之一：\r\n" \
                "    一、请移动和[或]修改文件名为：/CopyrightByGuoXianqiang/.../HugeCalc.exe；\r\n" \
                "    二、请在 HugeCalc.chm 中进行“(免费)注册”（一劳永逸）。\r\n\r\n" );

        return -1;
    }


    printf("Call %s\n", HugeCalc::GetVersion());
    while (TRUE)
    {
        printf("\nPlease input large number N, so that we could process for N!\n(if N is 0 then exit) N=");
        char buffer[1032] = { 0 };
        scanf("%s",buffer);
        CHugeInt N(buffer);
        if ( !N )
            return 0;

        if(N.GetDigitals()>1024){
            printf("Current only process number whose length no more than 1000\n");
            return -1;
        }

        int L;
        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<4)L=4;

        printf("Calcuate last %d non-zero digits of %s\n", L, N.ConvertToStr(FS_NORMAL));

        HugeCalc::ResetTimer();
        HugeCalc::EnableTimer();
        calc(N,L);
        printf("timer:%dms\n", HugeCalc::GetTimer() );
    }

    system("pause");
    return 0;
}

在 P4 1.7GHz/256MB 机器下：
N=123456789, T(18)=622499063305863168, timer:189ms
N=123456789, T(27)=922577745622499063305863168, timer:599ms

N=123456789123456789, T(18)=372919194014449664, timer:217ms
N=123456789123456789, T(27)=539805094372919194014449664, timer:650ms

注意：用楼主的程序计算：N=123456789123456789, T(18)=515278436455153664，可以肯定楼主结果不正确。

近期我将设计一个图形化界面的程序，并有针对性的优化算法，最终速度将可能比楼主的更高，敬请关注！

上面的程序还有个地方可以优化,特别是在N远远大于L时,
程序快结束时调用了:
mpz_powm(twomL,twomL,count_5,fiveL); 需要使用O(logX)次L位数乘法.
在count_5比较大时 (数量级上同X类似,再计算X!时)还是比较花费时间的
可以先计算
count_5'=count_5 %(5^(L-1)*4)
然后用count_5'代替count_5,
这里使用了 phi(5^L)=5^(L-1)*4的性质,这样这里就只需要O(L)次L位数乘法

对不起，请将函数 Ft4(..) 最后一行代码更正为：“( result += F(s, 0)) %= five_L;”
运行结果如下：

N=123456789, T(18)=626613857392263168, timer:188ms
N=123456789, T(27)=537575787626613857392263168, timer:591ms

N=123456789123456789, T(18)=515278436455153664, timer:216ms
N=123456789123456789, T(27)=213398051515278436455153664, timer:634ms

与楼主结果一致。

to mathe() : 
    在我上面公布的代码中已作了上述优化，咱们真是“英雄所见略同”啊。。。:)

我试了一下,在我的Linux下

Please input large number N, so that we could process for N!
Please input number of digits to get (no more than 100):Calcuate last 18 non-zero digits of 123456789
626613857392263168

real    0m0.011s
user    0m0.008s
sys     0m0.004s

Please input large number N, so that we could process for N!
Please input number of digits to get (no more than 100):Calcuate last 18 non-zero digits of 123456789123456789
515278436455153664

real    0m0.012s
user    0m0.008s
sys     0m0.004s
结果同楼主的Web版本相同的.

gxqcn(GxQ) 还是检查一下你的代码吧,可能哪里弄错了

呵呵,你的机器好像比我的慢很多倍呀:)

to mathe() : 
    我刚刚更正过来了，请见上上贴。

    注意到两者时间上的消耗差异比较大，这是因为没有充分发挥 HugeCalc 的优势，将代码略作改进，效率将会大幅提升（以上公布的代码更多是“直译”你发布的源代码）。

呵呵,发现算法还可以优化,主要部分可以提高L倍 ,也就是时间复杂度改为:
O(L^3 log(L)^2 + L*N Log(L) Log(N))

现在还是使用gmp来做运算,
输入是一个10000位的数,计算最后18位,原先算法1.880秒,现在只用0.749秒
输入是一个10000位的数,计算最后100位,原先算法34.668秒,现在只用2.993秒
对于位数比较低时,速度提高不是很明显(因为这时L^3 log(L)^2部分比重比较高)

这个算法就先不说了,让大家想一想:)

我试着在Windows下Cygwin环境下面用gcc编译了gmp,结果运行速度要比Linux下面慢一倍左右(而且我Windows机器比Linux要快).
比如10000位的数字,计算最后100位,要6.265秒, 计算18位0.969s
如果有人想要编译后的Windows程序,我可以发Email给他参考:)不过这个程序链接了一个cygwin1.dll,我不知道直接带上这个DLL是否就可以运行了

楼上的都有强人啊！
做个记号

我给出了个Windows下面的程序,大家可以试一试看:
http://218.1.231.240/iqbbs/viewFile.asp?Boardid=9&ID=3363

mathe()是个强人啊 ;)。

我的算法的起点也是计算 1*2*3*4*6*7*8*9*...*N mod (5^18)，不过有些地方还是稍有差异。最近比较忙，一直没有上来。今天刚上来就发现了强人 :)。现在只是粗看了一下以上各贴，启发多，容我抽时间再研究以上的内容。手头还有其他的活在忙着呢，抽不出身来啊。

利用午饭的时间仔细看了几贴，目前只看到了 mathe() 的前三贴。
其中使用的算法正是我程序中使用的算法，不过我稍有优化。计算的时候不超过 5 ^ 5 项，而不是 5 ^ 6 项。

回复人： gxqcn(GxQ) ( ) 信誉：100  2005-6-13 12:34:52  得分: 0  
   
to mathe() : 
    我刚刚更正过来了，请见上上贴。

    注意到两者时间上的消耗差异比较大，这是因为没有充分发挥 HugeCalc 的优势，将代码略作改进，效率将会大幅提升（以上公布的代码更多是“直译”你发布的源代码）。
------------------------------

    GMP 的优化就是针对基本运算，而且是相当好的。同样 mathe 的代码，用 HugeCalc 和 GMP ，运行时间差异很大，还是说明你的基本运算的优化是不够的。

    另外，象“library 的实现细节、哪些地方该如何做优化”这类的问题，用户知道的越少越好。说白了，User Interface 越简单越好，而且这些 Interface 是足够的，优化是到位的。GMP 一直秉承着这个原则。

前三贴内容还是比较一般的,再下面一个才是关键,可以推广到计算任意多位而不仅仅18位了

我经过简单优化，现在效率已提高数十倍。
但如果要一味追求效率，则可能需要为此重新定制大数运算，正如 xstring(麻雀) 那样。

其实如果能够有5进制运算的大数运算程序,上面的算法可以快很多的.
而且上面有部分代码应该可以通过FFT的方法来事现,比如PolyMult函数,这些也可以提高速度

to mathe():
    你我又想到一块了！（意指5进制运算的大数运算程序；但如需输出末尾0的数目，对于大数来说，进制转换将可能成为瓶颈）

仔细看完其余的贴子了.

mathe的算法果然高明，佩服之至啊。2000 分非你莫属了，；）

	
>>利用午饭的时间仔细看了几贴，目前只看到了 mathe() 的前三贴。
>>其中使用的算法正是我程序中使用的算法，不过我稍有优化。计算的时候不超过 5 ^ 5 项，而不是 5 ^ 6 项。
看了一下,对于F(*)在,L=18时,F(k,m)在k>=6时只有第一项F(k,0)非0,也就是说,
对于长度超过5^6项时,结果同起点没有关系了:)
理论可以证明,在k>=(L+1)/3时,F(k,m)中只有第一项F(k,0)非0 :) 这个可以用来提高一些速度.

我后来写的程序使用了另外一项优化,在对于N>>L时比较有用(N表示X的长度)
我们知道,上面的计算过程中先求出了函数
F(k)[x]=F(k,0)+F(k,1)x+F(k,2)x^2+...+F(k,L-1)*x^(L-1) (mod 5^L)
而且必然有5^s|F(k,s),可以发现F(k,L-1)总是0
而且对于k>=L-1, F(k)[x]=-1.

然后我们记G(k)[x]=F(k)[x-1]*F(k)[x-2]*...*F(k)[x-u] (where u=x%5)
如果x%5==0,那么G(k)[x]=1.  (F(0)[x]需要定义成x+1,这样公式就统一了)

对于N位大整数X,我们记
X=x0+x1*5+...+x(N-1)*5^(N-1)
X0=X=x0+x1*5+...+x(N-1)*5^(N-1)
X1=[X0/5]=x1+x2*5+...+x(N-1)*5^(N-2)
X2=[X1/5]
...
X(N-1)=[X(N-2)/5]=x(N-1)
我们可以看出
程序Get_Mult_No_5(Xi)
计算的就是
G(N-1-i)[Xi]*G(N-2-i)[X(i+1)]*G(N-3-i)[X(i+2)]*...G(0)[X(N-1)]
将上式对所有的i相乘,变成
G(0)[X(N-1)]G(1)[X(N-1)]...G(N-1)[X(N-1)]*
G(0)[X(N-2)]G(1)[X(N-2)]...G(N-2)[X(N-2)]*
...
G(0)[X(0)]

然后我们可以定义
H(0)[X]=F(0)[X]
H(1)[X]=H(0)[X]*F(1)[X]
...
H(N-1)[X]=H(N-2)[X]*F(N-1)[X]
J(k)[x]=H(k)[x-1]*H(k)[x-2]*...*H(k)[x-u] (where u=x%5)
如果上面的函数H(k)都可以事先算好,那么最终的乘积就可以表示成
N个J(k)函数的乘积,可以直接通过函数H算出.
而对于函数H(*)[X],我们同样只需要计算出前面L-1项,后面每一项都只是前面一项
乘上-1就可以了.
通过这种方法,我们可以将最后一步需要计算O(L*N)个多项式的值的过程降低为计算O(N)个多项式的值的过程

to   xstring (麻雀)
你的程序BUG！！！！！！！！！！！！！！！！！！！！！！！！！

在http://www.ysiw.info/T18/T18.html上下的“版本3”程序

运行时，输入某些数字出错 ！！！！！！！！！

经过运行，发现出错的数字 N 都有以下规律：

    N = a (mod 10^18) ,

按错键盘，没写完就发了出去 *^_^*

N = a (mod 10^18) , 1<=a<=8

即当N的后18位是17个0加一个1至8之间的数字时，程序运行错误。

我的机子的CPU是Celeron(R)2.4G

另外，为了测试随机超大整数的表现，
我选择了 N=Fibonacci(200000) (有 41,798 位)，计算最后 18 位，
用mathe()的程序用了 25.437 秒；而xstring(麻雀)因没有提供图形化界面的程序不便对比测试。

估计效率的瓶颈是在进制转换上，以前mathe曾与我有过交流讨论。

//========================================================

Q: 在 N>>L 时，如何快速计算末尾“0”的数目？

我测试了一下，N=Fibonacci(200000) 如果按普通的 [N/5]! 需花费 17.341 秒，
而如果改用如下算法则仅需 0.068 秒，下面给大家介绍这种快速算法：

1、将 N 转化为 5 进制的数N{5}（该过程有O(nlogn)算法），
2、将 N{5} 各个数位上的数累加，即为 sum，
3、计算 ( N - sum ) / 4 即可。

注：HugeCalc V5.x 的 ConvertToCarry() 函数就是经过优化的，时间开销为O(nlogn)算法.

to mathe():
    你的帖子写得都比较专业而易懂，大家从中受益非浅！
    但最近那篇因涉及了很多符号读起来比较费力，也许因为问题本身比较复杂，算法具有很高的技巧性。可否请你将这段算法的源代码也贴出来（因为阅读源代码有时比公式来得更直接些）? 谢谢！

呵呵,还是我写的东西的易读性不好.
修改1:
对于F(0,*)初始化的修改
    for(i=0;i<L;i++)mpz_set_si(F(0,i),0);
    mpz_set_si(F(0,0),1);mpz_set_si(F(0,1),1);
修改2:
    在F(i,j)计算都完成后,添加:
    for(i=1;i<L;i++){
	PolyMul(fun+i*L,fun+(i-1)*L, tmp4, L, fiveL);
	for(j=0;j<L;j++)mpz_set(F(i,j),tmp4[j]);
    }
在这一轮计算后,F(i,j)里面保存的值就变成上面帖子中提到的H函数了.
现在我们可以使用不同的Get_Mult_No_5(Xi)函数了:
void Get_Mult_No_5_2(mpz_t N,
		int L,
		mpz_t *fun,
		mpz_t five_L,
		mpz_t *five_i,
		mpz_t result,
		int step)
{
    mpz_t t1,t2,t3;
    int i,u;
    mpz_init(t1);mpz_init(t2);mpz_init(t3);
    mpz_fdiv_qr(t1,t2,N,five_L);
    u=mpz_mod_ui(t3,t2,5);
    if(step>=L&&(u&1)&&((step-L+1)&1)){
	mpz_sub_ui(result,five_L,1);
    }else{
        mpz_set_ui(result,1);
    }
    for(i=0;i<u;i++){
	mpz_sub_ui(t2,t2,1);
	Ft4(fun,t2,step>=L?L-1:step,five_L,L,t1);
	mpz_mul(t3,result,t1);
	mpz_mod(result,t3,five_L);
    }
    mpz_clear(t3);
    mpz_clear(t2);
    mpz_clear(t1);
}

请问新的 Get_Mult_No_5_2() 最后一个入参“step”的意义？在 Get_m5L() 函数中如何计算该变量？

j=0;
    while(mpz_sgn(N)>0){
	Get_Mult_No_5_2(N,L,fun,fiveL,tmp1,tmp,j++);
	mpz_mul(tmpp,m5L,tmp);
	mpz_mod(m5L,tmpp,fiveL);
	mpz_fdiv_q_ui(N,N,5);
	mpz_add(count_5,count_5,N);
    }

我的程序 FactorialTail.exe 终于完成了！

算  法：由 math() 提供，在其基础上再稍作优化改进；
算法库：直接由 HugeCalc.dll 提供支持（未自定制特殊版本的大数运算）；
输入项：阶乘数N、计算非零尾数位数T
输出项：规则化N、非零尾数L、尾部“0”串长度Z，以及它们值的各自10进制位数，计算总耗时微秒数；
范  围：0 <= N <= 10^100000000；1 <= T <= 255
文件大小：25 KB


经测试，计算N!最后18位非零尾数，
当N的位数在9位左右，基本上是楼主程序的 2~3 倍；
如果N更大，速率比估计将会更大，

比如（P4 1.7GHz / 256M RAM / WinXP）,
T18(N=  100!)= 377666050264399872, timer =    5452 us (N 的位数=158)
T18(N= 1000!)= 594893728159891456, timer =   92728 us (N 的位数=2568)
T18(N=10000!)= 534459594029924352, timer = 1642684 us (N 的位数=35660)

T18(10^10000-1)=979819424592166912, timer=143343us (N 的位数=10000)
T18(N=Fibonacci(200000))=539132392670232576, timer = 2480775 us (N 的位数=41798)

将 k 个 123456789 连续数字组成的数记作 N9(k)，如 N9(3) = 123456789123456789123456789;
下面测试计算最后99位非零数字的耗时：
T99(N9(1)), timer =  2720 us
T99(N9(2)), timer =  4746 us
T99(N9(3)), timer =  7766 us
T99(N9(4)), timer = 11689 us
T99(N9(5)), timer = 14893 us

近期我对 HugeCalc 作了一些小的升级，版本为 5.0.0.1，预定在下周正式发布，
本程序 FactorialTail.exe 将一同发布（也许同时发布修正后的源代码），敬请关注！

呵呵,看起来同我的程序相比较快了10倍:)
不过直接同楼主的程序相比意义不是很大,因为如果限定只计算18位,应该还有很多优化机会,
基本上可以不使用大数运算了

这个擂台摆得非常好，让大家从中学到了不少东西，受益匪浅！

楼主 xstring (麻雀) 与 mathe() 对本问题均有非常深入的研究，且均为算法界的高人，
这一点可从先前类似的帖子可以证明。

我是半途插进来的，目的只是为了实地检验一下 HugeCalc 的效率与可操作性，
好在是 mathe() 公开了全部算法（在此深表感谢），
所以我基本上就“望文生义”的“直译”过来，然后再作算法和代码上的适当优化而已。

楼主的大数算法库先由 GMP 提供，转而自已实现，两者之间可以有针对性优化的地方肯定非常多，
现在我用通用的 HugeCalc.dll 即可赶超楼主的程序，说明楼主的程序在算法上还有优化的空间。

本人优化的大致历程为：
1、确定选择十进制类 CHugeInt 还是 CHugeIntX？因本问题本身即为十进制的问题，所以选择 CHugeInt；
2、直译mathe()公开的源代码（原来的是Linux下，并调用GMP高精度算法库）；
3、优化代码：将循环体内的计算尽可能提到循环体外面；有指针代替数组；
4、优化算法：使用高效的进制转换算法、充分运用数论知识简化计算；
5、优化过程：将部分变量修改为静态，避免重复计算；
6、将控制台程序转化为图形化界面。

请大家下载如下文件：http://maths.diy.myrice.com/download/HugeCalcV5001.rar （414 KB）
该压缩包为 HugeCalc V5.0.0.1 精简尝鲜版（仅3天内有效），
压缩包内含我优化后的程序和全部源代码，欢迎大家指正。

另，to xstring (麻雀) 与 mathe():
    我计划将二位的程序一并发布在新版中，不知妥否？

That's NO PROBLEM