具有 mpz/mpfr 值的 numpy 数组答案

【问题标题】：numpy array with mpz/mpfr values具有 mpz/mpfr 值的 numpy 数组
【发布时间】：2013-03-09 05:34:15
【问题描述】：

我想要一个带有 mpz/mpfr 值的 numpy 数组。因为我的代码：

import numpy as np
import gmpy2
A=np.ones((5,5));
print A/gmpy2.mpfr(1);

生成：

RuntimeWarning: invalid value encountered in divide
  print A/gmpy2.mpfr(1);
[[1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]]

据我所知，不可能将 gmpy mpfr 转换为 numpy float64。那么我怎样才能首先获得一个带有 mpfr 值的 numpy 数组呢？

谢谢。

【问题讨论】：

标签： python numpy type-conversion gmp

【解决方案1】：

您需要使用dtype=object 创建数组，然后您可以在数组中使用任何python 类型。我没有安装 gmpy2，但下面的例子应该显示它是如何工作的：

In [3]: a = np.ones((5, 5), dtype=object)

In [5]: import fractions

In [6]: a *= fractions.Fraction(3, 4)

In [7]: a
Out[7]: 
array([[3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4]], dtype=object)

拥有一个 dtype=object 的 numpy 数组可能会产生误导，因为强大的 numpy 机制使得使用标准 dtype 的操作超快，现在由默认对象的 python 运算符处理，这意味着速度将不再存在：

In [12]: b = np.ones((5, 5)) * 0.75

In [13]: %timeit np.sum(a)
1000 loops, best of 3: 1.25 ms per loop

In [14]: %timeit np.sum(b)
10000 loops, best of 3: 23.9 us per loop

【讨论】：

再一次，fractions.Fraction 并不是一个特别快的课程。我想知道原生 Numpy 数组和 mpfr 数组之间的速度差异是多少，因为 mpfr 是一个相对低级的 C 包装类。
@nneonneo 我认为问题不在于操作的速度，而在于每一个都涉及 Python 函数调用，而另一个不会发生这种情况numpy 数据类型。
是的，有 Python 函数调用，但对于用 C 实现的类，这些调用的开销可能非常小。 Fraction 是用纯 Python 实现的，所以每次调用都是很多字节码指令。
@nneonneo 我只是用一组 numpy 浮点数和另一个 dtype 对象重复计时，所以 Python 浮点数，你确实有一点。对于(5, 5) 数组，性能差异小于 5%，尽管 numpy 领先。但是使用(500, 500) 数组，numpy 的速度要快 25 倍以上。
如果您使用= 而不是*=，则无需指定dtype。该问题可能是两个库之一中的问题。

【解决方案2】：

我相信这是两个库之一中的错误。我也相信它是固定的。

输入：

import sys
import numpy as np
import gmpy2

print(sys.version)
print(np.__version__)
print(gmpy2.version)

A=np.ones((5,5));
print(A/gmpy2.mpfr(1))

输出：

3.4.2 (v3.4.2:ab2c023a9432, Oct  6 2014, 22:15:05) [MSC v.1600 32 bit (Intel)]
1.9.1
2.0.5
[[mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]]

Numpy 没有正确说明遇到未知类型时该怎么做，或者 gmpy2 没有指定如何除以某物 (__rdiv__)。

没有必要指定ndarray 的dtype，除非您打算改写它的元素。诸如乘法之类的运算会产生一个新的ndarray，Numpy 会找出要使用的dtype。

【讨论】：