Syntactic Confectionery Delight PerlMonks

### Re^2: Computing pi to multiple precision

by Athanasius (Chancellor)
 on Sep 10, 2012 at 16:36 UTC ( #992798=note: print w/replies, xml ) Need Help??

in reply to Re: Computing pi to multiple precision
in thread Computing pi to multiple precision

I found a nice Perl 5 implementation of a spigot algorithm to generate pi at http://rosettacode.org/wiki/Pi#Perl. 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:

```#! perl
use strict;
use warnings;
use Math::BigInt try => 'GMP';
use bigint;

{
local \$| = 1;

my  \$digits   = int(\$ARGV[0] // 0) || 60;
my (\$len, \$d) = rosetta();

printf "pi = %s.%s", substr(\$d, 0, 1), substr(\$d, 1);

for (my \$count = \$len; \$count < \$digits; \$count += \$len)
{
(\$len, \$d) = rosetta();
print \$d;
}

print ";\n";

reset_rosetta();
(\$len, \$d) = rosetta();
printf "pi = %s.%s", substr(\$d, 0, 1), substr(\$d, 1);
}

BEGIN
{
my (\$ds, \$ns, \$n5, \$d5, \$n2, \$d2, \$p2, \$pow, \$x);

reset_rosetta();

sub reset_rosetta
{
\$ds  =          1;
\$ns  =          0;
\$n5  =      1_184;
\$d5  =        375;
\$n2  =    685_448;
\$d2  = 40_955_757;
\$p2  =          5;
\$pow =          1;
\$x   =          5;
}

sub rosetta
{
my \$out;

for (my \$ppow = 1; \$ppow == 1; \$x += 4)
{
\$ns  = (\$ns * \$d5) + (\$n5 * \$pow * \$ds);
\$ds *=  \$d5;
\$n5  =  16 * (25 * (\$x + 2) - \$x);
next_term(\$d5, 5, \$x);

while (\$d5 > \$d2)
{
\$ns  = (\$ns * \$d2) - (\$n2 * \$pow * \$ds);
\$ds *=  \$d2;
\$n2  =  4 * (57_121 * (\$p2 + 2) - \$p2);
next_term(\$d2, 239, \$p2);
\$p2 +=  4;
}

my \$product1 = \$n5 * 625;
my \$product2 = \$n2 * \$n2 * 3_262_808_641;
\$ppow     = 1;

while (\$pow * \$product1 < \$d5 &&
\$pow * \$product2 < \$d2)
{
\$pow  *= 10;
\$ppow *= 10;
}

if (\$ppow > 1)
{
\$ns *=  \$ppow;
\$out =  \$ns / \$ds;
\$ns %=  \$ds;
\$out = ('0' x (length(\$ppow) - length(\$out) - 1)) . \$o
+ut;
}

if (\$p2 % 20 == 1)
{
my \$g = Math::BigInt::bgcd(\$ds, \$ns);
\$ds  /= \$g;
\$ns  /= \$g;
}
}

return (length \$out, \$out);
}
}

sub next_term
{
my (\$coef, \$p) = @_[1, 2];

\$_[0] /= (\$p - 4) * (\$p - 2);
\$_[0] *=  \$p * (\$p + 2) * \$coef ** 4;
}
[download]```

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

Log In?
 Username: Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://992798]
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others avoiding work at the Monastery: (5)
As of 2017-06-28 06:13 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
How many monitors do you use while coding?

Results (625 votes). Check out past polls.