スターリングの近似

スターリングの公式の厳密な証明にはオイラーの和公式、あるいは鞍点法といった複素解析の技法などを用いられることが多いが、初等的に導くことも可能である。まず階乗の対数を積分で近似する。logが凹関数であることから k-1<x<k (k=2,3,...) に対して

${\displaystyle \log k-(k-x)\{\log k-\log(k-1)\}\;<\;\log x\;<\;\log k-{\frac {1}{k}}(k-x)}$

これを k-1 から k まで積分して

${\displaystyle {\begin{aligned}&\log k-{\frac {1}{2}}\{\log k-\log(k-1)\}\;<\;\int _{k-1}^{k}\log x\,dx\;<\;\log k-{\frac {1}{2k}}\\&\int _{k-1}^{k}\log x\,dx+{\frac {1}{2k}}\;<\;\log k\;<\;\int _{k-1}^{k}\log x\,dx+{\frac {1}{2}}\{\log k-\log(k-1)\}\end{aligned}}}$

k=m+1,m+2,...,n に対して足し合わせると

${\displaystyle {\begin{aligned}\log {\frac {n!}{m!}}&=\sum _{k=m+1}^{n}\log k\\&>\int _{m}^{n}\log x\,dx+\sum _{k=m+1}^{n}{\frac {1}{2k}}\\&>\int _{m}^{n}\log x\,dx+{\frac {1}{2n}}-{\frac {1}{2m}}+\int _{m}^{n}{\frac {1}{2x}}\,dx\\&=(n+1/2)\log n-n+1/(2n)-(m+1/2)\log m+m-1/(2m)\\\log {\frac {n!}{m!}}&=\sum _{k=m+1}^{n}\log k\\&<\int _{m}^{n}\log x\,dx+{\frac {1}{2}}\{\log n-\log m\}\\&=(n+1/2)\log n-n-(m+1/2)\log m+m\end{aligned}}}$

${\displaystyle {\frac {n^{n+1/2}e^{-n+1/(2n)}}{m^{m+1/2}e^{-m+1/(2m)}}}<{\frac {n!}{m!}}<{\frac {n^{n+1/2}e^{-n}}{m^{m+1/2}e^{-m}}}}$

ここで ${\displaystyle \;a_{n}={\frac {n!}{n^{n+1/2}e^{-n}}}\;}$ と定めると

${\displaystyle {\frac {e^{1/(2n)}}{e^{1/(2m)}}}<{\frac {a_{n}}{a_{m}}}<1}$

m,n→∞のとき最左辺は1に収束するから、特に n=2m のとき

${\displaystyle \lim _{m\to \infty }{\frac {a_{2m}}{a_{m}}}=1}$

これとウォリスの公式の系 ${\displaystyle \;\lim _{n\to \infty }{\frac {4^{n}(n!)^{2}}{{\sqrt {n}}(2n)!}}={\sqrt {\pi }}\;}$ と比較すると、

${\displaystyle {\frac {4^{n}(n!)^{2}}{{\sqrt {n}}(2n)!}}={\frac {4^{n}(a_{n}n^{n+1/2}e^{-n})^{2}}{{\sqrt {n}}a_{2n}(2n)^{2n+1/2}e^{-2n}}}={\frac {a_{n}^{2}}{{\sqrt {2}}a_{2n}}}}$

${\displaystyle \lim _{n\to \infty }a_{n}=\lim _{n\to \infty }{\frac {4^{n}(n!)^{2}}{{\sqrt {n}}(2n)!}}{\frac {{\sqrt {2}}a_{2n}}{a_{n}}}={\sqrt {2\pi }}}$

を得る。

精度を改善するために a_n を評価する。

${\displaystyle {\begin{aligned}\log {\frac {a_{m}}{a_{m+1}}}&=\log m!-\left(m+{\frac {1}{2}}\right)\log m+m-\log(m+1)!+\left(m+{\frac {3}{2}}\right)\log(m+1)-m-1\\&=\left(m+{\frac {1}{2}}\right)\log \left(1+{\frac {1}{m}}\right)-1\\&=\left(m+{\frac {1}{2}}\right)\sum _{k=1}^{\infty }{\frac {(-1)^{k-1}}{k}}{\frac {1}{m^{k}}}-1\\&=\sum _{k=0}^{\infty }{\frac {(-1)^{k}}{k+1}}{\frac {1}{m^{k}}}+\sum _{k=1}^{\infty }{\frac {(-1)^{k-1}}{2k}}{\frac {1}{m^{k}}}-1\\&=\sum _{k=2}^{\infty }{\frac {(k-1)(-1)^{k}}{2k(k+1)}}{\frac {1}{m^{k}}}\end{aligned}}}$

${\displaystyle {\begin{aligned}\log {\frac {a_{n}}{\sqrt {2\pi }}}&=\sum _{m=n}^{\infty }\log {\frac {a_{m}}{a_{m+1}}}=\sum _{k=2}^{\infty }{\frac {(k-1)(-1)^{k}}{2k(k+1)}}\sum _{m=n}^{\infty }{\frac {1}{m^{k}}}\\&=\sum _{k=2}^{\infty }{\frac {(-1)^{k}}{2k(k+1)}}\left\{{\frac {1}{n^{k-1}}}+O\left({\frac {1}{n^{k}}}\right)\right\}\\&={\frac {1}{12n}}+O\left({\frac {1}{n^{2}}}\right)\end{aligned}}}$

従って

${\displaystyle n!\sim {\sqrt {2\pi n}}\left({\frac {n}{e}}\right)^{n}\exp \left({\frac {1}{12n}}\right)}$

オイラーの乗積表示によるガンマ関数の定義の対数をとり

${\displaystyle {\begin{aligned}\log \Gamma (z)&=\log {\big \{}(z-1)\Gamma (z-1){\big \}}\\&=\lim _{N\to \infty }(z-1)\log N+\sum _{n=1}^{N}{\big \{}\log n-\log(n+z-1){\big \}}\end{aligned}}}$

${\displaystyle f(n)=\log n-\log(n+z-1)}$ にオイラーの和公式を適用すれば

${\displaystyle {\begin{aligned}\log \Gamma (z)&=\lim _{N\to \infty }(z-1)\log N+\int _{n=1}^{N}f(n)\,dn+{\frac {1}{2}}{\big (}f(N)+f(1){\big )}\\&\qquad +\sum _{k=1}^{m}{\frac {B_{2k}}{(2k)!}}\left(f^{(2k-1)}(N)-f^{(2k-1)}(1)\right)+\int _{n=1}^{N}{\frac {B_{2m+1}(n-\lfloor n\rfloor )}{(2m+1)!}}f^{(2m+1)}(n)\,dn\\&=\lim _{N\to \infty }(z-1)\log N+{\bigg [}n\log n-n-(n+z-1)\log(n+z-1)+(n+z-1){\bigg ]}_{n=1}^{N}+{\frac {1}{2}}{\big \{}\log N-\log(N+z-1)-\log z{\big \}}\\&\qquad +\sum _{k=1}^{m}{\frac {B_{2k}}{(2k)(2k-1)}}\left\{{\frac {1}{N^{2k-1}}}-{\frac {1}{(N+z-1)^{2k-1}}}-1+{\frac {1}{z^{2k-1}}}\right\}\\&\qquad +\int _{n=1}^{N}{\frac {B_{2m+1}(n-\lfloor n\rfloor )}{2m+1}}\left\{{\frac {1}{n^{2m+1}}}-{\frac {1}{(n+z-1)^{2m+1}}}\right\}dn\\&=\lim _{N\to \infty }\left(N+z-{\frac {1}{2}}\right){\big \{}\log N-\log(N+z-1){\big \}}+\left(z-{\frac {1}{2}}\right)\log z\\&\qquad +\sum _{k=1}^{m}{\frac {B_{2k}}{(2k)(2k-1)}}\left\{{\frac {1}{N^{2k-1}}}-{\frac {1}{(N+z-1)^{2k-1}}}-1+{\frac {1}{z^{2k-1}}}\right\}\\&\qquad +\int _{n=1}^{N}{\frac {B_{2m+1}(n-\lfloor n\rfloor )}{2m+1}}\left\{{\frac {1}{n^{2m+1}}}-{\frac {1}{(n+z-1)^{2m+1}}}\right\}dn\\&=-z+1+\left(z-{\frac {1}{2}}\right)\log z-\sum _{k=1}^{m}{\frac {B_{2k}}{(2k)(2k-1)}}\left(1-{\frac {1}{z^{2k-1}}}\right)\\&\qquad +\int _{n=1}^{N}{\frac {B_{2m+1}(n-\lfloor n\rfloor )}{2m+1}}\left\{{\frac {1}{n^{2m+1}}}-{\frac {1}{(n+z-1)^{2m+1}}}\right\}dn\end{aligned}}}$

となる。右辺の定数を集めて

${\displaystyle C=1-\sum _{k=1}^{m}{\frac {B_{2k}}{(2k)(2k-1)}}+\int _{n=1}^{N}{\frac {B_{2m+1}(n-\lfloor n\rfloor )dn}{(2m+1)n^{2m+1}}}}$

とすれば

${\displaystyle \log \Gamma (z)=C-z+\left(z-{\frac {1}{2}}\right)\log z+\sum _{k=1}^{m}{\frac {B_{2k}}{(2k)(2k-1)z^{2k-1}}}-\int _{n=1}^{N}{\frac {B_{2m+1}(n-\lfloor n\rfloor )dn}{(2m+1)(n+z-1)^{2m+1}}}}$

となり、この主要部をガンマ関数の相補公式に代入して ${\displaystyle z\to i\infty }$ とすれば

${\displaystyle \Gamma (z)\Gamma (1-z)=-z\Gamma (z)\Gamma (-z)={\frac {\pi }{\sin \pi z}}}$

${\displaystyle \pi i+\log z+\log \Gamma (z)+\log \Gamma (-z)-\log \sin \pi z-\log \pi =0}$

${\displaystyle \pi i+\log z+C-z+\left(z-{\frac {1}{2}}\right)\log z+C+z+\left(-z-{\frac {1}{2}}\right)\left(\pi i+\log z\right)-\log \sin \pi z-\log \pi =0}$

${\displaystyle 2C+{\frac {\pi i}{2}}-\pi iz-\log \sin \pi z-\log \pi =0}$

となるが

${\displaystyle \log \sin \pi z=\log \left(e^{\pi iz}-e^{-\pi iz}\right)-{\frac {\pi i}{2}}-\log 2\sim \pi iz-{\frac {\pi i}{2}}-\log 2}$

であるから

${\displaystyle C={\frac {\log 2\pi }{2}}}$

を得る。剰余項については

${\displaystyle \alpha ={\frac {2}{1+\cos \arg z}}}$

として

${\displaystyle {\begin{aligned}\left|\int _{n=1}^{N}{\frac {B_{2m+1}(n-\lfloor n\rfloor )dn}{(2m+1)(n+z-1)^{2m+1}}}\right|&\leq {\frac {\left|B_{2m}\right|}{2m}}\int _{n=1}^{N}{\frac {dn}{\left|n+z-1\right|^{2m+1}}}\\&\leq {\frac {\left|B_{2m}\right|}{2m\alpha ^{2m+1}}}\int _{n=1}^{N}{\frac {dn}{(n+\left|z\right|-1)^{2m+1}}}={\frac {\left|B_{2m}\right|}{\alpha ^{2m+1}\left|z\right|^{2m}}}=O\left(z^{-2m}\right)\\\end{aligned}}}$

である。故に

を得る。最初の数項を書き下せば

${\displaystyle \log \Gamma (z)\sim \log {\sqrt {2\pi }}-z+\left(z-{\frac {1}{2}}\right)\log z+{\frac {1}{12z}}-{\frac {1}{360z^{3}}}+{\frac {1}{1260z^{5}}}-{\frac {1}{1680z^{7}}}+{\frac {1}{1188z^{9}}}}$

${\displaystyle \Gamma (z)\sim {\sqrt {\frac {2\pi }{z}}}\left({\frac {z}{e}}\right)^{z}\exp \left({\frac {1}{12z}}-{\frac {1}{360z^{3}}}+{\frac {1}{1260z^{5}}}-{\frac {1}{1680z^{7}}}+{\frac {1}{1188z^{9}}}\right)}$

とやり、指数関数のテイラー展開により

${\displaystyle \Gamma (z)\sim {\sqrt {\frac {2\pi }{z}}}\left({\frac {z}{e}}\right)^{z}\left(1+{\frac {1}{12z}}+{\frac {1}{288z^{2}}}-{\frac {139}{51840z^{3}}}\right)}$

となる。

スターリングの公式は鞍点法の好適例とされることが多いが、実際に複素平面全体（負の実数を除く）で漸近近似が成立することを鞍点法によって示すのは困難であるから、ここでは ${\displaystyle z}$ を正の実数に限定する。ガンマ関数は ${\displaystyle t=z(1+u)}$ の置換により

${\displaystyle {\begin{aligned}\Gamma (z+1)&=\int _{0}^{\infty }t^{z}e^{-t}dt\\&=\int _{-1}^{\infty }z^{z}(1+u)^{z}e^{-z-zu}zdu\\&=z^{z+1}e^{-z}\int _{-1}^{\infty }e^{-z\left\{u-\log(1+u)\right\}}du\\&=z^{z+1}e^{-z}\left(\int _{-1}^{-\varepsilon }e^{-z\left\{u-\log(1+u)\right\}}du+\int _{-\varepsilon }^{\varepsilon }{e^{-z\left\{u-\log(1+u)\right\}}du}+\int _{\varepsilon }^{\infty }e^{-z\left\{u-\log(1+u)\right\}}du\right)\qquad (\varepsilon \ll 1)\\\end{aligned}}}$

となるが、${\displaystyle z}$ が十分に大きければ ${\displaystyle u=0}$ の付近が支配的であるから

${\displaystyle \Gamma (z+1)\approx z^{z+1}e^{-z}\int _{-\varepsilon }^{\varepsilon }e^{-z\left\{u-\log(1+u)\right\}}du\approx z^{z+1}e^{-z}\int _{-\varepsilon }^{\varepsilon }e^{-{\tfrac {zu^{2}}{2}}}du\approx z^{z+1}e^{-z}\int _{-\infty }^{\infty }e^{-{\tfrac {zu^{2}}{2}}}du}$

という近似が許され、ガウス積分により

${\displaystyle \Gamma (z+1)\sim z^{z+1}e^{-z}{\sqrt {\frac {2\pi }{z}}}={\sqrt {2\pi z}}\left({\frac {z}{e}}\right)^{z}}$

を得る。${\displaystyle \varepsilon =z^{-{\tfrac {1}{3}}}}$ として、近似の誤差は

${\displaystyle {\begin{aligned}\left|\int _{-1}^{-\varepsilon }e^{-z\left\{u-\log(1+u)\right\}}du\right|&\leq {\frac {1}{\varepsilon }}\int _{-1}^{-\varepsilon }\left|{\frac {-u}{(1+u)}}e^{-z\left\{u-\log(1+u)\right\}}\right|du={\frac {1}{\varepsilon z}}\left[e^{-z\left\{u-\log(1+u)\right\}}\right]_{-1}^{-\varepsilon }\\&={\frac {1}{\varepsilon z}}\left\{e^{-z\left({\tfrac {\varepsilon ^{2}}{2}}+O(\varepsilon ^{3})\right)}\right\}-0\approx z^{\tfrac {2}{3}}e^{-{\tfrac {z^{1/3}}{2}}}\end{aligned}}}$

${\displaystyle {\begin{aligned}\left|\int _{\varepsilon }^{\infty }e^{-z\left\{u-\log(1+u)\right\}}du\right|&\leq {\frac {2}{\varepsilon }}\int _{\varepsilon }^{\infty }\left|{\frac {2u}{(1+u)\varepsilon }}e^{-z\left\{u-\log(1+u)\right\}}\right|du={\frac {2}{\varepsilon z}}\left[-e^{-z\left\{u-\log(1+u)\right\}}\right]_{\varepsilon }^{\infty }\\&=0-{\frac {2}{\varepsilon z}}\left\{e^{-z\left({\tfrac {\varepsilon ^{2}}{2}}+O(\varepsilon ^{3})\right)}\right\}\approx 2z^{\tfrac {2}{3}}e^{-{\tfrac {z^{1/3}}{2}}}\end{aligned}}}$

${\displaystyle {\begin{aligned}\int _{-\varepsilon }^{\varepsilon }e^{-z\left\{u-\log(1+u)\right\}}du&=\int _{-\varepsilon }^{\varepsilon }e^{-z\left({\tfrac {u^{2}}{2}}-{\tfrac {u^{3}}{3}}+O(u^{4})\right)}du=\int _{-\varepsilon }^{\varepsilon }e^{-{\tfrac {zu^{2}}{2}}}e^{z\left({\tfrac {u^{3}}{3}}+O(u^{4})\right)}du\\&=\int _{-\varepsilon }^{\varepsilon }e^{-{\tfrac {zu^{2}}{2}}}\left(1+{\frac {zu^{3}}{3}}+zO(u^{4})\right)du\\&=\int _{-\varepsilon }^{\varepsilon }e^{-{\tfrac {zu^{2}}{2}}}\left(1+zO(u^{4})\right)du\\\end{aligned}}}$

${\displaystyle \left|\int _{-\varepsilon }^{\varepsilon }e^{-z\left\{u-\log(1+u)\right\}}du-\int _{-\varepsilon }^{\varepsilon }e^{-{\tfrac {zu^{2}}{2}}}du\right|\leq \left|zO(\varepsilon ^{4})\int _{-\varepsilon }^{\varepsilon }e^{-{\tfrac {zu^{2}}{2}}}du\right|=\left|O(z^{-{\tfrac {1}{3}}})\int _{-\varepsilon }^{\varepsilon }e^{-{\tfrac {zu^{2}}{2}}}du\right|}$

であり

${\displaystyle z^{\tfrac {2}{3}}e^{-{\tfrac {z^{1/3}}{2}}}\ll z^{-{\tfrac {1}{2}}}\qquad (z\to \infty )}$

であるから

${\displaystyle \Gamma (z+1)={\sqrt {2\pi z}}\left({\frac {z}{e}}\right)^{z}\left(1+O(z^{-{\tfrac {1}{3}}})\right)}$

を得る。これは

${\displaystyle \lim _{z\to \infty }{\frac {\Gamma (z+1)}{{\sqrt {2\pi z}}\left({\frac {z}{e}}\right)^{z}}}=1\qquad (|\arg z|<\pi )}$

を示すに十分である。ただし、実際の誤差は ${\displaystyle O(z^{-1})}$ であるが、それを鞍点法で示すのは困難である。