#include <math.h>
#include <iostream>
using namespace std;

double lagrange(double x, double *x_i, double *y, int n)
{
	double fx=0.0, tmp=1.0;
	int i, j;
	for(i=0; i<=n; i++)
	{
		tmp=1.0;
		for(j=0; j<i; j++)
			tmp = tmp*(x-x_i[j])/(x_i[i]-x_i[j]);
		for(j=i+1; j<=n; j++)
			tmp = tmp*(x-x_i[j])/(x_i[i]-x_i[j]);
		fx = fx+tmp*y[i];
	}
	return fx;
}

double n_max(double *y_f, double *y_p)
{
	int i,j;
	double max=0.0;
	for(i=0; i<101; i++)
	{
		if(max <= fabs(y_f[i]-y_p[i]))
		{
			max = fabs(y_f[i]-y_p[i]);
			j=i;
		}
	}
	return max;
}

void main(void)
{
	double x[101], y_f[101], y_p1[101], y_p2[101], x1[41], x2[41], y1[41], y2[41];
	int n[4]={5, 10, 20, 40};
	int i, j;
	double const pi=3.14159265358979323846;

	for(i=0; i<101; i++)
	{
		x[i] = i/10.0-5.0;
		y_f[i] = 1.0/(1.0+x[i]*x[i]);
	}

	for(i=0; i<4; i++)
	{
		for(j=0; j<=n[i]; j++)
		{
			x1[j] = 5.0-10.0/n[i]*j;
			y1[j] = 1.0/(1.0+x1[j]*x1[j]);
			
			x2[j] = -5.0*cos((2.0*j+1.0)/(2.0*n[i]+2.0)*pi);
			y2[j] = 1.0/(1.0+x2[j]*x2[j]);
		}
		for(j=0; j<101; j++)
		{
			y_p1[j] = lagrange(x[j], x1, y1, n[i]);
			y_p2[j] = lagrange(x[j], x2, y2, n[i]);
		}
		printf("N = %d\nMax Error of grid (1) : %10.12f\nMax Error of grid (2) : %10.12f\n\n"
			, n[i], n_max(y_f, y_p1), n_max(y_f, y_p2));
	}
	system("PAUSE");
}