7

我正在尝试使用Ramanujan 算法来近似 pi :

pi 的近似值

它应该计算总和,直到最后一个总和小于1e-15

这应该是为了好玩,最多占用我半个小时的时间……但我的代码没有产生任何接近 pi 的东西,我不知道为什么。很可能我忽略了一些愚蠢的事情,但不确定!

请注意:我从$k1 开始,因为 0 打破了我的factorialsub 并且从我的计算 k=0 无论如何都会返回 0。

我意识到可以更有效地编写代码;我尽可能简单地把它写出来,看看我是否能理解我哪里出错了。任何帮助表示赞赏!

#!/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";
4

1 回答 1

13

更新

脚本有几个问题。

  • 如果k==0那么项目是1103。所以$k从 0 开始,而不是 1。为此你应该修改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.


    关于阶乘的一些话。由于 Mark 的评论,我尝试了不同的实现,以找到最快的解决方案。为不同的实现运行了 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 秒

  • 带有调用尾优化的递归(调用_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]];
    }
    
  • 于 2013-05-16T10:49:09.097 回答