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

void aitken(double* x, double *y, double* res, int n) {
	int i,j=0;
	res[0] = y[0];
	for (i=1; i<n; i++) {
		for (j=0; j<n-i; j++) {
			y[j] = (y[j+1] - y[j])/( x[j+i] - x[j]);
		}	
		res[i] = y[0];
	}
	printf("%f", res[0]);
	for ( i = 1; i<n; i++) {
		printf(" + %f ", res[i]);
		for(j=0; j<i; j++) {
			printf("*(x-%f)", x[j]);
		}
	}
	printf("\n");
}

void lagrange(double *x, double *y, double* res, int n) {
	int i,j;
	for (i=0; i<n; i++) {
		double sum =1;
		for (j=0;j<n; j++) {
			if (i!=j)
				sum *= x[i] - x[j];
		}
		printf("+ (%f/%f) ", y[i], sum);
		for (j=0;j<n; j++) {
			if (i!=j)
				printf("*(x-%f)",x[j]); 
		}

	}
	printf("\n");
}

void catmull(double *x, double *y, int n) {
	int i;
	double* diffy = malloc (sizeof(double) *n);
	if (diffy == NULL) {
		perror("malloc");
		exit(EXIT_FAILURE);
	}
	diffy[0] = (y[1] - y[0]) / (x[1] - x[0]);
	diffy[n-1] = (y[n-1] - y[n-2]) / (x[n-1] - x[n-2]);
	for(i=1; i<n-1; i++)  
		diffy[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1]);
	for(i=0; i<n-1; i++) { 
		printf("p_%d(x) = ", i);
		double a_0 = y[i] / pow((x[i+1] - x[i]), 3);
		double a_3 = y[i+1] / pow((x[i+1] - x[i]), 3);
		double a_1 = 3*a_0 + (diffy[i] / pow((x[i+1] - x[i]),2));
		double a_2 = 3*a_3 - (diffy[i+1] / pow((x[i+1] - x[i]),2));
		printf("%f * (%f - x)^3 + %f * (%f-x)^2 * (x-%f) + %f * (%f-x) * (x-%f)^2 + %f * (x-%f)^3\n", 
			a_0, x[i+1], a_1, x[i+1], x[i], a_2, x[i+1], x[i], a_3, x[i]);
	}
	
}


int main(void) {
	int i=0;
	double x[] = {1, 2, 4,7,8};
	double y[] = {1, -1, 2,0,3};
	int n = 5;
	double res[n];
	printf("Lagrange \n ====================== \n");
	lagrange(x, y, res, n);
	printf("\nCatmull-Rom\n ======================== \n");
	catmull(x,y,n);
	printf("\nAitken-Neville \n ====================== \n");
	aitken(x, y, res, n);

}