Re: Polynomial fitting routines?

Steven G. Kargl wrote:

In article <quGdnWnhlu9pbhbeRVn-gA@xxxxxxxxxxx>,
	glen herrmannsfeldt <gah@xxxxxxxxxxxxxxxx> writes:

Steven G. Kargl wrote: 

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).


If I remember, it is not usual to write a polynomial in x, but
in (x-1) or some other rational function of x such as (x-1)/(x+1).

Indeed. After revisiting Abramowitz and Stegun's 4.1.26 and converting pfit.f to double precision, I can find a polynomial
over the range I'm interested. Unfortunately, it only produces
about O(1e-12) relative error on a few test values. This is
much too large for my purposes.

I think it would be usual to have a little more precision in the fit calculation than in the results of the fit...

The IBM OS/360 Fortran library, which is in the public domain
(I think it was before software copyright was allowed) uses:


 COMPUTE LOG((1+Z)/(1-Z)) BY MINIMAX
   APPROXIMATION OF THE FORM
     W+C1*W**3(W**2+C2+C3/
       (W**2+C4+C5/(W**2+C6)))

 with

C6     DC      X'C158FA4E0E40C0A5'    -0.5561109595943017E+1
C5     DC      X'C12A017578F548D1'    -0.2625356171124214E+1
C4     DC      X'C16F2A64DDFCC1FD'    -0.6947850100648906E+1
C3     DC      X'C38E5A1C55CEB1C4'    -0.2277631917769813E+4
C2     DC      X'422FC604E13C20FE'     0.4777351196020117E+2
C1     DC      X'3DDABB6C9F18C6DD'     0.2085992109128247E-3

If I understand the argument reduction, they get it down to between
sqrt(2)/2 and sqrt(2).

This supposed to be good for 56 bits.

-- glen




.