#include "nr3.h"
#include "poly.h"
#include "interp_1d.h"
#include "ludcmp.h"
#include "pade.h"
#include <fstream>

using namespace std;

int main() {
    ofstream out("out.txt");
    for (int N=4; N<17; N+=3) {
	const double PI=acos(-1); // define PI to machine accuracy
	VecDoub coef(N);  // the power series coefficients
	VecDoub x(81); // the interpolation grid points, x-abcissa
	VecDoub y(81); // the interpolation grid points, function value
	double rms[5]; // for the RMS
	for (int i=0; i<81; i++) { // fill the x,y vectors
		x[i]=0.25*i-10;
		y[i]=log(fabs(x[i]+PI));
	}
	Poly_interp pol(x,y,N); // pol does Nth-order interpolation
	Spline_interp splin(x,y); // the natural cubic spline
	coef[0]=log(PI); // function value at x=0
	for (int i=1; i<N; i++) {
		coef[i]=-pow(-PI,-i)/i; // power series, x>-PI
	}
	Ratfn rt=pade(coef); // rt = pade function
	for (int i=0; i<4; i++) rms[i]=0;
	for (double xx=-10; xx<10.01; xx+=0.01) {
		double yy = xx;
		if (xx<-PI) yy=-2*PI-xx;
		double power=0;
		for (int i=0; i<N; i++) power+=coef[i]*pow(yy,i);
		double calc=log(fabs(xx+PI));
		out << setw(6) << xx << " " <<N << " " << setw(12) << calc;
		out << " " << setw(12) << splin.interp(xx)-calc;
		out << " " << setw(12) << pol.interp(xx)-calc;
		out << " " << setw(12) << rt(yy)-calc;
		out << " " << setw(12) << power-calc ;
		out << endl;
		rms[0]+=pow(calc-splin.interp(xx),2);
		rms[1]+=pow(calc-pol.interp(xx),2);
		rms[2]+=pow(calc-rt(yy),2);
		rms[3]+=pow(calc-power,2);
	}
	cout << " N " << N << endl;
	cout << "spline , RMS " << sqrt(rms[0]/2001) << endl; 
	//cout << "3-rd order polynomial, RMS " << sqrt(rms[1]/2001) << endl; 
	cout << "N-th order polynomial, RMS " << sqrt(rms[1]/2001) << endl; 
	cout << "pade , RMS " << sqrt(rms[2]/2001) << endl; 
	cout << "power series , RMS " << sqrt(rms[3]/2001) << endl; 
	cout << endl;
   }
}