数学与计算机碰撞的地方

计算阶乘的 18 位非零尾数

摘要

xstring于2005年在csdn挑战 ,对于一个给定的大整数,如何快速计算这个整数阶乘去掉尾零后的最低18位

最后mathe给出了一种时间复杂度为O(L^3 \log(L)^2 +L^3 \log(L)T +T^2)的有效算法,其中L代表计算的非零尾数位数, T代表输入整数的位数。
gxqcn还利用他的HugeCalc 对算法重新实现并进行性能优化,比mathe用gmp 实现的代码更快了一个数量级。
贴子内容是mathe对https://blog.csdn.net/mathe/article/details/1132404 的备份。

详细内容

xstring首先发问: 阶乘不用解释吧,阶乘的 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
比如下图中结果
factail

mathe给出了如下分析:

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

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

对于N!,我们可以将5的倍数和不是5的倍数的数分开处理:
所有的5的倍数构成5^{[\frac{N}5]}\times[\frac{N}5]!
如果我们已经得到所有不是5的倍数的乘积关于5^{18}的余数,记为g(N),
所求的结果记为f(N),那么就有f(N)=g(N)f([\frac{N}5]) (\mod 5^{18}).
所以关键是求g(N),也就是1到N之间所有不是5的倍数的数的乘积关于5^{18}的余数,
也就是最前面列出的。
后面只讨论如何计算:
g(N)=1\times 2\times 3\times 4\times 6\times 7\times 8\times 9\times 11\times 12\times \dots\times 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) \times T+T\times\log(T)),
使用的空间复杂度为O(L^3+L^2 \log(L)+T)
计算过程中,我们总假设一个长度为L的数占用的空间为O(L),两个长度为L的数相乘花费的时间为O(L\times\log(L)),而相加的时间复杂度为O(L)。

方法如下:
i)计算杨辉三角形 关于5^L的余数,并保存,记为C_U^V (\mod 5^L) 其中1\le U\le L, 0\le V\le U
空间复杂度为O(L^3 (L^2个长度为L的整数),需要花费O(L^2)次加法,每次加法需要O(L)的时间,所以时间复杂度也是O(L^3).
ii)记F(k,x)=(x\times 5^k+1)(x\times5^k+2)(x\times5^k+3)(x\times5^k+4)(x\times5^k+6)\dots(x\times5^k+5^k-1) (\mod 5^L)
其中连乘中所有5的倍数全部被去掉
那么F(1,x)=24+50\times 5\times x+35\times 5^2\times x^2+10\times 5^3\times x^3+5^4\times x^4 (\mod 5^L).
假设F(n,x)=a_0+a_1\times 5\times x+a_2\times5^2\times x^2+\dots+a_{L-1}\times 5^{L-1}\times x^{L-1} (\mod 5^L),
那么F(n+1,x)=F(n,5x)\times F(n,5x+1)\times F(n,5\times x+2)\times F(n,5x+3)\times F(n,5x+4) (\mod 5^L).
首先将每个F(n,5x+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\le k\le L.
而保存所有的F(k,x),需要花费的空间是O(L^2 \log(L)) (L\log(L)个长度为L的数).
而对于F(k,x), k\gt L,我们知道F(k,x)=1,不需要再计算了.
iii)对于f(N,L),f(N,L) (\mod 5^L)=g(N,L)f([\frac{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([\frac{N}5],L),同样的,如果我们把N看成T位5进制数(事先转化一个数为5进制只需要T\log(T)的时间),那么\frac{N}5是T-1位5进制数,所以通过同样算法,可以再花费O(L^3 \log(L))的时间递归到T-2位数,…,这样总共经过T步后就会得出最终结果.
这个递归过程中总共花费时间最多为 O(T\times 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(超过\frac{10L}3肯定够了,而对于这么小的N,直接计算乘积就可以了),f(N,L) (\mod 2^L)总是0,所以通过中国剩余定理就可以计算出f(N,L) (\mod 10^L).
对于现在的计算机基本上在L\le 100而且T\le 1000 (也就是N\le 10^{1000})时不会有问题.
附件程序 使用了GMP库,
但是由于下面的程序中,对于输入的T位整数N,没有通过事先转化为5进制的方法,所以实际上花费的时间复杂度为 O(L^3 \log(L)^2 +L^3 \log(L)T +T^2).

Next Post

Previous Post

Leave a Reply

© 2020 数学研发论坛

Theme by Anders Norén