Perl6 Pi and continued fractions

I was looking through the code of the Rakudo implementation of Perl6 where I noticed that it defines pi as my constant pi = 3.14159_26535_89793_238e0; ( with an alias my constant π := pi; )
I immediately remembered that for a subset of fractional numbers Perl6 has a type that stores them without the loss of precision that generally accompanies floating point math. That type is of course the Rat (and FatRat) type. So of course I type pi.Rat into the REPL, it then prints 3.141593 which is obviously nowhere near as precise as the result of just typing pi into the REPL 3.14159265358979.
I wanted to see the numerator and denominator values that Perl chose to use so I typed pi.Rat.perl and got <355/113>, which is nowhere near as precise as the Rat type is capable of handling. That was just the entrance to the rabbit hole.

So I went searching online and found the Wikipedia article on approximations of pi. Which lead to the article of continued fractions, that has the sequence [3;7,15,1,292,1,1,1,2,1,3,1,…]. I won’t go into too much detail about continued fractions in this post, for more information I would recommend the Wikipedia article.

That bit about ... had me curious so I looked for a longer list, which lead to a pdf of the first 250 elements.


I think that is enough backstory, time to look at the code that I wrote to play around with this new bit of information.

I started by creating an array of the elements of the continued fraction.

#!/usr/bin/env perl6
use v6;

my @continued-fraction = <
  3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2 1 84 2 1 1 15 3 13 1 4 2 6 6 99 1
  2 2 6 3 5 1 1 6 8 1 7 1 2 3 7 1 2 1 1 12 1 1 1 3 1 1 8 1 1 2 1 6 1 1 5 2 2 3
  1 2 4 4 16 1 161 45 1 22 1 2 2 1 4 1 2 24 1 2 1 3 1 2 1 1 10 2 5 4 1 2 2 8 1
  5 2 2 26 1 4 1 1 8 2 42 2 1 7 3 3 1 1 7 2 4 9 7 2 3 1 57 1 18 1 9 19 1 2 18 1
  3 7 30 1 1 1 3 3 3 1 2 8 1 1 2 1 15 1 2 13 1 2 1 4 1 12 1 1 3 3 28 1 10 3 2
  20 1 1 1 1 4 1 1 1 5 3 2 1 6 1 4 1 120 2 1 1 3 1 23 1 15 1 3 7 1 16 1 2 1 21
  2 1 1 2 9 1 6 4 127 14 5 1 3 13 7 9 1 1 1 1 1 5 4 1 1 3 1 1 29 3 1 1 2 2 1 3 1
>\ >>.Int; # converts each element into an integer

To make the rest of the code simpler I created a helper function to get the reciprocal of a number.

multi sub reciprocal ( Rational ::R $rat ){ # Rat or FatRat
  R.new( $rat.denominator, $rat.numerator );
}
multi sub reciprocal ( Int $integer ){
  FatRat.new( 1, $integer );
}

Now that I think about it I could have created a new operator by naming it prefix:<1/>, but that would be a new rabbit hole into precedence and associativity of operators.

Anyway the basic formula is start from the right adding the reciprocal of the last number to the second to last, then adding the reciprocal of the previous result to the third to last, etc until you have added the first number.

@continued-fraction.reverse.reduce( ->$a,$b{ reciprocal($a) + $b } );

The result of which is:

11470607307161132716588859789437271840000798468219508525085448960062402821241970285202775001996784319885772780138414238670402064758
/
3651207706401417758577438077383177726074790491037819709434096621864877255825019944180015629313589782838738505101104441533461516175

A rather difficult to comprehend number. It is also way more precise than is necessary. Its accurate to more than 100 decimal places.

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955

To play around with the precision it would be nice if that algorithm was packaged up in a sub that took a count of elements as the only parameter. ( zero based )

multi sub pi-fraction ( 0 ){ Rat.new( 3 ) }
multi sub pi-fraction ( Int $index ){
  @continued-fraction[0..$index].reverse.reduce(->$a,$b{ reciprocal($a) + $b });
}

Let’s print out the first four fractions

for ^4 { pi-fraction($_).Rat.perl.say }
3.0
<22/7>
<333/106>
<355/113>

At this point it would be a good idea to create some MAIN subroutines. I also want to control how the numbers are printed, so I add some helper subs.

sub as-fraction ( Rational $rat ){
  $rat.numerator ~ ' / ' ~ $rat.denominator
}

sub display ( Rational $number ){
  my $as-fraction = as-fraction $number;
  my $display-width = $as-fraction.chars; # good enough for the chars we are using
  my $padding = ' ' x 8 - ($display-width % 4);
  $padding ~= ' ' x 4 if $display-width < 12;
  $padding ~= ' ' x 4 if $display-width < 8;
  return $as-fraction ~ $padding ~ $number;
}

#| list the first <count> fractions
multi sub MAIN ( 'list', $count = 0 ){
  my $total-count = @continued-fraction.elems;
  my @index = do given $count {
    when 0 { ^$total-count }
    when $_ < 0 { ($total-count + $count) ..^ $total-count }
    default { ^$count }
  }
  for @index -> $index {
    my $pi = pi-fraction $index;
    say display $pi;
  }
}

So now we can run it like this ./pi.p6 list 4

3 / 1               3
22 / 7              3.142857
333 / 106           3.141509
355 / 113           3.141593

I wanted to know what the largest one that could be stored as a Rat so I added this:

#| Display the most precise Rat
multi sub MAIN ( 'max', 'Rat' ) {
  my $max-rat;
  my $max-index;

  for @continued-fraction.keys -> $index {
    my $current = pi-fraction $index;

    last unless $current.Rat; # stop if it can't be represented as a Rat

    $max-rat = $current;
    $max-index = $index;
  }
  say '#', $max-index;
  say display $max-rat;
}
$ ./pi.p6 max Rat
#32
2646693125139304345 / 842468587426513207        3.1415926535897932385

If you are new to Perl6 you might not know why I put comments before the MAIN subs that start with #|. The reason should become apparent if you attempt to run it with the wrong arguments, or in this case, no arguments.

$ ./pi.p6
Usage:
  ./pi.p6 list [<count>] -- list the first <count> fractions
  ./pi.p6 max Rat -- Display the most precise Rat

Without those comments it will just print:

$ ./pi.p6
Usage:
  ./pi.p6 list [<count>]
  ./pi.p6 max Rat

I will probably continue to play around with this some more, as programming in Perl6 is a lot of fun.

2 Comments

This is tangentially related as it isn't looking at rationals.

I've seen a few places that use something like: Pi(n) that gives n digits of Pi. I added that to the Perl5 ntheory library a couple months ago, though I'm sure a clean Perl6 implementation wouldn't be as fractured. The performance difference is minor at a few digits, but at larger sizes it's quite dramatic (e.g. for 1M digits it can be a couple seconds vs. hours vs. weeks).

  • Chudnovsky / Ramanujan binary splitting. Looks like 4-5x faster than AGM, though more complicated. Pari uses this (they have AGM commented out after doing both and comparing). This is also what is used in the example on the GMP library site. I've tried it with Math::BigFloat and it was slower than AGM for me (probably due to more but smaller operations, which kills you with overhead).

  • AGM (Gauss-Brent-Salamin). Fast and easy (~15 lines of Perl 5), not too many lines with GMP. Obviously the Perl5 implementation is oodles slower than C+GMP, but it still has the good growth rate. MPFR uses this, as does ntheory's GMP and Perl5 code. There's a patch for Math::BigFloat to use this for larger inputs, but it hasn't been accepted. There's a Perl6 implementation on RosettaCode, but pretty slow.

  • Machin type formulas. Math::BigFloat uses this. It's pretty good for small amounts but AGM kills it for many digits because the growth rate is better. With the default Calc backend there are faster ways at all sizes.

  • Spigot style. The nice thing is that it doesn't need bigints and can be easily done in standard C or Perl. There is something similar for Perl6 on RosettaCode. Terrible growth rate vs. the others, but it's pretty fast for small sizes: at 2000 digits the pure Perl version is over 10x faster than either Machin or AGM with Math::BigFloat's Calc backend.

As for a rational, it would be nice to see something that produced a Rational or FatRat of the desired accuracy. Rationals are limited to 64-bit denominators (S02) so the fixed continuous fractions would seem fine. I think using the Chudnosky formula plus the RC trick for square roots could do this for FatRats.

As an aside, I'm happy that Perl6 is using arbitrary size ints, but it seems Rat vs. FatRat is akin to the native vs. bigint that drives me nuts with Perl5 and Perl6 gets rid of. For a library that expects exact answers, now I have to constantly deal with input and output conversion, and performance vs. correctness tradeoffs.

> "In the end I'm not sure he explained why pi.Rat is 355/113."

The comments in the Rat method in class Num reveal that .Rat returns an approximation of its invocant. It does so by finding "convergents of the continued fraction" and stopping when the result "has less error than any Rational with a smaller denominator".

The target maximum error defaults to 1.0e-6 if you don't pass an arg corresponding to the optional first parameter, as can be seen from the method's signature:

    method Rat(Num:D: Real $epsilon = 1.0e-6, :$fat)

(Both parameters are optional. If a value is specified for a parameter then that parameter is optional. And named arguments are optional unless marked in the signature as required.)

Thus:

    say pi.Rat(1.0e-1).perl  # <22/7>
    say pi.Rat(1.0e-4).perl  # <333/106>
    say pi.Rat(1.0e-5).perl  # <355/113>
    say pi.Rat(1.0e-6).perl  # <355/113>
    say pi.Rat().perl        # <355/113>
    say pi.Rat.perl          # <355/113>
    say pi.Rat(1.0e-7).perl  # <103993/33102>
    say pi.Rat(1.0e-15).perl # <80143857/25510582>
    say pi.Rat(1.0e-16).perl # <245850922/78256779>
    say pi.Rat(1.0e-99).perl # <245850922/78256779>

..........................................................................................................................................

The constant pi is defined in Num.pm as:

    my constant pi = 3.14159_26535_89793_238e0;

The appended E_notation makes the value be a Num literal.

pi could be specified as a Rat literal:

    my constant pi = 3.14159_26535_89793_238;
    # the stored value is a Rat
    say pi.Rat.perl        # 3.141592653589793238 
    say pi.perl            # 3.141592653589793238

When a Rat can be represented exactly as a decimal, it is. Thus, in contrast to the approximations of pi within epsilon 0.1 and less (see above), one gets:

    my constant pi = 3.14159_26535_89793_238e0;
    say pi.Rat(1.0e0).perl # 3.0

But of course, while 3.141592653589793238 is accurate to 18 decimal places, it's still an approximation of an irrational number so a floating point Num is much more appropriate than a Rat for pi.

..........................................................................................................................................

Perhaps things ought to be tweaked. Maybe it would be possible to have the exponent bit of the $epsilon error value (the '6' in '1.0e-6') default to the number of significant digits in a Num without that being a performance disaster. I plan to ask about this on #perl6. Hope to see you there. Maybe you'll beat me to it for asking the question. :)

Leave a comment

About Brad Gilbert

user-pic I blog about Perl.