On Thu, Apr 1, 2010 at 9:01 PM, Charles R Harris
<charlesr.harris@gmail.com> wrote:
On Thu, Apr 1, 2010 at 8:16 PM, Charles R Harris
<charlesr.harris@gmail.com> wrote:
On Thu, Apr 1, 2010 at 7:59 PM, David Cournapeau
<david@silveregg.co.jp> wrote:
Anne Archibald wrote:
>
> First I guess we should check which systems don't have log1p
This is already done - we do use the system log1p on linux (but note
that log2_1p is not standard AFAIK). I would guess few systems outside
windows don't have log1p, given that msun has an implementation,
I see that msun uses the same series I came to. However, my rational (Pade) approximation is a lot better than their polynomial.
In fact I get better than 119 bits using the same range as sun and the ratio of two 7'th degree polynomials. I suspect it's better than that, but I only have mpmath set to that precision. Sun got 58 bits with a 14'th degree polynomial. So we can definitely improve on sun.
A few notes. The Sun routine is pretty good, but results of almost equal quality can be obtained using barycentric interpolation at the Chebyshev II points, i.e., one doesn't need the slight improvement that comes from the Reme(z) algorithm. To get extended precision requires going from degree 14 to degree 16, and quad precision needs degree 30. Note that the function in question is even, so the effective degree is half of those quoted. Next up: using Chebyshev points to interpolate the polynomials from the Pade representation. The Pade version doesn't actually gain any precision over the power series of equivalent degree, but it does allow one to use polynomials of half the degree in num/den.
Short term, however, I'm going to take Anne's advice and simply use log1p. That should fix the problem for most users.
It's amusing that the Reme[z] typo has persisted in the Sun documentation all these years ;)
Chuck