#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>

using namespace std;

void fill_lhs(int ndata,int npoints,vector <double> x,int i, int j, double &res)
{
	res = 0.e0;
	double expo=double(i+j);
	for(int m=0;m<npoints;m++)
		res += pow(x[m],expo);
}

void fill_rhs(int ndata,int npoints, vector <double> x, vector <double> y, int i, double &res)
{
	res = 0.e0;
	for(int m=0;m<npoints;m++)
	{
		res=res+y[m]*pow(x[m],double(i));
	}
}

void gauss(int kdim, int n, vector < vector <double> > a, vector <double> b, vector <double> &x)
{
	double eps = 1.e-8;
	
	for(int k=0;k<n-1;k++)
	{
		if(abs(a[k][k]<eps))
		{
			cout<<"Divide by Zero !!!!"<<endl;
			break;
		}
		for(int i=k+1;i<n;i++)
		{
			double factor=a[i][k]/a[k][k];
			for(int j=k+1;j<n;j++)
				a[i][j]-=factor*a[k][j];
			b[i]-= factor*b[k];
		}
	}
	
	x[n-1]= b[n-1]/a[n-1][n-1];
	for(int i=n-2;i>=0;i--)
	{
		double sum=0.e0;
		for(int j=i+1;j<n;j++)
			sum+=a[i][j]*x[j];
		x[i]=(b[i]-sum)/a[i][i];
	}
}

int main(int argc, char** argv)
{
	string line;
	int ndata=100,npol=10;
	vector < vector <double> > a(npol, vector <double> (npol));
	vector <double> x(ndata),y(ndata),b(npol),coef(npol);
	ifstream in("data");
	int n=0;
	while(!in.eof())
	{
		in>>x[n]>>y[n];
		getline(in,line);
		n++;
	}
	n--;
	int npoints = n;
	int mpol;
	cout<<"Enter polynomial degree "<<endl;
	cin>>mpol;
	double res;
	for(int i=0;i<mpol+1;i++)
	{
		for(int j=i;j<mpol+1;j++)
		{
			fill_lhs(ndata,npoints,x,i,j,res);
			a[i][j] = res;
			a[j][i] = res;
		}
		fill_rhs(ndata,npoints,x,y,i,b[i]);
	}

	gauss(npol,mpol+1,a,b,coef);

	
	ofstream out("approx");
	for(int i=0;i<npoints;i++)
	{
		double app=0.e0;
		for(int j=0;j<mpol+1;j++)
			app=app+coef[j]*pow(x[i],j);
		double error = 100.e0*(app-y[i])/y[i];
		out<<x[i]<<"\t"<<app<<"\t"<<error<<endl;
	}
		
}