如何用gamma函数求π的阶乘?

snowblink

如题

战巡

直接套进去呗
\[\pi!=\Gamma(\pi+1)=\int_0^{+\infty}x^{\pi}e^{-x}dx\]
至于精确值就拉倒吧，顶多求个近似

hbghlyj

Here's a simple approximation for the gamma function, proposed in 2007 by Gergő Nemes. (See the wikipedia page on Stirling's approximation).\[\Gamma (z)\approx {\sqrt {\frac {2\pi }{z}}}\left({\frac {1}{e}}\left(z+{\frac {1}{12z-{\frac {1}{10z}}}}\right)\right)^{z},\]
In Javascript:

1. const π = Math.PI;
2. function gamma(z) {
3.   return Math.sqrt(2*π/z) * Math.pow((1 / Math.E) * (z + 1 / (12 * z - 1 / (10 * z))), z);
4. }
5. gamma(π+1)                // Output: 7.18807906494771
6. gamma(π+2)/(π+1)          // Output: 7.188081465246312
7. gamma(π+3)/(π+2)/(π+1)    // Output: 7.18808220428182

复制代码

hbghlyj

Lanczos approximation
numericana.com/answer/info/godfrey.htm
The Lanczos approximation consists of the formula
\[\Gamma(z+1) = \sqrt{2\pi} {\left( z + g + \tfrac12 \right)}^{z + 1/2 } e^{-(z+g+1/2)} A_g(z)\]
for the gamma function, with
\[A_g(z) = \frac12p_0(g) + p_1(g) \frac{z}{z+1} + p_2(g) \frac{z(z-1)}{(z+1)(z+2)} + \cdots\]
Here $g$ is a real constant that may be chosen arbitrarily subject to the restriction that $\text{Re}(z+g+\frac12)>0$.
The coefficients are given by
\[p_k(g) = \frac{\sqrt{2\,}}{\pi} \sum_{\ell=0}^k C_{2k+1,\,2\ell+1} \left(\ell - \tfrac{1}{2} \right)! {\left(\ell + g + \tfrac{1}{2} \right)}^{-(\ell+1/2)} e^{\ell + g + 1/2 }\]
where $C_{n,m}$ represents the $(n,m)$th element of the matrix of coefficients for the Chebyshev polynomials, which can be calculated recursively from these identities:
$$\begin{aligned}
    C_{1,\,1} &= 1 \\[5px]
    C_{2,\,2} &= 1 \\[5px]
  C_{n+1,\,1} &= -\,C_{n-1,\,1}  & \text{ for  } n &= 2, 3, 4\, \dots \\[5px]
C_{n+1,\,n+1} &= 2\,C_{n,\,n}    & \text{ for  } n &= 2, 3, 4\, \dots \\[5px]
C_{n+1,\,m+1} &= 2\,C_{n,\,m} - C_{n-1,\,m+1} & \text{ for  } n & > m = 1, 2, 3\, \dots
\end{aligned}$$

hbghlyj

Lanczos approximation in Javascript

1. const π = Math.PI;
2. function gamma(z) {
3.     const g = 7, C = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
4.     if (z < 0.5) return π / (Math.sin(π * z) * gamma(1 - z));
5.     else {
6.         z -= 1;
7.         var x = C[0];
8.         for (var i = 1; i < g + 2; i++)
9.         x += C[i] / (z + i);
10.         var t = z + g + 0.5;
11.         return Math.sqrt(2 * π) * Math.pow(t, (z + 0.5)) * Math.exp(-t) * x;
12.     }
13. }
14. gamma(π+1)                // Output: 7.188082728976035
15. gamma(π+2)/(π+1)          // Output: 7.188082728976032
16. gamma(π+3)/(π+2)/(π+1)    // Output: 7.188082728976043

复制代码

stackoverflow.com/questions/15454183