##########################
# Poly interp (Newton's Divided Difference given the x and y values. 
# The y values are replaced by the differences (locally).
##########################
InterpNewton <- function(knot.x, knot.y){
   n <- length(knot.x)
   for (k in 1:(n-1)){
      knot.y[(k+1):n] <- (knot.y[(k+1):n] - knot.y[k])/(knot.x[(k+1):n] - knot.x[k])
   }
   knot.y
}
##########################
# Use Horner's rule z*(c[n]+z*(c[n-1]+z*(...))) to evaluate a polynomial given 
#	the coefficients and knots and point(s)
##########################
HornerN <- function (coef, knot.x, z) {
   n <- length(knot.x)
   polyv <- coef[n]*rep(1,length(z))    
   for (k in (n-1):1){
      polyv <- (z - knot.x[k])*polyv +coef[k]
   }
   polyv
}