# 第2種不完全ガンマ関数 Incomplete Gamma Function Of The Second Kind # 引数 値 値 ($A, $X) # 戻り値 第2種不完全ガンマ関数 ($IncompleteGammaFunction) sub COMPLEMENTARYERRORFUNCTIONSECOND{ my ($A, $X) = @_; my $IncompleteGammaFunction = 0; my $GammaFunction = 0; my $IncompleteGammaFunctionFirst = 0; if($X == 0){ $GammaFunction = &GAMMAFUNCTION($A); $IncompleteGammaFunctionFirst = 0; }else { $GammaFunction = &GAMMAFUNCTION($A); $IncompleteGammaFunctionFirst = &COMPLEMENTARYERRORFUNCTIONFIRST($A, $X); } # 第2種不完全ガンマ関数 Incomplete Gamma Function Of The Second Kind $IncompleteGammaFunction = $GammaFunction - $IncompleteGammaFunctionFirst; return $IncompleteGammaFunction; } # ガンマ関数(近似) 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; # ログガンマ関数(近似) Log Gamma Function Approximation $LogGammaFunction = ((1 / 2) * log(2 * $Pi)) - log($Power1) - $X + (($X - (1 / 2)) * log($X)) + $Sum; return $LogGammaFunction; } # 第1種不完全ガンマ関数 Incomplete Gamma Function Of The First Kind # 引数 値 値 ($A, $X) # 戻り値 第1種不完全ガンマ関数 ($IncompleteGammaFunction) sub COMPLEMENTARYERRORFUNCTIONFIRST{ my ($A, $X) = @_; my $IncompleteGammaFunction = 0; my @NewtonsMethod = (); my $LegendrePolynomial = 0; my $a = 0; my $b = $X; my $x = 0; my $N = 10; my $Sum = 0; my $Weight = 0; my $AB2 = ($a + $b) / 2; my $BA2 = ($b - $a) / 2; # 値の確認 if($X < 0){ return 0; } # ガウス-ルジャンドル数値積分 Legendre-Gauss Quadrature for(my $i = 1; $i <= $N; $i++){ my $tmp = 0; # ニュートン法 Newton's Method @NewtonsMethod = &NEWTONSMETHOD($N, $i); $x = $NewtonsMethod[0]; $LegendrePolynomial = $NewtonsMethod[1]; # 重み $Weight = (2 * (1 - ($x * $x))) / (($N * $LegendrePolynomial) * ($N * $LegendrePolynomial)); $x = $AB2 + ($BA2 * $x); $tmp = &INTEGRAND($A, $x); $Sum += $Weight * $tmp; } # 第1種不完全ガンマ関数 Incomplete Gamma Function Of The First Kind $IncompleteGammaFunction = $BA2 * $Sum; return $IncompleteGammaFunction; } # ニュートン法 Newton's Method # 引数 次数 位置 ($N, $K) # 戻り値 ニュートン法 ($NewtonsMethod) sub NEWTONSMETHOD{ my ($N, $K) = @_; my @NewtonsMethod = (); my @LegendrePolynomial = (); my $Pi = atan2(1, 1) * 4; my $X = cos(($Pi * ($K - 0.25)) / ($N + 0.5)); my $Num = 0; my $Den = 0; my $PrevX = 0; my $Limit = 100; my $Epsilon = 1.0e-20; # 計算 for(my $i = 0; $i < $Limit; $i++){ # ルジャンドル多項式 Legendre Polynomial @LegendrePolynomial = &LEGENDREPOLYNOMIAL($N, $X); # 分子 $Num = $LegendrePolynomial[0]; # 分母 $Den = (($N * ($LegendrePolynomial[1] - ($X * $LegendrePolynomial[0]))) / (1 - ($X * $X))); # 近似解 $PrevX = $X; $X = $X - ($Num / $Den); # ニュートン法 Newton's Method # [0] 近似解 [1] 一つ前のルジャンドル多項式の値 $NewtonsMethod[0] = $X; $NewtonsMethod[1] = $LegendrePolynomial[1]; # 収束判定 last if(abs($X - $PrevX) < $Epsilon); } return @NewtonsMethod; } # ルジャンドル多項式 Legendre Polynomial # 引数 次数 値 ($N, $x) # 戻り値 ルジャンドル多項式 (@LegendrePolynomial) sub LEGENDREPOLYNOMIAL{ my ($N, $x) = @_; my @LegendrePolynomial = ($x , 1); # 計算 for(my $i = 1; $i < $N; $i++){ my $tmp = ((((2 * $i) + 1) * $x * $LegendrePolynomial[0]) - ($i * $LegendrePolynomial[1])) / ($i + 1); # [0] 最新 [1] 一つ前 # ルジャンドル多項式 Legendre Polynomial $LegendrePolynomial[1] = $LegendrePolynomial[0]; $LegendrePolynomial[0] = $tmp; } return @LegendrePolynomial; } # 被積分関数 Integrand # 引数 変数 被積分関数 ($X, \@Integrand) # 戻り値 被積分関数 ($Integrand) sub INTEGRAND{ my ($A, $X) = @_; my $Integrand = ($X ** ($A - 1)) * exp(-$X); return $Integrand;