[SciPy-dev] RE: divided differences

[Chuck Harris]

Hi all,

I have written fortran routines for real and complex interpolation using 
divided differences. The routines use the Steffensen form of evaluation
which zig-zags up the divided difference table in such a way as to 
minimize errors. These routines are essentially as accurate and fast as
the barycentric form of Lagrange interpolation, but without the
singularities. They are *much* more efficient than the Neville form
recommended in Numerical Recipes, and I suggest that the relevant pages
be razored out of that darned book.

f2py is used to generate the interface. What an amazing bit of software!
The complex routines can take real data and work just fine. I have written
both double and complex for efficiency. I intend to write python interfaces 
to call the appropriate routines, but haven't yet done so. The interface also
needs to convert arrays with one element to scalars on return --- maybe: this
could limit its usability as a subroutine, comments?

The routines are called:

make_table_double(xx,yy) -- xx,yy vectors, returns array dd of divided differences
make_table_complex(xx,yy) -- xx,yy vectors, returns array dd of divided differences
interpolate_double(dd,x) -- dd array, x vector, returns vector of interpolated values
interpolate_complex(dd,x) -- dd array, x vector, returns vector of interpolated values

In any case, is there an interest? If so, I suggest that these routines be included
in the interpolation directory, perhaps under a subdirectory DividedDifferences.
Perhaps there should also be a Splines subdirectory?

Chuck