public class cspline{

double[] x,y,b,c,d;
int n => x.Length;

int binsearch(double z){
if(z<x[0] || z>x[n-1]) throw new System.ArgumentException();
int i=0, j=n-1;
while(j-i>1){ int mid=(i+j)/2; if(z>x[mid]) i=mid; else j=mid;}
return i;
}

public cspline(double[] xs,double[] ys){
x=(double[])xs.Clone();
y=(double[])ys.Clone();
if(ys.Length<n) throw new System.ArgumentException("cspline: wrong y");
b=new double[n];
c=new double[n-1];
d=new double[n-1];
var D = new double[n];
var Q = new double[n-1];
var B = new double[n];
var h = new double[n-1];
var p = new double[n-1];
for(int i=0;i<n-1;i++) {
	h[i]=x[i+1]-x[i];
	if(h[i]<=0) throw new System.ArgumentException("cspline: wrong x");
	}
for(int i=0;i<n-1;i++) { p[i]=(y[i+1]-y[i])/h[i]; }

D[0]=2; Q[0]=1; B[0]=3*p[0]; D[n-1]=2; B[n-1]=3*p[n-2];
for(int i=0;i<n-2;i++){
	D[i+1]=2*h[i]/h[i+1]+2;
	Q[i+1]=h[i]/h[i+1];
	B[i+1]=3*(p[i]+p[i+1]*h[i]/h[i+1]);
	}

for(int i=1;i<n;i++){
	D[i]-=Q[i-1]/D[i-1];
	B[i]-=B[i-1]/D[i-1];
	}

b[n-1]=B[n-1]/D[n-1];
for(int i=n-2;i>=0;i--)
	b[i]=(B[i]-Q[i]*b[i+1])/D[i];

for(int i=0;i<n-1;i++){
	c[i]=(-2*b[i]-b[i+1]+3*p[i])/h[i];
	d[i]=(b[i]+b[i+1]-2*p[i])/h[i]/h[i];
	}
}//construcntor

public double eval(double z){/* evaluation of the spline at point z */
int i=binsearch(z);
double dx=z-x[i];/* calculate the inerpolating spline : */
return y[i]+dx*(b[i]+dx*(c[i]+dx*d[i]));
}

public double deriv(double z){/* derivative of the spline at point z */
int i=binsearch(z);
double dx=z-x[i];
return b[i]+dx*(2*c[i]+dx*3*d[i]);
}

public double integ(double z){/* derivative of the spline at point z */
int iz=binsearch(z);
double sum=0,dx;
for(int i=0;i<iz;i++){
	dx=x[i+1]-x[i];
	sum+=dx*(y[i]+dx*(b[i]/2+dx*(c[i]/3+dx*d[i]/4)));
	}
dx=z-x[iz];
sum+=dx*(y[iz]+dx*(b[iz]/2+dx*(c[iz]/3+dx*d[iz]/4)));
return sum;
}

}//cspline