/* $Id: pi-bbp-openmp3.c 859 2006-06-23 08:08:01Z olau $ */

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

#ifdef WIN32
#include <windows.h>
#include "getopt/getopt.h"
#else
#include <getopt.h>
#endif

#ifdef _OPENMP
#include <omp.h>
#include "ompdefs.h"
#endif

#include "globaldefs.h"
#include "timer/timer.h"

#define DEFAULT_STARTPOS       1000000
#define NHX                         16

/*  This returns, in chx, the first nhx hex digits of the fraction of x. */
void ihex(double x, int nhx, char chx[])
{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs(x);
  for (i = 0; i < nhx; i++) {
    y = 16.0 * (y - floor(y));
    chx[i] = hx[(int) y];
  }
}

/*  expm = 16^p mod ak.  This routine uses the left-to-right binary 
exponentiation scheme.  It is valid for  ak <= 2^24. */
double expm(double p, double ak)
{
  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.0;
    for (i = 1; i < ntp; i++)
      tp[i] = 2.0 * tp[i-1];
  }

  if (ak == 1.0) return 0.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.0;

  /*  Perform binary exponentiation algorithm modulo ak. */
  for (j = 1; j <= i; j++){
    if (p1 >= pt) {
      r = 16.0 * r;
      r = r - (int) (r / ak) * ak;
      p1 = p1 - pt;
    }
    pt = 0.5 * pt;
    if (pt >= 1.0) {
      r = r * r;
      r = r - (int) (r / ak) * ak;
    }
  }

  return r;
}

/*  This routine evaluates the series  sum_k 16^(ic-k)/(8*k+m) 
using the modular exponentiation technique. */
double series(int m, int ic)
{
  int k;
  double ak, p, s, t;
  const double eps = 1e-17;

  s = 0.0;
  /*  Sum the series up to ic. */
#pragma omp parallel for shared(m) private(ak, p, t) reduction(+:s)
  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.0, (double) (ic - k)) / ak;
    if (t < eps)
      break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}


void usage(void)
{
  printf("Aufruf mit: pi-bbp [-p <Position der ersten Ziffer>]\n");
}


/* here we go ... */
int main(int argc, char *argv[])
{
  double pid, s1, s2, s3, s4;
  double t;
  int ic = DEFAULT_STARTPOS;
  char chx[NHX];
  int option;

  while ((option = getopt(argc, argv, "p:?")) != -1) {
    switch (option) {
      case 'p':
        if (optarg)
          ic = atoi(optarg);
        break;
      case '?':
        usage();
        break;
      default:
        /* ignore */
        break;
    }
  }

  t = omp_get_wtime();
  /*  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.0;
  ihex(pid, NHX, chx);

  printf("Startposition = %i\n"
    "hexadezimale Nachkommastellen = %10.10s\n",
    ic, chx);
  printf("Ausfuehrungsdauer: %.2lf ms\n\n", 1000 * (omp_get_wtime() - t));

  return 0;
}

