10

I'm having a hard time understanding the Spigot algorithm for π (pi) found here at the bottom of the page.

I'm getting lost at the bottom of part 2 "Put A into regular form", I'm not exactly sure how to implement this in C (or any language really)

Shafik Yaghmour
  • 154,301
  • 39
  • 440
  • 740
Devan Buggay
  • 2,717
  • 4
  • 26
  • 35
  • 1
    You may be interested in this [question](http://stackoverflow.com/questions/20287513/can-anyone-make-heads-or-tales-of-this-spigot-algorithm-code-pitiny-c). – Shafik Yaghmour Dec 02 '13 at 14:07

3 Answers3

9
  #include <math.h>
  #include <stdio.h>
  #define N 100

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

  for(int i = 0; i < len; ++i) {
    A[i] = 2;
  }

  int nines    = 0;
  int predigit = 0;

  for(int j = 1; j < N + 1; ++j) {        
    int q = 0;

    for(int i = len; i > 0; --i) {
      int 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 (9 == q) {
      ++nines;
    }
    else if (10 == q) {
      printf("%d", predigit + 1);

      for (int k = 0; k < nines; ++k) {
        printf("%d", 0);
      }
      predigit, nines = 0;
    }
    else {
      printf("%d", predigit);
      predigit = q;

      if (0 != nines) {    
        for (int k = 0; k < nines; ++k) {
          printf("%d", 9);
        }

        nines = 0;
      }
    }
  }
  printf("%d", predigit);
Pedro Silva
  • 4,672
  • 1
  • 19
  • 31
  • Why are initializing in for loops? and why are your if statements reversed? – Devan Buggay Nov 03 '10 at 18:56
  • 1
    What do you mean 'initializing in for loops'? The `if` statements are reversed because that helps me catch inadvertent assignments (eg. `if(q = 0)` instead of `if(q == 0)`). The compiler will catch my version (`if(0 = q)`) but not `if(q = 0)`. – Pedro Silva Nov 04 '10 at 00:03
  • 2
    Oh, you mean doing `for (int i = 0 ...)`? Clear up clutter, I guess. No temporary initializations outside the actual scope of the loops. In theory, this should ensure that the index variables are lexically scoped inside the loop, but I think I read somewhere that they actually remain bound outside the loop. – Pedro Silva Nov 04 '10 at 00:08
1
// Spigot program for pi to NDIGITS decimals.
// 4 digits per loop.
// Thanks to Dik T. Winter and Achim Flammenkamp who originally made a compressed version of this. 
 
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
 
#define NDIGITS 15536          //max digits to compute
#define LEN (NDIGITS/4+1)*14   //nec. array length
 
long a[LEN];                   //array of 4-digit-decimals
long b;                        //nominator prev. base
long c = LEN;                  //index
long d;                        //accumulator and carry
long e = 0;                    //save previous 4 digits
long f = 10000;                //new base, 4 decimal digits
long g;                        //denom previous base
long h = 0;                    //init switch
 
int main(void) {
        for(; (b=c-=14) > 0 ;){    //outer loop: 4 digits per loop
                for(; --b > 0 ;){      //inner loop: radix conv
                        d *= b;            //acc *= nom prev base
                        if( h == 0 )
                                d += 2000*f;   //first outer loop
                        else
                                d += a[b]*f;   //non-first outer loop
                        g=b+b-1;           //denom prev base
                        a[b] = d % g;
                        d /= g;            //save carry
                }
                h = printf("%.4d",e+d/f);//print prev 4 digits
                                         // %.4d to add leading zeros
                d = e = d % f;         //save currently 4 digits
                                       //assure a small enough d
        }
        getchar();
        return 0;
}
dquadros
  • 43
  • 6
Progo
  • 3,452
  • 5
  • 27
  • 44
1

I am seeing a difference in the o/p digits of the above spigot pi program when compared with http://www.numberworld.org/misc_runs/pi-10t/details.html

Correct Value of 50 digits of Pi : http://www.numberworld.org/misc_runs/pi-10t/details.html

3.

1415926535 8979323846 2643383279 5028841971 6939937510

Above spigot pi :

3.

1415926535 8979323846 2643383279 5**(0)**28841971 6939937510

                                     ^^^ zero missing

Changed the 4 digits per loop to 8 digits by modifying long f = 100000000;

produced the correct result.

Hari
  • 11
  • 1
  • 1
    Seems like this is a known issue in Spigot algorithm which need corrective step - refer to haenel's implementation that fixes this issue and produce correct result. See the correct C implementation and description at the beging of http://www.jjj.de/hfloat/spigot.haenel.txt – Hari Feb 03 '15 at 02:57
  • Haenel's implementation of Pi Spigot for 32372 digits. Verified for correct digit sequence with Chudnovsky's algorithm o/p. 2. Here is a version which seems to be correct in this respect and which is hopefully without new bugs: (It is also faster and shorter.) Version "Fast": (approx. 25% faster) `/* Calculation of pi to 32372 decimal digits */ /* Size of program: 152 characters */ /* After Dik T. Winter, CWI Amsterdam */ unsigned a=1e4,b,c=113316,d,e,f[113316],g,h,i; main(){for(;b=c,c-=14;i=printf("%04d",e+d/a),e=d%a) while(g=--b*2)d=h*b+a*(i?f[b]:a/5),h=d/--g,f[b]=d-g*h;}` – Hari Feb 03 '15 at 03:00