在 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 打破了我的 factorial
sub 并且根据我的计算,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 = 1
到 infinity
的总和”
【参考方案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 < 2;
。在factorial
中传递总数大概是为了允许尾递归优化。
@MarkReed:您可能是对的,但无论如何都不需要factorial
函数... :)
@MarkReed:如果图片中有“尾递归优化”,我会推荐一个具有一个参数的代理函数,然后调用一个有两个参数(偶数函数)的内部函数。我认为它会更受用户欢迎。
这是一个比我预期的更好的响应!感谢您的时间,非常有见地,给了我很多实验。 :)
@bladepanthera:好的。在循环的第一次运行$k = 0; $P2 = 1103
。在第二个循环$k = 1; $P2 = 1103 + 26390 * 1
。第三个:$k = 2; $P2 = 1103 + 26390 * 2
等。在第 N 个循环中:$k = <N-1>; $P2 = 1103 + 26390 * <N-1>
我用在每个循环中添加 26390
替换了乘法,这更“便宜”。类似的乘法。我不是总是计算fact($k)
,而是将一个变量($F1
)与$k
相乘,因此每个循环中的ist 值将是fact($k)
。同样,乘法要便宜得多,然后在每个循环中计算fact($k)
。以上是关于在 Perl 中近似 pi - 我做错了啥?的主要内容,如果未能解决你的问题,请参考以下文章