polynomial interpolation with error estimate

A simple transcription, using hashed arrays, of polint
from numerical recipes fortran 77, 2nd. ed, 1992.
http://www.nrbook.com/a/bookfpdf.php

This approach is based on Neville's algorithm.

/* interpolate.mac */


/* polint transcribed from fortran 77 2nd. ed.
   of numerical recipes, ch. 3, p. 103.


   xa and ya are hashed arrays of length n
   c and d are hashed arrays of length n.
   For example, if n = 3, then we are given
   values of xa[j] and ya[j] for j = 1,2,3,
   and create c[j] and d[j] for j = 1,2,3,
   and the three data points (x,y) provided
   are used to create a quadratic which goes
   through the given three points, and then
   the quadratic is used to estimate y(x).

   polint(xa,ya,n,x) returns a list [y,dy],
   where y is an interpolated value for y(x),
   and dy is the estimated relative error.
                                */


polint(xa,ya,n,x) :=
block ([y,dy,i,m,ns:1,den,dif,dift,ho,hp,
        w,fails:false],

   dif : abs(x-xa[1]),
   for i thru n do
     (dift : abs(x-xa[i]),
     if dift < dif then
       (ns : i,
        dif : dift),
     c[i] : ya[i],
     d[i] : ya[i]),

   y : ya[ns],
   ns : ns - 1,
   for m thru n-1 do
     (for i thru n-m do
         (ho : xa[i] - x,
          hp : xa[i+m] - x,
          w : c[i+1] - d[i],
          den : ho - hp,
          if den = 0.0 then
            (fails:true,
             print ("  failure in polint "),
             return()),
          den : w/den,
          d[i] : hp*den,
          c[i] : ho*den),

      if fails then return(),

      if (2*ns < (n-m)) then
         dy : c[ns+1]
      else
         (dy : d[ns],
          ns : ns - 1),

      y : y + dy),

   remarray (c,d),

   if fails then [false,false]
   else [y,dy])$


/************ xpolint1  driver for polint  ***************/

/* transcribed from numerical recipes example book (fortran)
    2nd edition,p. 31. */

 xpolint1(nn) :=
  block ([ii,xv,yv,dyv  ],
    if nn >= 10 then
      (print (" choose nn < 10 "),
        return(false)),

    disp ("  sin(x)  for 0 < x <= %pi "),

    for ii thru nn do
      (xx[ii] : float(ii*%pi/nn),
       yy[ii] : sin(xx[ii])),

    /* based on the nn points, estimate the
      value of y(xv) at 10 points xv. */

    /* print("      xv   sin(xv)    yv      dyv  "), */
    printf(true,"  ~12a  ~12a  ~12a  ~12a  ~%","xv","sin(xv)","yv","dyv"),

    for ii thru 10 do
      ( xv : float((-0.05 + ii/10)*%pi),
        [yv,dyv] : polint(xx,yy,nn,xv),
        /* print("  ",xv,"  ",sin(xv),"  ",yv,"  ",dyv) */
        printf(true,"  ~12e  ~12e  ~12e  ~12e  ~%",xv,sin(xv),yv,dyv)),

    remarray(xx,yy),
    done)$

/*

(%i1) load(interpolate);

(%o1) "c:/work2/interpolate.mac"
(%i2) xpolint1(5);

"  sin(x)  for 0 < x <= %pi "

  xv                       sin(xv)               yv 
dyv
  1.5707963E-1  1.5643447E-1  1.5518926E-1  7.8254225E-2
   4.712389E-1   4.539905E-1  4.5341805E-1  1.3211752E-2
  7.8539816E-1  7.0710678E-1   7.074428E-1  -5.216948E-3
  1.0995574E+0  8.9100652E-1  8.9124652E-1  -3.048866E-3
  1.4137167E+0  9.8768834E-1  9.8748454E-1  2.3713402E-3
   1.727876E+0  9.8768834E-1  9.8748454E-1  2.3713402E-3
  2.0420352E+0  8.9100652E-1  8.9124652E-1  -3.048866E-3
  2.3561945E+0  7.0710678E-1   7.074428E-1  -5.216948E-3
  2.6703538E+0   4.539905E-1  4.5341805E-1  1.3211752E-2
   2.984513E+0  1.5643447E-1  1.5518926E-1  -5.216948E-3
(%o2) done
(%i3) xpolint1(9);

"  sin(x)  for 0 < x <= %pi "

   xv                      sin(xv)                yv                     dyv
  1.5707963E-1  1.5643447E-1   1.564342E-1   4.981823E-5
   4.712389E-1   4.539905E-1  4.5399055E-1  -3.279139E-6
  7.8539816E-1  7.0710678E-1  7.0710677E-1  5.9718307E-7
  1.0995574E+0  8.9100652E-1  8.9100653E-1  -1.604347E-7
  1.4137167E+0  9.8768834E-1  9.8768834E-1  3.6628832E-8
   1.727876E+0  9.8768834E-1  9.8768834E-1  3.6628832E-8
  2.0420352E+0  8.9100652E-1  8.9100653E-1  -1.604347E-7
  2.3561945E+0  7.0710678E-1  7.0710677E-1  5.9718307E-7
  2.6703538E+0   4.539905E-1  4.5399055E-1  -3.279139E-6
   2.984513E+0  1.5643447E-1   1.564342E-1  -2.969298E-6
(%o3) done



*/

Original fortran77 code:

      SUBROUTINE polint(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)pause 'failure in polint'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END
---------------------------------

Ted Woollett