[Numpy-discussion] Unexpected float96 precision loss

[David Cournapeau]

On Thu, Sep 2, 2010 at 10:03 AM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
>
>
> On Wed, Sep 1, 2010 at 3:30 PM, Michael Gilbert
> <michael.s.gilbert at gmail.com> wrote:
>>
>> On Wed, 1 Sep 2010 21:15:22 +0000 (UTC), Pauli Virtanen wrote:
>> > Wed, 01 Sep 2010 16:26:59 -0400, Michael Gilbert wrote:
>> > > I've been using numpy's float96 class lately, and I've run into some
>> > > strange precision errors.
>> > [clip]
>> > >   >>> x = numpy.array( [0.01] , numpy.float96 )
>> > [clip]
>> > > I would expect the float96 calculation to also produce 0.0 exactly as
>> > > found in the float32 and float64 examples.  Why isn't this the case?
>> >
>> > (i) It is not possible to write long double literals in Python.
>> >     "float96(0.0001)" means in fact "float96(float64(0.0001))"
>> >
>> > (ii) It is not possible to represent numbers 10^-r, r > 1 exactly
>> >      in base-2 floating point.
>> >
>> > So if you write "float96(0.0001)", the result is not the float96 number
>> > closest to 0.0001, but the 96-bit representation of the 64-bit number
>> > closest to 0.0001. Indeed,
>> >
>> > >>> float96(0.0001), float96(1.0)/1000
>> > (0.00010000000000000000479, 0.00099999999999999999996)
>>
>> Interesting.  float96( '0.0001' ) also seems to evaluate to the first
>> result. I assume that it also does a float64( '0.0001' ) conversion
>> first. I understand that you can't change how python passes in floats,
>> but wouldn't it be better to exactly handle strings since those can be
>> converted exactly, which is what the user wants/expects?
>>
>
> Well, yes. But then we would need to write our own routines for the
> conversions...

I think that it is needed at some point, though. There are quite a few
bugs related to this kind of issues in NumPy,

cheers,

David