revised　25/09/19　　updated 21/07/01



ガンマ関数は階乗を複素数に拡張したもので、名前は良く知られています。しかし、実際の計算方法についての解説が意外に少ないので、どういうことだろうと疑問に思っていたのですが、積分を使う面倒な計算をしなければいけないようで、誰でもすぐに判るものではなさそうです。しかし、私の発見した計算式を使うとガンマ関数をすぐに計算できます。積分を使わないので中学生でも理解できる内容となっています。そのことをこのＨＰをとおして広く知ってもらいたいと願っています。



■　　　　　　実数用の階乗公式を複素数に拡張します。

先のページで説明した実数用の階乗公式はそのまま複素数にも適用できます。ｎをｓに変えるだけです。ただし、このままではガンマー関数になりません。なぜなら、Γ(s) = (s-1)!　ということで、複素数階乗から１を引かなければなりません。もしくは複素数階乗に１を加えると言うべきでしょうか。とにかく、この計算は簡単なので、複素数階乗を少し修正すればガンマー関数になります。

実数階乗の公式のｎをｓに変えると次の式になります。

s!　=　(s/e)^s * √(2*pi*s) * e^ ｛ B(2)/2s + B(4)/12s^3 + B(6)/30s^5 + B(8)/56s^7 .... 　｝


計算は面倒ですが、プログラム化してしまえば簡単に答えが出ます。



■　　　　　　複素数用の階乗公式もあります。

複素数用の公式は実数も計算できますが、実数用よりもやや計算に時間がかかります。複素数計算の場合はこちらの公式のほうが安定しているように思いますが、そのあたりは現在検討中です。

複素数用の階乗公式にはふたつの表記方法があります。A型
s! = k(1-s) * s * (2pi)^s / 2 / k(s) / cos(pi*s/2)

もしくは、B型
s! = - k(-s) * (2pi)^(s+1) / 2 / k(s+1) / sin(pi*s/2)
です。

A型でも、B型でも同じ答えになりますが、k(s)　k(s+1) で割るところがあるので、k(s)=0 　 k(s+1)=0　のときは他の型を使います。



■　　　　　　ガンマ関数の公式　（複素数用を変更したもの）

ガンマ関数は　Γ(s+1) = s!　、Γ(s) = (s-1)!　という関係なので、１を引いて複素数階乗の計算をすれば、そのままガンマ関数の値になります。つまり、階乗の公式を　ｓ　で割るとガンマ関数の公式になります。

Γ(s) = 1/s * (s/e)^s * √(2*pi*s) * e^ ｛ B(2)/2s + B(4)/12s^3 + B(6)/30s^5 + B(8)/56s^7 .... ｝


複素数用の公式をｓで割ると以下のようになります。

A型
Γ(s) = k(1-s) * (2pi)^s / 2 / k(s) / cos(pi*s/2)

B型
Γ(s) = - k(-s) * (2pi)^(s+1) / 2s / k(s+1) / sin(pi*s/2)


(s-1) を代入してもガンマー関数公式になります。

cos(pi*(s-1)/2) = sin(pi*s/2)
sin(pi*(s-1)/2) = -cos(pi*s/2) なので、整理すると、

A型
Γ(s) = k(2-s) * (s-1) * (2pi)^s / 4pi / k(s-1) / sin(pi*s/2)

B型
Γ(s) = k(1-s) * (2pi)^s / 2 / k(s) / cos(pi*s/2)


1/s　のＡ型と　(s-1)　を代入したＢ型が同じになっているのは、当たり前かもしれませんが、面白いところです。


s=0+i　として　Γ(s) を計算したところ、A形も、B型も　(-1.54949E-1, -4.98015E-1)　となりました。

s=-1+i　として　s! を計算したところ同じ答えとなりました。式は間違っていません。


k(s) は、すでに他のページでも説明してありますが、再掲すると、

k(s) = lim_[n→∞] { zt(s,n) + ber(s,n) }

zt(s,n) = 1 + 1/2^s + 1/3^s + . . . + 1/n^s
ber(s,n) = Σ_[r=0,∞] B(r)*(s-2+r)!/r!/(s-1)!/n^(s-1)/n^r

ということです。

zt(s) と ber(s) の無限式をさらに無限に計算する結果なので、大変そうですが、収束は早いので、途中で打ち切れます。すでにjava を使った計算プログラムがあるので、それを使えば簡単に値が求められます。複素階乗も計算プログラムを作ってあります。ですから、答えはすぐに出ます。以下のところから入手できます。
Java による長桁計算



■　　　　　　副産物としてベルヌーイ数の新しい計算式が見つかったので記録しておきます。

階乗の公式に (n-1) を代入して、計算結果を纏めている過程で、意外にもベルヌーイ数の新しい計算式が見つかりました。ベルヌーイ数は各種の関数計算で利用するので大きな数のベルヌーイ数まで７０００桁まで計算しておかなければなりません。しかし、これにはかなりの時間がかかるので、今まで　B(1000)　あたりまでで終わらせていました。しかし、この新しい計算式を使うと時間がかなり短縮できるので便利です。

さて、その求め方ですが、　Γ(n) = (n-1)!　なので、階乗の公式に　(n-1) を代入します。一方、Γ(n) = 1/n * (n/e)^n * √2*pi*n) * ...なので、両者をイクオールで結びます。

代入式は　((n-1)/e)^(n-1) * √(2*pi*(n-1)) * e^{ B(2)/2(n-1) + B(4)/12(n-1)^3 + ... }

(n-1)^(n-1), √(n-1), 1/(n-2)^r　は「参考となる公式」HPに載っています。

(n-1)^(n-1) = n^(n-1) * e^{ -1 + 1/2n + 1/6n^2 +1/12n^3 +1/20n^4 + ...}

√(n-1) = √n * e^{ -1/2n - 1/4n^2 - 1/6n^3 - .... }

1/(n-1) = 1/n + 1/n^2 + 1/n^3 + 1/n^4 + ....

1/(n-1)^2 = 1/n^2 + 2/n^3 + 3/n^4 + 4/n^5 + ....

1/(n-1)^3 = 1/n^3 + 3/n^4 + 6/n^5 + 10/n^6 + ....

以下略

代入式を纏めると、ガンマ公式と同じところが現れるので、それを消して整理します。すると　e^f(n) のところだけが残ります。その部分をさらに整理すると、

B(2)は　　　0 = 1/6 -1/4 + B(2)/2

B(4)は　　　0 = 1/20 - 1/8 + B(2)/2 + B(4)*3/12

B(6)は　　　0 = 1/42 - 1/12 + B(2)/2 + B(4)*10/12 + B(6)*C(5,4)/30

B(8)は　　　0 = 1/72 - 1/16 + B(2)/2 + B(4)*4/12 + B(6)*C(5,4)/30 * B(8)*C(7,6)/56


これを一般的に表記すると

B(r)*C(r-1,r-2)/r(r-1) = - 1/r(r+1) + 1/2r - B(2)/2 - B(4)*C(r-1,2)/12 - B(6)*C(r-1,4)/30 - ... B(r-2)*C(r-1,r-3)/(r-2)(r-3)

という式になります。これで　B(r) を計算できます。

実際にやってみましょう。

B(2)/2 = - 1/6 + 1/4　　　なので、　B(2)=1/6
B(4)*3/12 = -1/20 + 1/8 - 1/12　なので　B(4)=-1/30
B(6)*5/30 = -1/42 + 1/12 -1/12 + 10/30*12　　なので　　B(6)=1/42

以下略

これで計算式が正しいことが確認できました。

これをもとにプログラムを作り、ベルヌーイ数を７０００桁計算したところ、n=2000 まで２時間３７分ほどで計算が終わりました。以前、n=1000 まで計算したときは二日ほどかかった記憶がありますが、今回パソコンのCPUが3.6GHzあるのと、プログラムが高速になっているので、このような短時間で計算できたのだろうと思います。

このままどこまでも計算してみたいという思いが湧いてきましたが、値はどんどん大きくなるし、計算結果を使うことはないだろうと思って断念しかけました。しかし、データのためではなく、増加率が判れば公式が見つかるのではないかと思い直し、調べてみるとことにしました。

B(n) ではあまりにも値が大きくなりすぎるので、B(n)/n! で計算してみることにしました。その結果、意外と綺麗な公式が見つかりました。これはホームランかもしれません。驚くとともに、とても喜んでいます。

B(n)/n! = 2/(2pi)^n * { 1 + 1/2^n + 1/3^n + 1/4^n + ....


B(n) のプラス・マイナスは無視しています。

B(n)の符号も考慮した式を求めるのは簡単で、上記の式に適切な値となる　sin関数を掛ければよいのです。sinがどういう関数になるを調べたところ以下のようになりました。

B(n)/n! = - cos(n*pi/2) * 2 / (2pi)^n * { 1 + 1/2^n + 1/3^n + 1/4^n + ..... }

これが　n 実数のときも成り立つことは、B(n)=-n*k(1-n)　および、階乗実数用公式を用いて検証できます。

n=7.4 としてみましょう。B(7.4)=-1.6732348E-2　　7.4!=1.1405887E4　、B(n)/n! = -1.4669922E-6
これに対して、公式を使って計算した結果は　-1.4669922E-6　でぴったり一致しています。

n=1 とすると無限大＊零　となり、計算できませんが、n→1　で計算すると　0.5　に近づいているように見えます。ただし、ここで問題なのは　B(1)=-1/2 であって　B(1)=1/2 ではないことです。この問題はすでにどこかで説明したように思いますが、本来のベルヌーイ数では　B(1)=1/2 であることです。ゼータ関数で使う場合、B(1)=-1/2 のほうが便利なのでこのように定義されました。また、B(1)=-1/2 でないと困ることもあるので、-1/2 としなければなりません。しかし、公式を使って計算すると　B(1)=1/2 になることは知っておく必要があります。


さて、B(n)/n! の公式の中に　( 1 + 1/2^n + 1/3^n + ... ) というゼータ関数が現れました。これがそのまま k(n) になるかどうかは微妙な問題ですが、強引に　= k(n) としてみると、B(n)/n! = - 2cos(n*pi/2) / (2pi)^n * k(n) となります。B(n) = - n*k(1-n) なので、代入して整理すると　k(1-n)/k(n) = 2cos(n*pi/2) * (n-1)! / (2pi)^n　　となります。これは以前紹介した　kv(n) = n! * 2 * cos(pi*n/2) / n * (2pi)^n　という　kv(n) の公式そのものです。そして、kv(n) 式は複素数に拡張されて　kv(s) となるので、B(n)/n! の公式も複素数に拡張することが出来そうです。もちろん、この場合の公式は k(s) を使わなければなりません。

B(s)/s! = - cos(s*pi/2) * 2 * k(s) / (2pi)^s