Numpy 用Cython求阶乘的近似值|极客教程

用Cython求阶乘的近似值，我们将使用两种求阶乘近似值的方法。首先使用的是斯特灵近似方法（详见http://en.wikipedia.org/wiki/Stirling%27s_approximation），其计算公式是：

然后将使用Ramanujan提出的近似计算方法，其计算公式为：

用Cython求阶乘的近似值

具体步骤

下面将具体介绍怎样用Cython计算阶乘的近似值。我们将用到类型声明。你也许还记得，在Cython中，类型声明是可选的。从理论上讲，声明静态类型会提高代码的执行速度。声明静态类型也带来了编写常规的Python代码时可能没有遇到过的有趣的挑战。但不要担心，我们将尽量让问题简单。

编写Cython代码。

我们将编写的Cython代码看上去和常规的Python代码差不多，只是把函数参数和一个局部变量声明为ndarray类型。为了使用静态类型，需要使用cimport引入NumPy，并且在声明局部变量时需要同时使用cdef关键字。

import numpy
cimport numpy

def ramanujan_factorial(numpy.ndarray n):
    sqrt_pi = numpy.sqrt(numpy.pi, dtype=numpy.float64)
    cdef numpy.ndarray root = (8 * n + 4) * n + 1 
    root = root * n + 1/30.
    root = root ** (1/6.)
    return sqrt_pi * calc_eton(n) * root

def stirling_factorial(numpy.ndarray n):
    return numpy.sqrt(2 * numpy.pi * n) * calc_eton(n)

def calc_eton(numpy.ndarray n):
    return (n/numpy.e) ** n

构建代码。

正如前面几个攻略中介绍过的，我们需要创建一个setup.py文件。在这个文件中，需要通过调用get_include函数，把NumPy相关的目录包括进来。修改之后的setup.py文件的内容如下。

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy

ext_modules = [Extension("factorial", ["factorial.pyx"], include_dirs = [numpy.get_include()])] 

setup(
    name = 'Factorial app',
    cmdclass = {'build_ext': build_ext},
    ext_modules = ext_modules
    )

绘制相对误差。

相对于用NumPy的cumprod函数计算出来的阶乘，我们将分别计算这两种近似计算方法的相对误差，并绘制出结果。

from factorial import ramanujan_factorial
from factorial import stirling_factorial
import numpy
import matplotlib.pyplot

N = 50
numbers = numpy.arange(1, N)
factorials = numpy.cumprod(numbers, dtype=float)

def error(approximations):
    return (factorials - approximations)/factorials

matplotlib.pyplot.plot(error(ramanujan_factorial(numbers)), 'b-')
matplotlib.pyplot.plot(error(stirling_factorial(numbers)), 'ro')
matplotlib.pyplot.show()

两种近似计算方法的相对误差图示如下。连续直线代表Ramanujan近似方法，点状线代表斯特灵近似方法。如你所见，Ramanuja近似方法比斯特灵近似方法更精确。

用Cython求阶乘的近似值