Approximating pi in Perl - what am I doing wrong?

Question

I am trying to approximate pi using the Ramanujan algorithm:

It should compute the sums until the last sum is less than 1e-15.

This was just supposed to be for fun and occupy at most half an hour of my time... but my code is not producing anything close to pi and I'm not sure why. Quite possibly something silly that I've overlooked, but not sure!

Just a note: I'm starting $k at 1 because 0 breaks my factorialsub and from my calculations k=0 would return 0 anyway.

I realise the code can be more efficiently written; I wrote it out as simply as possible to see if I could understand where I was going wrong. 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 the item is 1103. So start $k from 0, not 1. For this You should modify factorial:

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

There is no need to pass 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 interesting comment about the possible tail-call optimization issue in the original script. More about it at the end of this answer.)

Not the $sum values should be less then a threshold, but the kth difference item. So 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 anyway always recursively calling the factorial and calculating the 396 power of 4k is not effective, so they can be just left 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;
}

The result is:

my pi is: 3.14159265358979

I did some measures. Approx_pi function was run 1 000 000 times. The fixed original version takes 24 seconds, the other one 5 sec. For me it is somehow interesting that, $F1**4 is faster than $F1*$F1*$F1*$F1.

---

Some words about the factorial. Because of Mark's comment I tried different implementations, to find the fastest solution. 5 000 000 loops were run for different implementations to calculate 15!:

Recursive version

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

46.39 sec

Simple loop version

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

16.29 sec

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

Recursion with storing intermediate values

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

3.87 sec

Loop with storing intermediate values. The speed is the same as previous as only the first time is has to be run.

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