#include "stdio.h"
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define MAXN 101
// MAXN - rozmiar tablic węzłów i wartości funkcji w węzłach interpolacji
//        zakładamy, że n<MAXN
	int n; // maksymalny indeks węzła interpolacji, (n+1) - liczba węzłów
	int mn; // liczba pośrednich punktów między węzłami interpolacji w drugim pliku
	int opt; // opt=0 - węzły równoodległe, opt=1 - węzły dobrane optymalnie
	int init=0;
	double a, b; //granice przedziału
	double x[MAXN]; //tablica węzłów
	double f[MAXN]; //tablica wartości funkcji w węzłach interpolacji
	char *plkNazwaWezly, *plkNazwa; //nazwy pierwszego i drugiego pliku
// inne tabele i zmienne

double aa[MAXN];
double bb[MAXN];
double cc[MAXN];
double dd[MAXN];




void solveMatrix(double pa[], double pb[], double pc[], double pd[], double xsol[], int n);
//	funkcja na podstawie której wyznaczane są wartości w węzłach interpolacji





	double fun(double t)
{

	return 0.5 + fabs(t*t - t);
}



void ABCD()
{
  double uu[MAXN];
  double ww[MAXN];
  double vv[MAXN];
  double p2[MAXN];
  double hh[MAXN];
  for(int i = 0;i<=n;i++) p2[i] = 2;


	for(int i = 0 ;i<=n-1;i++)
	{
		hh[i] = x[i+1] - x[i];
	}

	for(int i = 0 ;i<n-1;i++)
	{
			uu[i+1] = (hh[i])/(hh[i] + hh[i+1]);
			ww[i+1] = (hh[i+1])/(hh[i] + hh[i+1]);
			vv[i+1] = ((f[i+2] - f[i+1])/(hh[i+1])-((f[i+1] - f[i])/(hh[i])))/(hh[i] + hh[i+1]);
	}

  solveMatrix(uu, p2, ww, vv, cc, n-1);
  cc[0] = cc[n] = 0;
  for(int i = 1;i<=n-1;i++)
  {
    cc[i] = 3*cc[i];
  }
  for(int i = 0;i<=n-1;i++)
  {
    aa[i] = f[i];
    bb[i] = ((f[i+1]-f[i]) / hh[i]) - ((hh[i]/3)*(cc[i+1] + 2*cc[i]));
    dd[i] = (cc[i+1] - cc[i])/(3*hh[i]);
  }

}

void solveMatrix(double pa[], double pb[], double pc[], double pd[], double xsol[], int n)
{
	double alfa[MAXN];
	double beta[MAXN];

	alfa[1] = - pc[1]/pb[1];
	beta[1] = pd[1]/ pb[1];
	for(int i =2; i<=n;i++)
	{
		alfa[i] = - pc[i]/(pb[i] + pa[i]*alfa[i-1]);
		beta[i] = (pd[i] - pa[i]*beta[i-1])/(pb[i] + pa[i]*alfa[i-1]);
	}

	xsol[n] = beta[n];
	for(int i = n-1; i>=1; i--)
	{
		xsol[i] = alfa[i]*xsol[i+1] + beta[i];
	}
}

//	funkcja interpolująca
	double inter(double t)
{

	int j=n-1;
	for(int i=0;i<=n-1;i++)
	{
		if(t<x[i+1])
		{
			j=i;
			break;
		}
	}
	double delta=t-x[j];



	return aa[j] + delta*(bb[j] + delta*(cc[j] + delta*dd[j]));
}

// Wyznaczenie węzłów interpolacji i wartości funkcji w węzłach
void generWezly()
{
	double delta;
	double pi=4*atan(1.0);
	if (opt==0)
	{
	// węzły równoodległe
		delta=(b-a)/n;
		for (int i=0; i<=n; i++)
		{
			x[i]=a+i*delta;
			f[i]=fun(x[i]);
			printf("%f     %f \n", x[i], f[i]);
		}
	}
	else
	{
	// węzły optymalnie dobrane
		delta=(b-a)/2;
		for (int i=0; i<=n; i++)
		{
			x[i]=(a+b)/2 - delta*cos(pi*(2*i+1)/(2*n+2));
			f[i]=fun(x[i]);
			printf("%f     %f \n", x[i], f[i]);
		}
	}
}

// Wyprowadzenie danych do pierwszego pliku
void generWezlyPlik()
{
	FILE *plk;
	plk = fopen( plkNazwaWezly, "w" );
// wydruk nagłówka dla pliku .rex
	fprintf(plk, "xi\t f_xi\n");
// wydruk par węzeł i wartość funkcji w węźle interpolacji
	for (int i=0; i<=n; i++)
	{
		fprintf(plk, "%f\t %f\n", x[i], f[i]);
	}
// wydruk stopki dla pliku .rex
//	fprintf(plk, "*InfKon*\nV 0 0 0 %d",n+1 );
	fclose(plk);
}

//Wyprowadzenie wiersza do drugiego pliku
void drukWiersz(double t, FILE *plk)
{
	double eps=1.0e-20;
	double f, wyn;
	f=fun(t);
	wyn=inter(t);
	fprintf(plk, "%lf\t %lf\t %lf\t %lf\t %lf\n", t, f, wyn, f-wyn, (f-wyn)/(fabs(f)+eps));

}

// Wyprowadzenie danych do drugiego pliku
void generWynikiPlik()
{
	double delta;
	int i, j, lbwier;
	FILE *plk;
	plk = fopen( plkNazwa, "w" );
// wydruk nagłówka dla pliku .rex
	fprintf(plk, "x\t f_x\t w_x\t blb_x\t blw_x\n");
	if (!(opt==0))
	{
	 // dla węzłów optymalnie dobranych wydruk punktów pośrednich w przedziale <a, x[0])
		delta=(x[0]-a)/mn;
		for(j=0; j<mn; j++)
		{
			drukWiersz(a+delta*j, plk);
		}
	}
	 // wydruk punktów w przedziale <x[0], x[n])
	for (i=0; i<n; i++)
	{
		delta=(x[i+1]-x[i])/mn;
		for(j=0; j<mn; j++)
		{
			drukWiersz(x[i]+delta*j, plk);
		}
	}
	if (!(opt==0))
	{
	 // dla węzłów optymalnie dobranych wydruk punktów pośrednich w przedziale <x[n], b)
		delta=(b-x[n])/mn;
		for(j=0; j<mn; j++)
		{
			drukWiersz(x[n]+delta*j, plk);
		};
		lbwier=(n+2)*mn+1;
	}
	else
		lbwier=n*mn+1;
	 // wydruk wydruk wartości dla x=b
	drukWiersz(b, plk);
// wydruk stopki dla pliku .rex
//	fprintf(plk, "*InfKon*\nV 0 0 0 %d", lbwier);
	fclose(plk);
}


int main(int argc, char* argv[])
{
// Odczyt danych z linii polecenia
	a=atof(argv[1]); //lewa granica przedziału
	b=atof(argv[2]); //prawa granica przedziału
	n=atoi(argv[3]); //maksymalny indeks węzła interpolacji
	mn=atoi(argv[4]); //liczba pośrednich punktów między węzłami interpolacji w drugim pliku
	opt=atoi(argv[5]); // opt=0 - węzły równoodległe, opt=1 - węzły dobrane optymalnie
	plkNazwa=argv[6];  //nazwa drugiego pliku (pliku z wynikami działania funkcji interpolującej)
	plkNazwaWezly=argv[7]; //nazwa pierwszego pliku (pliku z węzłami)
// Wyprowadzenie na ekran danych z linii polecenia
	printf("a= %f\n", a);
	printf("b= %f\n", b);
	printf("n= %d\n", n);
	printf("mn= %d\n", mn);
	printf("opt= %d\n", opt);
// Wyznaczenie węzłów interpolacji i wartości funkcji w węzłach
	generWezly();
// Wyprowadzenie danych do pierwszego pliku
	generWezlyPlik();
//Wywołania funkcji generujących parammetry pomocnicze



ABCD();




// Wyprowadzenie danych do drugiego pliku
	generWynikiPlik();



	printf("..............................................");
	return 0;
}