每日算法 — 使用java实现伪随机数生成器:线性同余法与梅森旋转算法

伪随机数生成器(Pseudorandom Number Generator, PRNG)是计算机科学中最基础也是最重要的算法之一。从游戏开发中的道具掉落、机器学习中的权重初始化,到密码学中的密钥生成和蒙特卡洛模拟中的采样,几乎所有领域都需要高质量的随机数。然而,计算机是确定性机器,真正的随机性必须通过外部物理熵源获取。因此,计算机中使用的大多是”伪随机数”——通过确定性算法产生、但在统计特性上近似于真随机数的数列。本文将用Java实现两种经典的PRNG算法:线性同余法(LCG)和梅森旋转算法(MT19937),并分析它们的周期、统计特性和适用场景。

线性同余法(LCG)

线性同余法由Lehmer于1951年提出,是最简单、最经典的PRNG之一。其递推公式为:

$$X_{n+1} = (a \cdot X_n + c) \mod m$$

其中:
– $m$ 为模数(modulus),决定随机数范围
– $a$ 为乘数(multiplier),$0 < a < m$
– $c$ 为增量(increment),$0 \leq c < m$
– $X_0$ 为种子(seed)

当参数满足Hull-Dobell定理的条件时,LCG可以达到最大周期 $m$。

梅森旋转算法(MT19937)

梅森旋转算法由松本和西山于1997年提出,是当今最广泛使用的PRNG之一,被Python、Ruby、PHP、MATLAB等主流语言采用。MT19937的名称来源于其周期为梅森素数 $2^{19937} – 1$,约为 $10^{6001}$,对于任何实际应用而言都是”无限”的。

算法核心是一个624维的状态向量,通过位运算(异或、与、移位)进行状态转换,产生高质量的32位随机整数。其统计特性通过了严格的Diehard和TestU01测试套件。

Java实现

import java.util.*;

/**
 * 伪随机数生成器实现
 * 包含线性同余法(LCG)和梅森旋转算法(MT19937)
 */
public class RandomNumberGenerator {

    // ==================== 线性同余法 ====================

    /**
     * 线性同余生成器(LCG)
     * 使用 glibc 的标准参数:a=1103515245, c=12345, m=2^31
     */
    public static class LCG {
        private long state;
        private static final long A = 1103515245L;
        private static final long C = 12345L;
        private static final long M = 1L << 31;  // 2^31

        public LCG(long seed) {
            this.state = seed;
        }

        /**
         * 生成下一个随机整数 [0, 2^31)
         */
        public int nextInt() {
            state = (A * state + C) % M;
            return (int) state;
        }

        /**
         * 生成 [0, bound) 范围内的随机整数
         */
        public int nextInt(int bound) {
            return Math.abs(nextInt()) % bound;
        }

        /**
         * 生成 [0.0, 1.0) 范围内的随机浮点数
         */
        public double nextDouble() {
            return nextInt() / (double) M;
        }
    }

    // ==================== 梅森旋转算法 MT19937 ====================

    /**
     * MT19937实现
     * 周期为 2^19937 - 1 的高质量伪随机数生成器
     */
    public static class MT19937 {
        private static final int N = 624;           // 状态向量长度
        private static final int M = 397;           //  Twist偏移量
        private static final long MATRIX_A = 0x9908b0dfL;  // 常量矩阵
        private static final long UPPER_MASK = 0x80000000L; // 最高位掩码
        private static final long LOWER_MASK = 0x7fffffffL; // 低31位掩码

        private long[] mt = new long[N];  // 状态向量
        private int mti = N + 1;          // 当前状态索引,N+1表示未初始化

        public MT19937(long seed) {
            init(seed);
        }

        /**
         * 初始化状态向量
         */
        private void init(long seed) {
            mt[0] = seed & 0xffffffffL;
            for (mti = 1; mti < N; mti++) {
                // 递推初始化:mt[i] = (1812433253 * (mt[i-1] ^ (mt[i-1] >> 30)) + i) & 0xffffffff
                mt[mti] = (1812433253L * (mt[mti - 1] ^ (mt[mti - 1] >>> 30)) + mti) & 0xffffffffL;
            }
        }

        /**
         * 生成下一个32位随机整数
         */
        public int nextInt() {
            long y;
            long[] mag01 = {0x0L, MATRIX_A};

            // Twist:当所有状态用完时重新生成
            if (mti >= N) {
                int kk;
                for (kk = 0; kk < N - M; kk++) {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
                    mt[kk] = mt[kk + M] ^ (y >>> 1) ^ mag01[(int)(y & 0x1L)];
                }
                for (; kk < N - 1; kk++) {
                    y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
                    mt[kk] = mt[kk + (M - N)] ^ (y >>> 1) ^ mag01[(int)(y & 0x1L)];
                }
                y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ mag01[(int)(y & 0x1L)];
                mti = 0;
            }

            y = mt[mti++];

            // Tempering:通过位运算提高随机性质量
            y ^= (y >>> 11);
            y ^= (y << 7) & 0x9d2c5680L;
            y ^= (y << 15) & 0xefc60000L;
            y ^= (y >>> 18);

            return (int) y;
        }

        /**
         * 生成 [0, bound) 范围内的随机整数
         */
        public int nextInt(int bound) {
            if (bound <= 0) throw new IllegalArgumentException("bound必须为正数");
            // 使用取模法,注意避免边界偏差
            int r = nextInt();
            return Math.abs(r) % bound;
        }

        /**
         * 生成 [0.0, 1.0) 范围内的随机浮点数
         */
        public double nextDouble() {
            // 使用53位精度(双精度浮点数的尾数位数)
            long a = nextInt() >>> 5;
            long b = nextInt() >>> 6;
            return ((a << 26) + b) / (double)(1L << 53);
        }
    }

    // ==================== 统计测试 ====================

    /**
     * 卡方检验:测试随机数的均匀分布性
     */
    public static double chiSquareTest(int[] frequencies, int totalSamples) {
        double expected = (double) totalSamples / frequencies.length;
        double chiSquare = 0.0;
        for (int freq : frequencies) {
            double diff = freq - expected;
            chiSquare += (diff * diff) / expected;
        }
        return chiSquare;
    }

    /**
     * 运行统计测试
     */
    public static void runStatisticalTests(String name, java.util.function.IntSupplier rng, int samples) {
        int buckets = 100;
        int[] freq = new int[buckets];

        for (int i = 0; i < samples; i++) {
            int val = Math.abs(rng.getAsInt()) % buckets;
            freq[val]++;
        }

        double chiSquare = chiSquareTest(freq, samples);
        System.out.println("\n=== " + name + " 统计测试 ===");
        System.out.println("样本数: " + samples);
        System.out.println("卡方统计量: " + String.format("%.4f", chiSquare));
        System.out.println("期望频数: " + (samples / buckets));
        System.out.print("前20桶频数: ");
        for (int i = 0; i < 20; i++) {
            System.out.print(freq[i] + " ");
        }
        System.out.println();
    }

    public static void main(String[] args) {
        int samples = 100000;
        long seed = 12345L;

        // 测试LCG
        LCG lcg = new LCG(seed);
        runStatisticalTests("LCG", lcg::nextInt, samples);

        // 测试MT19937
        MT19937 mt = new MT19937(seed);
        runStatisticalTests("MT19937", mt::nextInt, samples);

        // 演示:用MT19937生成随机浮点数
        System.out.println("\n=== MT19937随机浮点数示例 ===");
        MT19937 demo = new MT19937(42);
        for (int i = 0; i < 10; i++) {
            System.out.printf("%.6f%n", demo.nextDouble());
        }

        // 演示:洗牌算法
        System.out.println("\n=== 洗牌算法示例 ===");
        List<Integer> cards = new ArrayList<>();
        for (int i = 1; i <= 10; i++) cards.add(i);
        MT19937 shuffleRng = new MT19937(999);
        // Fisher-Yates洗牌
        for (int i = cards.size() - 1; i > 0; i--) {
            int j = shuffleRng.nextInt(i + 1);
            Collections.swap(cards, i, j);
        }
        System.out.println("洗牌结果: " + cards);
    }
}

算法特性对比

特性 LCG MT19937
周期 $2^{31}$(约20亿) $2^{19937}-1$(约 $10^{6001}$)
速度 极快(单次模乘) 快(批量生成后缓存)
内存占用 4字节 2.5KB
统计质量 一般(低位相关性高) 优秀(通过严格统计测试)
种子敏感性 一般 优秀
密码学安全 不安全 不安全(可预测)
适用场景 简单模拟、游戏 科学计算、机器学习

密码学安全随机数

需要特别强调的是,LCG和MT19937都不适用于密码学场景,因为它们的内部状态可以从输出序列中反推。密码学安全随机数生成器(CSPRNG)需要使用专门的算法,如Java的 SecureRandom(底层使用SHA-1或NativePRNG),或密码学标准算法如Yarrow、Fortuna等。

梅森旋转的数学原理

MT19937的核心是定义在二元域 $\mathbb{F}_2$ 上的线性递推关系。状态向量的更新可以看作一个 $19937 \times 19937$ 的矩阵乘法,周期等于矩阵的阶。由于 $19937$ 是梅森素数,矩阵被精心设计为使周期恰好为 $2^{19937}-1$,这是梅森素数在PRNG设计中的美妙应用。

理解伪随机数生成器的原理,不仅能帮助开发者选择合适的随机数源,还能避免因随机数质量问题导致的程序缺陷——例如使用低质量PRNG进行蒙特卡洛积分时可能产生严重偏差。