数学公式（２） 　階乗の公式

公式というものは、ただ並べればよいというものではありません。有益でなければならないと同時に、美しくなければなりません。私が特に重視したいのは美しさです。有益であるかどうかは、あとの話です。

昔、虚数が発見されたとき、その発見者自身がその有益さには気が付きませんでした。せいぜい便利と思った程度でした。しかし、その後、時代が変わった後、虚数は便利どころではなく、数学の根幹に関わる基本的数であることが理解されてきました。虚数なしには、今日の物理学は成り立たないでしょう。それほど有益であると判ったのは随分後のことでした。

ですから、有益かどうかは第二のことにしておきましょう。まず美しいかどうかです。

私が発見した公式も美しいものばかりですが、美しいかどうかは客観的に判断されることではありません。「自分の子は皆可愛い。」と言われるように、自分で発見したとなると可愛く見えるものです。しかし、可愛いことは事実ですから、せめて自分くらいは大事にしてゆきたいと思っています。

ご紹介できるほどのものになると、少し自信が無くなってきますが、すでに周知の公式である可能性もありますので、大きな声では言えません。小さい声で、一応、私が他の数学の書籍を見ることなく、独自に発見した式です。

それらは互いに関連していますので、説明すると長くなりますが、少しずつ説明していきましょう。


１）　　　　　　Σ(1/n)　を展開する

みなさん、Σ(1/n) の計算はしたことがあるでしょうか？ごく単純な式ですので、プログラムを書く練習とか、計算機が正常に動いているかどうかをチェックするのに便利です。実用性は、ほとんどありませんが、私にとっては愛すべき式となっています。

以前、プログラム電卓を購入した時、プログラムの練習と言うこともあり、また電卓機能の確認のため、Σ(1/n) を計算したことがあります。　f(n)=Σ(1/n) としましょう。　つまり、
f(n)=1+1/2+1/3+1/4+1/5+1/6+.....
と続くことになります。
f(1)=1　となり、　f(2)=1.5、　f(3)=1.833333...　となります。
f(10)=2.928968254 、　f(100)=5.187377518　、　　 f(1000)=7.485470859

f(1000)には、計算するのに５分もかかっています。今から考えると鈍いですが、当時はその電卓のスピードにビックリしたものです。
このまま計算を続けると無限になってしまいますので、それでは面白くありません。そこで、 f(n) がどの式に近づくかを調べてみました。つまり、f(n)が　n　に近づくのか、n^2　に近づくかということです。これはやらなくてもだいたい見当は付きますが、1000で 7.48ですから、そんなに大きくなりません。ですから、f(n)=log(n)に見当を付けることになります。

f(n)=a*log(n) とおいて、 a が何かの値に収束するのではとの予想のもとに計算してみました。これには結構時間がかかりましたが、 a=2.302585 という数字を得まして、これが　1/log(e)　であることが判ったときは天にも昇る気持ちでした。まさか、綺麗な定数になるとは思っていなかったからです。

計算すると、こういう良いことがあるのかと味をしめまして、次に　f(n) = 2.3025 * log(n) * { 1 + a/log(n) + b/log(n)^2 ..... } 　という式を想定し、この　a、b　を求めてみることにしました。これには手こずりました。電卓の誤差が出てくるのですが、当初その誤差のことが判らず、突然値が狂い始めるのですが、それが本来の式の値と誤解してしまいました。しかし、これも慣れてきますと、どれが誤差かはだいたい見当が付きます。

どうも　a=0.25　あたりらしいと言うことまで判りました。　b　についてはまったく手つかずの状態に終わりました。電卓ではこのあたりで終わるのもやむを得ません。しばらくここで止まったままでした。

その後、どう進んだか、自分でも記憶が薄れてきましたが、ある時、本屋で数学書を立ち読みしているとき、オイラー定数という数を発見しました。それによると、Σ(1/n) = 2.302log(n) + 0.577 ... となっていて、これはすでにオイラーが発見しているとのことでした。　0.577...の定数のことをオイラー定数と言うとのことでした。さすがオイラーさんは大したものだと感心した次第です。私の計算結果の 0.25　と 2.3025 を掛け合わせると、確かに　0.5756 になって　0.577 に近いことが判りました。私の計算も決して間違ってはいなかったようです。（オイラー定数を千桁まで計算しておきました。最初の-1は指数部を示しています。つまり、0.577...ということです。先頭が０の時は左詰に表示されているようなのでご注意ください。）

しかし、私が発見しようとしていた数をオイラーに先取りされたことは少し悔しい気がしました。やや無気になって、その続きの数を発見しようと決意したのです。つまり、Σ(1/n) をすべて　log(n) と　n　で表現してみようと言うことです。a、b、c　などという一つ二つの項数ではなく、すべての項数を発見するのですから、これは簡単ではありませんでした。当然のことながら電卓では壁にぶつかりました。誤差を少なくするうまい工夫を思いついたのですが、それでも第６項あたりになるともういけません。とても凡ての項数を発見するとまでは行きません。

その時、推定した式は、 f(n)=2.3025log(n)+0.577215+1/2n-1/12n^2+1/120n^4-1/288n^6+1/504n^8+...　というものでした。

ここで立ち止まってしまい、しばらく時が過ぎましたが、その後、別の公式を発見する過程で、この公式の項数の展開法則を見つけました。ラッキーとはこのことです。非常に綺麗な法則になっています。上記の推定式は、近似的には正しいのですが、やはり違っていました。

電卓の中にはΣ(1/n)　を計算できるようになっているのもありますので、私の発見した式はすでに知られている公式かもしれませんが、今のところ、数学公式集には書かれていません。私のオリジナルだと嬉しいのですが、そうでなくても、自分で意地張って見つけたのですから、満足しなければなりません。

その見つけた結果は、もったいぶっているようですが、次の述べる公式とも関係していますので、それを説明した後に書くことにします。



２）　　　　n!（階乗）　もしくは　Σlog(n)　の公式

n!　とΣlog(n)　は似ているというか、結局は同じ式ですから、一緒に扱います。

n!はあちこちで使われる重要関数で、基本的にはガンマ関数と同じことです。　n!=n*(n-1)*(n-2)*(n-3)....3*2*1　となります。　0!=1　　1!=1　　2!=2　　3!=6　　4!=24　　5!=120　　6!=720　　7!=5040　と、どんどん数が大きくなります。
これがどの式と同じレベルの無限であるかは、電卓でもなんとか見つけることが出来ます。つまり、n!=(n/e)^n*f(n)　あたりになることは判っていました。しかし、これはすでに発見されているようで、数学公式集にスターリングの公式として以下の式が載っていました。

スターリングの公式　　 n! ≒ (n/e)^n * √(2*pi*n) * { 1+1/12n+1/288n^2... }

これは近似値ですから、小さな数ですと誤差がでます。（注　　スターリングの公式については、さらに詳しい説明を別ページでしています。第６の広場の「スターリングの公式」を参照してください。）

1/288の続きの項数はどうなるか気になりましたので、n^3の項数を電卓で計算してみましたが、さすがに誤差ばかりで、法則が見えるような数には確定できませんでした。

この研究もここで止まってしまったのですが、まもなく、電卓で計算するだけが数学ではないことを教えられることになります。


３）　　　　　　(n+1)^(n+1)の展開

(n+1)^(n+1)の無限は n^(n+1) に近いことは容易に想定できますが、計算機を使う限り、　(n+1)^(n+1)≒e*n^(n+1) 程度のことしか判りません。しかし、別のやり方を使うと綺麗に求めることが出来ます。

f(n) = (n+1)^(n+1) とします。二項係数を c(n,r) で表現する事にします。

すると、
f(n) = n^(n+1) + c(n+1,1)*n^n + c(n+1,2)*n^(n-1) + c(n+1,3)*n^(n-2) + ..... と続くことが判ります。ただ展開しただけのことです。これをひとつひとつ整理します。

n^(n+1) はこのままですが、 c(n+1,1)*n^n は　(n+1)*n^n ですから、　n^(n+1)+n^n となります。同じように、c(n+1,2)*n^(n-1) は　{(n+1)*n}/2*n^(n-1) となります。つまり、 n^(n+1)/2+n^n/2 ということです。　c(n+1,3)*n^(n-2) は　{(n+1)*n*(n-1)/6}*n^(n-2) ですから、n^(n+1)/6-n^(n-1)/6 となるわけです。これを凡ての項数について行います。お判かりになるでしょうか。いずれ別紙に書いて張り付けて、もう少し見やすい形にする予定です。

今はこのまま計算したとして、その結果を書き出します。

これらの式を n^(n+1)、n^n、n^(n-1).... と整理してゆきますと、n^(n+1)の項数は　1+1+1/2+1/3!+1/4!+1/5!.... ですから　合計は　e となります。
n^n　の項数は　　1+1/2+1/6!-2/4!-5/5!-9/6!....　 ですから、合計は　e/2 となります。
n^(n-1)　の項数は　-1/3!-1/4!+5/5!+25/6!+70/7!+...　 ですから、合計は　-e/24 となります。
計算は面倒とは言え、結構簡単です。

ここで、上記の式を整理ておきましょう。

f(n) = e*n^(n+1)+e*n^n/2-e*n^(n-1)/24+e*n^(n-2)/48+....

となります。これを判りやすくするために、項数を表記する方法としてsekiと名付けた数列を使うことにします。
seki(0)=1 seki(1)=1/2 seki(2)=-1/24 seki(3)=1/48 となります。

もう一度書いておきましょう。
f(n) = seki(0)*e*n^(n+1) + seki(1)*e*n^n + seki(2)*e*n^(n-1) + seki(3)*e*n^(n-2) + .....
何を書いているかは理解していただけることと思います。
これを整理して書き直すと、 (n+1)^(n+1) = e*n^(n+1) * { seki(0) + seki(1)/n + seki(2)/n^2 + seki(3)/n^3 + ... } となります。

ですから、このseki(r)を明らかにすれば、(n+1)^(n+1)　を展開したと同じことになるわけです。

これについては計算機を使いました。誤差が出ないように工夫して計算すると、seki(4)=-73/48/120 であることがわかりました。　seki(5)=73*99/48/120/146 となりました。ここで、項数に現れる数字に２６の差があることに気が付きました。次のseki(6)を計算すると、　seki(6)=-73*99*125/48/120/146/172 となり、やはり２６の差となっています。その次は有効桁はとっくに過ぎているので、推測ですが　seki(7)=73*99*125*151/48/120/146/172/198 となります。これをもとに試しに計算してみますと、問題ない範囲であることが判りました。

あまり、綺麗ではありませんが、これをきっかけに別の公式を発見することになりますので、ここはこれで良しとしておきましょう。


seki(0)=1
seki(1)=1/2
seki(2)=-1/24
seki(3)=1/48
seki(4)=-73/48/120
seki(5)=73*99/48/120/146
seki(6)=-73*99*125/48/120/146/172　（注　　これ違ってますね。（＾＾；；；汗　修正点については付記を参照してください。）
seki(7)=73*99*125*151/48/120/146/172/198
seki(8)=? （この表を見れば計算することなく書くことが出来ますね。）

（注　　当時は電卓で計算したので間違えたのですが、それでも結論は間違ってないのは不思議ですね。seki(6), seki(7) などは以下の計算で使っていないのが原因のようです。このページは過去の歴史なのでそのまま修正せずに載せておきます。）

付記　　　正しい結果を載せておきます。
seki(6)=-3625/580608
seki(7)=5525/1161216
seki(8)=-5233001/1393459200
この項数法則を見つける意味が無くなったので分析していませんが、(n+1)^(n+1) = n^(n+1) * e^ { 1 + 1/2n - 1/6n^2 + 1/12n^3 - 1/20n^4 ... } であることは別のページで説明してあります。こちらは今のところ正しいみたいです。



４）　　　　n!（階乗）　の続き

一応、 (n+1)^(n+1) = e*n^(n+1) * { seki(0) + seki(1)/n + seki(2)/n^2 + seki(3)/n^3 + ... }　となることにしましょう。ここから n! （階乗の公式）を求める作業に戻ります。

スターリングの公式により　n! = (n/e)^n√(2*pi*n) * { 1+1/12/n+1/288/n^2 ... } であることは判っています。このとき、｛ ｝の中の項数をrui(r)を用いて表記する事にします。つまり、 1 + 1/12n + 1/288n^2... .= 1 + rui(1)/n + rui(2)/n^2 + rui(3)/n^3 + rui(4)/n^4... とします。
rui(1)=1/12　　rui(2)=1/288　　rui(3)=?　　　rui(4)=?　　ということです。

さて、ここでスターリング公式の応用として、(n+1)! = ((n+1)/e)^(n+1) * √(2*pi*(n+1)) * { 1 + rui(1)/(n+1) + rui(2)/(n+1)^2 ... } であることはすぐにおわかりいただけると思います。

また、一方で　(n+1)! = (n+1)*n! ですから、n! にスターリング公式を当てはめると、
(n+1)*n! = (n+1) * (n/e)^n * √(2*pi*n) * { 1 + rui(1)/n + rui(2)/n^2 ... }　となります。

どちらも同じ式ですから、右辺＝左辺とおいて、両辺にe^nを掛け、また、√(2*pi) /nを掛けます。これで少し整理されてきます。

左辺 = (n+1)^(n+1)/e * √((n+1)/n) * { 1 + rui(1)/(n+1) + rui(2)/(n+1)^2 .. .}
右辺 = (n+1)*n^n * { 1 + rui(1)/n + rui(2)/n^2 ... }

ここで左辺に (n+1)^(n+1)　があるので、３）のsekiを使った式に書き換えます。
つまり、(n+1)^(n+1) = e*n^(n+1) * { seki(0) + seki(1)/n + seki(2)/n^2 + seki(3)/n^3 + ... }　 ですから、これを代入するとeが消えます。

さらに、n^nが両辺にあるので、これを消して、その後で、ruiについては右辺に纏め、それ以外の変数を左辺に持ってきます。
左辺 = { √((n+1)/n) } * n/(n+1) * { 1 + seki(1)/n + seki(2)/n^2 .... }
右辺 = { 1 + rui(1)/n + rui(2)/n^2 ... } / { 1 + rui(1)/(n+1) + rui(2)/(n+1)^2 + .... }

この式をよく眺めると判るのですが、左辺は既知の数ばかりですから、計算すれば項数は凡て確定します。右辺はruiの数列だけなので、うまく纏まれば、ruiを確定させることが出来ることが判ります。そして、後で示すように、これは大成功でした。ruiは確定するだけでなく、ある法則のもとに並んでいたのです。

それでは、まず面倒な計算部分を処理してしまいましょう。この計算は、難しいわけではありませんが、膨大な計算量ですので、なるべく計算間違いをしないやり方で計算しなければなりません。そのために、いくつかの段階に区分けして計算することにしました。また、いくつかの工夫も導入しました。


５）　　　　1/f(n) a(r)、b(r) という数列

これは以下の式の中で定義されます。

f(n) = 1+a(1)/n+a(2)/n^2+a(3)/n^3+a(4)/n^4... 　　1/f(n) = 1 + b(1)/n + b(2)/n^2 + b(3)/n^3 + b(4)/n^4 ....

この式を計算しますと、a(r)　と　b(r)　の関係が明らかになります。

b(1)=-a(1)　　b(2)=a(1)^2-a(2) 　　b(3)=-a(1)^3+2*a(1)*a(2)-a(3)

b(4) = a(1)^4 - 3*a(1)^2*a(2) + 2*a(1)*a(3) + a(2)^2 - a(4)
b(5) = -a(1)^5 + 4*a(1)^3*a(2) - 3*a(1)^2*a(3) - 3*a(1)*a(2)^2 + 2*a(1)*a(4) + 2*a(2)*a(3) - a(5)

b(6) = a(1)^6 - 5*a(1)^4*a(2) + 4*a(1)^3*a(3) - 3*a(1)^2*a(4) + 6*a(1)^2*a(2)^2 - a(2)^3 + 2*a(2)*a(4) + 2*a(1)*a(5) - 6*a(1)*a(2)*a(3) + a(3)^2 - a(6)

b(7) = -a(1)^7 + 6*a(1)^5*a(2) - 5*a(1)^4*a(3) + 4*a(1)^3*a(4) - 3*a(1)^2*a(5) + 2*a(19*a(6) + 4*a(1)*a(2)^3 + 12*a(1)^2*a(2)*a(3) - 6*a(1)*a(2)*a(4) - 10*a(1)^3*a(2)^2 + 2*a(2)*a(5) + 2*a(3)*a(4) - 3*a(2)^2*a(3) - 3a(1)*a(3)^2 - a(7)

たとえば、f(n) = 1 - 1/2n + 1/3!n^2 - 1/4!n^3 + 1/5!n^4 .... という式の場合、 1/f(n) = 1 + 1/2n + 1/12n^2 + 0 - 1/720n^4 + ....　となります。a(r)　を上記の式に代入すれば良いと言うことです。



６）　　　　√(1+1/n) の展開

４）に登場した左辺に　n/(n+1)*√((n+1)/n)　という式が見られます。これを整理して√の中に入れると、√(n/(n+1))となります。これは　√(1+1/n)　の逆数になっています。ですから、√(n/(n+1))=1/√(1+1/n) です。

そこで、まずf(n)=√(1+1/n) を計算します。これは (1+1/n)^(0.5) ですから、二項定理を使って展開することが出来ます。 f(n) = 1 + c(0.5,1)/n + c(0.5,2)/n^2 + c(0.5,3)/n^3 + .... となります。　c(n,r) = n*(n-1)*(n-2)*.../r! ですから、　c(0.5,1)=1/2　　c(0.5,2)=0.5*(-0.5)/2!=-1/8　　c(0.5,3)=0.5*(-0.5)*(-1.5)/3!=1/16　　c(0.5,4)=0.5*(-0.5)*(-1.5)*(-2.5)/4!=-5/128　　 となります。
計算は面倒ですが、項数はすべて確定します。

a(1)=1/2　　
a(2)=-1/8　　
a(3)=1/16　　
a(4)=-5/128　　
a(5)=7/256
a(6)=-21/1024

√(n/(n+1)) = 1/f(n)　ですから、５）の式に代入すると、
√(n/(n+1)) = 1 - 1/2/n + 3/8/n^2 - 5/16/n^3 + 35/128/n^4 - 63/256/n^5 + .....

最後に、f(n) = 1 + seki(1)/n + seki(2)/n^2 .... を掛けて、数値部分の計算は終わります。
つまり、Ans(n) = 1 + 1/n * { -1/2 + seki(1) } + 1/n^2 * { 3/8 - seki(1)/2 + seki(2) } + 1/n^3 * { -5/16 + 3/8*seki(1) - seki(2)/2 + seki(3) } + 1/n^4 * { 35/128 - 5/16*seki(1) + 3/8*seki(2) - seki(3)/2 + seki(4) }

Ans(n) = 1 + kou(1)/n + kou(2)/n^2 + kou(3)/n^3 + ..... 　となるように　kou(r)　を定義しておくと、
kou(1)=0　kou(2)=1/12　kou(3)=-1/12　kou(4)=113/1440　　kou(5)=?

となります。



７）　　　　n!　の続き

４）で計算した右辺を纏めてみます。

右辺 = { 1 + rui(1)/n + rui(2)/n^2 ... } * { 1 + rui(1)/(n+1) + rui(2)/(n+1)^2 + .... }

このうち、1+rui(1)/n+rui(2)/n^2...　を [7.1]式、
1+rui(1)/(n+1)+rui(2)/(n+1)^2+....　を[7.2]式、と名付けます。


まず、[7.2]　1+rui(1)/(n+1)+rui(2)/(n+1)^2+....　部分を計算しておきましょう。
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-......

と並んでいますから、[7.2]式は　f(n) = 1 + den(1)/n + den(2)/n^2 + den(3)/n^3 + ..... となるように　den(r)　を定義すると、
den(1) = rui(1)
den(2) = -rui(1)+rui(2)
den(3) = rui(1)-2*rui(2)+rui(3)
den(4) = -rui(1)+3*rui(2)-3*rui(3)+rui(4)
den(5) = rui(1)-4*rui(2)+6*rui(3)-4*rui(4)+rui(5)
den(6) = -rui(1)+5*rui(2)-10*rui(3)+10*rui(4)-5*rui(5)+rui(6)

g(n) = 1/f(n) = 1 + sug(1)/n + sug(2)/n^2 + sug(3)/n^3 + sug(4)/n^4 ...　となるように　sug(r) を定義します。

すると、５）の公式により、sug(1) = -den(1) 　sug(2) = den(1)^2 - den(2) 　　sug(3) = -den(1)^3 + 2*den(1)*den(2) - den(3)　　 sug(4) = den(1)^4 - 3*den(1)^2*den(2) + 2*den(1)*den(3) + den(2)^2 - den(4)

となります。ここで den(r) に rui(r) を代入しておきます。
すると、sug(1) = -rui(1) 　sug(2) = rui(1)^2+rui(1)-rui(2) 　　sug(3) = -rui(1)^3 - 2*rui(1)^2 + 2*rui(1)*rui(2) - rui(1) + 2*rui(2) - rui(3)
sug(4) = rui(1)^4 + 3*rui(1)^3 + 3*rui(1)^2 + rui(1)^3*rui(1)^2*rui(2) - 6*rui(1)*rui(2) + 2*rui(1)*rui(3) + rui(2)^2 - 3*rui(2) + 3*rui(3) - rui(4)


これと[7.1]式を掛けます。[7.1]*g(n) = 1+a(1)/n+a(2)/n^2+a(3)/n^3....　とすると、
a(1) = rui(1) + sug(1)、　a(2) = rui(2) + rui(1)*sug(1) + sug(2)、 　a(3) = rui(3) + rui(2)*sug(1) + rui(1)*sug(2) + sug(3)、　a(4) = rui(4) + rui(3)*sug(1) + rui(2)*sug(2) + rui(1)*sug(3) + sug(4)

sug(r) を代入します。
a(1) = 0、 　a(2) = rui(1)、　a(3) = -rui(1)^2-rui(1)+2*rui(2)、 　a(4) = rui(1)^3 + 2*rui(1)^2 + rui(1) - 3*rui(1)*rui(2) - 3*rui(2) + 3*rui(3)


７）はすべて後ろから計算する形になりますが、６）の対応する(1/n)^r の項数は一致しますので、
kou(1)=a(1) kou(2)=a(2) kou(3)=a(3).... の形になります。
ですから、これをそのまま計算すればよいわけです。

すると
rui(1)=1/12
rui(2)=1/288
rui(3)=-139/51840
rui(4)=-571/2488320

これでスターリング公式より少し詳しい近似式を求められたことになります。



８）　　　　Σlog(n)

７）の表ではrui(r)の法則性が見えないので、さらに別の角度から分析してみることにします。

そこで、注目するのは　Σlog(n)　です。この式は、n! と基本的には同じもので、log(n!) のことです。ですから、Σlog(n)　を調べると　n! の特色も判ってきます。

n! を利用して、その対数形を作りますと、次のようになります。
Σlog(n) = n*log(n) - n*log(e) + 1/2*log(n) + 1/2*log(2*pi) + log(1+1/12/n+1/288/n^2+....)

となります。



９）　　　　log(f(n))の展開

８）で、log(1+1/12/n+1/288/n^2+...) という式が出てきます。これは n^r の式に変換できます。

log{ 1+a(1)/n+a(2)/n^2+a(3)/n^3+a(4)/n^4..... } = log(e)* {b(1)/n+b(2)/n^2+b(3)/^3.... } となるように b(r) を定義します。

すると、
b(1)=a(1)
b(2)=1/2{-a(1)^2+2a(2)}
b(3)=1/3{a(1)^3-3*a(1)*a(2)+3*a(3)}
b(4)=1/4{-a(1)^4+4*a(1)^2*a(2)-4*a(1)*a(3)-2*a(2)^2+4*a(4)}
b(5) = 1/5 { a(1)^5 - 5*a(1)^3*a(2) + 5*a(1)^2*a(3) + 5*a(1)*a(2)^2 - 5*a(1)*a(4) - 5*a(2)*a(3) + 5*a(5)}

たとえば、f(n)=1+1/n+1/n^2+1/n^3+....　という式を取り上げましょう。　log(f(n)) を展開しようと思うなら、上記の公式に代入して、b(r) を求めればよいことになります。やってみましょう。

log(f(n))=log(e)*(1/n+1/2n^2+1/3n^3+1/4n^4...)
となります。

この展開式を　log(1+1/12n+1/288n^2+...) に適用します。

すると、　log(f(n)) = log(e) * ( 1/12n + 0 - 1/360/n^3 + 0 +...) となります。これは結構綺麗です。n^2 の項目が 0 になりますが、これは検算の役目も果たしてくれます。

つまり、
Σlog(n) = n*log(n) - n*log(e) + 1/2*log(n) + 1/2*log(2*pi) + log(e)*(1/12/n-1/360/n^3+....)　と書くことができると言うことです。

ここで、rui(r) の計算はひどく面倒なので、さらに簡略化する方法を導入します。これには「発見した公式（１）」で見つけたsu(r) を使います。

1*rui(1)=su(2)*1!*rui(0)
2*rui(2)=su(2)*1!*rui(1)+su(3)*2!*rui(0)
3*rui(3)=su(2)*1!*rui(2)+su(3)*2!*rui(1)+su(4)*3!*rui(0)
4*rui(4)=su(2)*1!*rui(3)+su(3)*2!*rui(2)+su(4)*3!*rui(1)+su(5)*4!*rui(0)
5*rui(5)=su(2)*1!*rui(4)+su(3)*2!*rui(3)+su(4)*3!*rui(2)+su(5)*4!*rui(1)+su(6)*5!*rui(0)

これは綺麗に並んでいるので、間違えようがありません。r が奇数の時、su(r)=0 ですから、この式の半分は計算する必要はありません。並んでいることを確認するために書いただけです。また、rui(0)=1 であるとしておきます。

これを用いて　rui(5) を計算しますと、
rui(5)=163879/209018880
となります。

これを用いてΣlog(n) の続きを求めると、次の項数は　1/1260 という綺麗な数になります。これならば、何らかの法則があるだろうという気がしてきます。

そこで、これを静かに眺めていますと、偶数はすべて　0 になっていることに注目しました。これは su(r) の表にも現れる現象です。しかも　1/12 の項数もあります。それゆえ、試しにこの式をsu(r)を使って表現してみましょう。
Σlog(n) = ....+ log(e) * (su(2)/n + su(4)*2/n^3 + su(6)*4!/n^5...
となりました。うまい具合にマイナス符号まで消えてくれました。これですと、次の n^7 がどうなるか疑問の余地はありません。つまり、 su(8)*6!/n^7 と並ぶわけです。これを計算すると、-1/1680 となります。

さて、さすがに少し疲れてきました。あと、これを検証する作業をしなければなりませんが、それはいずれと言うことにして、ここで、これを結論と言うことにさせていただきます。


　　　　　　　結論---階乗(factorial)の公式---

Σlog(n) = n*log(n) - n*log(e) + 1/2*log(n) + 1/2*log(2*pi) + log(e) * ( su(2)/n + su(4)*2/n^3 + su(6)*4!/n^5 + su(8)*6!/n^7 + su(10)*8!/n^9 + .... )　と書くことができると言うことです。

これを前提にすると、 n! も綺麗に表現できます。
n! = (n/e)^n * √(2*pi*n) * e^( 1/12n - 1/360n^3 + 1/1260n^5 ... .)

つまり、
n! = (n/e)^n * √(2*pi*n) * e^ ｛ su(2)/n + su(4)*2/n^3 + su(6)*4!/n^5 + su(8)*6!/n^7 .... ｝
となります。

この公式を使うと、階乗計算が非常に楽になります。


（追記）　　　この公式はベルヌーイ数を用いても表現することができます。
それについては、あとのページで説明してありますので、そちらも参照してください。

（追記その２）　　　再計算が面倒なので、何とかこのまま検証する方法はないかと考えていたのですが、(-0.5)! = √pi を使って、計算結果が合うかどうかを調べてみました。上記の式を使って (9.5-10)! = 9.5!/9.5/8.5/7.5/..../0.5 ですから、まずは 9.5!を計算して、それから (-0.5)!を求めます。その結果 ans が√pi になっていれば良いと言うことです。
有効桁 ak=100 で、ベルヌーイ数30まで計算したところ、ans/√pi = 0.99999999999999999999....... と２２個の９が並びました。
(55.5-56)!=55.5!/55.5/54.5/....../0.5 を使って、有効桁 ak=100 で、ベルヌーイ数40まで計算したところ、ans/√pi = 1.000000000000000000000000....... と、５４個０が並びました。これで充分ではないでしょうか。


（追記その３）　　　もうひとつ別の階乗の公式を見つけました。
x! = e^{ -x*eu + x^2*k(2)/2 - x^3*k(3)/3 + x^4*k(4)/4 - x^5*k(5)/5 + x^6*k(6)/6 . . . }　
xの絶対値が小さいときは収束して値を持ちます。詳しい説明は「実数値での階乗計算方法」の最後にある「追加」を参照してください。