/*  This program employs the recently discovered digit extraction scheme

    to produce hex digits of pi.  This code is valid up to ic = 2^24 on 

    systems with IEEE arithmetic. */



/*  David H. Bailey     2000-03-28 */



#include <stdio.h>

#include <math.h>



main()

{

  double pid, s1, s2, s3, s4;

  double series (int m, int n);

  void ihex (double x, int m, char c[]);

  int ic = 1000000;

#define NHX 16

  char chx[NHX];



/*  ic is the hex digit start position. */



  s1 = series (1, ic - 1);

  s2 = series (4, ic - 1);

  s3 = series (5, ic - 1);

  s4 = series (6, ic - 1);

  pid = 4. * s1 - 2. * s2 - s3 - s4;

  pid = pid - (int) pid + 1.;

  ihex (pid, NHX, chx);

  printf (" start position = %i\n hex digits =  %10.10s\n", ic, chx);

}



void ihex (double x, int nhx, char chx[])



/*  This returns, in chx, the first nhx hex digits of the fraction of x. */



{

  int i;

  double y;

  char hx[] = "0123456789ABCDEF";



  y = fabs (x);



  for (i = 0; i < nhx; i++){

    y = 16. * (y - floor (y));

    chx = hx[(int) y];

  }

}



double series (int m, int ic)



/*  This routine evaluates the series  sum_k 16^(ic-k)/(8*k+m) 

    using the modular exponentiation technique. */



{

  int k;

  double ak, eps, p, s, t;

  double expm (double x, double y);

#define eps 1e-17



  s = 0.;



/*  Sum the series up to ic. */



  for (k = 0; k < ic; k++){

    ak = 8 * k + m;

    p = ic - k;

    t = expm (p, ak);

    s = s + t / ak;

    s = s - (int) s;

  }



/*  Compute a few terms where k >= ic. */



  for (k = ic; k <= ic + 100; k++){

    ak = 8 * k + m;

    t = pow (16., (double) (ic - k)) / ak;

    if (t < eps) break;

    s = s + t;

    s = s - (int) s;

  }

  return s;

}



double expm (double p, double ak)



/*  expm = 16^p mod ak.  This routine uses the left-to-right binary 

    exponentiation scheme.  It is valid for  ak <= 2^24. */



{

  int i, j;

  double p1, pt, r;

#define ntp 25

  static double tp[ntp];

  static int tp1 = 0;



/*  If this is the first call to expm, fill the power of two table tp. */



  if (tp1 == 0) {

    tp1 = 1;

    tp[0] = 1.;



    for (i = 1; i < ntp; i++) tp = 2. * tp[i-1];

  }



  if (ak == 1.) return 0.;



/*  Find the greatest power of two less than or equal to p. */



  for (i = 0; i < ntp; i++) if (tp > p) break;



  pt = tp[i-1];

  p1 = p;

  r = 1.;



/*  Perform binary exponentiation algorithm modulo ak. */



  for (j = 1; j <= i; j++){

    if (p1 >= pt){

      r = 16. * r;

      r = r - (int) (r / ak) * ak;

      p1 = p1 - pt;

    }

    pt = 0.5 * pt;

    if (pt >= 1.){

      r = r * r;

      r = r - (int) (r / ak) * ak;

    }

  }



  return r;

}