#include <iostream>
#include <iomanip>
using namespace std;
#include <math.h>
#include <stdlib.h>
#include "Rational.h"

#if defined(__KCC_33)
// Work-around version of operator<< for long long.
// Templatizing this should be straight-forward.
// This version presumes decimal output and that (long double) can
// hold a (long long) without loss of precision.
ostream& operator<<( ostream& o, long long x ) {
    // Set precision wide enough to avoid exponential notation.
    streamsize save_precision = o.precision( 20 );

    // Turn off fixed-point or scientific notation.
    ios_base::fmtflags save_fmt = o.flags();
    o.unsetf( ios_base::floatfield );

    // Print the number
    o << (long double)x;

    // Restore the format information 
    o.precision( save_precision );
    o.setf( save_fmt );

    // Return the stream.
    return o;
}
#endif


#ifndef LONGONLY
typedef long long base;
const int HNUM = 32;
const int BNUM = 32;
const int PREC = 12;
#else
typedef long base;
const int HNUM = 25;
const int BNUM = 23;
const int PREC = 8;
#endif

template class Rational<base>;

Rational<base> harmonic(int n) {
  Rational<base> r = 1;
  for (int i=2; i<=n; ++i) r += Rational<base>(1,i);
  return r;
}

Rational<base> binomial(int n, int k) {
  if (n < 0) { cerr << "1st index out of range" << endl; exit(1); }
  if (k < 0 || k > n) { cerr << "2nd index out of range" << endl; exit(1); }
  if (k > n-k) k = n-k;
  if (k == 0) return 1;
  return Rational<base>(n-k+1, k) * binomial(n, k-1);
}

Rational<base> bernoulli(int n) {
  if (n < 0) { cerr << "index out of range" << endl; exit(1); }
  else if (n == 0) return 1;
  else if (n == 1) return Rational<base>(-1,2);
  if (n % 2) return 0;
  Rational<base> r = 0;
  for (int k=0; k<n; ++k) {
    r -= binomial(n+1, k) * bernoulli(k);
  }
  r /= Rational<base>(n+1);
  return r;
}

int main() {
  int n;

  cout.precision(PREC);
  cout << "n\tEuler Approx.\tHarmonic(n)" << endl;
  cout << "========================================" << endl;
  for (n=1; n<HNUM; ++n) {
    Rational<base> r = harmonic(n);
    double g = toDouble(r) - log((double)n);
    g -= (1.0/(2*n)) * (1.0 - (1.0/(6*n)));
    cout << n << '\t' << g << '\t' << r << endl;
  }

  cout << endl << "n\tBernoulli(n)" << endl;
  cout << "========================================" << endl;
  Rational<base> b;
  for (n=0; n<BNUM; ++n) {
    if ((b=bernoulli(n)) != Rational<base>(0,1)) cout << n << '\t' << b << endl;
  }
}