阶乘有多少位？

wayne

一个数的阶乘有多少位，有好的算法吗 比如四百万的阶乘有 24671066位

KeyTo9_Fans

这是什么公式来着？ $n! = \sqrt(2\pi(n+1/6+1/(72n))) * (n/e)^n$ 其中$1/6$和$1/(72n)$这两项是我自己凑的。 据测试，多了那两项可以修正结果，以至于有效数字比原来多了2倍。 两边取对数，得 $log(n!) = log(2\pi(n+1/6+1/(72n)))/2 + n*(log(n)-log(e))$ 把$n=4000000$代进去算，得到 $log(4000000!) = 24671065.737818781141945821796116$ 大概这几位数字 $24671065.7378187811419458217961$ 是有效的吧，后面的$...16$可能就不精确了。 即便如此，我们已经可以得到精确值的位数，以及精确值的前面若干位有效数字。

medie2005

Stirling's formula

数学星空

斯特林级数公式: http://mathworld.wolfram.com/StirlingsSeries.html

wayne

哦，谢谢各位了

wayne

这是什么公式来着？ $n! = \sqrt(2\pi(n+1/6+1/(72n))) * (n/e)^n$ 其中$1/6$和$1/(72n)$这两项是我自己凑的。 $log ... KeyTo9_Fans 发表于 2010-1-20 00:01

这个KeyTo9_Fans 公式虽然在逼近n!上还可以有改进的空间， 但用来求阶乘的位数，恰到好处，已经是简洁到不能再简洁的了！！！

wayne

哦，用来求大于1的数的阶乘的位数， Stirling's formula 已经足够了 $\lceil n(\log _{10}n-\log _{10}\e)+\frac{1}{2}\log _{10}(2\pi n)\rceil $ 而逼近n!,更好的是： $n! = \sqrt(2\pi(n+1/6+1/(36.4+70.3n))) * (n/e)^n$

KeyTo9_Fans

关于根号里第3项的分母，楼上是怎么得到$36.4+70.3n$的？ 我通过计算1000000000的阶乘的准确值(仅保留100位有效数字)，然后倒过来推根号里面的多项式，可以得到 $n! =sqrt(2\pi p)*(n/e)^n$ 其中 $p=n+1/6+a/n+b/n^2+c/n^3+d/n^4+e/n^5+f/n^6+g/n^7$ $a = 1/72$ $b = -31/6480$ $c = -139/155520$ $d = 9871/6531840$ $e = 324179/1175731200$ $f = -8225671/7054387200$ $g = -69685339/338610585600$ 到这里，已经把100位有效数字全部用完了，无法再继续考察下一项了。 要得到更多的项，需要计算比1000000000更大的阶乘的准确值，并保留更多的有效数字。 把$n=16$代进去算，得到的结果为 $20922789887999.748608502448646388$ 四舍五入到整数，与真实值 $16! =20922789888000$ 完全相同。 因为根号里的多项式有9项，所以上述公式的大致趋势是$n$每增大10倍，$n!$的有效数字增加9个。

数学星空

wayne

用C打印大数的阶乘，终于完成了，刚出锅，还是热的，高手们看看还有没有可以改进的地方

2. #include<stdio.h>
3. #include<stdlib.h>
4. #include<math.h>
6. int main(){
7. int i,j,*f,tmp,c=0;
8. long int n,bits;
9. const double PI=2*asin(1.0),E=exp(1.0);
11. scanf("%ld",&n);
12. bits=(long)ceil(n*(log10(n)-log10(E))+log10(2*PI*n)/2);
13. printf("there are %ld digits\n%ld!=",bits,n);
14. f=(int*)calloc(bits,sizeof(int));
15. f[0]=1;
16. for(i=2;i<=n;i++){
17. for(j=0;j<bits;j++){
18. tmp=f[j]*i+c;
19. c=tmp/10;
20. f[j]=tmp%10;
21. } }
23. for(i=bits-1;i>=0;i--) printf("%d",f[i]);
24. printf("\n");
25. return 0;
26. }

复制代码