Numpy float128 overflows on ppc64le but works on x86_64

[npanpaliya]

While working on a tensorflow's test case, I came across a problem on ppc64le Linux system that boils down to numpy. I've also written a small test program to show case the behavior on both Power8 and x86_64. This is with numpy version 1.13.3.
Test Program

import numpy
t = numpy.float128(1.7976931e+308)
res = numpy.sqrt(t**2 + 1)
print("Result of sq root of square of 1.7976931e+308:  ", res)

Output on Power with numpy version 1.13.3:

# python test_numpyfloat128.py
test_numpyfloat128.py:3: RuntimeWarning: overflow encountered in longdouble_scalars
  res = numpy.sqrt(t**2 + 1)
('Result of sq root of square of 1.7976931e+308:  ', inf)

Min/max range of float128 on Power:

>>> numpy.finfo(numpy.float128)
finfo(resolution=1e-31, min=-1.79769313486e+308, max=1.79769313486e+308, dtype=float128)

Output on x86 with numpy version 1.13.3:

$ python test_numpyfloat128.py
('Result of sq root of square of 1.7976931e+308:  ', 1.7976931000000000494e+308)

Min/max range of float128 on x86:

>>> numpy.finfo(numpy.float128)
finfo(resolution=1e-18, min=-1.18973149536e+4932, max=1.18973149536e+4932, dtype=float128)

I've read many documents/links about float128 on Power and in general. float128 being platform dependent, is long double on x86 and some other architectures but on Power, it is double double[2]. On Power, I checked long double type is available and it is of 128-bit through a small C/C++ program. Some of the below links say that __float128 type is also available on Power7 and Power8 and we need to have VSX enabled[3], but unfortunately I'm unable to test it in a C/C++ program.

1. https://www.ibm.com/support/knowledgecenter/en/ssw_aix_71/com.ibm.aix.genprogc/128bit_long_double_floating-point_datatype.htm
2. https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic
3. https://gcc.gnu.org/wiki/Ieee128PowerPC#A1.3_Sizes

I have also come across many issues in numpy on Power around this float128. Also as per npy_common.h, it seems that float128 is not available on systems where size of double is 64bits, instead float64 will be used. But I'm still not sure whether Power has an unsupported feature or it's a bug in numpy on PPC.
Kindly provide me some pointers on this. Or if there is any way, with which above test can give me same results as that of x86_64.

[eric-wieser]

Min/max range of float128 on Power

I'm sure you noticed this, but that looks like the range of double to me.

What does np.longdouble is np.float128 give?

[npanpaliya]

Yes, that's the range of double. np.longdouble == np.float128 gives True on Power.

[eric-wieser]

Can you try with is as well as ==, just to be sure?

[npanpaliya]

Okay. I tried -

>>> np.longdouble == np.float128
True
>>> np.longdouble is np.float128
True

[eric-wieser]

Your last link explains everything:

The PowerPC GCC compiler for Linux, AIX, and Darwin do not use the IEEE 754R 128-bit floating point type for the long double data type. Instead the PowerPC compiler uses a software type known as IBM extended double type. The IBM extended double type is a pair of double values that give the user more bits of precision, but no additional range for the mantissa.

The only remaining question is why numpy calls this type float128, for which the answer is in the size section:

IBM extended double | 16 bytes

So at best, the bug report here is that float128 is a misleading name that can be confused with the " IEEE 754R" spec. If you want to open a new issue about naming, feel free.

[eric-wieser]

From the same page, compiling numpy with -mabi=ieeelongdouble will solve your problem.

[npanpaliya]

Thanks @eric-wieser for your prompt replies and concluding answers. I'll try building numpy with -mabi=ieeelongdouble. I also agree that naming of float128 should be appropriate on PPC, may be some error should be given while using it.
Thanks again!

[npanpaliya]

I also saw a numpy's fork on PPC64 repo. Probably it might have addressed this. I will check that too.

[charris]

The float128 name may become problematic at some point. Currently, it is either IEEE extended precision, IBM extended double, or IEEE quad precision (was on SUN). NumPy does recognize the different types and they should show up in finfo for the type. At some point I expect quad precision to make a comeback -- gcc has support -- and then we should probably come up with a better name than float128.

[charris]

Maybe we should also add the "official" name of the type to finfo?

[eric-wieser]

And for iinfo while we're at it.

[ahaldane]

Hmm I didn't specifically account for IBM extended double in the 1.14 Dragon4 float-printing updates. I didn't realize they were used on ppc. That means they probably print incorrectly in 1.14.

I think I have a ppc system somewhere, I'll investigate.

[charris]

Hmm I didn't specifically account for IBM extended double in the 1.14 Dragon4 float-printing updates.

Oops, I missed that. Might not be simple, the precision is a bit odd ...

[eric-wieser]

Doesn't this mean that arrays saved with dtype f16 are inherently non-portable? Do we need a different type code for each type of long double, so that we can at least raise an error when loading pickle files from another machine, rather than doing the wrong thing?

We already solve this for integer types by recording the number of bytes rather than the C type name.