Numpy longdouble算术似乎没有转换的长双倍

Question

我一直在玩C99的四倍精度长双.据我所知,(平台特定的)numpy支持long double和128bit float.

我遇到了一些我无法解释的事情.

鉴于:

>>> import numpy as np

计算需要超过64位但小于128位的数字表示为整数:

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

如果我在C99 128位长双倍中计算相同的数字,我得到18446744073709551618.000000

现在,如果我使用numpy long double:

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

这些不正确的结果怎么样:

>>> 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

但是,这有效:

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

问题:numpy是否存在将值正确转换为长双打的问题？有什么我做错了吗？

Answer 1

您正尝试在非直接可转换类型之间执行类型转换.看看堆栈:

#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

正如您所看到的,Python long(无限精度整数)2**64 + 2正在被转换为float(即64位双精度),这会失去精度; 然后使用float初始化long double,但精度已经丢失.

问题是128位双精度不是本机Python类型,所以long没有本机转换,只有64位双精度.NumPy可能有可能检测到这种情况并使用longC API执行自己的转换,但可能相对较少的好处相当复杂(您可以np.longdouble从一开始就进行算术运算).