12

I have been playing C99's quad precision long double. It is my understanding that (platform specific) numpy supports long double and 128bit floats.

I have run across something I cannot explain however.

Given:

>>> import numpy as np

Calculate a number that will require more than 64 bits but less than 128 bits to represent as an integer:

>>> 2**64+2
18446744073709551618          # note the '8' at the end
>>> int(2**64+2)
18446744073709551618          # same obviously

If I calculate the same number in C99 128 bit long double, I get 18446744073709551618.000000

Now, if I use numpy long double:

>>> a=np.longdouble(2)
>>> b=np.longdouble(64)
>>> a**b+a
18446744073709551618.0              # all good...

What about these incorrect results:

>>> np.longdouble(2**64+2)
18446744073709551616.0             # Note '6'; appears 2**64 not done in long double
>>> np.longdouble(int(2**64+2))
18446744073709551616.0             # can't force the use of a Python long
>>> n=int(2**64+2)
>>> np.longdouble(n)
18446744073709551616.0
>>> np.longdouble(18446744073709551618)
18446744073709551616.0             # It really does not want to do '8' at the end

But, this works:

>>> np.longdouble(2**64)+2
18446744073709551618.0

Question: Does numpy have issues converting values correctly into long doubles? Is there something I am doing incorrect?

dawg
  • 98,345
  • 23
  • 131
  • 206
  • Did you try to do the same in C? The code you linked does *not* do what you are trying to do with numpy. Try to first compute `2**64 + 2` as an integer and *then* assign it to a `long double`. (I mean something like: `long long val = (1 << 64) + 2; long double res = (long double)val`) – Bakuriu Aug 30 '13 at 15:47
  • @Bakuriu: Yes, I tried in C. Yes, I tried `n=int(2**64+2); np.longdouble(n)` which still does not add the 2. – dawg Aug 30 '13 at 15:51
  • 1
    Using the expression `n=(2**x+2); np.longdouble(n)==n` will return false for `x>53`. Seems like a strange break point. – Daniel Aug 30 '13 at 16:12
  • 1
    @Ophion: A [very interesting breaking point](http://en.wikipedia.org/wiki/Double-precision_floating-point_format) since 53 would lead me to believe that the conversion then is being done at 64 bits vs 128 bits. – dawg Aug 30 '13 at 16:16
  • @drewk yes, the conversion is being done via 64-bit double (52 fraction bits). See below. – ecatmur Aug 30 '13 at 16:36
  • 1
    @drewk: What platform are you on? Is numpy.longdouble really the 128-bit IEEE 754 quadruple precision format on your system, or is it just the usual 80-bit x87 format in disguise? C99 Annex F allows either of those possibilities for long double, but the latter is *far* more common. – Mark Dickinson Aug 30 '13 at 17:01
  • @MarkDickinson: It is 64 bit OS X and it is true 128-bit long double floats IF numpy used all the correct compiler pragmas. In this particular case tho, 80-bit vs 128-bit would be a correct result, correct? It is the conversion to the 64-bit float that is the issue. – dawg Aug 30 '13 at 17:05
  • 1
    Ah, nice; I didn't realise it was possible to get real 128-bit floats in NumPy. Agreed that for this question it doesn't make a difference. I guess you could confirm 80-bit versus 128-bit by trying something like `longdouble(2**50 + 1) * longdouble(2**50 - 1) - longdouble(2 ** 100)`, which should give a result of `-1` for quad precision. – Mark Dickinson Aug 30 '13 at 17:12
  • Quite simply, yes numpy has an issue there (I think it goes via floats, maybe python floats). If you care about it, maybe you can fix it ;)? Btw. longdouble (and the equivalent float96/128) are **not** quad precision, and the implementation depends on the system. – seberg Aug 30 '13 at 18:50

2 Answers2

10

You're trying to perform a type conversion between non-directly-convertible types. Take a look at the stack:

#0  0x00002aaaaab243a0 in PyLong_AsDouble ()
   from libpython2.7.so.1.0
#1  0x00002aaaaab2447a in ?? ()
   from libpython2.7.so.1.0
#2  0x00002aaaaaaf8357 in PyNumber_Float ()
   from libpython2.7.so.1.0
#3  0x00002aaaae71acdc in MyPyFloat_AsDouble (obj=0x2aaaaae93c00)
    at numpy/core/src/multiarray/arraytypes.c.src:40
#4  0x00002aaaae71adfc in LONGDOUBLE_setitem (op=0x2aaaaae93c00, 
    ov=0xc157b0 "", ap=0xbf6ca0)
    at numpy/core/src/multiarray/arraytypes.c.src:278
#5  0x00002aaaae705c82 in PyArray_FromAny (op=0x2aaaaae93c00, 
    newtype=0x2aaaae995960, min_depth=<value optimized out>, max_depth=0, 
    flags=0, context=<value optimized out>)
    at numpy/core/src/multiarray/ctors.c:1664
#6  0x00002aaaae7300ad in longdouble_arrtype_new (type=0x2aaaae9938a0, 
    args=<value optimized out>, __NPY_UNUSED_TAGGEDkwds=<value optimized out>)
    at numpy/core/src/multiarray/scalartypes.c.src:2545

As you can see, the Python long (unlimited-precision integer) 2**64 + 2 is being converted to float (i.e. 64-bit double), which loses precision; the float is then used to initialise the long double but the precision has already been lost.

The problem is that 128-bit double is not a native Python type, so long doesn't have a native conversion to it, only to 64-bit double. It probably would be possible for NumPy to detect this situation and perform its own conversion using the long C API, but might be fairly complicated for relatively little benefit (you can just do arithmetic in np.longdouble from the start).

ecatmur
  • 152,476
  • 27
  • 293
  • 366
  • 4
    It is a potentially sneaky bug, no? – dawg Aug 30 '13 at 16:52
  • 2
    @dawg that's why [numpy docs](http://docs.scipy.org/doc/numpy-dev/user/basics.types.html#extended-precision) suggest to _"test your code with the value `1 + np.finfo(np.longdouble).eps`."_ – ivan_pozdeev Sep 30 '16 at 09:04
2

NumPy does not provide quad precision on x86 machines. It does provide access to the C long double type (as provided by the compilation environment; with MSVC this may be 64 bits, with GCC it is normally 80 bits) as np.longdouble. The types np.float96 and np.float128 are simply long doubles padded to 96 or 128 bits (for aligned memory access). See the numpy docs. To get quad precision in numpy you need to use a hardware platform and compiler where long double is actual quad precision.

While it would be possible for numpy to support quad precision using compiler support (GCC's float128) or external libraries, this has not been implemented. It would also be possible to write a third-party importable module that made one of these available, but that has not been done either.

Note too that even when using np.longdouble, it is easy to lose precision: for example, the % operator forces numpy to pass its numbers through python floats, throwing away any extra precision.

user2475529
  • 260
  • 2
  • 5