FC2ブログ

# ガンマ関数(近似) Gamma Function Approximation
# 引数 値 ($X)
# 戻り値 ガンマ関数(近似) ($GammaFunction)
sub GAMMAFUNCTION{
my ($X) = @_;
my $GammaFunction = 0;
my $Factorial = 1;
my $Diff = abs($X - int($X));
my $Pi = atan2(1, 1) * 4;

# 値の確認(近似しない)
if(((0 < $X) && ($X != 0.5) && ($X < 1)) || (($X < 0) && ($Diff != 0.5))){
return 0;
}

if($Diff == 0){
for(my $i = $X - 1; $i >= 2; $i--){
$Factorial *= $i;
}

# X-1の階乗 ガンマ関数 Gamma Function
$GammaFunction = $Factorial;
}elsif($Diff == 0.5){
my $Start = ($X > 0 ? ((2 * int($X)) - 1) : ((2 * int(abs($X) + 1)) - 1));
my $Temp = 0;

for(my $i = $Start; $i > 1; $i -= 2){
$Factorial *= $i;
}

if($X > 0){
$Temp = ($Factorial / (2 ** int($X))) * sqrt($Pi);
}else {
$Temp = ((((-1) ** int(abs($X) + 1)) * (2 ** int(abs($X) + 1))) / $Factorial) * sqrt($Pi);
}

# 半整数 ガンマ関数 Gamma Function
$GammaFunction = $Temp
}else {
my $X1 = 1 / (12 * $X);
my $X2 = 1 / (288 * ($X * $X));
my $X3 = 139 / (51840 * ($X * $X * $X));
my $X4 = 571 / (2488320 * ($X * $X * $X * $X));
my $X5 = 163879 / (209018880 * ($X * $X * $X * $X * $X));
my $X6 = 5246819 / (75246796800 * ($X * $X * $X * $X * $X * $X));
my $X7 = 534703531 / (902961561600 * ($X * $X * $X * $X * $X * $X * $X));
my $X8 = 4483131259 / (86684309913600 * ($X * $X * $X * $X * $X * $X * $X * $X));
my $X9 = 432261921612371 / (514904800886784000 * ($X * $X * $X * $X * $X * $X * $X * $X));

# スターリングの近似 Stirling's Approximation
# ガンマ関数(近似) Gamma Function Approximation
$GammaFunction = sqrt(2 * $Pi) * ($X ** ($X - (1 / 2))) * exp(-$X) * (1 + $X1 + $X2 - $X3 - $X4 + $X5 + $X6 - $X7 - $X8 + $X9);
}

return $GammaFunction;
}


ログガンマ関数を使用 1
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/GammaFunction3.txt

C言語による最新アルゴリズム事典のログガンマ関数を使用
ログガンマ関数の計算式が少し違う

# ガンマ関数(近似) Gamma Function Approximation
# 引数 値 ($X)
# 戻り値 ガンマ関数(近似) ($GammaFunction)
sub GAMMAFUNCTION{
my ($X) = @_;
my $GammaFunction = 0;
my $Factorial = 1;
my $Diff = abs($X - int($X));
my $Pi = atan2(1, 1) * 4;
my $Temp = 0;

# 値の確認
if(($X <= 0) && ($Diff == 0)){
return 0;
}

if($Diff == 0){
for(my $i = $X - 1; $i >= 2; $i--){
$Factorial *= $i;
}

# X-1の階乗 ガンマ関数 Gamma Function
$GammaFunction = $Factorial;
}elsif($Diff == 0.5){
my $Start = ($X > 0 ? ((2 * int($X)) - 1) : ((2 * int(abs($X) + 1)) - 1));

for(my $i = $Start; $i > 1; $i -= 2){
$Factorial *= $i;
}

if($X > 0){
$Temp = ($Factorial / (2 ** int($X))) * sqrt($Pi);
}else {
$Temp = ((((-1) ** int(abs($X) + 1)) * (2 ** int(abs($X) + 1))) / $Factorial) * sqrt($Pi);
}

# 半整数 ガンマ関数 Gamma Function
$GammaFunction = $Temp
}else {
if($X > 0){
$Temp = exp(&LOGGAMMAFUNCTION($X));
}else {
$Temp = $Pi / (sin($Pi * $X) * exp(&LOGGAMMAFUNCTION(1 - $X)));
}

# ガンマ関数(近似) Gamma Function Approximation
$GammaFunction = $Temp
}

return $GammaFunction;
}

# ログガンマ関数(近似) Log Gamma Function Approximation
# 引数 値 ($x)
# 戻り値 ログガンマ関数(近似) ($LogGammaFunction)
sub LOGGAMMAFUNCTION{
my ($x) = @_;
my $LogGammaFunction = 0;
my @Bernoulli = ((1 / 12), (1 / 360), (1 / 1260), (1 / 1680), (1 / 1188), (691 / 360360), (1 / 156), (3617 / 122400), (43867 / 244188), (174611 / 125400), (77683 / 5796), (236364091 / 1506960));
my $Pi = atan2(1, 1) * 4;
my $X = $x;
my $Sum = 0;
my $Sign = 0;
my $Power1 = 1;
my $Power2 = 1;
my $Count = @Bernoulli - 1;

for($X = $x; $X <= $Count; $X++){
$Power1 *= $X;
}
$Power2 = 1 / ($X * $X);

for(my $i = $Count; $i >= 0; $i--){
# 符号
$Sign = (($i % 2) == 0 ? 1: -1);

$Sum = $Sum + ($Sign * $Bernoulli[$i]);
$Sum = $Sum * $Power2 if($i != 0);
}
$Sum = $Sum / $X;

# スターリングの近似 Stirling's Approximation
# ログガンマ関数(近似) Log Gamma Function Approximation
$LogGammaFunction = ((1 / 2) * log(2 * $Pi)) - log($Power1) - $X + (($X - (1 / 2)) * log($X)) + $Sum;

return $LogGammaFunction;
}


参考URL
Gamma Function -- from Wolfram MathWorld
ガンマ関数 - 高精度計算サイト

一言
ログガンマ関数の方は、0付近だと値が近似しない
『C言語による最新アルゴリズム事典』にあるログガンマ関数だと近似する
ちがいがわからん
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

Author:雨宮
Firefoxを使用しているので気づかなかったけど、IE6でソースコードを上手くコピーできない

5/3
携帯用ならIE6でもソースコードをコピーできる
携帯用

検索フォーム


あわせて読みたいブログパーツ
一寸先は闇 RSS