FC2ブログ

# ポリガンマ関数(近似) Poly Gamma Function Approximation
# 引数 微分階数 値 ($n, $x)
# 戻り値 ポリガンマ関数(近似) ($PolyGammaFunction)
sub POLYGGAMMAFUNCTION{
my ($n, $x) = @_;
my $PolyGammaFunction = 0;
my $PolyGamma = 0;
my @Bernoulli = (0.166666666666667, -0.0333333333333333, 0.0238095238095237, -0.0333333333333323, 0.0757575757575658, -0.253113553113418, 1.16666666666417, -7.09215686268424, 54.971177942975, -529.12424235158, 6192.12318500442, -86580.2529233161);
my $N = int($n);
my $X = $x;
my $Sum = 0;
my $Temp1 = 1;
my $Temp2 = 0;
my $Power = 1;
my $Number = 2 * @Bernoulli;
my $Count = @Bernoulli - 1;

# 微分階数の確認
if($N < -1){
return 0;
}

if($N == -1){
# ログガンマ関数(近似) Log Gamma Function Approximation
$PolyGamma = &LOGGAMMAFUNCTION($X);
}elsif($N == 0){
# ディガンマ関数(近似) Digamma Function Approximation
$PolyGamma = &DIGGAMMAFUNCTION($X);
}else {
for(my $i = 1 - $N; $i < 0; $i++){
$Temp1 *= $i;
}

for($X = $x; $X <= $Count; $X++){
$Sum += 1 / ($X ** ($N + 1));
}
$Power = ($X * $X);

$Temp2 = $Bernoulli[$Count];
for(my $i = $Count - 1; $i >= 0; $i--){
$Temp2 = ($Temp2 * ((($N + $Number - 1) * ($N + $Number - 2)) / ($Number * ($Number - 1) * $Power))) + $Bernoulli[$i];
$Number = $Number - 2;
}
$Temp2 = ($Temp2 * ((($N + $Number - 1) * ($N + $Number - 2)) / ($Number * ($Number - 1) * $Power))) + (0.5 * ($N / $X)) + 1;

# ポリガンマ関数(近似) Poly Gamma Function Approximation
$PolyGamma = $Temp1 * (($Temp2 / ($X ** $N)) + ($N * $Sum));
}

# ポリガンマ関数(近似) Poly Gamma Function Approximation
$PolyGammaFunction = $PolyGamma;

return $PolyGammaFunction;
}

# ディガンマ関数(近似) Digamma Function Approximation
# 引数 値 ($x)
# 戻り値 ディガンマ関数(近似) ($DigammaFunction)
sub DIGGAMMAFUNCTION{
my ($x) = @_;
my $DigammaFunction = 0;
my @Bernoulli = ((1 / 12), (1 / 120), (1 / 252), (1 / 240), (5 / 660), (691 / 32760), (1 / 12), (3617 / 8160), (43867 / 14364), (174611 / 660), (77683 / 276), (236364091 / 65520));
my $X = $x;
my $Sum1 = 0;
my $Sum2 = 0;
my $Sign = 0;
my $Power = 1;
my $Count = @Bernoulli - 1;

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

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

$Sum2 = ($Sum2 + ($Sign * $Bernoulli[$i])) * $Power;
}
$Sum1 += $Sum2 + (0.5 / $X);

# ディガンマ関数(近似) Digamma Function Approximation
$DigammaFunction = log($X) - $Sum1;

return $DigammaFunction;
}

# ログガンマ関数(近似) 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;

# ログガンマ関数(近似) Log Gamma Function Approximation
$LogGammaFunction = ((1 / 2) * log(2 * $Pi)) - log($Power1) - $X + (($X - (1 / 2)) * log($X)) + $Sum;

return $LogGammaFunction;
}


参考URL
『C言語による最新アルゴリズム事典』
ポリガンマ関数 - 高精度計算サイト
Polygamma function - Wikipedia, the free encyclopedia
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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