相关疑难解决方法(0)

我正在编写一个 cython 应用程序，我需要在紧密的嵌套循环中动态生成高斯随机变量。我想在不引入任何额外依赖项（例如 GSL）的情况下执行此操作。

对于我目前能够使用统一随机数即时执行此操作的最小版本：

from libc.stdlib cimport rand, RAND_MAX
import numpy as np

cdef double random_uniform():
    cdef double r = rand()
    return r/RAND_MAX

def my_function(int n):
    cdef int i
    cdef double[:] result = np.zeros(n, dtype='f8', order='C')
    for i in range(n):
        result[i] = random_uniform()
    return result

上面的代码在功能上等同于 numpy.random.rand(n)，并且可以使用以下最小设置文件进行编译：

from distutils.core import setup
from Cython.Build import cythonize
import numpy as np

setup(ext_modules=cythonize("example.pyx"), include_dirs=[np.get_include()])

# compile instructions:
# python setup.py build_ext --inplace

为了回答这个问题，我正在寻找的是与 np.random.randn(n) 功能等效的相同类型的最小解决方案，同样理想的是出于可移植性的原因直接从 libc.stdlib 导入任何依赖项。

维基百科条目上有Box-Muller 算法的 …

生成伪统一随机数（[0，1中的双精度数]）的最佳方法是：

1. 跨平台（理想情况下具有相同的样品序列）
2. 线程安全（显式传递prng的突变状态或在内部使用线程局部状态）
3. 没有GIL锁
4. 可轻松包装在Cython中

3年前也有类似的文章，但是很多答案并不符合所有标准。例如，drand48特定于POSIX。

我知道的（但不确定）满足所有某些条件的唯一方法是：

from libc.stdlib cimport rand, RAND_MAX

random = rand() / (RAND_MAX + 1.0)

注意@ogrisel 大约在3年前问了同样的问题。

编辑

调用rand不是线程安全的。感谢您指出@DavidW。