我试图使用
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那么项目是1103.所以从0开始$k,而不是1.为此你应该修改阶乘:
- sub factorial {
- my ($i,$total) = @_;
- return $total if $i <= 1;
>没有必要在factorial中传递产品.它可能是这样的:
- 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++;
- }
>但无论如何总是递归地调用阶乘并且计算4k的396次幂是无效的,所以它们可以被忽略.
- 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秒
>使用call-tail优化的递归(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]];
- }