Perl-Sensitive Sunglasses 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; }

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

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 romping around the Monastery: (3)
As of 2018-05-21 23:34 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
World peace can best be achieved by:

Results (163 votes). Check out past polls.

Notices?