Computing pi to multiple precision

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'; } [download] 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). Perl 6 - the future is here, just unevenly distributed	[reply] [d/l] [select]
Re^2: Computing pi to multiple precision by ambrus (Abbot) on Sep 09, 2012 at 16:30 UTC
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.)	[reply]
Re^3: Computing pi to multiple precision by moritz (Cardinal) on Sep 09, 2012 at 19:46 UTC
The point is exactly that you don't need a full bigint or bigfloat implementation to compute a few hundred digits of pi. Perl 6 - the future is here, just unevenly distributed	[reply]
Re^2: Computing pi to multiple precision by Athanasius (Archbishop) on Sep 10, 2012 at 16:36 UTC
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: Read more... (3 kB) 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	[reply] [d/l]
Re^2: Computing pi to multiple precision by ambrus (Abbot) on Sep 09, 2012 at 16:00 UTC
As for algorithms like that, Hexadecimal representation of pi also implements one, only that doesn't try to convert to decimal.	[reply]
Re^2: Computing pi to multiple precision by martin (Friar) on Sep 12, 2012 at 13:21 UTC
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.	[reply]
Re: Computing pi to multiple precision by marioroy (Prior) on Jul 28, 2022 at 22:36 UTC
Back in 2018, I came across the article Parallel GMP-Chudnovsky using OpenMP with factorization and thought how cool it would be to do something similarly involving Perl, MCE, and Inline::C. The examples demonstrate nested parallelization using OpenMP and pthreads. https://github.com/marioroy/Chudnovsky-Pi The code supports GMP and MPIR. Moreover, I created a complementary extra folder, my attempt at making GMP/MPIR's mpn_get_str parallel via divide-and-conquer. I used prime_test.cpp by Anthony Hay for testing. mpn_get_str $ cd Chudnovsky-Pi/src/extra $ g++ -O3 -DTEST6 -fopenmp prime_test.cpp -l gmp -o prime6_test -Wno-a +ttributes $ time ./prime6_test Calculating 2^n-1 for n=1257787... Calculating 2^n-1 for n=1398269... Calculating 2^n-1 for n=2976221... Calculating 2^n-1 for n=3021377... Calculating 2^n-1 for n=6972593... Calculating 2^n-1 for n=13466917... Calculating 2^n-1 for n=20996011... Calculating 2^n-1 for n=24036583... Calculating 2^n-1 for n=25964951... Calculating 2^n-1 for n=30402457... Calculating 2^n-1 for n=32582657... Calculating 2^n-1 for n=37156667... Calculating 2^n-1 for n=42643801... Calculating 2^n-1 for n=43112609... Calculating 2^n-1 for n=57885161... total failures 0 real 0m11.722s user 0m11.626s sys 0m0.087s [download] parallel mpn_get_str Run top in another window and have it refresh every 0.1 seconds. You will see the parallel demonstration consume many threads not exceeding 8 threads max. $ g++ -O3 -DPARALLEL -DTEST6 -fopenmp prime_test.cpp -l gmp -o prime6_ +test -Wno-attributes $ time ./prime6_test Calculating 2^n-1 for n=1257787... Calculating 2^n-1 for n=1398269... Calculating 2^n-1 for n=2976221... Calculating 2^n-1 for n=3021377... Calculating 2^n-1 for n=6972593... Calculating 2^n-1 for n=13466917... Calculating 2^n-1 for n=20996011... Calculating 2^n-1 for n=24036583... Calculating 2^n-1 for n=25964951... Calculating 2^n-1 for n=30402457... Calculating 2^n-1 for n=32582657... Calculating 2^n-1 for n=37156667... Calculating 2^n-1 for n=42643801... Calculating 2^n-1 for n=43112609... Calculating 2^n-1 for n=57885161... total failures 0 real 0m4.368s user 0m11.966s sys 0m0.134s [download] Chudnovsky Pi demonstration `$ cd Chudnovsky-Pi/src $ make $ make pi-gmp # requires GMP $ make pi-mpir # requires MPIR $ cd ../bin` [download] pi-hobo.pl Once the computation is completed (upon seeing end date...), mpf_get_str is called. It will consume 1 thread for a while, then 2, 4, max 8 threads. $ perl pi-hobo.pl 100000000 1 auto \| md5sum # start date = Thu Jul 28 17:50:15 2022 # terms = 7051366, depth = 24, threads = 64, logical cores = 64 sieve cputime = 0.27s wallclock = 0.27s factor = 1. +0 bs cputime = 54.74s wallclock = 1.00s factor = 54. +5 sum cputime = 29.47s wallclock = 3.89s factor = 7. +6 div/sqrt cputime = 8.21s wallclock = 5.02s factor = 1. +6 mul cputime = 2.09s wallclock = 2.09s factor = 1. +0 total cputime = 94.78s wallclock = 12.26s factor = 7. +7 1.58m 0.20m # P size = 158218296 digits (1.582183) # Q size = 158218289 digits (1.582183) # end date = Thu Jul 28 17:50:27 2022 969bfe295b67da45b68086eb05a8b031 - [download] pi-gmp $ ./pi-gmp 100000000 1 auto \| md5sum # start date = Thu Jul 28 17:54:00 2022 # terms = 7051366, depth = 24, threads = 64, logical cores = 64 sieve cputime = 0.28s wallclock = 0.28s factor = 1. +0 bs cputime = 55.93s wallclock = 0.95s factor = 58. +7 sum cputime = 28.54s wallclock = 3.87s factor = 7. +4 div/sqrt cputime = 8.18s wallclock = 5.03s factor = 1. +6 mul cputime = 2.13s wallclock = 2.13s factor = 1. +0 total cputime = 95.07s wallclock = 12.27s factor = 7. +8 1.58m 0.20m # P size = 158218296 digits (1.582183) # Q size = 158218289 digits (1.582183) # end date = Thu Jul 28 17:54:12 2022 969bfe295b67da45b68086eb05a8b031 - [download]	[reply] [d/l] [select]
Re: Computing pi to multiple precision by ambrus (Abbot) on Jul 24, 2022 at 16:11 UTC
John Savard's page on �Squaring the circle� talks about formulas similar to the arctangent one that I'm using here.	[reply]
Re: Computing pi to multiple precision by Discipulus (Canon) on Oct 23, 2017 at 20:12 UTC
thanks ambrus for this very interesting and useful post, even if I'm a bit late. With the help of haukex and other kind monks, I've ported the C code used above to Perl as you can read at Re^2: porting C code to Perl -- solved. Probably the used algorithm has room for improvement. L* There are no rules, there are no thumbs.. Reinvent the wheel, then learn The Wheel; may be one day you reinvent one of THE WHEELS.	[reply] [d/l]