如何用Java生成和Python Numpy一样的随机数序列？

一、问题提出

蒙特卡洛模拟在金融产品定价、估值、风险管理等方面有着广泛的应用。由于种种原因，这些定价、估值、风险管理方面的系统往往是由不同团队使用不同的语言开发的。因此，在实践中，有时会需要对各个系统的计算结果进行核对验证，以确保各系统的结果是一致的。但是，由于这些计算结果都是在随机模拟的基础上生成的，而不同开发语言产生的随机数序列又不相同，因此当结果出现差异时，很难分辨出是由于程序算法方面的原因造成的，还是由于随机数序列的不同造成的。

我最近就遇到了这样的问题，其中一个系统是Java开发的，另一个系统是Python开发的，使用了numpy中的随机数模块。在网上搜索发现，在stackoverflow上有人提出过类似的问题，但没有得到解答。因此我分析了一下相关的源代码，做了一个解决方案，以供有需要的同行参考。 

二、Python numpy源代码分析

   经查阅numpy文档，可知numpy生成随机数采用的是MersenneTwister算法，这是随机数生成的经典方法，来自于1998年ACM Transactions on Modeling and Computer Simulation中的一篇论文，在维基上有详细介绍（https://en.wikipedia.org/wiki/Mersenne_Twister），论文也可在此处下载。其基本思路是，根据给定的随机数种子，先通过一系列运算，生成624个整数key；再以这624个key为基础，生成一系列满足均匀分布的、在0到2^32-1之间的伪随机整数。其他程序可以对这些伪随机整数进行再加工，生成符合其他概率分布形态的随机数。

经阅读numpy源码，其生成随机数的核心程序为numpy/random/mtrand目录下的randomkit.c，其中几个主要的函数如下：

1、rk_seed(unsigned long seed, rk_state *state)：基于Mersenne Twister算法，初始化624个整数key，放在state结构体中。

2、rk_random(rk_state *state)：基于Mersenne Twister算法，生成一个32位的、服从均匀分布的随机数。

3、rk_double(rk_state *state)：生成一个0到1之间的、服从均匀分布的随机数。基本原理是：首先调用rk_random()函数生成两个32个随机数a、b，然后将a右移5位、保留27位，将b右移6位、保留26位，再将a、b按位串接在一起，形成一个53位的整数，再除以2^53。源代码如下：

    long a = rk_random(state)>> 5, b = rk_random(state) >> 6;
    return (a * 67108864.0 +b) / 9007199254740992.0;

   4、rk_gauss(rk_state*state)：生成一个服从正态分布N(0, 1)的随机数。这一部分numpy利用了Box–Muller变换公式，它是一种将均匀分布转换成正态分布的经典方法，在维基上有详细介绍（https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform）。Box-Muller变换有基本公式和极坐标公式两种形态，numpy采用的是极坐标公式，即：假设u, v是两个服从均匀分布的、0到1之间的随机数，则以下z0、z1服从正态分布N(0, 1)：

源代码如下：

       do {
           x1 = 2.0*rk_double(state) - 1.0;
           x2 = 2.0*rk_double(state) - 1.0;
           r2 = x1*x1 + x2*x2;
       } while (r2 >= 1.0 || r2 == 0.0);
       f = sqrt(-2.0*log(r2)/r2);
       state->gauss = f*x1;
       return f*x2;

三、Java apache math3源代码分析

经搜索，java社区的apache math3已经提供了一个Mersenne Twister算法的开源实现，于是下载源代码分析了一下。Apache math3的Mersenne Twister算法实现类为org.apache.commons.math3.random.MersenneTwister，其中主要方法有：

1、setSeed(intseed)：初始化624个整数key；

2、next(intbits)，生成不超过32位的伪随机整数。

经阅读源码，发现上述算法的具体实现虽然与numpy的C语言版本略有差异，但结果应该是一致的。经测试验证，果然如此，在给定相同seed情况下，两者生成的32位随机整数序列是相同的。

因此，numpy和math3随机数序列的差异，主要是由于对Mersenne Twister算法结果的后续加工处理方法不同。Apache math3生成其他形态随机数的方法在MersenneTwister的父类BitsStreamGenerator中定义，其中主要方法有：

1、nextDouble()：生成一个0到1之间的、服从均匀分布的随机数。基本原理是：首先调用MersenneTwister的next()函数生成两个26位的随机整数a、b，再将a、b按位串接在一起，形成一个52位的整数，再除以2^52。源代码如下：

        final long high = ((long) next(26))<< 26;
        final int  low  =next(26);
        return (high | low) * 0x1.0p-52d;

由此可见，math3和numpy的差异在于，math3使用的是52位随机数，而numpy使用的是53位随机数。

2、nextGaussian()：生成一个服从正态分布N(0,1)的随机数。这一部分math3也是利用了Box–Muller变换公式，但是与numpy不同的是，它没有采用Box–Muller极坐标公式，而是使用Box–Muller基本公式，即：假设x, y是两个服从均匀分布的、0到1之间的随机数，则以下z0、z1服从正态分布N(0, 1)：

源代码如下：

            final double x = nextDouble();
            final double y = nextDouble();
            final double alpha = 2 *FastMath.PI * x;
            final double r = FastMath.sqrt(-2 * FastMath.log(y));
            random = r * FastMath.cos(alpha);
            nextGaussian = r * FastMath.sin(alpha);

四、解决方案

基于以上分析，我们可以在自己的Java代码中继承MersenneTwister类，覆盖其中的nextDouble、nextGaussian等方法，将其改为与numpy一致的算法即可。源代码如下：

public class NumpyRandom extends MersenneTwister {

	private boolean has_gauss = false;
	private double gauss = 0;

	public double nextDouble() {
		int a = next(32) >>> 5;
		int b = next(32) >>> 6;
		return (a * 67108864.0 + b) / 9007199254740992.0;
	}

	public double nextGaussian() {
		if (has_gauss) {
			double tmp = gauss;
			gauss = 0;
			has_gauss = false;
			return tmp;
		} else {
			double x1, x2, s;

			do {
				x1 = 2.0 * nextDouble() - 1.0;
				x2 = 2.0 * nextDouble() - 1.0;
				s = x1 * x1 + x2 * x2;
			} while (s >= 1.0 || s == 0.0);

			double f = Math.sqrt(-2.0 * Math.log(s) / s);
			gauss = f * x1;
			has_gauss = true;
			return f * x2;
		}
	}

	public static void main(String[] args) {
		NumpyRandom rand = new NumpyRandom();
		rand.setSeed(1000);
		for (int i = 0; i < 10; i++)
			System.out.println(" " + rand.nextGaussian());
	}
}

经测试，上述程序产生的正态分布随机数与python numpy的结果完全相同。