; 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
; ==============================
; 
;   In the following example, we use Thiele's interpolation formula
;   for computing the inverse cos function cos^(-1) from 32 points
;   given by the cos function itself.
; 
;   We create two vectors:
;     * the vector 'a' contains regularly spaced points in [0..1.55]
;     * the vector 'b' contains cos(x) for each x in vector 'a'
; 
;   We set the keepfloat flag to true in order to prevent Maxima
;   to replace the values with rational approximations.
; 
;   Then we compute the inverse function by using thiele(b,a);
;   we evaluate this interpolation at 0.5 expecting some value
;   close to pi/3.
; 
; (%i1) a:makelist(0.05*i,i,0,31)$
; (%i2) b:makelist(float(cos(0.05*i)),i,0,31)$
; (%i3) keepfloat:true$
; (%i4) subst(x=0.5, thiele(b,a))*3;
; (%o4)                         3.14159265357928
;
;    The function may be used with exact values, but it is highly sensitive
;    to input values, and equally spaced abscissæ often yield some division
;    by zero error.
;
; (%i5) thiele([1,2,5,6], [10,12,11,13]);
;                                   2
;                                9 x  + 29 x - 238
; (%o5)                          -----------------
;                                    8 x - 28


(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))))))