Re: Polynomial fitting routines?

In article <5s5mo1laoh1dlnfjv3aobnogic83q0371k@xxxxxxx>,
Herman D. Knoble <SkipKnobleLESS@xxxxxxxxxxxxxxx> writes:
>
> On Sun, 27 Nov 2005 23:05:58 +0000 (UTC), kargl@xxxxxxxxxxxxxxxxxxxxxxxxxxxx (Steven G.
> Kargl) wrote:
>
> -|I'm looking for a to polynomial fitting routine at accepts a
> -|dataset and returns a fit to fairly high precision. Anyone
> -|have a pointer to a favorite routine?
>
> Do you mean Splines? See http://www.netlib.org/cgi-bin/search.pl
> and do a Search for: Splines.
>
> If you want do straight polynomial fit,
> search for: Polynomail least squares
>
> Also see:
> http://users.bigpond.net.au/amiller/lsq/spline5.f90
> or
> http://users.bigpond.net.au/amiller/lsq/fit_poly.f90
>
> That is search the site:
> http://users.bigpond.net.au/amiller/

I'm looking for an algorithm that will fit a polynomial to
log(x) for 0.25 <= x <= 2 with at least 64-bit precision.
I investigated NSWC Math LIb's pfit.f last night. pfit.f
worked well on a test on sin(x) and -pi/2 <= x <= pi/2, but
I got horrible results for log(x).

I've written an implementation of C99's logl(x) for FreeBSD
(and any other project that needs it and does not object
to a BSD license) where x is a C long double with a 64-bit
significand. My logl(x) works well over the entire range
of x except in the interval (0.36,2) where catastrophic
cancellation causes problems. FWIW, the algorithm filters
out exceptional cases, then splits x into f*2**n and I
write log(x) = n * log(2) + log(f) where f is in [1/2,1).
log(f) is evaluated via an 8th order minimax polynomial
approximation.

--
Steve
http://troutmask.apl.washington.edu/~kargl/
.