/*
 * Computation of the n'th decimal digit of pi with very little memory.
 * Written by Fabrice Bellard on February 26, 1997.
 * 
 * We use a slightly modified version of the method described by Simon
 * Plouffe in "On the Computation of the n'th decimal digit of various
 * transcendental numbers" (November 1996). We have modified the algorithm
 * to get a running time of O(n^2) instead of O(n^3log(n)^3).
 *
 * This program uses a variation of the formula found by Gosper in 1974 :
 * 
 * pi = sum( (25*n-3)/(binomial(3*n,n)*2^(n-1)), n=0..infinity);
 *
 * This program uses mostly integer arithmetic. It may be slow on some
 * hardwares where integer multiplications and divisons must be done by
 * software. We have supposed that 'int' has a size of at least 32 bits. If
 * your compiler supports 'long long' integers of 64 bits, you may use the
 * integer version of 'mul_mod' (see HAS_LONG_LONG).  
 */

/* this is the HACKED around version ready for calculators
 * -- hugh@voidware.com 
 */

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

int mul_mod(int a, int b, int m)
{
    /* find a*b mod m
     * 
     * this version assumes there will be no 10 digit overflow.
     */

#if 0
    /* check for overflow incase we are wrong!
     */
    double t = (double)a * (double)b;
    if (t > 1e10) {
        printf("a*b = %f\n", t);
    }
#endif
    return fmod( (double) a * (double) b, m);
}

int mul_mod2_slow(int a, int b, int m)
{
    /* calculate a*b mod 2^m
     *
     * this is the version of logic used by the calculator
     * execept that 2^m and 2^(floor((m+1)/2)) are in already in
     * registers.
     */
    int a1, a0, b1, b0;
    int B;
    int t;

    B = 1<<(m/2);
    if (m & 1) B <<= 1;
    
    a1 = a/B;
    a0 = a - a1*B;

    b1 = b/B;
    b0 = b - b1*B;

    t = (a1*b0 + a0*b1) % B;
    t = (t*B + a0*b0) % (1<<m);
    
    return t;
}

int mul_mod2(int a, int b, int m)
{
    /* calculate a*b mod 2^m
     *
     * nice and easy for binary computer. slow version above
     * for calculators.
     */
    int v =  (a*b) & ((1<<m)-1);
    return v;
}

int inv_mod(int x,int y) {
  int q,u,v,a,c,t;

  /* return the inverse of x mod y */
  
  u=x;
  v=y;
  c=1;
  a=0;
  do {
    q=v/u;
    
    t=c;
    c=a-q*c;
    a=t;
    
    t=u;
    u=v-q*u;
    v=t;
  } while (u!=0);
  a=a%y;
  if (a<0) a=y+a;
  return a;
}

int pow_mod(int a,int b,int m)
{
    /* calculate a^b mod m */

    int r,aa;
   
    r=1;
    aa=a;
    for (;;) {
        if (b&1) r=mul_mod(r,aa,m);
        b=b>>1;
        if (b == 0) break;
        aa=mul_mod(aa,aa,m);
    }
    return r;
}

int pow_mod2(int a,int b,int m)
{
    /* calculate a^b mod 2^m 
     *
     * use the special 2 power versions.
     */
    int r,aa;
    int c;
   
    r=1;
    aa=a;
    for (;;) {
        c = b/2;
        if (b != c*2) r=mul_mod2(r,aa,m);
        b=c;
        if (b == 0) break;
        aa=mul_mod2(aa,aa,m);
    }
    return r;
}

int is_prime(int n)
{
    /* return true if odd n is prime */
    
   int r,i;

   r=(int)(sqrt(n));
   for(i=3;i<=r;i+=2) if ((n % i) == 0) return 0;
   return 1;
}

int next_prime(int n)
{
    /* return the next odd prime */
    do {
        n += 2;
    } while (!is_prime(n));
    return n;
}

int DIVN(int* t,int a,int* kq)
{	
    /* factor deflation subroutine */
    int c = 0;
    if (*kq >= a) {			
        do { *kq-=a; } while (*kq>=a);	
        if (*kq == 0) {			
            do {				
                (*t)/=a;
                ++c;
            } while ((*t % a) == 0);		
        }					
    }		
    return c;
}

#define DIGLIMIT 8.846944397064768934339805266660e-1

void doPi(int n)
{
    int av,a,vmax,N,num,den,k,kq1,kq2,kq3,kq4,t,v,s,i,t1;
    double sum;

    N=(int)((n+20)*DIGLIMIT);
    sum = 0;

    /* here we start the case a=2.
     * this is unrolled here and simplified at the start.
     *
     * furthermore, use the special 2 power subroutine versions.
     * a crucial speed up for calculators.
     */
    vmax=(int)(log(3*N)/log(2));
    vmax=vmax+(N-n);
    if (vmax<=0) {
        goto l1;
    }

    /* av = 2^vmax */
    av=2;
    for(i=1;i<vmax;i++) av=av*2;
    
    s=0;
    den=1;
    num=1;
    v=-n; 

    for(k=1;k<=N;k++) {

      t=k;
      while ((t & 1) == 0) { --v; t /=2; }
      t *= 2*k-1;
      num=mul_mod2(num,t,vmax);

      t=3*(3*k-1);
      while ((t & 1) == 0) { ++v; t /=2; }
      den=mul_mod2(den,t,vmax);

      t=(3*k-2);
      while ((t & 1) == 0) { ++v; t /=2; }
      den=mul_mod2(den,t,vmax);
      
      if (v > 0) {
          t=inv_mod(den,av);
          t=mul_mod2(t,num,vmax);
          
          t1 = 1 << (vmax-v);
          t1 *= (25*k-3);
          t = mul_mod2(t, t1, vmax);

          s+=t;
          if (s>=av) s-=av;
      }
    }

    t=pow_mod2(5,n-1,vmax);
    s=mul_mod2(s,t,vmax);
    sum=fmod((double) s/ (double) av,1.0);

l1:

    /* now walk the odd primes. use the subroutine versions
     * that expect no overflow.
     */

    for(a=3;a<=(3*N);a=next_prime(a)) {
        vmax=(int)(log(3*N)/log(a));
        av=a;
        for(i=1;i<vmax;i++) av=av*a;

        s=0;
        den=1;

        num=pow_mod(2,n,av);
        v=0;

        kq1=0;
        kq2=-1;
        kq3=-3;
        kq4=-2;

        for(k=1;k<=N;k++) {

            t=2*k;
            kq1 += 2;
            v -= DIVN(&t,a,&kq1);
            num=mul_mod(num,t,av);
      
            t=2*k-1;
            kq2 += 2;
            v -= DIVN(&t,a,&kq2);
            num=mul_mod(num,t,av);

            t=3*(3*k-1);
            kq3 += 9;
            v += DIVN(&t,a,&kq3);
            den=mul_mod(den,t,av);

            t=(3*k-2);
            kq4 += 3;
            v += DIVN(&t,a,&kq4);
            t=t*2;
      
            den=mul_mod(den,t,av);
      
            if (v > 0) {
                t1 = pow_mod(a, vmax-v, av);
                t1 = mul_mod(t1, 25*k-3, av);

                t=inv_mod(den,av);
                t=mul_mod(t,num,av);
                t=mul_mod(t,t1,av);

                s+=t;
                if (s>=av) s-=av;
            }
        }

        t=pow_mod(5,n-1,av);
        s=mul_mod(s,t,av);
        sum=fmod(sum+(double) s/ (double) av,1.0);
    }

    /* using doubles we get 10 (11 maybe) correct digits.
     * just change 1e6 to 1e10 and skip 10 below. 
     * this is because we need at least 4 to protect against 
     * rounding error accumulation.
     *
     * i think 7 might work on a 10 digit calculator, but after
     * seeing 12 not work in binary, im playing safe with 6.
     */
    printf("(%d) %06d\n", n, (int)(sum*1e6));
}
   
int main(int argc,char *argv[])
{
    int i;
    for (i = 1; i <= 1000; i += 6) {
        doPi(i);
    }
    return 0;
}