计算阶乘的 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

代码以及预编译过的可执行程序可以在下面链接下载到:
http://download.csdn.net/source/270547

下面介绍原理：

分析这个问题,可以发现主要在于计算
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 <stdio.h>
#include <stdlib.h>
#include <gmp.h>

#define TRI