The problem with Pi (ii)

To print out precise numbers, we need another approach. One such approach to π was offered in 1990, when  Rabinowitz and Wagon discovered a “spigot” algorithm for π [1]. This algorithm allows for successive digits of π to be calculated using a simple recursive algorithm.

Here is the Pascal program from the original paper:

program pi_spigot;
const n = 1000;
      len = 10*n div 3;
var i, j, k, q, x, nines, predigit : integer;
    a : array[1..len] of longint;
begin
   for j := 1 to len do a[j] := 2;
   nines := 0;
   predigit := 0;
   for j := 1 to n do
   begin
      q := 0;
      for i := len downto 1 do
      begin
         x := 10*a[i] + q*i;
         a[i] := x mod (2*i-1);
         q := x div (2*i-1);
      end;
      a[1] := q mod 10; 
      q := q div 10;
      if q = 9 then nines := nines + 1
      else if q = 10 then
      begin
         write(predigit+1);
         for k := 1 to nines do write(0);
         predigit := 0;
         nines := 0;
      end
      else begin
         write(predigit);
         predigit := q;
         if nines <> 0 then
         begin
            for k := 1 to nines do write(9);
            nines := 0;
         end
      end
   end;
   writeln(predigit);
end.

The reader is left to browse the original paper for algorithm details. When running the code (I used the Free Pascal Compiler), it gives the following output:

03141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923459099207118856806516745194252854340964271825877594867237227818052636476240001253966769588037865014608282022810293731532093125719608182508805545474738567199841858637460939653266515205973230200470327367748014549418069976867613243271975438847522856914469354541235658607836752109748167637744042444357591344102954987267348532155345226074874607261272294756681440117171146164262329464280245194165531036682264959997624747949649719425809626723943536390912536760246664166117527686634968595119533790744923019073345333321686533168484290724097901531077274052175927055788110268222469183375946146833893295768948223454910618813165693638675020776147309271039650546608291414446931319129743860775232869818000001373465107828774083538328579350606183962012466

Only the first 265 digits are correct. This is probably due to overflow errors due to datatypes being too small. Convert the program to C, and it works as it should:


#include <stdio.h>
#include <stdlib.h>

void spigot();

int main(void) {
   spigot();
   return 0;
}

void spigot() {
   int i, j, k, q, x;
   int len, nines=0, predigit=0;
   int N=1000;

   len = (10*N/3)+1;
   int a[len];

   // Initialize A to (2,2,2,2,2,...,2)
   for (i=0; i<len; i=i+1)
      a[i] = 2;

   // Repeat n times
   for (j=1; j<=1000; j=j+1) {
      q = 0;
      for (i=len; i>0; i=i-1) {
         x = 10 * a[i-1] + q*i;
         a[i-1] = x % (2*i-1);
         q = x / (2*i-1);
      }
      a[0] = q % 10;
      q = q / 10;
      if (q == 9)
         nines = nines + 1;
      else if (q == 10) {
         printf("%d", predigit+1);
         for (k=0; k<nines; k=k+1)
            printf("%d",0);
         predigit = 0;
         nines = 0;
      }
      else {
         printf("%d", predigit);
         predigit = q;
         if (nines != 0) {
            for (k=0; k<nines; k=k+1)
               printf("%d",9);
            nines = 0;
         }
     }
  }
  printf("%d\n", predigit);
}

[1] ] Rabinowitz, S.D., Wagon, S., ” A spigot algorithm for pi”. American Mathematical Monthly, 102, pp.195–203 (1995)

 

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.