// сплайн.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <math.h>
#include <fstream>
#define Pi 3.1415926535897932384626433832795
#define a 0.0
#define b 2.0
#define n 20
using namespace std;

//пишем функцию тесттовую
double f(double x)
{
	return  3*x*x*x+2*x*x+x+11;//sin(x)+cos(x);
}

int main()
{
	int i;
	//массив сетки неравномерной
	//объявление
	double *h, *x, *y, *g, *z, *d, *c, *e, h_max=0;
	//создаем
	h=new double [n+1];
	h[0]=0;//h[0] не будем использовать-для удобства нумерации
	x=new double [n+1];
	y=new double [n+1];
	z=new double [n];// массив правой части
	e=new double [n];// верхние диагональные элементы, 0й элемент не будем использовать для удобства нумерации
	//массив неизвестных коэффициентов c
	c=new double [n+1];
	g=new double [n+1];
	d=new double [n+1];


	//cчитываем сетку с файла
	ifstream ff("in.txt");
		for(i=0;i<=n;i++) ff>>x[i];
	ff.close();


	for(i=1;i<=n; i++)
		h[i]=x[i]-x[i-1];


/*	//задаем равномерную сетку
	h_max=(b-a)/(n+0.0);
	for(i=1;i<=n;i++) h[i]=h_max;
	for(i=0;i<=n;i++) x[i]=a+i*h_max;


*/

	for (i=0; i<=n; i++)
		y[i]=f(x[i]);

	//вычисляем вторую производную

	if(h[1]<h[2]) h_max=h[2];
	else h_max=h[1];
	c[0]=(y[0]-2*y[1]+y[2])/h_max;

	if(h[n-1]<h[n]) h_max=h[n];
	else h_max=h[n-1];
	c[n]=(y[n]-2*y[n-1]+y[n-2])/h_max;

	c[0]=4;
	c[n]=40;
	//найдем коэффициенты c с помощью метода прогонки

	//зададим правый столбец системы
	for (i=1; i<n; i++)
		z[i]=6.0*(y[i+1]-y[i])/h[i+1]-6.0*(y[i]-y[i-1])/h[i];

	z[1]-=h[1]*c[0];
	z[n-1]-=h[n]*c[n];

	// вычисление верхних диагональных элементов
	e[1]=h[2]/(2.0*(h[1]+h[2]));

	for (i=2; i<n; i++)
	e[i]=h[i+1]/(2.0*(h[i]+h[i+1])-h[i]*e[i-1]);

	//пересчитываем правый столбец
	z[1]=z[1]/(2.0*(h[1]+h[2]));
	for (i=2; i<n; i++)
		z[i]=(z[i]-h[i]*z[i-1])/(2.0*(h[i]+h[i+1])-h[i]*e[i-1]);

	//обратный ход
	c[n-1]=z[n-1];

	for(i=n-2; i>=1; i--)
		c[i]=z[i]-e[i]*c[i+1];

	//УРА!!!! нашли коэффициенты с, теперь находим коэффициенты d and g
	for(i=1;i<=n;i++)
	{
		d[i]=(c[i]-c[i-1])/h[i];
		g[i]=h[i]*c[i]/2.0-h[i]*h[i]*d[i]/6.0+(y[i]-y[i-1])/h[i];
	}
	double xo=0.2, s;//s - значение сплайна в точке хо

	//найдем интервал, которому принадлежит точка хо
	for(i=1;i<=n;i++)
		if(x[i]-xo>0) break;

	//s=y[i-1]*(x[i]-xo)/h[i]+y[i]*(xo-x[i-1])/h[i]+g[i-1]*(pow((x[i]-xo),3)-h[i]*h[i]*(x[i]-xo))/6*h[i]+g[i]*(pow((xo-x[i-1]),3)-h[i]*h[i]*(xo-x[i]))/6*h[i];
	s=y[i]+g[i]*(xo-x[i])+c[i]*(xo-x[i])*(xo-x[i])/2.0+d[i]*(xo-x[i])*(xo-x[i])*(xo-x[i])/6.0;
	cout<<s<<'\n';
	cout<<f(xo)<<'\n';

	//записать в файл для того, что бы затем построить график

	ofstream off("out.txt");
	//off<<"with(plots);"<<endl<<endl;
	off<<"p1:=plot(function, x="<<a<<".."<<b<<", color=blue):"<<endl<<endl;
	for(i=1; i<=n;i++)
		off<<"s"<<i<<":=plot("<<y[i]<<"+("<<g[i]<<")*(x-"<<x[i]<<")+("<<c[i]/2.0<<")*(x-"<<x[i]<<")^2+("<<d[i]/6.0<<")*(x-"<<x[i]<<")^3, x="<<x[i-1]<<".."<<x[i]<<"):"<<endl<<endl;
	off<<"display([ p1, ";
	for(i=1;i<=n;i++)
		off<<"s"<<i<<", ";
	off<<"]);";

	off.close();



	return 0;
}