#include <stdlib.h>
#include <assert.h>

typedef struct {int n; double *x, *y, *b, *c;} qspline;

qspline* qspline_make(int n,double* x,double* y){
	double *b  = (double*)malloc(n*sizeof(double));
	double *c  = (double*)malloc(n*sizeof(double));
	double *sx = (double*)malloc(n*sizeof(double));
	double *sy = (double*)malloc(n*sizeof(double));
	qspline* S = (qspline*)malloc(sizeof(qspline));

	for(int i=0;i<n;i++){sx[i]=x[i]; sy[i]=y[i];}
	S->n=n; S->x=sx; S->y=sy; S->b=b; S->c=c;

	double p[n-1], h[n-1]; //VLA

	for(int i=0;i<n-1;i++) {
		h[i]=x[i+1]-x[i];
		p[i]=(y[i+1]-y[i])/h[i]; }

	c[0]=0; //recursion up
	for(int i=0;i<n-2;i++) c[i+1]=(p[i+1]-p[i]-c[i]*h[i])/h[i+1];

	c[n-2]/=2; //recursion down
	for(int i=n-3;i>=0;i--) c[i]=(p[i+1]-p[i]-c[i+1]*h[i+1])/h[i];

	for(int i=0;i<n-1;i++) b[i]=p[i]-c[i]*h[i];
	return S;
}// end qspline

double qspline_eval(qspline *S, double z){
	assert(z>=S->x[0] && z<=S->x[S->n-1]);
	int i=0, k=S->n-1;
	while(k-i>1){int m=(i+k)/2; if(z>S->x[i]) i=m; else k=m;} //bisection:
	double h=z-S->x[i];
	return S->y[i]+h*(S->b[i]+h*S->c[i]); //inerpolating polynomial:
}

void qspline_free(qspline *S){
free(S->x);
free(S->y);
free(S->b);
free(S->c);
free(S);
}