不完整的 gamma 函数：此代码在 cython、C 或 Fortran 中能否变得更快？

Question

作为一大段代码的一部分，我需要计算不完整的伽马函数数组。例如，我需要一个函数 returns (log of) (gamma(z + m, a, inf)/m!) for m in [0, m_max], for various values of m_max（通常在 400 左右）、z 和 a。我需要快点做这件事。目前，这一步是我代码中最慢的大约 2 倍。然而，完整的代码需要大约一天的时间才能完成运行，因此将这一步的计算时间减少 2 会节省我很多时间。

我正在使用以下 cython 代码进行计算：

import numpy as np
cimport numpy as np
from mpmath import mp

sp_max = 5000 

def log_factorial(k):
    return np.sum(np.log(np.arange(1., k + 1., dtype=np.float)))

log_factorial_ary = np.vectorize(log_factorial)(np.arange(sp_max))

gamma_mem = mp.memoize(mp.gamma)
gammainc_mem = mp.memoize(mp.gammainc)

def gammainc_up_fct_ary_log(np.int m_max, np.float z, np.float a):
    cdef np.ndarray gi_list = np.zeros(m_max + 1, dtype=np.float)    
    gi_list[0] = np.float(gammainc_mem(z, a))
    cdef np.ndarray i_array = np.arange(1., m_max + 1., dtype=np.float)
    cdef Py_ssize_t i 
    for i in np.arange(1, m_max + 1):
        gi_list[i] = (i_array[i-1] - 1. + z)*gi_list[i-1]/i + np.exp((i_array[i-1] - 1. + z)*np.log(a) - a - log_factorial_ary[i])
    return gi_list

例如，当我调用 gammainc_up_fct_ary_log(400,-0.3,10.0) 时，大约需要 0.015-0.025 秒。我想将其速度至少提高 2 倍（或者，理想情况下，尽可能快）。

是否有明确的方法可以使用 cython 加速此计算？如果不是，C 或 Fortran 会快得多吗？如果是这样，用那种语言编写这个函数然后从 python 调用代码的最快方法是什么（我的其余代码是用 python/cython 编写的）。

提前致谢。

Answer 1

你的cython版本有几个大问题：

i_array 没用，你可以安全地用 i

i_array[i-1]

您没有充分利用 cython。如果您查看代码中 cython -a 的输出，您会发现 cython 只是生成对 C-API 的调用，而您需要调用 C 代码才能拥有它运行快.

这是您可以实现的示例（不完整，但加速已经很好）

import numpy as np
cimport numpy as np
cimport cython
from mpmath import mp

cdef extern from "math.h":
    double log(double x) nogil
    double exp(double x) nogil

sp_max = 5000 

def log_factorial(k):
    return np.sum(np.log(np.arange(1., k + 1., dtype=np.float)))

factorial_ary = np.array([np.float(mp.factorial(m)) for m in np.arange(sp_max)])
log_factorial_ary = np.vectorize(log_factorial)(np.arange(sp_max))

gamma_mem = mp.memoize(mp.gamma)
gammainc_mem = mp.memoize(mp.gammainc)

def gammainc_up_fct_ary_log(m_max, z, a):
    return gammainc_up_fct_ary_log_impl(m_max, z, a)

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef gammainc_up_fct_ary_log_impl(int m_max, double z, double a):
    cdef double[::1] gi_list = np.zeros(m_max + 1, dtype=np.float)
    gi_list[0] = gammainc_mem(z, a)
    cdef Py_ssize_t i
    for i in range(1, m_max + 1):
        t0 = (i - 1. + z)
        t1 = (i - 1. + z)*log(a) - a
        gi_list[i] = t0*gi_list[i-1]/i + exp(t1 - log_factorial_ary[i])
    return gi_list

运行宁此代码给我：

python -m timeit -s 'from ff import gammainc_up_fct_ary_log' 'gammainc_up_fct_ary_log(400,-0.3,10.0)'

10000 次循环，3 次循环中的最佳次数：每次循环 132 微秒

虽然你的版本几乎没有：

python -m timeit -s 'from ff import gammainc_up_fct_ary_log' 'gammainc_up_fct_ary_log(400,-0.3,10.0)'

100 次循环，3 次循环最佳：每次循环 2.44 毫秒

不完整的 gamma 函数：此代码在 cython、C 或 Fortran 中能否变得更快？

Incomplete gamma functions: can this code get any faster in cython, C, or Fortran?

c

python

performance

fortran

cython