如何使用 numpy longdouble dtype？

Question

我正在尝试在我的 Python 代码中使用 np.longdouble dtype，并正在尝试使用 NumPy 来处理我从使用 Cython 编译的 C 模块中获得的长双精度数。

假设我这样做：

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 转换为浮点数，其表示为inf。但是，有人可以提出解决方法吗？有没有办法创建不首先转换为浮点数的 NumPy longdoubles？

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

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

Answer 1

您可能需要考虑一些事项。

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

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

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

据我所知，细节是特定于实现的。因此，这是可行的，但需要格外小心，并且可能存在可移植性问题。