FC2ブログ
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

# 対数積分 li(x) Logarithmic Integral li
# 引数 変数 ($X)
# 戻り値 対数積分 En(x) ($LogarithmicIntegralli)
sub LOGARITHMICINTEGRALLI{
my ($X) = @_;
my $LogarithmicIntegralli = 0;

# 変数の確認
if($X < 0){
return "Error";
}

if($X == 0){
# 対数積分 li(x) Logarithmic Integral li
$LogarithmicIntegralli = 0;
}
elsif($X == 1){
# 対数積分 li(x) Logarithmic Integral li
$LogarithmicIntegralli = "Inf";
}else {
# 対数積分 li(x) Logarithmic Integral li
$LogarithmicIntegralli = &EXPONENTIALINTEGRALEI(log($X) / log(exp(1)));
}

return $LogarithmicIntegralli;
}

# 指数積分 Ei(x) Exponential Integral Ei
# 引数 変数 ($X)
# 戻り値 指数積分 En(x) ($ExponentialIntegralEi)
sub EXPONENTIALINTEGRALEI{
my ($X) = @_;
my $ExponentialIntegralEi = 0;
my $EulersConstant = 0.577215665;
my $SimpsonsRule = 0;
my $h = 0.1;
my $i = 1;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;
my $PrevSum = 0;
my $Factorial = 1;
my $Limit = 100;
my $Epsilon = 1.0e-20;

if($X == 0){
# 指数積分 En(x) Exponential Integral
$ExponentialIntegralEi = "-Inf";

return $ExponentialIntegralEi;
}

if(0 < $X){
for(my $i = 1; $i < $Limit; $i++){
# 階乗
$Factorial *= $i;

$Sum += (($X ** $i) / ($i * $Factorial));
}

# 指数積分 Ei(x) Exponential Integral Ei
$ExponentialIntegralEi = $EulersConstant + (log($X) / log(exp(1))) + $Sum;
}else {
do {
# 刻み幅 $h は適当
my $T = -$X + ($i * $h);
my $Number = ($i % 2 == 0 ? 2 : 4);
my $Temp = &INTEGRAND($T);

$PrevSum = $Sum;
$Sum += $Number * $Temp;
$i++;
} while(($Sum - $PrevSum) >= $Epsilon);

$SumA = &INTEGRAND(-$X);
$SumB = &INTEGRAND((($h * $i) - $X));

# シンプソンの公式 Simpson's Rule
$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);

# 指数積分 Ei(x) Exponential Integral Ei
$ExponentialIntegralEi = -$SimpsonsRule;
}

return $ExponentialIntegralEi;
}

# 被積分関数 Integrand
# 引数 変数 ($T)
# 戻り値 被積分関数 ($Integrand)
sub INTEGRAND{
my ($T) = @_;
my $Integrand = 0;

# 被積分関数 Integrand
$Integrand = exp(-$T) / $T;

return $Integrand;
}


参考URL
対数積分 - Wikipedia
対数積分 li(x) - 高精度計算サイト
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。