#include <iostream>
#include <iomanip>
using namespace std;
#include <math.h>
#include <stdlib.h>
#include "Rational.h"
#if defined(__KCC_33)
ostream& operator<<( ostream& o, long long x ) {
streamsize save_precision = o.precision( 20 );
ios_base::fmtflags save_fmt = o.flags();
o.unsetf( ios_base::floatfield );
o << (long double)x;
o.precision( save_precision );
o.setf( save_fmt );
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;
}
}