perl - Perl で pi を概算する - 何が間違っているのでしょうか?

Question

ラマヌジャン アルゴリズムを使用して pi を概算しようとしています。

最後の合計が 未満になるまで合計を計算する必要があり1e-15ます。

これは楽しみのためであり、せいぜい私の時間の 30 分を占めるはずでした... しかし、私のコードは pi に近いものを生成しておらず、その理由はわかりません。私が見落としていたばかげた何かかもしれませんが、確かではありません!

ただのメモ: $k0 は私のサブを壊し、factorial私の計算から k=0 はいずれにせよ 0 を返すので、私は 1 から始めています。

コードをより効率的に記述できることがわかりました。どこが間違っているのかを理解できるかどうかを確認するために、できるだけ簡単に書きました。どんな助けでも大歓迎です！

#!/usr/bin/perl
use warnings;
use strict;

sub approx_pi {
    my $const = (2 * sqrt(2)) / 9801;

    my $k = 1;
    my $sum = 0;
    while ($sum < 1e-15) {
        my $p1 = factorial((4 * $k), 1);
        my $p2 = 1103 + (26390 * $k);
        my $p3 = (factorial($k, 1))**4;
        my $p4 = 396**(4 * $k);

        $sum = $sum + ( ($p1 * $p2) / ($p3 * $p4) );

        $k++;
    }

    #print "Const: $const\nSum: $sum\n";
    return (1 / ($const * $sum));
}

sub factorial {
    my ($i, $total) = @_;
    return $total if $i == 1;

    $total = $total * $i;
    #print "i: $i total: $total\n";

    factorial($i-1, $total);
}

my $pi = approx_pi();
print "my pi is: $pi\n";

Answer 1

更新しました

スクリプトにはいくつかの問題があります。

その場合k==0、アイテムは1103です。したがって$k、1 ではなく 0 から開始します。このためには、以下を変更する必要がありますfactorial。

sub factorial {
    my ($i, $total) = @_;
    return $total if $i <= 1;

積を階乗で渡す必要はありません。それは次のようなものかもしれません：

sub fact {
  my $n = shift;
  return 1 if $n == 0 || $n ==1;
  return $n * fact($n -1);
}

(元のスクリプトで発生する可能性のある末尾呼び出し最適化の問題に関する Mark Reed の興味深いコメントを参照してください。詳細については、この回答の最後を参照してください。)

$sum値がしきい値未満であってはなりませんが、k 番目の差分項目です。したがって、approx_pi次のようなものを使用する必要があります。

my $Diff = 1;
while ($Diff > 1e-15) {
    my $p1 = factorial((4 * $k), 1);
    my $p2 = 1103 + (26390 * $k);
    my $p3 = (factorial($k, 1))**4;
    my $p4 = 396**(4 * $k);

    $Diff = ($p1 * $p2) / ($p3 * $p4);
    $sum += $Diff;

    $k++;
}

ただし、常に再帰的に を呼び出してfactorialを計算すること396 power of 4kは効果的ではないため、そのままにしておくことができます。

sub approx_pi {
    my $const = 4900.5 / sqrt(2);

    my $k = 0;
    my $k4 = 0;
    my $F1 = 1;
    my $F4 = 1;
    my $Pd = 396**4;
    my $P2 = 1103;
    my $P4 = 1;
    my $sum = 0;

    while (1) {
        my $Diff = ($F4 * $P2) / ($F1**4 * $P4);
        $sum += $Diff;
        last if $Diff < 1e-15;

        ++$k;
        $k4 += 4;
        $F1 *= $k;
        $F4 *= ($k4 - 3)*($k4 - 2)*($k4 - 1)*$k4;
        $P2 += 26390;
        $P4 *= $Pd;
    }

    return $const / $sum;
}

結果は次のとおりです。

my pi is: 3.14159265358979

いくつかの対策を行いました。Approx_pi関数は 1 000 000 回実行されました。修正された元のバージョンは 24 秒かかり、他のバージョンは 5 秒かかります。$F1**4私にとって、が よりも速いことは、なんとなく興味深いことです$F1*$F1*$F1*$F1。

---

階乗についてのいくつかの言葉。マークのコメントのため、最速の解決策を見つけるために、さまざまな実装を試しました。5 000 000 ループが実行され、さまざまな実装が計算されました15!。

再帰バージョン

sub rfact;
sub rfact($) {
    return 1 if $_[0] < 2;
    return $_[0] * rfact $_[0] - 1 ;
}

46.39秒

シンプルなループバージョン

sub lfact($) {
     my $f = 1;
     for(my $i = 2; $i <= $_[0]; ++$i) { $f *= $i }
     return $f;
}

16.29秒

コールテール最適化による再帰 (call _fact 15,1):

sub _fact($$) {
    return $_[1] if $_[0] < 2;
    @_ = ($_[0] - 1, $_[0] * $_[1]);
    goto &_fact;
}

108.15秒

中間値を格納する再帰

my @h = (1, 1);
sub hfact;
sub hfact($) {
    return $h[$_[0]] if $_[0] <= $#h;
    return $h[$_[0]] = $_[0] * hfact $_[0] - 1;
}

3.87秒

中間値を格納してループします。最初だけ実行する必要があるため、速度は以前と同じです。

my @h = (1, 1);
sub hlfact($) {
    if ($_[0] > $#h) {
      my $f = $h[-1];
      for (my $i = $#h + 1; $i <= $_[0]; ++$i) { $h[$i] = $f *= $i }
    }
    return $h[$_[0]];
}