计算器上阶乘如何实现?

阶乘一般写为 n!  ,如果为正整数则,n!=n*(n-1)*(n-2)...*3*2*1;
计算器上支持小数的阶乘运算,一般内部为gamma函数(Γ函数)实现, gamma(n+1)=n!
如果只支持到双精度浮点(1e-16)以内,那么一般用近似公式就可实现,如: http://www.rskey.org/CMS/index.php/the-library/11;
gamma函数,wiki上的介绍: http://en.wikipedia.org/wiki/Gamma_function

我现在想实现更高的精度(上千位),并支持更大的n输入值
当n为正整数数时,n!=n*(n-1)*(n-2)...*3*2*1 ; 可以用统计表达式质因数个数的方法(素数筛法,并快速获得因数个数),这样n在千万级别速度都不是太大问题;
( 当n<0.5时 , gamma(n)=PI/(gamma(1-n)*sin(PI*n)) ,这样,就只剩下>0.5的正实数的通用求解方法 (当输入为正整数+0.5 时, gamma(n+0.5)= (2n)!*PI^0.5/(n!*4^n); ) )

当n为正实数时,我借鉴了这篇文章,gamma函数表达为一个级数和一个连分数:
http://www.cws.net.cn/Journal/slxb/199912/13.html
其中,连分数的收敛性很差,运行慢;
我从该文章整理的计算公式:  Γ(x)=e^(x*ln(y)-(y))*[S(x,y)+U(x,y)]
          S(x,y)=y^0/x + y^1/(x*(x+1)) + y^2/(x*(x+1)*(x+2)) + ...
          U(x,y)=a0/(Q0*Q1)-a0*a1/(Q1*Q2)+a0*a1*a2/(Q2*Q3)- ...  
               a0=1; a(k)= ((k+1)>>1)-x ,k%2==1; a(k)=((k+1)>>1), k%2==0;
               b0=0; b(k)= y ,k%2==1; b(k)=1,k%2==0;
               Q(-1)=0; Q0=1; Q(k)=b(k)*Q(k-1)+a(k-1)*Q(k-2);

        其中(作者建议)y取,y=x+1

但Γ(x)当x==172时,计算结果不对(文章作者也应该只是在双精度内验证过),这时修改y为其他值,还可以继续计算,我寻找了更多的可能y值,但当x==244时,适应性最强的y=x+718也失败; 有没有可能得到其他合适的y的计算方法?(根据U的收敛性分析?)

有没有其他合适的计算公式推荐?

16 个解决方案

#1

记得以前在书店翻过一本算法的书好像也讲了相关的问题，书名常用算法程序集
楼主可以翻翻里面有没有相关的东西

#2

这个可以使用大数模拟。就是自己开一个数组一个一个模拟过去，这个算法我写过如果只是开32位的话12！就溢出了....大数模拟的话随便多大！

#3

引用 2 楼  的回复:

这个可以使用大数模拟。就是自己开一个数组一个一个模拟过去，这个算法我写过如果只是开32位的话12！就溢出了....大数模拟的话随便多大！

我已有一个任意精度的大数类,现在需要的是一个合适的n!计算公式

#4

找到另一个gamma参考资料: http://mathworld.wolfram.com/GammaFunction.html

#5

用连分数的方法算是为了不丢精度还是为了速度？

引用 4 楼  的回复:

找到另一个gamma参考资料: http://mathworld.wolfram.com/GammaFunction.html

#6

不知道lz的运算范围有多大？不大的话可以*10的幂，转为大整数来算。

#7

引用 5 楼  的回复:

用连分数的方法算是为了不丢精度还是为了速度？

引用 6 楼  的回复:

不知道lz的运算范围有多大？不大的话可以*10的幂，转为大整数来算。

为了精度

浮点数小数位很长,保留上千10进制精度

#8

那个式子用到的ln，sin，pi，e，是否都会有精度的问题？如果只是算1000以内，个人感觉转为大整数算应该可以的。

引用 7 楼  的回复:

为了精度
浮点数小数位很长,保留上千10进制精度

#9

找到问题了,那篇文章没有问题, 我的代码处理double下溢出数(10^-309..10^-324)的bug;
另外,找到了一个更好地计算公式,计算速度更快:  http://en.literateprograms.org/Gamma_function_with_Spouge's_formula_(Mathematica)   不过其中用来控制精度误差的a值的计算公式总感觉有问题,调试中 

#10

这个比较全: http://bh0.physics.ubc.ca/People/matt/Doc/ThesesOthers/Phd/pugh.pdf

#11

我用易语言编写了个求阶乘的计算器，很简单啊。有需要的话，源码发你。一起交流，QQ：869321271
以下是源码：
.版本 2

.子程序 求阶乘, 双精度小数型
.参数 n, 整数型

.如果 (n ＝ 0)
    返回 (1)
.否则
    返回 (n × 求阶乘 (n － 1))
.如果结束


.子程序 _按钮3_被单击

标签3.标题 ＝ 编辑框3.内容 ＋ “的阶乘等于：” ＋ 到文本 (求阶乘 (到整数 (编辑框3.内容)))

#12

楼主的意思是用公式来计算？

#13

如果存1个表，1对1的，常用的也够了吧 

#14

楼主，你好我有一个计算N!的java程序，N<=10000以的内阶乘在1s以内是可以计算的。再大没试过，不过我忘记是怎么求的了，只记得是在数论上看过分解N!的式子。

#15

回复11--14 楼:
    我要计算的是输入参数范围数是实数域,而不仅仅是整数;计算结果的10进制精度需要上万位;  1楼我的说明里已经解决了大部分的问题,只是一些小问题和想寻找一个更好用的计算公式; 9楼我的回复里已经说明我已经解决了那个小问题,并说明自己使用了一个更好一些的计算公式;

#16

我看到楼上易语 言了，莫非现在有人用易语言编程？