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

ガンマ関数を計算するには階乗の公式を使います。なぜなら、ガンマー関数と階乗は基本的に同じものだからです。Γ(n) = (n-1)!　という法則がありますが、これは定義とも言えることで、階乗の複素数化がガンマ関数であると理解してそれほど間違いではありません。階乗の公式の　ｎ（実数）　を　ｓ（複素数）　に変えれば良いだけの話なのです。

階乗の公式は、２種類あります。ひとつはおもに実数の計算に使うもので、もうひとつは複素数に使うものです。実数用の公式と言っても複素数を計算できないわけではありません。複素数用の公式は実数も計算できますが、実数用よりもやや計算に時間がかかります。


実数用の階乗の公式は
n! = (n/e)^n * √(2*pi*n) * e^ ｛ B(2)/2n + B(4)/12n^3 + B(6)/30n^5 + B(8)/56n^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)!　という関係なので、１を引いて複素数階乗の計算をすれば、そのままガンマ関数の値になります。

階乗の公式はまだ知られてないようで、私としては自慢できるのですが、具体的計算方法については前頁とか、hirokuroリーマン証明ver20 第一部などを参照していただければと思います。

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

式を再掲しておきます。

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




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


階乗の公式を変更して、ガンマ関数の公式を導いてみようと思い立ちました。計算が早くなるとか、便利なことがあるということではなく、単なる興味本位の思いつきでやったことですが、意外と重要な結果を得ることができました。

計算をするためだけならば
Γ(n) = 1/n * (n/e)^n * √(2*pi*n) * e^ ｛ B(2)/2n + B(4)/12n^3 + B(6)/30n^5 + B(8)/56n^7 .... ｝
で良いわけで、面白いところは何もありません。

これを複素数用にするには　n を　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 .... ｝




興味深いのは、階乗の公式に (n-1) を代入して、計算結果を纏めてみたところ、意外にもベルヌーイ数の計算式が見つかったことです。

Γ(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





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


階乗の公式（複素数用）をガンマ関数に変更してみます。

階乗の公式を　ｓ　で割るとガンマ関数の公式になります。

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! を計算したところ同じ答えとなりました。式は間違っていません。