/* bench.cpp -- Matt Mahoney, mmahoney@cs.fit.edu

This program is a benchmark for Integer and Real.  It computes pi by
two different methods at increasingly higher precisions and compares
the results.  Normally all but the last few digits will match due
to roundoff errors in the calculations.

The calculations are deliberately inefficient to test all
the operators.  It accepts one argument indicating the number of
comparisons to be done (default 10).  After each comparison it prints
the average of the two values of pi (3.1415926535897...), the number
of matching decimal digits, and the total run time so far.  For UNIX, it
is recommended that the time command be used to allow for CPU time used
by other processes.

The program assumes the unlimited range and precision numeric classes
Integer and Real as described at http://cs.fit.edu/~mmahoney/cse2050
are available in two files in the current directory.

  mmahoney.h   -- Class declarations and inlined code
  mmahoney.cpp -- Function definitions

To compile: g++ bench.cpp mmahoney.cpp
To run, e.g.: time a.out 10

*/

#include "mmahoney.h"  // Replace this with your header file
#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cmath>
using namespace std;

// Absolute value
Real abs(const Real& a) {
  if (a<Integer(0))
    return Integer(0)-a;  // 0-a might increase precision to 20
  else
    return a;
}

// Square root to prec decimal places by Newton's method.  Each iteration
// doubles the precision of the result.
Real sqrt(const Real& a, int prec) {
  Real r(1, prec);
  for (int i=1; i<prec;) {
    r=(r+a/r)/Integer(2);
    if (abs(r*r-a)<0.1)
      i*=2;
  }
  return r;
}

// Factorial, a >= 0
Integer factorial(const Integer& a) {
  Integer r=1;
  for (Integer i=1; i<=a; i+=1)
    r*=i;
  return r;
}

// Exponentiation, b >= 0,
template <class T>
T pow(const T& a, const Integer& b) {
  T r=1;
  for (Integer i=0; i<b; i+=1)
    r*=a;
  return r;
}

// Compare 2 values, return the number of matching decimal places up to prec
int compare(const Real& a, const Real& b, int prec) {
  int matches=0;
  Real diff=abs(a-b);
  for (int i=0; i<prec; ++i)
    if ((diff*=10) < 1)
      ++matches;
  return matches;
}

// Compute pi to prec places by Ramanujan's method.
// Each iteration increases the precision by 8 decimal places.
Real ramanujan(int prec) {
  Real sum(0, prec);
  for (int k=0; k*8<prec; k+=1) {
    sum += factorial(4*k)*(1103+k*Real(26390, prec)) /
           (pow(factorial(k), 4) * pow(Integer(396), 4*k));
  }
  return 9801/(sum*sqrt(8, prec));
}

// Compute pi to prec places by Borwein's quartic method.
// Each iteration quadruples the precision.
Real borwein(int prec) {
  Real y=sqrt(2, prec)-1;
  Real a=6-4*sqrt(2, prec);
  for (Integer k=0; k<0.3+log(prec)/log(4); k+=1) {
    Real y4=sqrt(sqrt(1-y*y*y*y, prec), prec);
    y=(1-y4)/(1+y4);
    a=a*pow(1+y, 4)-pow(2, 2*k+3)*y*(1+y+y*y);
  }
  return 1/a;
}

// Take one argument and do that number of comparisons (default 10).
// Start with precision 16 and increase by 50% each time.
int main(int argc, char** argv) {
  clock_t start_time=clock();  // Start timer
  int limit=argc>1 ? atoi(argv[1]) : 10;
  int prec=16;
  for (int i=0; i<limit; ++i) {
    Real pi1 = ramanujan(prec);
    Real pi2 = borwein(prec);
    cout << (pi1+pi2)/2 << endl;  // Average of two pi's
    cout << compare(pi1, pi2, prec) << " of " << prec
      << " digits match after "
      << double(clock()-start_time)/CLOCKS_PER_SEC << " seconds" << endl;
    prec=prec*3/2;
  }
  return 0;
}