/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*- */
/*
 * pi_hex_calc.c
 * Copyright (C) dmgualtieri 2012 <gualtieri**at**ieee.org>
 * 
 * pi_hex_calc is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * pi_calc is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License along
 * with this program.  If not, see <http://www.gnu.org/licenses/>.
 */


/*  Program to produce hex digits of pi.
    Adapted from David H. Bailey code of 960429 
    This code is valid up to ic = 2^24 on systems with IEEE arithmetic.
    First drdev version, August, 2002.  Revised March 13, 2012
*/

#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include <time.h>

//Prototype
void exit( int exit_code );

char fn1[64] = "pi_hex.txt";  // File of hex digits of pi.
FILE *outdata;

int main(void)
{
  double pid, s1, s2, s3, s4;
  double series (int m, int n);
  void ihex (double x, int m, char c[]);
  int ic=0;
  int j=0;
#define NHX 16
  char chx[NHX];
  time_t t1,t2;

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

printf("\nOutput file selected = %s\n",fn1);

if ((outdata = fopen(fn1,"w"))==NULL)
	{printf ("\nOutput file cannot be opened.\n");
	exit (1);}

	fprintf(outdata,"\nHex Digits of pi\n");
	printf("0:\t");
    fprintf(outdata,"0:\t");

  (void) time(&t1); /* get starting time of calculation */

for(ic = 0;ic<20001;ic=ic+10)
{
  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);
  printf ("%.10s",chx);
  fprintf(outdata,"%.10s",chx);
  j=j+1;
  if(j==8)
  {
   j=0;
   printf("\n%d:\t",ic+10);
   fprintf(outdata,"\n%d:\t",ic+10);
  }
}

(void) time(&t2); /* get ending time of calculation */
printf ("\ndigit = %d, time = %d seconds\n",ic,(int)(t2-t1));

fclose(outdata);
   printf("\n");
   return 1;
}

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, 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;
}