Stirlingの公式
Stirlingの公式： ln N! N ln N - N
これを数値計算によって気軽に確認するにはどうしたらよいだろう．例えばN=10^1, 10^2, ..., 10^10 についてそれぞれの値と相対誤差のグラフを書きたい（ことによれば対数スケールで）．
まぁ詳しい本にはもっと細かいStirlingの公式 ln N! (N+1/2) ln N - N + (ln2π)/2 が載っているのでそれで 糸冬 了 としてもいいんだが，試しにRubyスクリプトを書いてみた．（正規化のあたりが）あまり気軽ではないしグラフも書けないが，まぁいいや．Cでも書いたが，これは「Cは単純な数値計算ではRubyより速い」という事実の再確認に過ぎない．
いや，要するに私が言いたかったのは，IEEE 754 doubleを指数の幅で*1超えるような「小数」を気軽に扱えて，なおかつ超越数の演算が（せめてIEEE 754の浮動小数点数に対しては）使えて，欲を言えばグラフプロットができるような環境が欲しかったということ．次にその環境が入用になるのがいつかは知らないけどさ．
Mathematicaだと実装の負担と処理速度はどの程度なのだろう？
#! /usr/bin/ruby
include Math
def fac(n)
v = 1.0
e = 0
while n > 1
v *= n
if v > 1e100
v /= 1e100
e += 100
end
n -= 1
end
e0 = log10(v).to_i
v /= 10**e0
e += e0
return v, e
end
n = (ARGV[0] || 10).to_i
v, e = fac(n)
puts "#{n}! = #{v}E#{e}"
x = log(v) + e*log(10)
puts "ln #{n}! = #{x}"
d = log(2 * atan2(1,1)*4) / 2
y = (n+0.5) * log(n) - n + d
puts "(n+1/2)*ln n - n + d = #{y}"
puts "rel.err = #{(x-y)/y}"ただPenM 1.2GHz, colinux (64MB RAM)上のRuby 1.8.2でN=10^7だと20秒，N=10^8だと3分かかる．N=10^9はやりたくない．
#include <stdio.h>
/* #define __USE_BSD */
#include <math.h>
void fac(int n, double *vout, int *eout)
{
double v = 1.0;
int e = 0, e0;
while (n > 1) {
v *= n;
if (v > 1e100) {
v /= 1e100;
e += 100;
}
--n;
}
e0 = (int) floor(log(v)/log(10));
v /= pow(10, e0);
e += e0;
*vout = v;
*eout = e;
}
int main(int argc, char *argv[]) {
int n, e;
double v, x, y, d;
n = argc > 1 ? atoi(argv[1]) : 10;
fac(n, &v, &e);
printf("%d! = %fE%d\n", n, v, e);
x = log(v) + e*log(10.0);
printf("ln %d! = %30.16f\n", n, x);
d = log(2 * M_PI)/2;
y = (n+0.5) * log(n) - n + d;
printf("(n+1/2)*ln n - n + d = %30.16f\n", y);
printf("rel.err = %30.16f\n", (x-y)/y);
return 0;
} N=10^9で16秒，N=10^10で35秒．でもCは書きたくないのですよ...
