如何使用numpy的longdouble数据类型？

Question

我正在尝试在我的Python代码中使用np.longdouble这种数据类型，并且想用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'))

其实，这个问题和这个问题是同一个原因。我明白Python首先把1e+346转换成了浮点数，而这个浮点数的表示是inf（无穷大）。不过，有没有人能给我一些解决办法？有没有办法创建NumPy的长双精度数，而不先转换成浮点数？

我有一个C模块，可以输出长双精度数，我想把它们放进一个数据类型为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 字节的对齐方式。

据我所知，具体的细节是依赖于实现的。所以，这个方法是可行的，但需要额外小心，并可能会有可移植性的问题。