98 lines
4.1 KiB
Common Lisp
98 lines
4.1 KiB
Common Lisp
; 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))))))
|