#include <iostream.h>
#include <stdlib.h>
//#define NDEBUG
//uncomment the above to remove assert() checking
//when everything works correctly
#include <assert.h>

/***************************

$C^1$ splines

We imagine the interval $[a,b]$ partitioned into $n$ subintervals by
the points $t_i$, $0\le i \le n$. The intervals are $[t_i,t_{i+1}]$,
$0\le i \le n-1$, and so that $t_i = a + i(b-a)/n$ for $0 \le i \le
n$.

There's two kinds of $C^1$ cubic elements adapted to the partition.
One, the $\phi_i$ has

$$\phi_i(x) =
\begin{cases}
x*x*(3.0-2.0*x) & \mbox{$0 \le x \le 1$}\\
(2.0-x)*(2.0-x)*(2.0*x-1) & \mbox{$1 \le x \le 2}\\
0 & \mbox{ otherwise}
\end{cases}
$$

		  * *
		*     *
	 *         *
  *             *
**---------------**
$t_{i-1}$       $t_{i+1}$

The $\psi_i(x)$ look like  :

				* *
			  *   *
**-------*------**
	*    *
	 * *
$$\psi_i(x) =
\begin{cases}
x*x*(x-1.0) & \mbox{$0 \le x \le 1$ }\\
(2.0-x)*(2.0-x)(1.0-x) & \mbox{$1 \le x \le 2$}\\
0 & \mbox{ otherwise}
\end{cases}
$$


*/

/*
Here $[a,b]$ is the interval on which we want the spline
elements, and $n$ is the number of subintervals on which we
define the elements. Linear combinations of th elements
will be piecewise $C^1$ and adapted the partition of $[a,b]$ into
$n$ intervals of equal length. There are then $2n+2$ elements, so
eval(x,i) makes sense for $0 <= i <=   2n+1$
*/
class CubicCOneSpline
{
	double xlow,xhigh,stepsize;
	int steps;

	double phi(double x)
	{
		if(0.0 <= x && x <= 1.0) return x*x*(3.0-2.0*x);
		else if( 1.0 <=x && x <=2.0) return (2.0-x)*(2.0-x)*(2.0*x-1.0);
		else return 0.0;
	}
	double psi(double x)
	{
		if(0.0 <=x && x <= 1.0) return x*x*(x-1.0);
		else if(1.0 <= x && x<= 2.0) return (2.0-x)*(2.0-x)*(x-1.0);
		else return 0.0;
	}

	public:
	CubicCOneSpline(double low, double high, int n)
	{
		assert( high > low && n >=1);
		xlow=low;
		xhigh=high;
		steps=n;
		stepsize=(xhigh-xlow)/steps;
	}
	double value(double t, int i);

};
double CubicCOneSpline::value(double t, int i)
{
  if( 0 <= i && i <= steps)
	 {
		//use phi
	  if( 1 <=i && i < steps) return phi( (t- (xlow+(i-1)*stepsize))/stepsize);
		if(i==0)
	{
	  t= (t-xlow)/stepsize;
	  if(0.0 <= t && t <=1.0) return phi(1.0+t);
	  else return 0.0;
	}
		if(i==steps)
	{
	  t = (t - (xlow+(steps-1)*stepsize)) /stepsize;
	  if( 0.0 <=t && t <= 1.0) return phi(t);
	  else return 0.0;
	}
	 }
  if( steps+1 <=i && i <= 2*steps+1)
	 {
		i=i-(steps+1);
		//use psi
	  if( 1 <=i && i < steps) return psi( (t- (xlow+(i-1)*stepsize))/stepsize);
		if(i==0)
	{
	  t= (t-xlow)/stepsize;
	  if(0.0 <= t && t <=1.0) return psi(1.0+t);
	  else return 0.0;
	}
		if(i==steps)
	{
	  t = (t - (xlow+(steps-1)*stepsize)) /stepsize;
	  if( 0.0 <=t && t <= 1.0) return psi(t);
	  else return 0.0;
	}


	 }
  return 0.0;
}

/**********
Piecewise linear ($C^0$) spline elements
Same deal --- divide the interval $[a,b]$ into
$n$ subintervals and define the $i$ spline element
to be a "tent function" supported on two of the subintervals,
or "half tent functions" on the left and right sub intervals.
There's n+1 spline elements, so eval(i,x)
makes sense for $0 \le i \le n$

******************/
class LinearCZeroSpline
{
  double xlow,xhigh,stepsize;
  int steps;

  double phi(double x)
	 {
		if(0.0 <= x && x <= 1.0) return x;
		else if( 1.0 <=x && x <=2.0) return 2.0-x;
		else return 0.0;
	 }

public:
  LinearCZeroSpline(double low, double high, int n)
    {
      assert( high > low && n >=1);
      xlow=low;
      xhigh=high;
      steps=n;
      stepsize=(xhigh-xlow)/steps;
    }
  double value(double t, int i);

};
double LinearCZeroSpline::value(double t, int i)
{
  if( 0 <= i && i <= steps)
    {
      //use phi
	  if( 1 <=i && i < steps) return phi( (t- (xlow+(i-1)*stepsize))/stepsize);
      if(i==0)
	{
	  t= (t-xlow)/stepsize;
	  if(0.0 <= t && t <=1.0) return phi(1.0+t);
	  else return 0.0;
	}
      if(i==steps)
	{
	  t = (t - (xlow+(steps-1)*stepsize)) /stepsize;
	  if( 0.0 <=t && t <= 1.0) return phi(t);
	  else return 0.0;
	}
    }
  
  return 0.0;
}


/*
main()
{
	int i;
	double x;
	CubicCOneSpline spline(0.0,1.0,4);
	cout << "Step is okay" <<endl;
	for( x=0.0; x <= 1.0; x += 0.02)
	{
		cout << x << " " << spline.value(x, 9)<< endl;
	}

	for(i=0;i<10;i++) cout << (spline.value(0.5,i)) << endl;

	for( x=0.0; x <= 1.0; x += 0.02)
	{
		cout << x << " " << spline.value(x, 9)<< endl;
	}
}
*/