# April 22, 2000

#	Smoothest cubic spline interpolating { (t1,f1),...,(t_m+1,f_m+1) }

smo.spline _ function(fk,tk=seq(-1,1,length=length(fk))) {
    m _ length(fk)-1; if(length(tk) != m+1) {stop("tk wrong length")}
    hk_ diff(tk) ; if(any(hk) <= 0) {stop("mesh has ties or negative orientation")}
    hk2 _ hk**2; hk3 _ hk**3; z _ rep(0,m-1); o _ rep(1,m-1)
    if(m==1) {U _ matrix(c(4*hk,c(6,6)*hk^2,12*hk^3),2,2)}
    if(m>1) { U _ cbind(rbind( diag(4*hk),diag(6*hk2) ), rbind( diag(6*hk2),diag(12*hk3)))}
    if(m==1) warning("m = 1 not implemented")
    if(m==2) {V _ rbind(c(hk,c(2,1)*hk^2),c(1,-1,3*hk[1],0))}
    if(m>2) {
	Vt_ matrix(0,m-1,m); Vi_ ( row(Vt) == col(Vt)-1 )
	V11 _ cbind(diag(hk[-m]),z); V11[Vi] _ hk[-1]
	V12 _ cbind(diag(2*hk[-m]^2),z); V12[Vi] _ hk[-1]^2
	V21 _ cbind(diag(o),z); V21[Vi] _ -o
	V22 _ cbind(diag(3*hk[-m]),z)
	V _ cbind( rbind(V11,V12), rbind(V21,V22) ) }
    df _ ( fk[-1]-fk[-(m+1)] ) / hk
    w  _ c(diff(df),z); plot(range(tk),nicerange(fk,.2),type="n"); points(tk,fk,pch=16)
    x _ solve(U) %*% t(V) %*% solve( V %*% solve(U) %*% t(V) ) %*% w
    ck _ x[1:m]; dk _ x[-(1:m)]; bk _ df - hk*ck - hk^2*dk; ak _ fk[1:m]
    crt _ sum( 4*hk*ck^2 + 12*hk2*ck*dk + 12*hk^3*dk^2 ); tot_0
    for( k in 1:m ) { t_seq(tk[k],tk[k+1],l=101)-tk[k]
	y _ ak[k] + bk[k]*t + ck[k]*t^2 + dk[k]*t^3; lines(t+tk[k],y)
	tot _ tot + sum(diff(diff(y))^2)/diff(t[1:2])^3 }
    print(c(crt,tot))
      print("r=0"); lhs_ak+hk*bk+hk2*ck+hk3*dk; rhs_fk[-1]; print(cbind(lhs,rhs))
      print("r=1"); lhs_(bk+2*hk*ck+3*hk2*dk)[-m]; rhs_bk[-1]; print(cbind(lhs,rhs))
      print("r=2"); lhs_(2*ck+6*hk*dk)[-m]; rhs_2*ck[-1]; print(cbind(lhs,rhs))
browser()
}