[Numpy-discussion] Bug in logaddexp2.reduce
Charles R Harris
Sat Apr 3 12:23:05 CDT 2010
On Thu, Apr 1, 2010 at 9:01 PM, Charles R Harris
> On Thu, Apr 1, 2010 at 8:16 PM, Charles R Harris <
> email@example.com> wrote:
>> On Thu, Apr 1, 2010 at 7:59 PM, David Cournapeau <firstname.lastname@example.org>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 ;)
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion