diff --git a/README.md b/README.md
index 5d81658..91714d4 100644
--- a/README.md
+++ b/README.md
@@ -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.
diff --git a/carleman.lisp b/carleman.lisp
new file mode 100644
index 0000000..6e6543a
--- /dev/null
+++ b/carleman.lisp
@@ -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))))))))
+  
diff --git a/convolution.lisp b/convolution.lisp
new file mode 100644
index 0000000..28193b2
--- /dev/null
+++ b/convolution.lisp
@@ -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)))
diff --git a/pslq.mac b/pslq.mac
new file mode 100644
index 0000000..049f73e
--- /dev/null
+++ b/pslq.mac
@@ -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)
+)$
diff --git a/rinterp.mac b/rinterp.mac
new file mode 100644
index 0000000..57325aa
--- /dev/null
+++ b/rinterp.mac
@@ -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 ));
diff --git a/thiele.mac b/thiele.mac
new file mode 100644
index 0000000..6124647
--- /dev/null
+++ b/thiele.mac
@@ -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) ))$