我正在尝试使用 Ramanujan algorithm 来近似 pi :
它应该计算总和,直到最后一个总和小于 1e-15
.
这只是为了好玩,最多占用我半个小时的时间……但是我的代码没有产生任何接近 pi 的结果,我不知道为什么。很可能我忽略了一些愚蠢的事情,但不确定!
请注意:我正在开始 $k
在 1 因为 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
那么该项目是 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 关于原始脚本中可能的 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]];
}
关于perl - 在 Perl 中逼近 pi - 我做错了什么?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16584635/