readlib(bspline):			# to get the convert function
pairs := proc(x,y)
	local i;
	[x[i], y[i]] $ i=1..nops(x);
end:

param := proc(p,i)
	local X, s;

	X := indets(p);

	s := (i[2]-i[1])*X[1] - i[1];

	subs(t=s,p);
end:

#
#	Polysegment returns the CR cubic polynomial on four vertices.
#
CRsegment := proc(x,t)
	local i,CR;
	CR  :=  [-t^3/2 + t^2 - t/2, 3/2*t^3 - 5/2*t^2 + 1,
	  -3/2*t^3 + 2*t^2 + 1/2*t, 1/2*t^3 - 1/2*t^2]:

	sum(CR[i]*x[i], i=1..4);
end:

#
#	Compute piecewise cubic spline
#
CRsegments := proc(x,t)
	local segment, i, nosegs;

	if (nops(x) < 4) then
		ERROR(`I need at least four points`);
	fi;

	nosegs := nops(x) - 3;

	segment := [ [t < 0, 0], [t<=i,
		      param(CRsegment([x[i],x[i+1],x[i+2],x[i+3]],t),[i-1,i])]
		   $ i=1..nosegs, [t > nosegs, 0]];
	
	convert(segment,procedure);
end: