// interpolacja_newtona_11k2.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <cmath>
#include <iostream>

using namespace std;

double f(double x);
double W(double x1, double x2, double x3, double x, double a0, double a1, double a2, double a3);

int _tmain(int argc, _TCHAR* argv[])
{
	double x1, x2, x3, dx, er, ermax, erdop, xStart, xEnd, xp, a0, a1, a2, a3;
	int m = 3; // przedzialy
	double d[5][5];
	double x[5];

	xStart = 0.0; 
	xEnd = 2.0*asin(1.0);

	cout << endl << "Podaj blad="; 
	cin >> erdop;

	do
	{
		ermax = 0.0;
		xStart = 0.0;
		xEnd = 2.0*asin(1.0);
		dx = (xEnd - xStart) / m;

		x[1] = xStart;
		x[4] = xEnd / m;
		x[2] = x[1] + (x[4] - x[1]) / 3;
		x[3] = x[1] + (x[4] - x[1]) * 2 / 3;
		xp = x[1] + dx / 2.0;
		

		for (int i = 0; i<m; i++)
		{
			d[1][1] = f(x[1]);
			d[2][1] = f(x[2]);
			d[3][1] = f(x[3]);
			d[4][1] = f(x[4]);
			
			for (int i = 2; i < 5; i++)
			{
				for (int j = 2; j < 5; j++)
				{
					d[i][j] = (d[i][j - 1] - d[i - 1][j - 1]) / (x[i] - x[i - j + 1]);
				}
			}

			a0 = d[1][1];
			a1 = d[2][2];
			a2 = d[3][3];
			a3 = d[4][4];

			er = f(xp) - W(x[1], x[2], x[3], xp, a0, a1, a2, a3);

			if (er<0)
				er = -er;

			if (er>ermax)
				ermax = er;

			x[1] = x[4];
			x[4] = x[1] + dx;
			x[2] = x[1] + (x[4] - x[1]) / 3;
			x[3] = x[1] + (x[4] - x[1]) * 2 / 3;
			xp = x[1] + dx / 2.0;
		}
		m++;
	} while (ermax>erdop);

	printf("\nliczba przedzialow=%d blad=%f", m - 1, ermax);
	cout << endl;
	system("pause");
	return 0;
}

double f(double x)
{
	return sin(x);
}
//---------------------------------------
double W(double x1, double x2, double x3, double x, double a0, double a1, double a2, double a3)
{
	return  a0+ a1*(x - x1) + a2*(x - x1)*(x - x2) + a3*(x - x1)*(x - x2)*(x - x3);
}