; Thiele's interpolation formula
; ==============================

(defun $thiele (u v)
  (loop for rh on
    (loop for w on (cdr u) ; we build rho with axes transposed
          ; first, building column 1 and 2
          for r = (list
                    (loop for a on u
                          for b on v
                          while (cdr a)
                          collect (div (sub (car a) (cadr a))
                                       (sub (car b) (cadr b)))) v)
          ; then building starting from column 3
                then (cons
                       (loop for a1 in u
                             for a2 in w
                             for y1 on (car r)
                             for y2 in (cdr (cadr r))
                             collect (add (div (sub a1 a2)
                                               (sub (car y1) (cadr y1)))
                                        y2)) r)
          ; finally get the row 1
                ; finally (return (loop for i in r ; i is already a 1-list
                ;                       for j = i then (cons (car i) j)
                ;                       finally (return j))))
                finally (return (nreverse (mapcar #'car r))))
        for k in (cdr u)
        with a = 0
        do (setf a (add a (div (sub 'x k)
                               (sub (caddr rh) (car rh)))))
        finally (return (add (car v)
                             (div (sub 'x (car u))
                                  (sub (cadr rh) a))))))