/*  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     960429 */

// extended to show a pseudo-phase space of pi.
// Grant K. Sherman     991006

#include <stdio.h>
#include <iostream.h>
#include	<iomanip.h>
#include <math.h>
#include <fstream.h>

#define INTERVAL 16
int interval[INTERVAL][INTERVAL];

int nold = 3;
int nnew;

char hx[] = "0123456789ABCDEF";

void initialise_interval();
void display_interval();

main()
{
  double pid, s1, s2, s3, s4;
  double series (int m, int n);
  void ihex (double x, int m, char c[]);
  int iint (double x);
  void ialpha (double x, int m, char c[], char b[]);
  int ic;
  int count;
#define NHX 16
#define INTERVAL 16
  char chx[NHX];
  char calpha[NHX];
  char cbeta[NHX];
  char hx[] = "0123456789ABCDEF";

  ofstream output_file("pi500.txt");


  initialise_interval();


/*  ic is the hex digit position -- output begins at position ic + 1. */

for (ic=0; ic < 100; ic++)
{
  s1 = series (1, ic);
  s2 = series (4, ic);
  s3 = series (5, ic);
  s4 = series (6, ic);
  pid = 4. * s1 - 2. * s2 - s3 - s4;
  pid = pid - (int) pid + 1.;
  //ihex (pid, NHX, chx);
  //ialpha (pid, 1, calpha, cbeta);
  nnew = iint (pid);

  interval[nold][nnew]++;
  nold = nnew;

  //printf ("Pi hex digit computation\n");
  //printf ("position = %i + 1\n", ic);
  //printf ("position = %i + 1\n %20.15f\n %12.12s\n", ic, pid, chx);
  //printf ("position = %i + 1\n %20.15f\n %s\n", ic, pid, chx);

  cout << ic << endl;

  //cout << " 1;+" << calpha[0] << cbeta[0] << "; $30;" << endl;
  //cout << " 59;-" << calpha[0] << cbeta[0] << "; $30;" << endl;
  //output_file << " 1;+" << calpha[0] << cbeta[0] << "; $30;" << endl;
  //output_file << " 59;-" << calpha[0] << cbeta[0] << "; $30;" << endl;

  //	printf ("%s\n",calpha);
  // cout << "position = " << ic << "+ 1\n";
  // cout << pid << endl << chx << calpha << endl;
  }


  display_interval();

  output_file.close();
}


// extra function definitions - G.K.S. @ anti-copyright  1999

void display_interval()
{
int count1;
int count2;

for (count1 = 0; count1 < INTERVAL; count1++)
{
	cout << hx[count1];
	for (count2 = 0; count2 < INTERVAL; count2++)
	{
		cout << setw(3) << interval[count1][count2];
		}
	cout << endl;
	}


}

void initialise_interval()
{
int count1;
int count2;

for ( count1 = 0; count1 < INTERVAL; count1++)
{
	for ( count2 = 0; count2 < INTERVAL; count2++)
		interval[count1][count2]=0;
	}
}








// modified function definitions

int iint (double x)

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

{
  int i,value;
  double y;

  y = fabs (x);
  y = 16. * (y - floor (y));
  value = (int) y;
  return value;
}


void ialpha (double x, int nhx, char calpha[], char cbeta[])

/*  this function is part of the midi file program c2, d2, etc. are note and
octave values that are converted by Gunter Nagler's txt2midi program
*/

{
  int i,d;
  double y;
  char alpha[] = "c2d2e2f2g2a3b3c3d3e3f3g3a4b4c4d4";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
	 y = 16. * (y - floor (y));
	 d = (int) y;
	 calpha[i] = alpha[d*2];
	 cbeta[i] = alpha[d*2+1];
  }
}

//



// original function definitions



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[i] = 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[i] = 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[i] > 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;
}