[Numpy-discussion] OS X PPC problem with Numpy 1.3.0b1
Robert Pyle
rpyle@post.harvard....
Mon Mar 23 13:34:34 CDT 2009
Hi all,
This is a continuation of something I started last week, but with a
more appropriate subject line.
To recap, my machine is a dual G5 running OS X 10.5.6, my python is
Python 2.5.2 |EPD Py25 4.1.30101| (r252:60911, Dec 19 2008,
15:28:32)
and numpy 1.3.0b1 was installed from the source tarball in the
straightforward way with
sudo python setup.py install
On Mar 19, 2009, at 3:46 PM, Charles R Harris wrote:
> On Mar 19, 2009, at 1:38 PM, Pauli Virtanen wrote:
>
> Thanks for tracking this! Can you check what your platform gives for:
>
> > import numpy as np
> > info = np.finfo(np.longcomplex)
> > print "eps:", info.eps, info.eps.dtype
> > print "tiny:", info.tiny, info.tiny.dtype
> > print "log10:", np.log10(info.tiny), np.log10(info.tiny/info.eps)
>
> eps: 1.3817869701e-76 float128
> tiny: -1.08420217274e-19 float128
> log10: nan nan
>
> The log of a negative number is nan, so part of the problem is the
> value of tiny. The size of the values also look suspect to me. On my
> machine
>
> In [8]: finfo(longcomplex).eps
> Out[8]: 1.084202172485504434e-19
>
> In [9]: finfo(float128).tiny
> Out[9]: array(3.3621031431120935063e-4932, dtype=float128)
>
> So at a minimum eps and tiny are reversed.
>
> I started to look at the code for this but my eyes rolled up in my
> head and I passed out. It could use some improvements...
>
> Chuck
I have chased this a bit (or perhaps 128 bits) further.
The problem seems to be that float128 is screwed up in general. I
tracked the test error back to lines 95-107 in
/PyModules/numpy-1.3.0b1/build/lib.macosx-10.3-ppc-2.5/numpy/lib/
machar.py
Here is a short program built from these lines that demonstrates what
I believe to be at the root of the test failure.
######################################
#! /usr/bin/env python
import numpy as np
import binascii as b
def t(type="float"):
max_iterN = 10000
print "\ntesting %s" % type
a = np.array([1.0],type)
one = a
zero = one - one
for _ in xrange(max_iterN):
a = a + a
temp = a + one
temp1 = temp - a
print _+1, b.b2a_hex(temp[0]), temp1
if any(temp1 - one != zero):
break
return
if __name__ == '__main__':
t(np.float32)
t(np.float64)
t(np.float128)
######################################
This tries to find the number of bits in the significand by
calculating ((2.0**n)+1.0) for increasing n, and stopping when the sum
is indistinguishable from (2.0**n), that is, when the added 1.0 has
fallen off the bottom of the significand.
My print statement shows the power of 2.0, the hex representation of
((2.0**n)+1.0), and the difference ((2.0**n)+1.0) - (2.0**n), which
one expects to be 1.0 up to the point where the added 1.0 is lost.
Here are the last few lines printed for float32:
19 49000010 [ 1.]
20 49800008 [ 1.]
21 4a000004 [ 1.]
22 4a800002 [ 1.]
23 4b000001 [ 1.]
24 4b800000 [ 0.]
You can see the added 1.0 marching to the right and off the edge at 24
bits.
Similarly, for float64:
48 42f0000000000010 [ 1.]
49 4300000000000008 [ 1.]
50 4310000000000004 [ 1.]
51 4320000000000002 [ 1.]
52 4330000000000001 [ 1.]
53 4340000000000000 [ 0.]
There are 53 bits, just as IEEE 754 would lead us to hope.
However, for float128:
48 42f00000000000100000000000000000 [1.0]
49 43000000000000080000000000000000 [1.0]
50 43100000000000040000000000000000 [1.0]
51 43200000000000020000000000000000 [1.0]
52 43300000000000010000000000000000 [1.0]
53 43400000000000003ff0000000000000 [1.0]
54 43500000000000003ff0000000000000 [1.0]
Something weird happens as we pass 53 bits. I think lines 53 and 54
*should* be
53 43400000000000008000000000000000 [1.0]
54 43500000000000004000000000000000 [1.0]
etc., with the added 1.0 continuing to march rightwards to extinction,
as before.
The calculation eventually terminates with
1022 7fd00000000000003ff0000000000000 [1.0]
1023 7fe00000000000003ff0000000000000 [1.0]
1024 7ff00000000000000000000000000000 [NaN]
(7ff00000000000000000000000000000 == Inf).
This totally messes up the remaining parts of machar.py, and leaves us
with Infs and Nans that give the logs of negative numbers, etc. that
we saw last week.
But wait, there's more! I also have an Intel Mac (a MacBook Pro).
This passes numpy.test(), but when I look in detail with the above
code, I find that the float128 significand has only 64 bits, leading
me to suspect that it is really the so-called 80-bit "extended
precision".
Is this true? If so, should there be such large differences between
architectures for the same nominal precision? Is float128 just
whatever the underlying C compiler thinks a "long double" is?
I've spent way too much time on this, so I'm going to bow out here
(unless someone can suggest something for me to try that won't take
too much time).
Bob
More information about the Numpy-discussion
mailing list