This commit is contained in:
Thomas Baruchel 2022-09-14 17:54:10 +02:00
parent b334b9424b
commit b6f8b489bd
6 changed files with 680 additions and 1 deletions

View File

@ -1,2 +1,15 @@
# numerical-routines
# Numerical routines for Maxima
This repository contains various scripts to be used with the (https://maxima.sourceforge.io)[Maxima] computer algebra system. The files are written either in Lisp or in Maxima's own scripting language.
Each script contains some documentation in the initial comments at the top of the file.
Each file can be loaded with `load("myscript.mac");` or `load("myscript.lisp");`
You can also put some of these files in your `~/.maxima/` directory, which will make you able to load them from any current location (even if the required file is not in your current directory).
Files written in Lisp can also be compiled; just type:
:lisp (compile-file "myscript.lisp");
then look for the newly compiled file with a name ending in `.o` or `.fasl` (according to the Lisp implementation used by your version of Maxima) in your current directory and load it instead of the original one.

160
carleman.lisp Normal file
View File

@ -0,0 +1,160 @@
; 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("carleman.lisp")$
; * with compilation (must be compiled only once):
; :lisp (compile-file "carleman.lisp");
; look for the compiled file like "convolution.o" and from now on:
; load("carleman.o")$
; Usage
; -----
; Two functions are provided here:
; * carleman(v)
; * carleman_diag(v)
; The most important is the second one.
;
; Both functions compute the Carleman matrix of a function given as
; a list of coefficients from its Taylor expansion at 0.
; For instance: carleman([0,2,1,0,0,0,0,0]); for f(x)=x^2+2x
;
; The size of resulting matrices will match the length of the initial list.
;
; The second function returns the diagonalized Carleman matrix as a list
; of three matrices whose dot product is the Carleman matrix.
; Compute the Carleman matrix for series whose coefficients are in v.
; Return a list of lists.
; The vector v contains Maxima objects, but the car '(mlist simp) should
; be removed before calling the function.
(defun carleman (v)
(let ((n (list-length v)))
(loop repeat n
with w = (reverse v)
for u = (cons 1 (make-list (1- n) :initial-element 0))
then (loop for x on w
with s = NIL
do (setf s
(cons
(addn
(loop for a in x
for b in u
collect (mul a b))) s))
finally (return s))
collect u)))
; Formula (4.17) - but there seems to be some misprint in the PDF and
; the sign "-" has been added here
(defun carleman-diag-left (m d)
(apply #'mapcar #'list
(loop for w in m
for i from 1
for y = (list (cdr w)) then (cons (nthcdr i w) y)
for d1 = 1 then (meval (list 'mtimes d1 d))
for r = (cdar m) then (cdr r) ; exactly the required number of 0's !!!
collect (loop with z = (list 1)
for u in y
for x = z
then (cons
(div
(addn
(loop for e in x
for f in u
collect (mul e f)))
(sub d1 d2)) x)
for d2 = (div d1 d) then (div d2 d)
finally (setf (cdr z) r)
(return x)))))
(defun carleman-diag-middle (m d)
(loop for NIL in m
for z = (cdar m) then (cdr z)
for q = NIL then (cons 0 q)
for d1 = 1 then (mul d1 d)
collect (append q (cons d1 z))))
; Formula (4.16) in "Continuous time evolution form iterated maps and
; Carleman linearization" (Gralewicz and Kowalski)
(defun carleman-diag-right (m d)
(loop for NIL in m
for d1 = 1 then (mul d1 d)
; transpose matrix m to z and iterate on rows of z
for z = (cdr (apply #'mapcar #'list m)) then (mapcar #'cdr (cdr z))
for q = NIL then (cons 0 q)
collect (loop with x = (list 1)
for y = x then (cdr y)
for u in z
for d2 = (mul d1 d) then (mul d2 d)
do (push (div
(addn
(loop for e in x
for f in u
collect (mul e f)))
(meval (sub d1 d2))) (cdr y))
finally (return (append q x)))))
;;; MAXIMA interface
;;; ================
; Return the Carleman matrix of a function whose Taylor expansion is
; given as a list of coefficients.
(defun $carleman (v)
(simplifya (cons '($matrix)
(mapcar #'(lambda (x) (simplify (cons '(mlist) x)))
(carleman (cdr v))))))
; Let M be the Carleman matrix of a function having 0 as a fixed point
; (ie. f(0)=0) and f'(0) not in {0, 1} ; now, V(M) is such
; that M = V^(-1) . L . V with L a diagonal matrix of eigenvalues.
;
; Return the diagonalized Carleman matrix of a function whose Taylor
; expansion is given as a list of coefficients.
; The first coefficient MUST be 0 (since f(0)=0 is a fixed point).
; The second coefficient MUST be some positive value different from 1.
;
; The function returns a list of three matrices whose product is the
; Carleman matrix.
(defun $carleman_diag (v)
(let ((m (carleman (cdr v)))
(d (caddr v)))
(simplifya (list '(mlist)
; left part
(simplifya
(cons '($matrix) (mapcar #'(lambda (x) (simplify (cons '(mlist) x)))
(carleman-diag-left m d))))
; diagonal matrix
(simplifya
(cons '($matrix) (mapcar #'(lambda (x) (simplify (cons '(mlist) x)))
(carleman-diag-middle m d))))
; right part
(simplifya
(cons '($matrix) (mapcar #'(lambda (x) (simplify (cons '(mlist) x)))
(carleman-diag-right m d))))))))

228
convolution.lisp Normal file
View File

@ -0,0 +1,228 @@
; Copyright (c) 2016 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.
;
; Functions for Maxima related to convolutions of series
; ======================================================
;
; Installation
; ------------
; The functions can be used with or without compilation:
; * without compilation:
; load("convolution.lisp")$
; * with compilation (must be compiled only once):
; :lisp (compile-file "convolution.lisp");
; look for the compiled file like "convolution.o" and from now on:
; load("convolution.o")$
;
; Functions
; ---------
; * Function 'recvec': detecting a recurrence relation (return a list of
; rational coefficients);
; * Function 'recvecn': detecting a recurrence relation (return a list of
; integer coefficients);
; * Function 'ggf': compute the generating function of a sequence
;
; Examples
; --------
; (%i5) recvec(makelist( 1/2^i, i, 12));
; (%o5) [- 2, 1]
; (%i6) recvec(makelist( 1/2^i+1, i, 12));
; (%o6) [2, - 3, 1]
; (%i7)
; (%i6) recvec(makelist( fib(i), i, 12));
; (%o6) [- 1, 1, 1]
; (%i7) recvec(makelist( 1/2^(12-i), i, 12));
; (%o7) [- 1/2, 1]
; (%i10) recvecn(makelist( 1/2^(12-i), i, 12));
; (%o10) [- 1, 2]
; (%i64) ggf(makelist( fib(i), i, 12));
; 1
; (%o64) - -----------
; 2 1
; x + x - 1
; (%i65) ggf(makelist( fib(i), i, 12),y);
; 1
; (%o65) - -----------
; 2 1
; y + y - 1
;;; (proclaim '(optimize (speed 3) (safety 0) (debug 0)))
; Compute the smallest (integer) coefficient for converting (by multiplication)
; a list of rational numbers to a list of integers; this is the LCM of all
; denominators.
(defmacro coeff-normalize-list-fractions (v)
`(reduce #'lcm (mapcar #'denominator ,v)))
; Convolution between two series (lists of coefficients); the final size is the
; size of the shortest list
(defun convolution (a b)
(loop
for NIL in a
for y in b
for z = (list (car b)) then (cons y z)
collect (loop
for i in a
for j in z
sum (* i j))))
; Convolution between one series and one polynomial; the final size is the
; size of the longest list (first argument).
; The polynomial (second argument) MUST have less coefficients than the series.
(defun convolution-poly (a b)
(loop
for NIL in a
for y = b then (if (cdr y) (cdr y) '(0))
for z = (list (car b)) then (cons (car y) z)
collect (loop
for i in a
for j in z
sum (* i j))))
; Compute the reciprocal of a series (list of coefficient); the first coefficient
; MUST not be zero.
(defun convolution-reciprocal (l)
(loop
for NIL in l
for m = (list (/ 1 (car l))) then
(cons (/ (-
(loop
for i in (cdr l)
for j in m
sum (* i j)))
(car l)) m)
finally (return (nreverse m))))
; Compute the minimal recurrence vector; returned coefficients are rational
; (though integer coefficients may often been returned, non integer ones can
; also been returned in some cases).
; Warning: an empty sequence returns (1)
(defun recurrence-vector-raw (v)
(loop named main
with z = (floor (/ (length v) 2))
with l = v
with q1 = '(0)
for q2 = '(1)
then (multiple-value-bind (m qq1)
(loop
for a on l
for b = q1 then (cons 0 b)
finally (return (values NIL NIL))
do (if (/= 0 (car a)) (return (values a b))))
(if qq1
(loop named add-and-compute-size
for k2 = q2 then (if (cdr k2) (cdr k2) '(0))
for k1 = qq1 then (if (cdr k1) (cdr k1) '(0))
for o = (cons (+ (/ (car q2) (car m)) (car qq1)) o)
then (cons (+ (/ (car k2) (car m)) (car k1)) o)
for ss from 1
do (if (not (or (cdr k1) (cdr k2)))
(loop
for oo on o
for i downfrom ss
do (if (/= 0 (car oo))
(if (> i z)
(return-from main NIL)
(progn
(setf q1 (cons 0 q2))
(setf l (cdr (convolution-reciprocal m)))
(return-from add-and-compute-size
(nreverse oo))))))))
(return-from main q2)))))
; Compute the minimal recurrence vector; returned coefficients are integers.
(defun recurrence-vector (v)
(let ((l (recurrence-vector-raw v)))
(if l
(mapcar #'(lambda (a) (* (coeff-normalize-list-fractions l) a)) l)
NIL)))
; Compute the generating function for a sequence if g.f is a rational function.
(defun ggf (v)
(let ((l (recurrence-vector-raw v)))
(if l
(let* ((s (nreverse
(loop
for x on (nreverse (convolution-poly v l))
do (if (/= 0 (car x)) (return x)))))
(c (lcm
(coeff-normalize-list-fractions l)
(coeff-normalize-list-fractions s))))
(list
(mapcar #'(lambda (a) (* c a)) s)
(mapcar #'(lambda (a) (* c a)) l)))
NIL)))
;;; MAXIMA interface
;;; ================
(defun from-maxima-list (l)
(mapcar #'(lambda(r)(if (and (listp r)(eq (caar r) 'rat))
(/ (second r)(third r))
r))
(cdr l)))
(defun to-maxima-list (l)
(cons '(mlist)
(mapcar #'(lambda (r)
(if (rationalp r)
(list '(rat)(numerator r)(denominator r)) r)) l)))
(defun to-maxima-polynomial (v x)
(if (cdr v)
(cons '(mplus simp)
(cons (car v)
(loop
for c in (cdr v)
for n from 1
when (/= 0 c)
collect
(if (= 1 c)
(if (= 1 n) x (list '(mexpt simp) x n))
(if (= 1 n)
(list '(mtimes simp) c x)
(list '(mtimes simp) c (list '(mexpt simp) x n)))))))
(car v)))
(defun to-maxima-ratfrac (f x)
(list '(mtimes simp)
(list '(mexpt simp) (to-maxima-polynomial (cadr f) x) -1)
(to-maxima-polynomial (car f) x)))
(defun $recvec (v)
(let* ((l (from-maxima-list v))
(g (if l (recurrence-vector-raw l) NIL)))
(if g
(to-maxima-list g)
NIL)))
(defun $recvecn (v)
(let* ((l (from-maxima-list v))
(g (if l (recurrence-vector l) NIL)))
(if g (cons '(mlist simp) g) NIL)))
(defun $ggf (v &optional (x (quote $x)))
(let* ((l (from-maxima-list v))
(g (if l (ggf l) NIL)))
(if g
(to-maxima-ratfrac g x)
NIL)))

146
pslq.mac Normal file
View File

@ -0,0 +1,146 @@
/*
Copyright (c) 2015 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.
*/
/*
An integer relation algorithm for finding integer relations between constants
(bfloat type being the intended type to be used by the function).
See: https://en.wikipedia.org/wiki/Integer_relation_algorithm
Examples
========
(%i4) fpprec: 64$
(%i5) pslq([bfloat(zeta(2)), %pi^2], 1b-32);
(%o5) [- 6, 1]
(%i78) fpprec: 64$ phi: bfloat((1+sqrt(5))/2)$
(%i79) algdep(phi,2,1b-32);
(%o79) x^2-x-1
(%i65) fpprec: 64$ zerobern: false$
(%i67) b:makelist(bfloat(bern(k)), k, 10);
(%o67) [- 5.0b-1,
1.666666666666666666666666666666666666666666666666666666666666667b-1,
- 3.333333333333333333333333333333333333333333333333333333333333333b-2,
7.575757575757575757575757575757575757575757575757575757575757576b-2,
- 2.531135531135531135531135531135531135531135531135531135531135531b-1,
1.166666666666666666666666666666666666666666666666666666666666667b0,
- 7.092156862745098039215686274509803921568627450980392156862745098b0,
5.497117794486215538847117794486215538847117794486215538847117795b1,
- 5.291242424242424242424242424242424242424242424242424242424242424b2,
6.192123188405797101449275362318840579710144927536231884057971015b3]
(%i68) makelist( algdep( b[k], 1, 1b-32), k, 10);
(%o68) [2 x + 1, 1 - 6 x, 30 x + 1, 66 x - 5, - 2730 x - 691, 7 - 6 x,
510 x + 3617, 43867 - 798 x, 330 x + 174611, 138 x - 854513]
*/
define_variable(pslqMaxItr, 256, integer)$
/* TODO:
known bugs:
(1) fpprec: 64$
pslq([bfloat(1/90),1],1b-32); Works
pslq([1,bfloat(1/90)],1b-32); Fails
*/
pslq(x, prec) := block([
n:length(x),
teps:prec*1.6b1,
gam: 1.2b0,
t, y, s, h, a, b, m,
done: false
],
/* Compute the vector s */
s:makelist(0.0b0, n),
s[n]: bfloat(abs( x[n] )),
t: x[n]*x[n],
for i:(n-1) step -1 thru 1 do (
t: bfloat(t + x[i]*x[i]),
s[i]: sqrt(t) ),
/* Normalize the vector x, s */
t: s[1],
y: makelist( bfloat(x[i]/t), i, n),
s: map(lambda([z], z/t), s),
/* Construct matrix H, a, b */
h: zeromatrix(n, n-1),
a: ident(n),
b: ident(n),
for i:1 thru n do
for j:1 thru min(i, n-1) do
h[i,j]: if i=j then s[j+1]/s[j] else -y[i]*y[j]/(s[j]*s[j+1]),
/* Reduce matrix H */
for i:2 thru n do
for j:(i-1) step -1 thru 1 do (
t: round( h[i,j]/h[j,j] ),
y[j]: y[j]+t*y[i],
for k:1 thru j do h[i,k]: h[i,k] - t *h[j,k],
for k:1 thru n do (
a[i,k]: a[i,k] - t * a[j,k],
b[k,j]: b[k,j] + t * b[k,i] ) ),
for itr:1 while (itr <= pslqMaxItr) and (not done)
do block([ mval:-1, g:gam ],
for i:1 thru n-1 do (
t: abs(g * h[i,i]),
if t > mval then ( mval:t, m:i ),
g: g*gam),
t: y[m], y[m]: y[m+1], y[m+1]: t,
a:rowswap(a, m, m+1),
h:rowswap(h, m, m+1),
b:columnswap(b,m, m+1),
if m < n-1 then block([t0: sqrt( h[m,m]^2 + h[m,m+1]^2), t1, t2, t3, t4],
t1: h[m,m]/ t0,
t2: h[m,m+1]/ t0,
for i:m thru n do (
t3: h[i,m],
t4: h[i,m+1],
h[i,m]: t1*t3 + t2*t4,
h[i,m+1]: t1*t4 - t2*t3)
),
for i:m+1 thru n do
for j:min(i-1,m+1) step -1 thru 1 do (
t: round( h[i,j]/h[j,j] ),
y[j]: y[j] + t * y[i],
for k:1 thru j do h[i,k]: h[i,k] - t*h[j,k],
for k:1 thru n do (
a[i,k]: a[i,k] - t * a[j,k],
b[k,j]: b[k,j] + t * b[k,i] )),
mval: 0,
for j: 1 thru n-1 do (
t: abs( h[j,j] ),
if t > mval then mval: t),
for i:1 thru n do (
t: abs( y[i] ),
if t < teps then (done: true, m: i))
),
if done then makelist(b[i,m],i,n) else makelist(0, n)
)$
algdep(x, n, prec) := block([k, v],
v: makelist(1, n+1), for k:n thru 1 step -1 do v[k]:v[k+1]*x,
v: pslq(v, prec),
sum( v[i] * ('x)^(n+1-i), i, 1, n+1)
)$

68
rinterp.mac Normal file
View File

@ -0,0 +1,68 @@
/*
Copyright (c) 2015 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.
Adapted from: http://axiom-wiki.newsynthesis.org/RationalInterpolation
*/
/* version 1 */
/*
rinterp(xlist, ylist, m, k) := block([collist, i, j, x, l, ol],
if length(xlist) # length(ylist)
then error("Different number of points and values."),
if length(xlist) # m+k+1
then error("Wrong number of points."),
collist: genmatrix(lambda([i,j],1), m+k+1, m+k+2),
for j:1 thru max(m,k) do
for i: 1 thru m+k+1 do
collist[i,j+1]: collist[i,j] * xlist[i],
for j:1 thru k+1 do
for i:1 thru m+k+1 do
collist[i,m+k+3-j]: -collist[i,k+2-j]*ylist[i],
x: makelist( gensym(), j, m+k+2 ),
ol: linsolve_params, linsolve_params: true,
l: linsolve(makelist(sum( x[j]*collist[i,j], j, 1, m+k+2 ), i, m+k+1), x),
for j:1 thru length(%rnum_list) do l:subst(1, %rnum_list[j], l),
linsolve_params: ol,
sum( rhs(l[j])*('x)^(j-1), j, 1, m+1 ) /
sum( rhs(l[m+1+j])*('x)^(j-1), j, 1, k+1 ));
*/
/* version 2 */
rinterp(xlist, ylist, m, k) := block([collist, i, j, x, l, ol],
if length(xlist) # length(ylist)
then error("Different number of points and values."),
if length(xlist) # m+k+1
then error("Wrong number of points."),
collist: genmatrix(lambda([i,j],1), m+k+2, m+k+1),
for j:1 thru max(m,k) do
for i: 1 thru m+k+1 do
collist[j+1,i]: collist[j,i] * xlist[i],
for j:1 thru k+1 do
for i:1 thru m+k+1 do
collist[m+k+3-j,i]: -collist[k+2-j,i]*ylist[i],
x: genmatrix(lambda([i,j],gensym()), 1, m+k+2),
ol: linsolve_params, linsolve_params: true,
l:linsolve(list_matrix_entries(x.collist), list_matrix_entries(x)),
for j:1 thru length(%rnum_list) do l:subst(1, %rnum_list[j], l),
linsolve_params: ol,
sum( rhs(l[j])*('x)^(j-1), j, 1, m+1 ) /
sum( rhs(l[m+1+j])*('x)^(j-1), j, 1, k+1 ));

64
thiele.mac Normal file
View File

@ -0,0 +1,64 @@
/*
Copyright (c) 2020 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.
*/
/*
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.
(%i18) a:makelist(0.05*i,i,0,31)$
(%i19) b:makelist(float(cos(0.05*i)),i,0,31)$
(%i20) keepfloat:true$
(%i21) subst(x=0.5, thiele(b,a))*3;
(%o21) 3.14159265357928
*/
/* Thiele's interpolation formula */
thiele(u, v) := block([rho:makelist(
makelist(v[i], length(v)-i+1),
i, length(v)), a:0],
for i:1 thru length(rho)-1
do rho[i][2]: (u[i]-u[i+1]) / (rho[i][1] - rho[i+1][1]),
for i:3 thru length(rho)
do (for j:1 thru length(rho)-i+1
do rho[j][i]: (u[j]-u[j+i-1])
/ (rho[j][i-1]-rho[j+1][i-1])
+ rho[j+1][i-2]),
rho: rho[1],
for i:length(rho) thru 3 step -1
do a: ratsimp(( 'x - u[i-1])/(rho[i]-rho[i-2]+a)),
ratsimp( v[1] + ( 'x - u[1] ) / (rho[2] + a) ))$