[Numpy-discussion] Bug in logaddexp2.reduce

Anne Archibald peridot.faceted@gmail....
Thu Apr 1 01:04:52 CDT 2010


On 1 April 2010 01:59, Charles R Harris <charlesr.harris@gmail.com> wrote:
>
>
> On Wed, Mar 31, 2010 at 11:46 PM, Anne Archibald <peridot.faceted@gmail.com>
> wrote:
>>
>> On 1 April 2010 01:40, Charles R Harris <charlesr.harris@gmail.com> wrote:
>> >
>> >
>> > On Wed, Mar 31, 2010 at 11:25 PM, <josef.pktd@gmail.com> wrote:
>> >>
>> >> On Thu, Apr 1, 2010 at 1:22 AM,  <josef.pktd@gmail.com> wrote:
>> >> > On Thu, Apr 1, 2010 at 1:17 AM, Charles R Harris
>> >> > <charlesr.harris@gmail.com> wrote:
>> >> >>
>> >> >>
>> >> >> On Wed, Mar 31, 2010 at 6:08 PM, <josef.pktd@gmail.com> wrote:
>> >> >>>
>> >> >>> On Wed, Mar 31, 2010 at 7:37 PM, Warren Weckesser
>> >> >>> <warren.weckesser@enthought.com> wrote:
>> >> >>> > T J wrote:
>> >> >>> >> On Wed, Mar 31, 2010 at 1:21 PM, Charles R Harris
>> >> >>> >> <charlesr.harris@gmail.com> wrote:
>> >> >>> >>
>> >> >>> >>> Looks like roundoff error.
>> >> >>> >>>
>> >> >>> >>>
>> >> >>> >>
>> >> >>> >> So this is "expected" behavior?
>> >> >>> >>
>> >> >>> >> In [1]: np.logaddexp2(-1.5849625007211563, -53.584962500721154)
>> >> >>> >> Out[1]: -1.5849625007211561
>> >> >>> >>
>> >> >>> >> In [2]: np.logaddexp2(-0.5849625007211563, -53.584962500721154)
>> >> >>> >> Out[2]: nan
>> >> >>> >>
>> >> >>> >
>> >> >>> > Is any able to reproduce this?  I don't get 'nan' in either 1.4.0
>> >> >>> > or
>> >> >>> > 2.0.0.dev8313 (32 bit Mac OSX).  In an earlier email T J reported
>> >> >>> > using
>> >> >>> > 1.5.0.dev8106.
>> >> >>>
>> >> >>>
>> >> >>>
>> >> >>> >>> np.logaddexp2(-0.5849625007211563, -53.584962500721154)
>> >> >>> nan
>> >> >>> >>> np.logaddexp2(-1.5849625007211563, -53.584962500721154)
>> >> >>> -1.5849625007211561
>> >> >>>
>> >> >>> >>> np.version.version
>> >> >>> '1.4.0'
>> >> >>>
>> >> >>> WindowsXP 32
>> >> >>>
>> >> >>
>> >> >> What compiler? Mingw?
>> >> >
>> >> > yes, mingw 3.4.5. , official binaries release 1.4.0 by David
>> >>
>> >> sse2 Pentium M
>> >>
>> >
>> > Can you try the exp2/log2 functions with the problem data and see if
>> > something goes wrong?
>>
>> Works fine for me.
>>
>> If it helps clarify things, the difference between the two problem
>> input values is exactly 53 (and that's what logaddexp2 does an exp2
>> of); so I can provide a simpler example:
>>
>> In [23]: np.logaddexp2(0, -53)
>> Out[23]: nan
>>
>> Of course, for me it fails in both orders.
>>
>
> Ah, that's progress then ;) The effective number of bits in a double is 53
> (52 + implicit bit). That shouldn't cause problems but it sure looks
> suspicious.

Indeed, that's what led me to the totally wrong suspicion that
denormals have something to do with the problem. More data points:

In [38]: np.logaddexp2(63.999, 0)
Out[38]: nan

In [39]: np.logaddexp2(64, 0)
Out[39]: 64.0

In [42]: np.logaddexp2(52.999, 0)
Out[42]: 52.999000000000002

In [43]: np.logaddexp2(53, 0)
Out[43]: nan

It looks to me like perhaps the NaNs are appearing when the smaller
term affects only the "extra" bits provided by the FPU's internal
larger-than-double representation. Some such issue would explain why
the problem seems to be hardware- and compiler-dependent.

Anne


More information about the NumPy-Discussion mailing list