The Digits of Pi
February 20, 2009
A spigot algorithm spits out the digits of an irrational number, one at a time, like a leaky faucet. Jeremy Gibbons gives a spigot algorithm for pi at http://web.comlab.ox.ac.uk/people/Jeremy.Gibbons/publications/spigot.pdf:
(define (pi-spigot z)
(let loop ((z z) (ds '()) (q 1) (r 0) (t 1) (k 1) (n 3) (m 3))
(cond ((zero? z) (reverse ds))
((< (+ q q q q r (- t)) (* n t))
(loop (- z 1) (cons n ds) (* 10 q) (* 10 (- r (* n t)))
t k (- (quotient (* 10 (+ q q q r)) t) (* 10 n)) m))
(else (loop z ds (* q k) (* (+ q q r) m) (* t m) (+ k 1)
(quotient (+ (* q (+ k k k k k k k 2)) (* r m))
(* t m)) (+ m 2))))))
The answer, which you can see at http://programmingpraxis.codepad.org/qfgnpRbn, is 9:
> (pi-spigot 1001)
(3 1 4 1 5 9 ... 2 0 1 9 8 9)
Just in case you’re curious, the one-millionth digit of pi is 1, and the one-billionth digit of pi is 9.
Programming Praxis 
My solution, written in C:
http://codepad.org/WOEQnbTt
This version calculates the requested digit directly, based on the method described in “Computation of the nth decimal digit of pi with low memory”: http://numbers.computation.free.fr/Constants/Algorithms/nthdecimaldigit.pdf (and I must admit I also peeked a little at the implementation at http://numbers.computation.free.fr/Constants/Algorithms/nthdigit.html)
It lacks the conciseness and elegance of the spigot solution, but makes up for it in efficiency. The spigot version works well until it runs out of main memory and starts paging, at which point performance really hits a wall. On my computer that point is somewhere between the 20,000th digit (took 1:48) and the 50,000th (I killed it after 3.5 hours). For comparison, this solution returns the 50,000th digit in ~12 seconds.
Alexander J. Yee & Shigeru Kondo recently computed five trillion digits of pi.
I modified the fixed-point approach for arccot found
here
and used it with Wikipedia’s most efficient Machin-Like formula found
here
My submission
This seems like a really interesting exercise, however I cannot understand how this spigot algorithm works. I read the PDF, but I could not follow the maths in the diagram, is there an easier explanation for someone who doesn’t have a massive background in maths?
Cheers,
-Sam
It isn’t pretty but it works:
#!/usr/bin/perl
$| = 1;
push(@B, 0);
$counter = 3;
for($i = 0;$i<=33223;$i++){
push(@A,$i);
push(@remainder, 2);
push(@B, $counter);
$counter += 2;
}
for($x = 1; $x=1;$i--){
$m10[$i] = $remainder[$i] * 10;
$sum[$i] = $m10[$i] + $carried[$i];
$remainder[$i] = $sum[$i] % $B[$i];
$carried[$i - 1] = (($sum[$i] - $remainder[$i]) / $B[$i]) * $A[$i];
}
$m10[0] = $remainder[0] * 10;
$sum[0] = $m10[0] + $carried[0];
$remainder[0] = $sum[0]%10;
$predigit = int($sum[0]/10);
push(@predigits, $predigit);
if($predigit < 9){
for($i = 0; $i<((scalar @predigits)-1); $i++){
print $predigits[$i];
}
undef @predigits;
push(@predigits, $predigit);
}
if($predigit == 10){
for($i = 0; $i<((scalar @predigits)-1); $i++){
if($predigits[$i] == 9) { print 0; }
else {
$predigits[$i] = $predigits[$i] + 1;
print $predigits[$i];
}
}
undef @predigits;
push(@predigits, 0);
}
}
Python:
def pi_spigot(): (q, r, t, k, n, l) = (1, 0, 1, 1, 3, 3) while True: if 4*q+r-t < n*t: yield n (q, r, t, k, n, l) = (10*q, 10*(r-n*t), t, k, (10*(3*q+r))//t-10*n, l) else: (q, r, t, k, n, l) = (q*k, (2*q+r)*l, t*l, k+1, (q*(7*k+2)+r*l)//(t*l), l+2) def print_pi(n): digits = pi_spigot() print('pi =', str(next(digits)) + '.', end='') while n > 0: print(str(next(digits)), end='') n -= 1 print()Running:
For the curious, all first 1000 digits (it ran off last post there…)
https://github.com/ftt/programming-praxis/blob/master/20090220-the-digits-of-pi/pi.py