(defun rk-4 (f ynot start end num)
  "
;; Description: implements Runge-Kutta four 

	(rk-4 f ynot start end num)

;; Examples:
 (mprint 
   (bind-columns
     (rseq 0 2 11)
     (mapcar (lambda(x) (- (square (+ x 1)) (* .5 (exp x)))) (rseq 0 2 11)) 
     (rk-4 (lambda(time y) (- y (square time) -1)) .5 0 2 10)
     )
   )
"
  (let* (
	 (time 0)
	 (h (/ (- end start) num))
	 (half-h (/ h 2))
	 (yvalues (list ynot))
	 (w ynot)
	 )
    (dotimes
	(i num)
      (let* (
	     (time+ (+ time half-h))
	     (k1 (* h (funcall f time w)))
	     (k2 (* h (funcall f time+ (+ w (* .5 k1)))))
	     (k3 (* h (funcall f time+ (+ w (* .5 k2)))))
	     (time+ (+ time+ half-h))
	     (k4 (* h (funcall f time+ (+ w k3))))
	     (del_w (/ (+ k1 (* 2 (+ k2 k3)) k4) 6.0))
	     )
	(setq 
	 time (+ time h)
	 w (+ w del_w)
	 yvalues (append yvalues (list w))
	 )
	)
      )
    yvalues
    )
  )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; I adapt my routine 
;;
;; 			integrate-simpson-recursive
;;
;; so that it updates a global variable that I can use to graph (the boundaries
;; of each panel).
;;
;; As mentioned in the text (p. 216), the algorithm is a lot simpler if 
;; implemented in a recursive programming language, like lisp. Here it is!
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun integrate-simpson-recursive (f a b eps &key int-last verbose)
  "
;; Examples:
 (func f(x) 100/x^2*sin(10/x))
 (integrate-simpson-recursive #'f 1 8 1e-2)
 (integrate-simpson-recursive #'f 1 3 1e-4)
"
  (let* (
	 (xmean (mean (list a b)))
	 (int-l (simple-simpson f a xmean))
	 (int-r (simple-simpson f xmean b))
	 (int (+ int-l int-r))
	 )
    ;; here's the global....
    (setq xs (append xs (list a b)))

    ;; after the first trial, the "grosser" integral (int-last) will be passed
    ;; in (why recalculate?), however, if not, then we need to calculate it:
    (if (not int-last)
	(setq int-last (simple-simpson f a b))
      )
    
    ;; If verbose, printout the results so far:
    (if verbose (print (list int int-last)))

    ;; If our tolerance is met,
    (if (< (abs (- int int-last)) (* 15 eps)) 
	;; bail on that subinterval -- we're done!
	;; But when we do bail, perform a final extrapolation:
	(/ (- (* 16 int) int-last) 15)
      ;; otherwise, recurse, splitting the error in halves:
      (let ((new-eps (/ eps 2)))
	(+ (integrate-simpson-recursive 
	    f a xmean new-eps :int-last int-l :verbose verbose)
	   (integrate-simpson-recursive 
	    f xmean b new-eps :int-last int-r :verbose verbose)
	   )
	)
      )
    )
  )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; 			Simple-Simpson:
;;
;; Simpson's rule on a single panel
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun simple-simpson (f a b)
  (let (
	(f0 (funcall f a))
	(f1 (funcall f (mean (list a b))))
	(f2 (funcall f b))
	)
    (* 
     (/ (- b a) 6.0) ;; h/3 -- since 2h=b-a
     (+ f0 (* 4.0 f1) f2)
     )
    )
  )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; 			Euler's method
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun euler-method (f alpha a b n)
  "
;; Examples:
 (euler-method (lambda(time y) (- y (square time) -1)) .5 0 2 10)
 (euler-method (lambda(time y) (- (* time (exp (* 3 time))) (* 2 y))) 0 0 1 2)
 (euler-method (lambda(time y) (- (/ y time) (square (/ y time)))) 1 1 2 10)
"
  (let* (
	 (h (/ (- b a) n))
	 (time a)
	 (current alpha)
	 (result (list current))
	 )
    (dotimes
	(i n)
      (setq 
       current (+ current (* h (funcall f time current)))
       result (append result (list current))
       time (+ time h)
       )
      )
    result
    )
  )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;   Define the Taylor methods of orders 2 and 3
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun taylor_2(alpha a b n)
  "
Define function f before using this. It will compute the partials necessary to
arrive at the Taylor method of order 2.

;; Examples:
 (defun f(time y) (- y (square time) -1))
 (taylor_2 .5 0 2 10)
 (defun f(time y) (- (* time (exp (* 3 time))) (* 2 y)))
 (taylor_2  0 0 1 2)
 (defun f(time y) (- (/ y time) (square (/ y time))))
 (taylor_2 1 1 2 10)
"
  (let* (
	 (h (/ (- b a) n))
	 (time a)
	 (current alpha)
	 (result (list current))
	 )
    (derfunc fy f y)
    (derfunc ft f time)
    ;; (print (gle f))
    (derfunc fy f y)
    ;;     (print (gle fy))
    ;;     (print (gle ft))
    ;;     (print 'hi)
    (dotimes
	(i n)
      (let ((fval (f time current)))
	(setq 
	 current (+ current 
		    (* h 
		       (+ 
			fval
			(* (/ h 2) (+
				    (ft time current)
				    (* fval (fy time current))
				    )
			   )
			)
		       )
		    )
	 result (append result (list current))
	 time (+ time h)
	 )
	)
      )
    result
    )
  )
(defun taylor_3(alpha a b n)
  "
Define function f before using this. It will compute the partials necessary to
arrive at the Taylor method of order 3.

;; Examples:
 (defun f(time y) (- y (square time) -1))
 (taylor_3 .5 0 2 10)
 (defun f(time y) (- (* time (exp (* 3 time))) (* 2 y)))
 (taylor_3  0 0 1 2)
 (defun f(time y) (- (/ y time) (square (/ y time))))
 (taylor_3 1 1 2 10)
"
  (let* (
	 (h (/ (- b a) n))
	 (time a)
	 (current alpha)
	 (result (list current))
	 )
    (derfunc fy f y)
    (derfunc ft f time)
    (derfunc fyy fy y)
    (derfunc fyt fy time)
    (derfunc ftt ft time)
    (dotimes
	(i n)
      (let (
	    (fval (f time current))
	    (ftval (ft time current))
	    (fyval (fy time current))
	    (fttval (ftt time current))
	    (fytval (fyt time current))
	    (fyyval (fyy time current))
	    )
	(setq 
	 current (+ current 
		    (* h 
		       (+ 
			fval
			(* (/ h 2) 
			   (+
			    (ft time current)
			    (* fval (fy time current))
			    (* (/ h 3)
			       (+
				fttval
				(* fyval ftval)
				(* fval 
				   (+ (* 2 fytval)
				      (square fyval)
				      (* fval fyyval)
				      )
				   )
				)
			       )
			    )
			   )
			)
		       )
		    )
	 result (append result (list current))
	 time (+ time h)
	 )
	)
      )
    result
    )
  )

(defun taylor_4(alpha a b n)
  "
Define function f before using this. It will compute the partials necessary to
arrive at the Taylor method of order 4.

;; Examples:
 (defun f(time y) (- y (square time) -1))
 (mprint (bind-columns (rseq 0 2 11) (taylor_4 .5 0 2 10)))
 (defun f(time y) (- (* time (exp (* 3 time))) (* 2 y)))
 (taylor_4  0 0 1 2)
 (defun f(time y) (- (/ y time) (square (/ y time))))
 (taylor_4 1 1 2 10)
"
  (let* (
	 (h (/ (- b a) n))
	 (time a)
	 (current alpha)
	 (result (list current))
	 )
    (derfunc fy f y)
    (derfunc ft f time)
    (derfunc fyy fy y)
    (derfunc fyt fy time)
    (derfunc ftt ft time)
    (derfunc fyyy fyy y)
    (derfunc fyyt fyy time)
    (derfunc fytt fyt time)
    (derfunc fttt ftt time)
    (defun block1(time fval fyval fytval fyyval) 
      (+
       (square fyval)
       (* 2 fytval)
       (* fval fyyval)
       )
      )
    (defun block2(time y) 
      (+
       (* 2 (fy time y) (fyt time y))
       (* 2 (fytt  time y))
       (* (ft time y) (fyy time y))
       (* (f time y) (fyyt time y))
       )
      )
    (defun block3(time y) 
      (+
       (* 2 (fy time y) (fyy time y))
       (* 2 (fyyt  time y))
       (* (fy time y) (fyy time y))
       (* (f time y) (fyyy time y))
       )
      )
    (dotimes
	(i n)
      (let* (
	     (fval (f time current))
	     (ftval (ft time current))
	     (fyval (fy time current))
	     (fttval (ftt time current))
	     (fytval (fyt time current))
	     (fyyval (fyy time current))
	     (fyyyval (fyyy time current))
	     (fyytval (fyyt time current))
	     (fyttval (fytt time current))
	     (ftttval (fttt time current))
	     (block1 (block1 time fval fyval fytval fyyval) )
	     (block2 (block2 time current))
	     (block3 (block3 time current))
	     )
	(setq 
	 current (+ current 
		    (* h 
		       (+ 
			fval
			(* (/ h 2) 
			   (+
			    ftval
			    (* fval fyval)
			    (* (/ h 3)
			       (+
				(+
				 fttval
				 (* fyval ftval)
				 (* fval 
				    block1
				    )
				 )
				(* (/ h 4)
				   (+
				    ftttval
				    (+ (* fytval ftval) (* fyval fttval))
				    (* ftval block1)
				    (* fval block2)
				    (*
				     fval
				     (+
				      fyttval
				      (+ (* fyyval ftval) (* fyval fytval))
				      (* fyval block1)
				      (* fval block3)
				      )
				     )
				    )
				   )
				)
			       )
			    )
			   )
			)
		       )
		    )
	 result (append result (list current))
	 time (+ time h)
	 )
	)
      )
    result
    )
  )

(defun ext-round(round n factor)
  "
;; Examples:
 (mprint (ext-round 1 4 2))
 (mprint (ext-round 2 4 2))
 (mprint (ext-round 3 4 2))
"
  (let* (
	 (factor (^ factor round))
	 (n (- n round))
	 (mat (bind-columns 
	       (- (id n))
	       (if (> n 1) (combine (repeat 0 (- n 1)) factor)
		 (list factor))))
	 )
    (dotimes (i (- n 1))
      (setf (aref mat i (+ 1 i)) factor)
      )
    (/ mat (- factor 1))
    )
  )

(defun forward-difference(x y)
  "
;; Examples:
 (forward-difference (list .2 .4) 
		     (list .9798652 .917771))
"
  (/ 
   (apply #'- y)
   (apply #'- x)
   )
  )

(defun no-quotes(string)
  "
;; Examples:
 (no-quotes \"hi\")
"
  (format t "~a" string)
  )

(defun simpson-1-4-2-4-2-step-dance(n)
  "
;; Examples:
 (simpson-1-4-2-4-2-step-dance 3)
 (simpson-1-4-2-4-2-step-dance 5)
 (simpson-1-4-2-4-2-step-dance 4)
"
  (if (oddp n)
      (+
       (combine -1 (repeat 0 (- n 2)) -1)
       (mapcar (lambda(i) (^ -1 i)) (iseq 1 n))
       (repeat 3 n)
       )
    (error "Simpson's rule requires an odd number of points")
    )
  )

(defun trapezoid-1-2-2-step-dance(n)
  "
;; Examples:
 (trapezoid-1-2-2-step-dance 3)
 (trapezoid-1-2-2-step-dance 5)
 (trapezoid-1-2-2-step-dance 4)
"
  (+
   (combine -1 (repeat 0 (- n 2)) -1)
   (repeat 2 n)
   )
  )

(defun horner(x coefs &key (x0 0)) 
  "
;; Examples:
 (horner 1 (list 1))
 (horner 1 (list 1 1))
 (horner 1 (list 1 1 1))
 (horner 2 (list 1 1 1))
 (horner 2 (list 1 2 4))
 (horner 2 (list 1 2 4) :x0 2)
 (horner (list 2 3 4) (list 1 2 4) :x0 2)
"
  (let* (
	 (xdiff (- x x0))
	 (n (length coefs))
	 (result (cond
		  ((= 0 n)
		   (error "Sorry: you must have the coefficients of a polynomial!"))
		  ((= 1 n)
		   (first coefs)
		   )
		  (t
		   (final coefs)
		   )
		  )
		 )
	 )
    ;; use Horner's method to compute the value(s):
    (mapcar
     (lambda (fore)
       (setq result
	     (+ fore (* xdiff result)))
       )
     (reverse (butlast coefs))
     )
    result
    )
  )

(defun hermite(knots f derivs)
  "
;; Examples:
 ;; 3.3 Exercise 1a, 2a
 
 (func f(x) x*ln(x))
 (derfunc fx f x)
 (fx 8.3)
 (fx 8.6)
 (hermite (list 8.3 8.3 8.6 8.6) #'f (mapcar #'fx (list 8.3 8.6)))
 (hermite (list 1.3 1.3 1.6 1.6 1.9 1.9) (lambda(x) (bessel-j 0 x))
    (list -.5220232 -.5698959 -.5811571)
    )
"
  (let* (
	 (n (length knots))
	 (mat (zeros n n))
	 (fvals (if (functionp f) (mapcar f knots) f))
	 )
    (setf (select mat (iseq n) 0) (bind-columns fvals))
    (dolist
	(i (iseq 1 (- n 1)))
      (dolist
	  (j (iseq 1 i))
	(let* (
	       (diff (- (elt knots i) (elt knots (- i j))))
	       (fa (if (= diff 0) 
		       ;; use the derivative value given:
		       (pop derivs) 
		     ;; use the derivative approximation given by the finite difference:
		     (/ (-
			 (aref mat i (- j 1))
			 (aref mat (- i 1) (- j 1))
			 )
			diff)))
	       )
	  (setf (aref mat i j) fa)
	  )
	)
      )
    (mprint (bind-columns knots mat))
    mat
    )
  )

(defun divided-difference(x knots f)
  "
;; Examples:
 ;; 3.1 Example 5:
 (divided-difference 1.5 (list 1 1.3 1.6 1.9 2.2) (lambda(x) (bessel-j 0 x)))
 ;; 3.1 Example 6:
 (divided-difference 2.1 (list 2 2.2 2.3) (lambda(x) (ln x)))
 ;; 3.1 #1c
 (divided-difference .45 (list 0 .6 .9) (lambda(x) (ln (+ x 1))))
 (divided-difference -1/3 (list -.75 -.5 -.25 0) 
    (list -.0718125 -.02475 .3349375 1.101))
"
  (let* (
	 (n (length knots))
	 (mat (zeros n n))
	 (vals (if (listp f) f (mapcar f knots)))
	 )
    (setf (select mat (iseq n) 0) (bind-columns vals))
    (dolist
	(i (iseq 1 (- n 1)))
      (dolist
	  (j (iseq 1 i))
	(setf (aref mat i j)
	      (/ (-
		  (aref mat i (- j 1))
		  (aref mat (- i 1) (- j 1))
		  )
		 (- (elt knots i) (elt knots (- i j)))
		 )
	      )
	)
      )
    (mprint (bind-columns knots mat))
    mat
    )
  )

(defun bezier-interpolation(knots &key diffs (factor 3))
  "
;; Description:

	Function bezier-interpolation (knots &key diffs (factor 3))
creates a function which interpolates using Bezier polynomials to join a list
of knots which the user supplies. Factor is a factor of 3 that Bezier
polynomials use (an alternative is 1). If diffs is true, then the middle two
knots of a set of four is assumed to be the differences from the first and
fourth knots (see Burden and Faires, \"Numerical Analysis\"). 

	This is the algorithm as given in Burden and Faires, \"Numerical
Analysis\".

;; Examples:
 (def knots '((0 0) (.5 .5) (2 -1) (1 0)))
 (bezier-interpolation knots)
"
  (let (result)
    (do
     ;; start with a copy of knots, then shift by three to get the next foursome:
     ((z (copy-list knots) (rest (rest (rest z)))))
     
     ;; do this until z winds up with fewer than four knot (presumably zero!):
     ((< (length z) 4))

     (let* (
	    ;; I prep the values of the knots, both interpolating and control:
	    ;; control knots sometimes require a little preconditioning. 
	    ;; v1 and v2 are the knots to interpolate:
	    (v1 (first z))
	    (v2 (fourth z))

	    ;; Here are PRELIMINARY values of v+ and v-:
	    (v+ (second z))
	    (v- (third z))
	    ;; To get results that agree with Figure 3.16, p. 160 of Edition 7,
	    ;; I need to take the given control knots and divide by 3 (the
	    ;; default value of factor). BF discuss this in their text, but
	    ;; it's not clear that their algorithm would give the right plots
	    ;; as is....

	    ;; v+ is the control knot to v1:
	    (v+ (+ v1 (/ (if diffs v+ (- v+ v1)) factor)))
	    ;; v- is the control knot to v2:
	    (v- (+ v2 (/ (if diffs (- v-) (- v- v2)) factor)))

	    ;; Step 2 of BF's algorithm
	    (p1 v1)
	    (p2 (* 3 (- v+ v1)))
	    (p3 (* 3 (+ v- v1 (- (* 2 v+)))))
	    (p4 (- v2 v1 (* 3 (- v- v+))))

	    ;; output is v for each curve
	    (v (list p1 p2 p3 p4))
	    )
       ;; (print (list v+ v-))

       ;; v is tacked on to the result
       (setq result (append result v))
       (print v)
       )
     )

    ;; here's the long list of coefficients:
    result
    )
  )

(defun bezier-points(result n)
  "
;; Examples:
 (bezier-points (bezier-interpolation '((0 0) (-1 1) (2 1) (1 0))))
"
  (let* (
	 (steps (rseq 0 1 n))
	 (all (list (repeat 1 n) steps (^ steps 2) (^ steps 3)))
	 stuff 
	 coefs
	 xs
	 ys
	 xdata
	 ydata
	 )
    (while 
	result
      (setq 
       coefs (select result (iseq 4))
       xcoefs (mapcar #'first coefs)
       ycoefs (mapcar #'second coefs)
       result (rest (rest (rest (rest result))))
       )
      (setq 
       xs (apply #'+ (* xcoefs all))
       ys (apply #'+(* ycoefs all))
       xdata (if xdata (combine xdata xs) xs)
       ydata (if ydata (combine ydata ys) ys)
       )
      )
    (list xdata ydata)
    )
  )      

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                 clamped-spline
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun clamped-spline (xx yy &key m1 m2)
  "
;; Description: function clamped-spline

                (xx yy &key m1 m2)

Computes a clamped spline to data (xx,yy). Clamped-spline requires (ordinarily)
derivative values at the endpoints. If you don't have a clue, clamped-spline
will use the values of the first order divided-differences. This should give a
good fit provided points are dense at the ends, and fairly linear there.

         xx: a list of x values
         yy: a corresponding list of y values

Optional:
         m1: left endpoint slope
         m2: right endpoint slope

;; Note: following BF, Algorithm 3.5, p. 148 (7th edition)

;; Examples:
 (clamped-spline '(0 1) '(0 1) :m1 1 :m2 1) ;; should be y=x

 (print \"This should be the cubic x^3:\")
 (clamped-spline '(0 1 2 3) '(0 1 8 27) :m1 0 :m2 36 )

 (title \"A Duck?\")
 (def a '(.9 1.3 1.9 2.1 2.6 3 3.9 4.4 4.7 5 6 7 8 9.2 10.5
         11.3 11.6 12 12.6 13 13.3))
 (def b '(1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3
         2.25 1.95 1.4 .9 .7 .6 .5 .4 .25))
 (clamped-spline a b :m1 .54 :m2 -0.552)
"
  (let* (
	 (x (simul-sort xx yy))
	 (a (second x))
	 (x (first x))

	 ;; step 1
	 (h (- (rest x) (butlast x)))
	 (xdoublediff (if (rest h) (+ (rest h) (butlast h))))
	 (adiff (- (rest a) (butlast a)))

	 ;; if given, use the slopes; otherwise use a divided-difference approximation
	 (m1 (if m1 m1 (/ (- (second y) (first y)) (- (second x) (first x)))))
	 (m2 (if m2 m2 (/
			(- (second (reverse y)) (first (reverse y)))
			(- (second (reverse x)) (first (reverse x))))))

	 (n+1 (length x))
	 (n (1- n+1))
	 (n-1 (1- n))

	 (alpha (repeat 0 n+1))
	 (b (repeat 0 n))
	 (c (repeat 0 n+1))
	 (d (repeat 0 n))
	 (l (repeat 0 n+1))
	 (z (repeat 0 n+1))
	 (mu (repeat 0 n))
	 )

    ;; Step 2
    (setf 
     (elt alpha 0) (* 3 (- (/ (elt adiff 0) (elt h 0)) m1))
     (elt alpha n) (* 3 (- m2 (/ (elt adiff n-1) (elt h n-1))))

     ;; Step 4
     (elt l 0) (* 2 (elt h 0))
     (elt mu 0) .5
     (elt z 0) (/ (elt alpha 0) (elt l 0))
     )     

    ;; Step 3 and 5
    (mapcar
     (lambda(i)
       (let (
	     (i-1 (- i 1))
	     (i+1 (+ i 1))
	     )
	 (setf 
	  (elt alpha i) (* 3 (-
			      (/ (elt adiff i) (elt h i))
			      (/ (elt adiff i-1) (elt h i-1))
			      )
			   )
	  (elt l i) (- (* 2 (elt xdoublediff i-1))
		       (* (elt h i-1)
			  (elt mu i-1)
			  )
		       )
	  (elt mu i) (/ (elt h i) (elt l i))
	  (elt z i) (/ (- (elt alpha i)
			  (* (elt h i-1)
			     (elt z i-1)
			     )
			  )
		       (elt l i)
		       )
	  )
	 )
       )
     (if (> n-1 0) (iseq 1 n-1))
     )

    ;; Step 6
    (setf 
     (elt l  n) (* (elt h n-1)
		   (- 2 (elt mu n-1)))
     (elt z  n) (/ (- (elt alpha n)
		      (* (elt h n-1)
			 (elt z n-1)
			 )
		      )
		   (elt l n)
		   )
     (elt c  n) (elt z n)
     )

    ;; Step 7
    (mapcar
     (lambda(j)
       (let ((j+1 (+ j 1)))
	 (setf 
	  (elt c j) (- (elt z j) (* (elt mu j) (elt c j+1))) 
	  (elt b j) (- (/ (elt adiff j) (elt h j))
		       (* (elt h j)
			  (/ 
			   (+ (elt c j+1) 
			      (* 2 (elt c j))
			      )
			   3)
			  )
		       )
	  (elt d j) (/
		     (- (elt c j+1) 
			(elt c j)
			)
		     (* 3 (elt h j)))
	  )
	 )
       )
     (iseq n-1 0)
     )
    (list 
     (butlast x)
     (transpose (list (butlast a) b (butlast c) d))
     )
    )
  )  

(defun cubic-spline-points(x result)
  "
;; Examples:
 (cubic-spline-points 1.2 (clamped-spline '(0 1 2 3) '(0 1 8 27) :m1 0 :m2 36 ))
"
  (let* (
	 (xs (first result))
	 (order (where xs x #'<))
	 (n (if order (final order) 0))
	 (coefs (elt (second result) n))
	 (xdiff (- x (elt xs n)))
	 )
    ;; use Horner's method to compute the value:
    (+ (elt coefs 0)
       (*
	xdiff
	(+ (elt coefs 1)
	   (*
	    xdiff
	    (+ (elt coefs 2)
	       (*
		xdiff
		(elt coefs 3)
		)
	       )
	    )
	   )
	)
       )
    )
  )      

(defun cubic-spline-prime-points(x result)
  "
;; Examples:
 (cubic-spline-points 1.2 (clamped-spline '(0 1 2 3) '(0 1 8 27) :m1 0 :m2 36 ))
"
  (let* (
	 (xs (first result))
	 (order (where xs x #'<))
	 (n (if order (final order) 0))
	 (coefs (elt (second result) n))
	 (xdiff (- x (elt xs n)))
	 )
    ;; use Horner's method to compute the value:
    (+ (elt coefs 1)
       (*
	xdiff
	(+ (* 
	    2 
	    (elt coefs 2))
	   (*
	    3
	    xdiff
	    (elt coefs 3)
	    )
	   )
	)
       )
    )
  ) 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                 natural-spline
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun natural-spline (xx yy)
  "
;; Description: function natural-spline

                (xx yy)

Computes a natural spline to data (xx,yy). Natural-spline requires (ordinarily)
derivative values at the endpoints. If you don't have a clue, natural-spline
will use the values of the first order divided-differences. This should give a
good fit provided points are dense at the ends, and fairly linear there.

         xx: a list of x values
         yy: a corresponding list of y values

;; Note: following BF, Algorithm 3.4, p. 146 (7th edition)

;; Examples:
 (natural-spline '(0 1) '(0 1)) ;; should be y=x

 (print \"This should be the cubic x^3:\")
 (natural-spline '(0 1 2 3) '(0 1 8 27))

 (title \"A Duck?\")
 (def a '(.9 1.3 1.9 2.1 2.6 3 3.9 4.4 4.7 5 6 7 8 9.2 10.5
         11.3 11.6 12 12.6 13 13.3))
 (def b '(1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3
         2.25 1.95 1.4 .9 .7 .6 .5 .4 .25))
 (natural-spline a b)
"
  (let* (
	 (x (simul-sort xx yy))
	 (a (second x))
	 (x (first x))
	 
	 ;; step 1
	 (h (- (rest x) (butlast x)))
	 (xdoublediff (if (rest h) (+ (rest h) (butlast h))))
	 (adiff (- (rest a) (butlast a)))
	 
	 (n+1 (length x))
	 (n (1- n+1))
	 (n-1 (1- n))
	 
	 (alpha (repeat 0 n+1))
	 (b (repeat 0 n))
	 (c (repeat 0 n+1))
	 (d (repeat 0 n))
	 (l (repeat 0 n+1))
	 (z (repeat 0 n+1))
	 (mu (repeat 0 n))
	 )
    
    ;; Step 3
    (setf 
     (elt l 0) 1
     (elt mu 0) 0
     (elt z 0) 0
     )     
    
    ;; Step 2 and 4
    (mapcar
     (lambda(i)
       (let (
	     (i-1 (- i 1))
	     (i+1 (+ i 1))
	     )
	 (setf 
	  (elt alpha i) (* 3 (-
			      (/ (elt adiff i) (elt h i))
			      (/ (elt adiff i-1) (elt h i-1))
			      )
			   )
	  (elt l i) (- (* 2 (elt xdoublediff i-1))
		       (* (elt h i-1)
			  (elt mu i-1)
			  )
		       )
	  (elt mu i) (/ (elt h i) (elt l i))
	  (elt z i) (/ (- (elt alpha i)
			  (* (elt h i-1)
			     (elt z i-1)
			     )
			  )
		       (elt l i)
		       )
	  )
	 )
       )
     (if (> n-1 0) (iseq 1 n-1))
     )
    
    ;; Step 6
    (setf 
     (elt l  n) 1
     (elt z  n) 0
     (elt c  n) 0
     )
    
    ;; Step 7
    (mapcar
     (lambda(j)
       (let ((j+1 (+ j 1)))
	 (setf 
	  (elt c j) (- (elt z j) (* (elt mu j) (elt c j+1))) 
	  (elt b j) (- (/ (elt adiff j) (elt h j))
		       (* (elt h j)
			  (/ 
			   (+ (elt c j+1) 
			      (* 2 (elt c j))
			      )
			   3)
			  )
		       )
	  (elt d j) (/
		     (- (elt c j+1) 
			(elt c j)
			)
		     (* 3 (elt h j)))
	  )
	 )
       )
     (iseq n-1 0)
     )
    (list 
     (butlast x)
     (transpose (list (butlast a) b (butlast c) d))
     )
    )
  )  

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Deriving the combination of Taylor series for the points 
;;    f(x0), f(x0+h), f(x0+2h), f(x0+3h), f(x0+4h)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Let's make it into a function, that will find the coefficients of any 
;; five-point formula:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun five-point(which &key (target (list 0 1 0 0 0)) (int 12) (error 5))
  "
;; Examples:

 ;; derives the coefficients of the centered-difference approximation to the
 ;; first derivative:
 (five-point (iseq -2 2))

 ;; finds the five-point approximation to the second derivative:
 (five-point (iseq -2 2) :target (list 0 0 1 0 0) :int 24 :error 6)
"
  (let* (
	 (A (matrix '(5 5)
		    (combine
		     (mapcar 
		      (lambda(n)
			(mapcar 
			 (lambda(m) (^ n m))
			 (iseq 5)
			 )
			)
		      which
		      )
		     )
		    )
	    )
	 (B (transpose A))

	 ;; solve for the linear combination of the Taylor series that gives
	 ;; only the first derivative term:
	 (soln (solve B target))
	 )
    (mprint B)
    ;; Here's the solution:
    (format t "The coefficients: ~{~,5f ~}~%" soln)
    
    ;; looks better if we round:
    (format 
     t 
     "rounded integer coefficients (assuming division by ~dh): ~{~d ~}~%" 
     int (round (* soln int)))

    ;; Now what's left of the error term?
    (format t "coefficient of the error (~dth derivative term): ~f~%" 
	    error
	    (-
	     (/ 
	      (dot 
	       soln 
	       (mapcar (lambda(x) (^ x error)) which)
	       )
	      (factorial error)
	      )
	     )
	    )
    soln
    )
  )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Let's do the same with the three-point formula:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun three-point(which &key (target (list 0 1 0)) (int 2) (error 3))
  "
;; Examples:

 ;; derives the coefficients of the forward-difference approximation to the
 ;; first derivative:
 (three-point (iseq 3))
 ;; backwards:
 (three-point (iseq -3))
 ;; centered difference
 (three-point (iseq -1 1))
"
  (let* (
	 (A (matrix '(3 3)
		    (combine
		     (mapcar 
		      (lambda(n)
			(mapcar 
			 (lambda(m) (^ n m))
			 (iseq 3)
			 )
			)
		      which
		      )
		     )
		    )
	    )
	 (B (transpose A))

	 ;; solve for the linear combination of the Taylor series that gives
	 ;; only the first derivative term:
	 (soln (solve B target))
	 )
    (mprint B)
    ;; Here's the solution:
    (format t "The coefficients: ~{~,5f ~}~%" soln)
    
    ;; looks better if we round:
    (format 
     t 
     "rounded integer coefficients (assuming division by ~dh): ~{~d ~}~%" 
     int (round (* soln int)))

    ;; Now what's left of the error term?
    (format t "coefficient of the error (~dth derivative term): ~f~%" 
	    error
	    (-
	     (/ 
	      (dot 
	       soln 
	       (mapcar (lambda(x) (^ x error)) which)
	       )
	      (factorial error)
	      )
	     )
	    )
    soln
    )
  )