; Copyright (c) 2022 Thomas Baruchel
; 
; Permission is hereby granted, free of charge, to any person obtaining a copy
; of this software and associated documentation files (the "Software"), to deal
; in the Software without restriction, including without limitation the rights
; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
; copies of the Software, and to permit persons to whom the Software is
; furnished to do so, subject to the following conditions:
; 
; The above copyright notice and this permission notice shall be included in
; all copies or substantial portions of the Software.
; 
; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
; SOFTWARE.
;
; Installation
; ------------
; The functions can be used with or without compilation:
;   * without compilation:
;         load("thiele.lisp")$
;   * with compilation (must be compiled only once):
;         :lisp  (compile-file "thiele.lisp");
;     look for the compiled file like "thiele.o" or "thiele.fasl"
;     and from now on:
;         load("thiele.fasl")$

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

(defun $thiele (u v)
  (loop for rh on
    (loop for w on (cdr u) ; we build rho with axes transposed (and reversed)
          ; first, building column 2 and 1
          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))))))