PI: the movie

Calculating digits of pi: The Ruby release has a remarkably short and simple Ruby program that calculates arbitrarily many decimal digits of pi:

Ruby program to calculate pi
#!/usr/local/bin/ruby
k, a, b, a1, b1 = 2, 4, 1, 12, 4
loop do
  p, q, k = k*k, 2*k+1, k+1
  a, b, a1, b1 = a1, b1, p*a+q*a1, p*b+q*b1
  d, d1 = a/b, a1/b1
  while d == d1
    print d
    $stdout.flush
    a, a1 = 10*(a%b), 10*(a1%b1)
    d, d1 = a/b, a1/b1
  end
end
Output of a run
% pi.rb
3141592653589793238462643383279502884 ...

Here is an expanded verion that arranges the digits by 10s and in rows of 100. Given an optional command-line argument, it calculates pi to any base B <= 36.

Ruby program to calculate pi to any base
#!/usr/bin/env ruby
k, a, b, a1, b1 = 2, 4, 1, 12, 4
h = %w{0 1 2 3 4 5 6 7 8 9 a b c d e f g h
       i j k l m n o p q r s t u v w x y z}
if ARGV[0] == nil
  B = 10
else 
  B = ARGV[0].to_i
end 
dig = -1
row = -1
loop do
  # Next approximation
  p, q, k = k*k, 2*k+1, k+1
  a, b, a1, b1 = a1, b1, p*a+q*a1, p*b+q*b1
  # Print common digits
  d = a / b
  d1 = a1 / b1
  while d == d1
    if dig == -1 then
       printf "     "
       print h[d]
       dig = dig+1
    else
       print h[d]
       dig = dig+1
    end
    if dig%10 == 0 then printf " " end
    if dig%100 == 0 then 
       row = row+1
       printf "\n%4i ", row
    end
    $stdout.flush
    a, a1 = B*(a%b), B*(a1%b1)
    d, d1 = a/b, a1/b1
  end
end
Output of runs
% piA.rb
     3
   0 1415926535 8979323846 2643383279 ...
   1 8214808651 3282306647 0938446095 ...
% piA.rb 16
     3 
   0 243f6a8885 a308d31319 8a2e037073 ...
   1 29b7c97c50 dd3f84d5b5 b547091792 ...
% piA.rb 36
     3
   0 53i5ab8p5f sa5jhk72i8 asc47wwzla ...
   1 067iy9wnzd z9n4jltedt iw2b65acrp ...

Here are outputs of the above program: pi to 10000 decimal digitspi to 10000 hex digitspi to 10000 base 36 digits. I checked that the answers are correct. (Well, spot checked the decimal version, and checked at the end. I also checked the first 55 hex digits, and checked that 9 hex digits starting at the 5000th and at the 10000th position are correct.)

The May 6, 1993 episode of The Simpsons has the character Apu boast "I can recite pi to 40,000 places. The last digit is one." See 40000th digit in red. (This fact taken from: The Life of Pi, by J.M. Borwein.)

This program even works for bases 2 and 3, although it gives an initial digit in both cases of 3, the subsequent base 2 or base 3 digits are correct. (I checked base 2 against the hex digits and the initial few digits of base 3.)

The most interesting decimal run in pi starts in position 762 (row 7, column 7), where 9999998 occurs. Another interesting run: pi base 2 is: 11.00100100001111110110..., so that 11.0010010001 = 3 + 1/8 + 1/64 + 1/1024 = 3.1416015625 is just 0.00000890891 bigger than pi -- another good approximation, though not as good as 355/113 = 3.14159292035398: just 0.000000266764 bigger than pi.

Actually, this isn't a very efficient way to calculate pi, but it was able to give 5000 digits in 15 seconds and 10000 digits in under 60 seconds on a 800 Mhz Pentium.

This program seems to be based on the continued fraction below (taken from Wolfram Functions Pages):

PI as continued fraction

If one solves for pi, the result is buried in the Ruby program below, which gives 16 digits of pi as a double. Of course the calculations above used far more terms of the continued fraction.

Ruby: pi as double
#!/usr/local/bin/ruby
     pi =         4/
       ( 1 + ( 1* 1)/
       ( 3 + ( 2* 2)/
       ( 5 + ( 3* 3)/
       ( 7 + ( 4* 4)/
       ( 9 + ( 5* 5)/
       (11 + ( 6* 6)/
       (13 + ( 7* 7)/
       (15 + ( 8* 8)/
       (17 + ( 9* 9)/
       (19 + (10*10)/
       (21 + (11*11)/
       (23 + (12*12)/
       (25 + (13*13)/
       (27 + (14*14)/
       (29 + (15*15)/
       (31 + (16*16)/
       (33 + (17*17)/
       (35 + (18*18)/
       (37 + (19*19)/
       (39 + (20*20)/
       (41 + (21*21)/
       (43 + (22*22.0)
    )) )))))  )))) ))))) ))))))
print pi
% pi_constant.rb
3.14159265358979

Here a Java version of the second Ruby program above (in a slightly simplified form), using BigInteger. It produces exactly the same output as the Ruby program, only slower.

Java program to calculate pi
import java.math.*; // for BigInteger
public class PiBig2 {

   public static void main(String[] args) {
      BigInteger TWO = new BigInteger("2");
      BigInteger TEN = new BigInteger("10");
      BigInteger k = new BigInteger("2");
      BigInteger a = new BigInteger("4");
      BigInteger b = new BigInteger("1");
      BigInteger a1 = new BigInteger("12");
      BigInteger b1 = new BigInteger("4");
      int dig = -1, row = -1;
      for (;;) {
         // Next approximation
         BigInteger p = k.multiply(k);
         BigInteger q = (TWO.multiply(k)).add(BigInteger.ONE);
         k = k.add(BigInteger.ONE);
         BigInteger tempa1 = a1;
         BigInteger tempb1 = b1;
         a1 = (p.multiply(a)).add(q.multiply(a1));
         b1 = (p.multiply(b)).add(q.multiply(b1));
         a = tempa1;
         b = tempb1;
         // Print common digits
         BigInteger d = a.divide(b);
         BigInteger d1 = a1.divide(b1);
         while (d.equals(d1)) {
            if (dig == -1) {
               System.out.print("     ");
               System.out.print(d);
               dig++;
            }
            else {
               System.out.print(d);
               dig++;
            }
            if (dig%10 == 0) System.out.print(" ");
            if (dig%100 == 0) {
               row++;
               System.out.println();
               if (row < 10) System.out.print("   " + row);
               else if (row < 100) System.out.print("  " + row);
               else if (row < 1000) System.out.print(" " + row);
               else System.out.print(row);
               System.out.print(" ");
            }
            a =  TEN.multiply(a.mod(b));
            a1 = TEN.multiply(a1.mod(b1));
            d = a.divide(b);
            d1 = a1.divide(b1);
         }
      }
   }
}

The formula above is a special case of a much more general formula:

PI as continued fraction

Plugging in z = 1 gives the result above. (See Wolfram Functions Pages):


Revision date: 2004-03-16. (Use ISO 8601, an International Standard.)