#include <fstream.h>
#include <iomanip.h>
#include <vector>


// Chapter 7.1, p.320
// Spline 1
template <class Type>
Type spline1(int n, std::vector<Type> t, std::vector<Type> y, Type x)
{
	int i;
	for (i = n - 1; i >= 0; i--)
		if (x - t[i] >= 0)
			break;
	return y[i] + (x - t[i]) * ((y[i+1] - y[i]) / (t[i+1] - t[i]));
}


// Chapter 7.2, p.338
// Spline 3
template <class Type>
void spline3_coef (int n, std::vector<Type> t, std::vector<Type> y, std::vector<Type>& z)
{
	std::vector<Type> h(n+1, 0.0);
	std::vector<Type> b(n+1, 0.0);
	std::vector<Type> u(n+1, 0.0);
	std::vector<Type> v(n+1, 0.0);
	int i;
  
	for (i = 0; i <= n - 1; i++)
	{
		h[i] = t[i+1] - t[i];
		b[i] = (y[i+1] - y[i])/h[i];
	}
	u[1] = 2*(h[0] + h[1]);
	v[1] = 6*(b[1] - b[0]);
	for (i = 2; i <= n - 1; i++)
	{
		u[i] = 2*(h[i] + h[i-1]) - (h[i-1]*h[i-1])/u[i-1];
		v[i] = 6*(b[i] - b[i-1]) - (h[i-1]*v[i-1])/u[i-1];
	}
	z[n] = 0;
	for (i = n - 1; i >= 0; i--)
		z[i] = (v[i] - h[i]*z[i+1])/u[i]; 
	z[0] = 0;
}

template <class Type>
Type spline3_eval(int n, std::vector<Type> t, std::vector<Type> y, std::vector<Type> z, Type x)
{
	int i;
	Type h, tmp;

	for (i = n -1; i >= 0; i--)
		if (x - t[i] >= 0)
			break;

	h = t[i+1] - t[i];
	tmp = (z[i]/2) + (x - t[i])*(z[i+1] - z[i])/(6*h);
	tmp = -(h/6)*(z[i+1] + 2*z[i]) + (y[i+1] - y[i])/h + (x - t[i])*tmp;
	return (y[i] + (x - t[i])*tmp);
}


// Chapter 7.4, p.360-361
// BSpline 2
template <class Type>
void bspline2_coef(int n, std::vector<Type> t, std::vector<Type> y, std::vector<Type>& a, std::vector<Type>& h)
{
	int i;
	Type delta, gamma, p, q, r;

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

	h[0] = h[1];
	h[n+1] = h[n];
	delta = -1.0;
	gamma = 2*y[0];
	p = delta*gamma;
	q = 2;
	for (i = 1; i <= n; i++)
	{
		r = h[i+1]/h[i];
		delta = -r*delta;
		gamma = -r*gamma + (r + 1)*y[i];
		p += gamma*delta;
		q += delta*delta;
	}
	a[0] = -(p/q);
	for( i = 1; i <= n + 1; i++)
		a[i] = ((h[i-1] + h[i])*y[i-1] - h[i]*a[i-1]) / h[i-1];   
}

template <class Type>
Type bspline2_eval(int n, std::vector<Type> t, std::vector<Type> a, std::vector<Type> h, Type x)
{
	int i;
	Type d,e;

	for (i = n - 1; i >= 0; i--)
		if(x - t[i] >= 0) 
			break;

	i += 1;
	d = ( a[i+1]*(x - t[i-1]) + a[i]*(t[i] - x + h[i+1]) ) / (h[i] + h[i+1]);
	e = ( a[i]*(x - t[i-1] + h[i-1]) + a[i-1]*(t[i-1] - x + h[i]) ) / (h[i-1] + h[i]);
	return (d*(x - t[i-1]) + e*(t[i] - x)) / h[i];
}


// Chapter 7.4, p.363-364
// Schoenberg
template <class Type>
void sch_coef (Type f(Type), Type a, Type b, int n, std::vector<Type>& d)
{
	Type h;
	int i;

	h = (b - a)/n;

	for (i = 2; i <= n + 2 ; i++)
		d[i] = f(a + (i-2)*h);

	d[1] = 2*d[2] - d[3];
	d[n+3] = 2*d[n+2] - d[n+1];
}

template <class Type>
Type sch_eval(Type a, Type b, int n, std::vector<Type> d, Type x)
{
	int k;
	Type c, e, h, p;

	h = (b - a)/n;
	k = (int) ((x - a)/h + 5/2);
	p = x - a - (k - 5/2)*h;
	c = (d[k+1]*p + d[k]*(2*h - p))/(2*h);
	e = (d[k]*(p + h) + d[k-1]*(h - p))/(2*h); 
	return (c*p + e*(h - p))/h;
}