如何使用numpy longdouble数据类型？

Question

我正在尝试使用Python代码中的np.longdoubledtype，并尝试使用NumPy操作从用Cython编译的C模块中获得的长double。

假设我这样做：

import numpy as np

print np.finfo(np.longdouble)
Machine parameters for float128
---------------------------------------------------------------------
precision= 18   resolution= 1e-18
machep=   -63   eps=        1.08420217249e-19
negep =   -64   epsneg=     5.42101086243e-20
minexp=-16382   tiny=       3.36210314311e-4932
maxexp= 16384   max=        1.18973149536e+4932
nexp  =    15   min=        -max
---------------------------------------------------------------------


a = np.longdouble(1e+346)

a
Out[4]: inf

b = np.longdouble(1e+347)

b
Out[6]: inf

c = a/b
/usr/lib/python2.7/site-packages/spyderlib/widgets/externalshell/start_ipython_kernel.py:1:
RuntimeWarning: invalid value encountered in longdouble_scalars
  # -*- coding: utf-8 -*-

c
Out[8]: nan

a.dtype, b.dtype, c.dtype
Out[9]: (dtype('float128'), dtype('float128'), dtype('float128'))

实际上，它与this问题中的问题是同一个问题，我理解Python首先将1e+346转换为一个float，其表示形式是inf。不过，有人能提出解决办法吗？有没有办法创建不首先转换为浮点的NumPy longdouples？

我有一个C模块，可以输出长双精度数，我想在数据类型为np.longdouble的numpy数组中使用它。

即使解决方案涉及重新编译Python/NumPy，我也愿意尝试。

Answer 1

有一些事情你可能需要考虑。

首先，这里一团糟。NumPy知道longdouble和float128。不幸的是，这些名称是误导性的，底层实现是C long double，它通常（但不一定总是）是一个80位浮点。（实际上，通过查看“精度”，您可以在这里看到它；18位大约是60位，80位浮点在尾数中有64位。如果使用真正的128位浮点运算，精度将在34位左右。）

可能没有任何直接的方法可以将长double作为参数传递给C函数，但是如果改为传递指针，则可以避免此问题。例如，可以将数组数据作为uint8（通过使用myarray.view(dtype='uint8')）传递，并将指向缓冲区的指针转换为C程序中的long double*。至少Python与类型转换无关。（很可能您不需要使用view，因为毕竟您只导出指向数组缓冲区的指针。）

请注意，这个技巧依赖于编译器在编译Python和C程序时具有相同类型的设置。除了精度差异外，还可能存在字节顺序差异（如果程序在同一台计算机上运行，则很少）和对齐差异。我的Python似乎将longdouble项以16字节的边界对齐（即每个元素总是有6字节的零），但是C编译器可能使用10/12/16字节的对齐方式。

据我所知，细节是具体实现的。因此，这是可行的，但需要一些额外的注意，可能存在可移植性问题。