用户：Peter Luschny/FactorialFunction

近似值
阶乘函数。

关键词：阶乘函数、伽马函数，斯特林公式，斯蒂尔特杰斯公式，连分数近似。

与序列有关：

公式	分子	分母
斯特林	A001163号	A001164号
德莫伊夫尔	A182935号	A144618号
斯蒂尔切斯	A005146号	A005147号
兰索士	A090674号	A090675号
尼姆（Nemes）	A181855号	A181856号
复仇女神G	A182912号	A182913号
韦迈尔	A182916号	A182917号
高斯珀	A182919号	A182920号
/新建/	A182914号	A182915号
/新建/	A277000型	A277001型
/新建/	A277002型	A277003型

介绍

n的阶乘定义为

${\显示样式n！=\int_{0}^{\infty}\x^{n}\e^{-x}}$.

此函数最著名的近似值是De Moivere-Stirling公式

$显示样式S（n）={\sqrt{2\pi}}\n^{n+1/2}\e^{-n}}$.

为了加速该公式的渐近收敛，提出了几个修正函数σ（n）。

$｛\displaystyle S^｛*｝（n）＝｛\sqrt｛2\pi｝｝\n^｛n+1/2｝\e^｛-n｝\\sigma（n）｝$.

例如，最简单的是σ（n）=（1+1/（12n））。

我们想在这里看到的公式是从一个不同于De Moivre-Stirling公式的基础结构开始的。基础结构更加均匀；它不混合n个和n个+ 1/2.

$显示样式T（n）={\sqrt{2\pi}}\（n+1/2）^{n+1/2}\e^{-n-1/2}}$.

同样，这个公式必须进行调整才能变得有用。这里我们提出了一个修正函数τ（n）它保持了半移位论点的基本结构。

${\显示样式T^{*}（n）={\ sqrt{2\pi}}\（n+1/2）^{n+1/2}\e^{-n-1/2}\\tau（n+1/2）^}$

下面将讨论此校正函数的详细信息。这个据我所知，公式是第一次引入的作者于2007年1月9日在其主页上发表[1]。事实证明，新公式的效率可与斯蒂尔特杰斯的著名公式媲美阶乘函数的连分式逼近。

斯蒂尔特杰斯近似

1894年，T.J.Stieltjes提出了一种有趣的阶乘函数方法。

${\显示样式n！=\exp\left（\log（2\pi）/2-n+\ left（n+{\frac{1}{2}}\ right）\ log（n）+p（n）\right）。}$

对于p（n），Stieltjes发展了一个连续分数

$｛\displaystyle p（n）＝｛\crac｛a_｛0｝｝｛n+｛\crac｛a_｛1｝｝｛n+｛\crac｛a_｛2｝｝｛n+｛\crac｛a_｛3｝｝｝｛n+\dots$

前几个系数a_n个在中A005146号和A005147号.

n个	0	1	2	三	4
一n个	1/12	1/30	53/210	195/371	22999/22737

首先我们必须计算所有涉及的常数。这远不是小事一桩。斯蒂尔特杰斯自己对此呻吟道：“计算是一个过程……这是一个非常复杂的过程”。使用Rutishauser和Akiyama-Tanigawa的算法，我们现在可以计算系数比斯蒂尔特杰斯容易得多。可以找到有关此的更多信息在这里. 一个简短的Maple程序是在这里.

新近似值

我们还将考虑一个连分数，这次是τ（n）

${\displaystyle\tau（n）=n\times{\cfrac{c{0}}{n+{\cfras{c{1}}{n+{\frac{c}2}}{n+{\cfrac{c{3}}{n\ddots}}}}{}}}$

前几个系数c_n个在中A182914号和A182915号.（顺便注意，45360是最小的整数，正好有100个除数(A163816号).)

n个	0	1	2	三	4
c（c）n个	1	1/24	3/80	18029/45360	6272051/14869008

现在的问题是如何计算这些系数。当然这里的数值给出了近似值的精确定义。

新公式的系数

c：=数组（0..20，[1/24]）：w：=“w”：近似值：=proc（deg，n）局部n，A，k，p；N：=N+1/2；A:=N；对于k从deg到-1到0的do A:=N+c[k]/A od；p:=N^2/A；ln（2*Pi）/2+N*（ln（p）-1）；exp（%）结束：coeffCF:=proc（度）全局w，c；局部R，k；R:=1..0美元；如果w=“w”，则R：=$1..度；w:=度elif deg>w，然后R：=$w+1..deg；w：=度fi；对于R do中的kc[k]：=a；转换（asympt（近似（k，m）-m！，m、 2*k+4），多项式）：c[k]：=求解（%=0，a）；打印（`coeff:`，k，c[k]）；od端：

假设系数已经计算完毕，使用函数“近似”，否则调用函数“coeffCF”第一。参数“deg”是所请求公式的阶数。例如：

系数CF（5）；近似值（3,10）：evalf（%，20）；3628800.0000007818727近似值（5,10）：evalf（%，20）；3628800.0000000002108

Wehmeier公式

最后一节中使用的技术显然也可以应用于其他阶乘近似的基本结构。我学会了这种方法2006年3月，来自Usenet group de.sci.matematik的Stefan Wehmeier数据源1. Wehmeier用它推导出Gosper近似到n！的推广！。数字平方米.他接着说：

W:=数组（0..20）：W:=“W”：wehmeier:=proc（deg，n）局部A，k；A:=n；对于k从deg到-1到0的do A:=A+W[k]/n^k od；sqrt（2*Pi*A）*n^n*exp（-n）结束：coeffW：=程序（deg）全局w，w；局部R，k；R：=$0.0；如果w='w'，则R:=$0..deg；w:=度elif deg>w，然后R：=$w.deg；w：=度fi；对于R do中的kW[k]：=a；转换（asympt（wehmeier（k，m）-m！，m、 k+2），多项式）：W[k]：=求解（%=0，a）；打印（`coeff:`，k，W[k]）；od端：

使用与上面相同的调用约定，我们得到

系数W（5）；韦梅尔（3,10）：评估值（%，20）；3628799.9727503853301wehmeier（5,10）：评价（%，20）；3628800.0002087858324

前几个系数W_n个在中A182916号和A182917号.

n个	0	1	2	三	4
W公司n个	1/6	1/72	-31 / 6480	-139 / 155520	9871/6531840

Nemes-G公式

本月（2011年3月），GergöNemes发布了一个新的近似公式，“伽玛函数的更精确近似”，《泰国数学杂志》，第9卷，（1），p。21-28. Nemes和Wehmeier一样，从Gosper的近似值开始。但随后他继续具有以下基本结构：A（n+1/4）sqrt（2Pi（n+1/6））（n/e）^n。同样，此扩展可以用同样的技术很容易导出以上使用；只需稍作调整。

H:=数组（0..20）：w:=“w”：nemesG：=程序（deg，n）局部A，n，k；A： =1；N：=N+1/4；对于k，从度乘-1到2 do A:=A+H[k]/N^k od；sqrt（2*Pi*（n+1/6））*n^n*exp（-n）*A端：系数H：=proc（deg）全局w，H；局部R，k；R:=2..0美元；如果w=“w”，则R：=$2.deg；w:=度elif deg>w，然后R:=$w.deg；w：=度fi；对于R do中的kH[k]：=a；转换（asympt（nemesG（k，m）-m！，m、 k+1），多项式）：H[k]：=求解（%=0，a）；print（`coeff:`，k，H[k]）；od端：

使用与上面相同的调用约定，我们得到

系数H（5）；nemesG（3,10）：评价（%，20）；3628800.4061837742630 nemesG（5,10）：评价（%，20）；3628799.9981087549165

前几个系数H_n个在中A182912号和A182913号.

n个	0	1	2	三	4
H（H）n个	1	0	1/144	-1 / 12960	-257 / 207360

高斯珀公式

甚至还有一种更简单的方法来推广高斯珀近似。基本公式是A（n）sqrt（2Pi（n+1/6））（n/e）^n。与Nemes公式中的Ansatz相同，但参数没有移位。

G:=数组（0..20）：w:=“w”：gosper:=proc（deg，n）局部A，k；A:=1；对于k从deg乘以-1到2的do A:=A+G[k]/n^k od；A*sqrt（2*Pi*（n+1/6））*（n/exp（1））^n结束：coeffG：=proc（deg）全局w，G；局部R，k；R：=$2.0；如果w=“w”，则R：=$2.deg；w:=度elif deg>w，然后R：=$w.deg；w：=度fi；对于R do中的kG[k]：=a；convert（asympt（gosper（k，m）-m！，m、 k+1），多项式）：G[k]：=求解（%=0，a）；打印（`coeff:`，k，G[k]）；od端：

使用与上面相同的调用约定，我们得到

系数G（8）；鹅（3,10）：评估值（%，20）；3628799.9289952224556 小鹅（5,10）：evalf（%，20）；3628800.0001794192645

前几个系数G_n个（系数（尚未）在OEIS数据库中）。

n个	0	1	2	三	4
G公司n个	1	0	1/144	-23 / 6480	5/41472

斯特林公式

斯特林近似是阶乘函数最著名的近似-出于历史原因。如今，没有理由将这些公式用于计算目的。性能更好的公式是已知的。

h:=proc（k）选项记忆；局部j`如果`（k=0,1，（h（k-1）/k-加（（h（k-j）*h（j））/（j+1），j=1..k-1））/（1+1/（k+1））结束：coeffStirling:=proc（n）选项记忆；h（2*n）*2^n*pochhammer（1/2，n）结束：

n个	0	1	2	三	4
S公司n个	1	1 / 12	1/288	-139 / 51840	-第571/2488320号

一个小基准

首先，我们以更紧凑的方式重写公式。

新：=proc（n）局部C，n，A；N：=N+1/2；C:=[1/24,3/8018029/453606272051/14869008]；A:=N^2/（N+C[1]/（N+C[2]/（N+C[3]/（N+C[4]/N）））；sqrt（2*Pi）*exp（N*（ln（A）-1））结束：wehmeier:=proc（n）局部W，n，A；N:=N+1/6；宽度：=[1/72，-31/6480，-139/1555209871/6531840]；A:=N+W[1]/N+W[2]/N^2+W[3]/N^3+W[4]/N^4；sqrt（2*Pi*A）*n^n*exp（-n）结束：nemes:=proc（n）局部n，H，A；N：=N+1/4；H:=[1/144，-1/12960，-257/207360，-53/2612736]；A:=1+H[1]/N^2+H[2]/N^3+H[3]/N^4+H[4]/N^5；A*sqrt（2*Pi*（n+1/6））*n^n*exp（-n）结束：gosper:=proc（n）局部A，G；G:=[1/144，-23/6480,5/414724939/6531840]；A:=1+G[1]/n^2+G[2]/n^3+G[3]/n^4+G[4]/n^5；A*sqrt（2*Pi*（n+1/6））*n^n*exp（-n）结束：斯特林：=proc（n）局部A，S；S：=[1/12,1/288，-139/51840，-571/2488320]；A:=1+S[1]/n+S[2]/n^2+S[3]/n^3+S[4]/n^4；A*sqrt（2*Pi）*n^（n+1/2）*exp（-n）结束：

为了衡量公式的性能，我们计算了精确的十进制数字。

edd:=程序（a，t）评估（-log[10]（abs（1-a/t）），60）结束：对于[100100010000]中的V，执行VF:=V！：printf（“斯特林%.5d！%6.1f\n”，V，edd（斯特林（V），VF））；printf（“Nemes-G%.5d！%6.1f\n”，V，edd（Nemes（V），VF））；printf（“Wehmeier%.5d！%6.1f\n”，V，edd（Wehmeier（V），VF））；printf（“Gosper%.5d！%6.1f\n”，V，edd（Gosper（V），VF））；printf（“新%.5d！%6.1f\n”，V，edd（新（V），VF））；日期：

精确小数 数字	100!	1000!	10000！
斯特林	13.1	18.1	23.1
复仇女神-G	15.2	21.2	27.2
韦迈尔	15.9	21.9	27.9
高斯珀	17.5	23.1	29.1
新建	21.5	30.5	39.5

在更大的范围内比较Stirling、Gosper和新公式的精确小数位数是如下图所示。

蓝色涡流函数中的峰值是由于以下事实高斯珀近似精确到x=67.0033148435486248。。。

另请参阅

[1] 有关阶乘公式的更多信息，请访问我的主页.

[2] （2016年9月23日增补）关于阶乘（和伽玛函数）到今天是：王伟，伽马函数近似的统一方法,J.数论（2016）