1

我有一些涉及非常快爆炸的阶乘的计算，所以我决定使用任意精度库mpmath。

我的代码如下所示：

import numpy as np
import mpmath as mp
import time

a    = np.linspace( 0, 100e-2, 100 )
b    = np.linspace( 0, np.pi )
c    = np.arange( 30 )

t    = time.time()
M    = np.ones( [ len(a), len(b), len(c) ] )
A, B = np.meshgrid( a, b, indexing = 'ij' )
temp = A**2 + B
temp = np.reshape( temp, [ len(a), len(b), 1 ] )
temp = np.repeat( temp, len(c), axis = 2 )
M   *= temp
print 'part1:      ', time.time() - t
t    = time.time()

temp = np.array( [ mp.fac(x) for x in c ] )
temp = np.reshape( temp, [ 1, 1, len(c) ] )
temp = np.repeat(  temp, len(a), axis = 0 )
temp = np.repeat(  temp, len(b), axis = 1 )
print 'part2 so far:', time.time() - t
M   *= temp
print 'part2 finally', time.time() - t
t    = time.time()

似乎花费最多时间的是最后一行，我怀疑这是因为M有一堆floats 和temp一堆mp.mpfs。M我尝试使用s 进行初始化，mp.mpf但随后一切都变慢了。

这是我得到的输出：

part1:        0.00429606437683
part2 so far: 0.00184297561646
part2 finally 1.9477159977

有什么想法可以加快速度吗？

4

1 回答 1

3

gmpy2比mpmath这种类型的计算要快得多。以下代码在我的机器上运行速度提高了大约 12 倍。

import numpy as np
import gmpy2 as mp
import time

a = np.linspace(0, 100e-2, 100)
b = np.linspace(0, np.pi)
c = np.arange(30)

t = time.time()
M = np.ones([len(a), len(b), len(c)])
A, B = np.meshgrid( a, b, indexing = 'ij' )
temp = A**2+B
temp = np.reshape(temp, [len(a), len(b), 1])
temp = np.repeat(temp, len(c), axis=2)
M *= temp
print 'part1:', time.time() - t
t = time.time()

temp = np.array([mp.factorial(x) for x in c])
temp = np.reshape(temp, [1, 1, len(c)])
temp = np.repeat(temp, len(a), axis=0)
temp = np.repeat(temp, len(b), axis=1)
print 'part2 so far:', time.time() - t
M *= temp
print 'part2:', time.time() - t
t = time.time()

mpmath是用 Python 编写的，通常使用 Python 的本机整数进行计算。如果gmpy2可用，它将使用由 提供的更快的整数类型gmpy2。如果您只需要 直接提供的功能之一gmpy2，那么gmpy2直接使用通常会更快。

更新

我进行了一些实验。实际发生的事情可能不是你所期望的。计算temp时，值可以是整数（math.factorial、gmpy.fac或gmpy2.fac）或浮点值（gmpy2.factorial、mpmath.fac）。当numpycomputes时M *= temp，所有的值temp都被转换为 64 位浮点数。如果该值是整数，则转换会引发溢出错误。如果该值是浮点数，则转换返回无穷大。您可以通过更改c为np.arange(300)并在最后打印来看到这一点M。如果你使用gmpy.facor math.factorial，你会得到 and OverflowError。如果您使用mpmath.factorialor gmpy2.factorial，您将不会得到 anOverflowError但结果M将包含无穷大。

如果您试图避免OverflowError，则需要temp使用浮点值进行计算，以便转换为 64 位浮点数将导致无穷大。

如果您没有遇到OverflowError，那么math.factorial是最快的选择。

如果您试图避免两者OverflowError和无穷大，那么您将需要在整个过程中使用 thempmath.mpf或gmpy2.mpfr浮点类型。（不要尝试使用gmpy.mpf。）

更新#2

这是一个使用gmpy2.mpfr200 位精度的示例。使用c=np.arange(30)，它比您的原始示例快约 5 倍。我使用它来展示它，c = np.arange(300)因为它会生成一个OverflowError或无穷大。较大范围的总运行时间与您的原始代码大致相同。

import numpy as np
import gmpy2
import time

from gmpy2 import mpfr

gmpy2.get_context().precision = 200

a = np.linspace(mpfr(0), mpfr(1), 100)
b = np.linspace(mpfr(0), gmpy2.const_pi())
c = np.arange(300)

t = time.time()
M = np.ones([len(a), len(b), len(c)], dtype=object)
A, B = np.meshgrid( a, b, indexing = 'ij' )
temp = A**2+B
temp = np.reshape(temp, [len(a), len(b), 1])
temp = np.repeat(temp, len(c), axis=2)
M *= temp
print 'part1:', time.time() - t
t = time.time()

temp = np.array([gmpy2.factorial(x) for x in c], dtype=object)
temp = np.reshape(temp, [1, 1, len(c)])
temp = np.repeat(temp, len(a), axis=0)
temp = np.repeat(temp, len(b), axis=1)
print 'part2 so far:', time.time() - t
M *= temp
print 'part2:', time.time() - t
t = time.time()

免责声明：我坚持gmpy2.

于 2014-10-23T23:49:34.283 回答