在 Perl 中近似 pi - 我做错了啥?

Posted

技术标签:

【中文标题】在 Perl 中近似 pi - 我做错了啥?【英文标题】:Approximating pi in Perl - what am I doing wrong?在 Perl 中近似 pi - 我做错了什么? 【发布时间】:2013-05-11 04:20:46 【问题描述】:

我正在尝试使用Ramanujan algorithm 近似 pi:

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

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

请注意:我从 1 开始 $k,因为 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";

【问题讨论】:

“根据我的计算,k=0 无论如何都会返回 0。” - 跟我们讨论一下? 将 k = 0 放入等式中,您至少会得到一个分子零。我不确定k!,我读到如果k = 0,那么k!= 1(不知何故)但如果分子为0,则结果为零。除非我的基础数学不正确...... 哎呀,什么是k!如果 k = 0 ?? ....报废那个,这就是问题....谢谢@AakashM。一个糟糕的假设 没问题。请记住,数学家非常简洁,如果不需要 k = 0 术语,那么公式肯定是“从 k = 1infinity 的总和” 【参考方案1】:

更新

脚本有几个问题。

如果k==0,则项目为1103。所以从0开始$k,而不是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 关于原始脚本中可能存在的 tail-call optimization 问题的有趣评论。有关此问题的更多信息,请参见此答案的末尾。)

不是$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]];

【讨论】:

或只是return $total if $i &lt; 2;。在factorial 中传递总数大概是为了允许尾递归优化。 @MarkReed:您可能是对的,但无论如何都不需要factorial 函数... :) @MarkReed:如果图片中有“尾递归优化”,我会推荐一个具有一个参数的代理函数,然后调用一个有两个参数(偶数函数)的内部函数。我认为它会更受用户欢迎。 这是一个比我预期的更好的响应!感谢您的时间,非常有见地,给了我很多实验。 :) @bladepanthera:好的。在循环的第一次运行$k = 0; $P2 = 1103。在第二个循环$k = 1; $P2 = 1103 + 26390 * 1。第三个:$k = 2; $P2 = 1103 + 26390 * 2 等。在第 N 个循环中:$k = &lt;N-1&gt;; $P2 = 1103 + 26390 * &lt;N-1&gt; 我用在每个循环中添加 26390 替换了乘法,这更“便宜”。类似的乘法。我不是总是计算fact($k),而是将一个变量($F1)与$k 相乘,因此每个循环中的ist 值将是fact($k)。同样,乘法要便宜得多,然后在每个循环中计算fact($k)

以上是关于在 Perl 中近似 pi - 我做错了啥?的主要内容,如果未能解决你的问题,请参考以下文章

单对象 json 解析 - 我做错了啥?

在 Popover 中显示 ActionSheet - 我做错了啥?

django - Ajax 实时搜索,我做错了啥?

JSON.parse,我做错了啥?

核心数据多线程——我做错了啥

sectionName TableView - 我做错了啥?