Per approximation in Perl - what am I doing wrong?

Question

I am trying to approximate pi using the Ramanujan algorithm:

He must calculate the amounts until the last amount is less than 1e-15 .

It was just meant to be fun and take no more than half an hour of my time ... but my code doesn't produce anything close to pi, and I'm not sure why. Maybe something stupid that I forgot, but not sure!

Just a note: I start $k at 1 because 0 breaks my factorial sub and from my calculations k = 0 will return 0 anyway.

I understand that code can be written more efficiently; I wrote this as simple as possible to understand if I can understand what I'm wrong about. Any help appreciated!

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

Answer 1

UPDATED

There are several problems with the script.

if k==0 , then element 1103 . So, run $k from 0, not 1. To do this, you need to change factorial :

 sub factorial { my ($i, $total) = @_; return $total if $i <= 1; 

There is no need to transfer the product in factorial. It could be something like this:

 sub fact { my $n = shift; return 1 if $n == 0 || $n ==1; return $n * fact($n -1); } 

(See Mark Reed’s wonderful comment about possible tail optimization problems in the original script. For more on this, see the answer at the end.)

The values ​​of $sum should not be lower than the threshold value, but of the k-th difference element. Therefore, in approx_pi you should use something like this:

 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++; } 

But in any case, the factorial call is always recursive and the calculation of 396 power of 4k inefficient, so you can simply turn them off.

 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; } 

Result:

 my pi is: 3.14159265358979 

I did some action. The approx_pi function has approx_pi run 1,000,000 times. The fixed original version takes 24 seconds, the other takes 5 seconds. For me it’s somehow interesting that $F1**4 faster than $F1*$F1*$F1*$F1 .

---

A few words about factorial. Due to Mark's comment, I tried different implementations to find the fastest solution. 5,000,000 cycles were performed for various implementations to calculate 15! :

Recursive version

 sub rfact; sub rfact($) { return 1 if $_[0] < 2; return $_[0] * rfact $_[0] - 1 ; } 

46.39 s

Simple version of the loop

 sub lfact($) { my $f = 1; for(my $i = 2; $i <= $_[0]; ++$i) { $f *= $i } return $f; } 

16.29 s

Recursion with call-tail optimization (call _fact 15,1 ):

 sub _fact($$) { return $_[1] if $_[0] < 2; @_ = ($_[0] - 1, $_[0] * $_[1]); goto &_fact; } 

108.15 sec

Intermediate Recursion

 my @h = (1, 1); sub hfact; sub hfact($) { return $h[$_[0]] if $_[0] <= $#h; return $h[$_[0]] = $_[0] * hfact $_[0] - 1; } 

3.87 s

Loop preserving intermediate values. The speed is the same as the previous one, since it should be performed only for the first time.

 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]]; }