Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl Monk, Perl Meditation

Computing pi to multiple precision

by ambrus (Abbot)
on Sep 09, 2012 at 12:16 UTC ( #992580=perlmeditation: print w/replies, xml ) Need Help??

This post shows two different ways to compute pi to multiple precision efficiently using elementary arithmetic and square root. We're using the Math::BigFloat library for the multiple precision arithmetic.

Math::BigFloat already has an implementation of pi computation built in of course, but I'd like to show how to write our own.

#!perl use warnings; use strict; use 5.010; use Math::BigFloat try => "GMP"; my $digits = int($ARGV[0] // 0) || 60; { say "Newton iteration, using Taylor series of sine and cosine."; my $b = Math::BigFloat->new(0, undef, -$digits); my $x = $b + 0.5; my %x; while (!$x{$x}++) { my $m = 1; my $s = 0; my $c = 1; my $k = 0; while (0 != $m) { $s += $m = $m * $x / ++$k; $c += $m = -$m * $x / ++$k; } $x += (0.5 - $s) / $c; } say "pi = ", 6*$x, ";"; } { say "Using Taylor series of arctangent."; my $b = Math::BigFloat->new(0, undef, -$digits); my $x = 2 - sqrt($b + 3); my $f = -$x * $x; my $m = $x * ($b + 12); my $a = $m; my $k = 1; while (0 != $m) { $m *= $f; $k += 2; $a += $m / $k; } say "pi = ", $a, ";"; } __END__

Example output:

Newton iteration, using Taylor series of sine and cosine. pi = 3.141592653589793238462643383279502884197169399375105820974942; Using Taylor series of arctangent. pi = 3.141592653589793238462643383279502884197169399375105820974945;

Update: have I mentioned that this post is partly based on this mailing list post I wrote?

Replies are listed 'Best First'.
Re: Computing pi to multiple precision
by moritz (Cardinal) on Sep 09, 2012 at 14:45 UTC

    Ever since I've read about them, I've been fascinated by spigot algorithms for producing digits of pi.

    The basic idea for those algorithms is that most "interesting" transcendental numbers (like pi, e, ln(2)) have a pretty simple representation if the base is allowed (in a regular pattern) for each digit. Then the task of computing the first $N first digits is just that of a base conversion. And the fascinating part is that you can work with integers only-

    To stay a bit on topic, I've ported this C implementation of a spigot algorithm for pi to Perl 6:

    sub pi-spigot(Int $digits) { my $len = 1 + floor 10 * $digits / 3; my @a = 2 xx $len; my Int $nines = 0; my Int $predigit = 0; join '', gather for 1 .. ($digits + 1) -> $j { my Int $q = 0; loop (my int $i = $len; $i > 0; $i = $i - 1) { my int $x = 10 * @a[$i - 1] + $q * $i; @a[$i - 1] = $x % ( 2 * $i - 1); $q = $x div (2 * $i - 1); } @a[0] = $q % 10; $q div= 10; if $q == 9 { ++$nines; } elsif $q == 10 { take $predigit + 1, 0 xx $nines; $nines = 0; $predigit = 0; } else { take $predigit; $predigit = $q; take 9 xx $nines; $nines = 0; } } } multi MAIN($n = 100) { say pi-spigot($n.Int); } multi MAIN('test') { use Test; plan 1; is pi-spigot(100), '03141592653589793238462643383279502884197169399375105820974944 +59230781640628620899862803482534211706', 'it works'; }

    Since Rakudo is still pretty slow for this kind of stuff, I've traded a bit of readabilty for speed by using a native int in the inner loop, which means that Rakudo can inline most operators, but means I have to write $i = $i - 1 instead of $i-- (because native ints are value types, and you cannot (yet?) pass them as writable values to routines, so the -- operator cannot work on them).

      But wait, if you print the value all at once, then why do you need the complicated incremental radix conversion routine? Why not just generate the whole binary representation at once, then convert it to decimal all at once, using the built in bigfloat or bigint operations? (Sorry, I won't write that in code now.)

      I found a nice Perl 5 implementation of a spigot algorithm to generate pi at With a bit of tweaking, I was able to convert this code into a true spigot: a function which returns a further sequence of digits on each successive call:

      Some advantages of this approach:

      • The output is in decimal.
      • The output can be displayed progressively, so that, for a large number of digits, the user can ‘see’ that progress is being made.
      • With use of the GMP library, performance is surprisingly fast (10,000 digits in under a minute).
      • Flexibility: the caller is free to display or otherwise use the data returned as required.

      Let me emphasize, the code is not mine, I have only massaged it into a (hopefully) more useful form. Perhaps it will prove interesting or helpful to others exploring this topic.

      Athanasius <°(((><contra mundum

      This article has a couple of variations of the algorithm with example code in Haskell.

      As the author points out, representation change algorithms like these are showcase examples for lazy evaluation. This in turn makes Perl 6 a very good language to implement them, besides Haskell.

      With the techniques explained in the article you can replace the large fixed-size array @a for the state by two FatRats or four Integers and actually return an infinite series of decimal digits (limited only by memory resources as the state numbers grow).

      These algorithms are not particularly fast, but not half bad either. Perfect, if you want to incrementally increase precision as you need it. In fact, they could be used in a framework for arbitrarily precise real arithmetic. I would love Perl 6 to support high precision math beyond rational arithmetic.

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: perlmeditation [id://992580]
Front-paged by Corion
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others imbibing at the Monastery: (4)
As of 2017-09-23 16:38 GMT
Find Nodes?
    Voting Booth?
    During the recent solar eclipse, I:

    Results (272 votes). Check out past polls.