尧图网站开发 
简介这套资源提供五种随机数发生器的C与MATLAB实现适合需要学习伪随机数算法、进行模拟或数据分析的开发者与学生。涵盖平方取中法、乘积取中法、Mersenne Twister、ISAAC和PCG并配有平方取中、线性同余、组合LCG等具体代码文件。压缩包共12个文件包含6个m脚本、1个cpp源文件、编译依赖文件、可执行程序及Code::Blocks工程文件整体约320KB便于直接查看算法实现与运行效果。已有228人学习下载。通过对比这两种语言下的实现读者不仅可以掌握不同随机数生成器的原理与参数选择还能了解C标准库 和MATLAB内置随机函数之外的自主构造方法为编写可靠模拟程序和进一步研究随机算法打下基础。示例中展示了多种生成器的核心结构与中间结果代码文件分类清晰适合入门到进阶循序渐进地实践。1. 五种随机数发生器为什么一套模拟里要塞五种算法很多人写模拟时拿rand()一把梭种子设好循环跑完结果复现不出来就去怀疑编译器版本。等到把 C 的rand()换成库里的梅森旋转再和 MATLAB 对拍时又会发现同一颗种子在两边的输出完全不同模拟结果又翻了一轮。其实随机数发生器是个有脾气的组件周期、内存、可复现性、并发下的独立性都不一样选错一个后面所有统计结论都要跟着遭殃。下面这套方案不是只给五种能跑的代码而是把 C 侧的模板类直接放进业务代码、MATLAB 侧的 classdef 用来做快速原型和数据复用两种语言各自实现同一套算法方便交叉验证。五种算法覆盖了最快的、最省内存的、最稳定的以及最不该在生产里用的教学案例。看完后你会愿意花十分钟把项目里的随机数发生器显式替换掉而不是继续依赖某个编译器自带的全局状态。2. 五种随机数发生器的 C 与 MATLAB 实现从线性同余到 MT199372.1 LCG线性同余发生器用 3 行代码建立随机数发生器线性同余发生器是所有伪随机数发生器里最朴素的一种递推式只有一行x_{n1} (a * x_n c) mod ma是乘数c是增量m是模数。当c 0时它退化为乘法发生器所以 LCG 也是后面 Park-Miller 的基础。它的特点是状态量只有一个整数保存现场、并行分片都非常方便缺点是周期最多只有m而且连续随机点的散布图里存在明显的格点结构。下面这份 C 代码直接使用 32 位无符号整数的自然回绕把m 2^32的取模运算省掉// LCG32.h #pragma once #include cstdint class LCG32 { uint32_t x; // 当前状态也是唯一的种子状态 public: explicit LCG32(uint32_t seed 1) : x(seed) {} uint32_t next() { x 1664525u * x 1013904223u; // 溢出回绕等价于 mod 2^32 return x; } double unit() { // 取高 24 位转 [0,1)避开低位的短周期问题 return (next() 8) * (1.0 / 16777216.0); } };逻辑说明无符号整数乘法溢出后自动截断mod 2^32由编译器完成不需要手写%。unit()刻意右移 8 位而不是直接除以 2^32原因在于 LCG 的低位比特周期较短比如最低位的周期是 2取高位才能保留均匀性。同样的行为写成 MATLAB 类需要在每个赋值处显式转uint32语法上比 C 啰嗦但思路一致% LCG32.m classdef LCG32 handle properties x uint32 end methods function obj LCG32(seed) if nargin 1 seed uint32(1); end obj.x uint32(seed); end function r next(obj) obj.x uint32(1664525) * obj.x uint32(1013904223); r obj.x; end function u unit(obj) r obj.next(); u double(bitshift(r, -8)) / 16777216.0; end end end参数说明a 1664525, c 1013904223是 Numerical Recipes 里推荐的参数组合对模数2^32满足完整周期条件即c与m互质、a-1能被 4 整除。如果自己随意换成别的数值周期会从2^32缩水到几万甚至几百均匀性也会跟着崩掉。这两组参数在 C 与 MATLAB 两侧是同一份后文做跨语言对拍时以它作为起点最稳。2.2 Park-Miller 乘式发生器避开大素数求模的溢出陷阱Park-Miller 是 LCG 的一个著名特例c 0, m 2147483647递推式为x_{n1} 16807 * x_n mod 2147483647。由于模数m是梅森素数2^31 - 1它天然拥有m - 1的满周期前提是种子不能落在0上。这个算法在上世纪八十年代是 Park 和 Miller 推荐的“最低标准”发生器很多老系统的rand()都沿用它的变体。直接写16807 * x % 2147483647看似简单但在 32 位平台上有隐患16807 * x的最大值接近16807 * 2147483646 ≈ 3.6 × 10^13远超 32 位整数范围。现代 PC 上用uint64_t可以一了百了// ParkMiller.h #include cstdint class ParkMiller { uint64_t x; public: explicit ParkMiller(uint64_t seed 1) : x(seed ? seed : 1) {} uint32_t next() { x (x * 16807ull) % 2147483647ull; return static_castuint32_t(x); } double unit() { return (next() 8) * (1.0 / 16777216.0); } };如果目标平台连 64 位乘法都嫌贵可以用 Schrage 分解把一次大数乘法拆成小整数运算避免中间结果溢出uint32_t schrage_next(uint32_t x) { const uint32_t a 16807, m 2147483647; const uint32_t q m / a; // 127773 const uint32_t r m % a; // 2836 int32_t t static_castint32_t(a * (x % q) - r * (x / q)); if (t 0) t m; return static_castuint32_t(t); }逻辑说明Schrage 方法利用m a * q r的恒等式把x拆成x q * (x / q) (x % q)从而让每一步乘法都落在 31 位整数范围内。它解决的正是 32 位平台上做伪随机数发生器最容易踩的溢出坑结果对了但中间乘积先溢出了。MATLAB 侧不需要这么小心因为uint64是内置类型可以照搬第一版 C 代码state uint64(1); state mod(state * uint64(16807), uint64(2147483647));2.3 平方取中法演示原理可以生产环境别碰平方取中法由冯·诺依曼在 1940 年代提出思路很直观取当前数的平方再截取平方结果中间的几位作为下一个状态。比如四位种子1234平方是1522756补零成八位01522756取中间四位得到5227再平方、再取中间如此循环。用 C 写一个十进制演示版本// middle_square_demo.cpp #include cstdint #include cstdio uint32_t middle_square(uint32_t x, int digits 4) { uint64_t sq static_castuint64_t(x) * x; uint64_t mask 1; for (int i 0; i digits; i) mask * 10; // 平方补零后取中间 digits 位等价于先除以 10 的幂再取模 return static_castuint32_t((sq / (mask / 10)) % mask); } int main() { uint32_t x 1234; for (int i 0; i 20; i) { x middle_square(x); std::printf(%04u\n, x); } return 0; }这里用了十进制取位方便肉眼观察序列退化。实际工程如果真要实现建议改用二进制位操作把 32 位状态平方后取中间的 16 位或 32 位。无论哪种写法它都有一个致命问题序列很容易进入循环甚至收敛到常数。种子换成0000后会立刻变成全零后面每一次平方取中间都是零。保留这个算法的意义在于教学它能很直观地解释“为什么随机数发生器需要精心设计状态转移”以及“为什么不能拿上一轮结果简单算一算就假装随机”。部署到线上代码前请先把这一节划掉。2.4 xorshift位运算型高性能随机数发生器xorshift 是 George Marsaglia 在 2003 年提出的一族位运算发生器。名字来源于它的核心操作异或xor加移位shift。没有乘法没有取模只靠位运算所以在嵌入式平台和游戏引擎里非常流行。32 位版本的状态转移写成 C 只有三行// Xorshift32.h #include cstdint class Xorshift32 { uint32_t x; public: explicit Xorshift32(uint32_t seed 1) : x(seed ? seed : 1) {} uint32_t next() { x ^ x 13; x ^ x 17; x ^ x 5; return x; } double unit() { return (next() 8) * (1.0 / 16777216.0); } };MATLAB 里没有运算符重载要用bitxor和bitshift显式表达% Xorshift32.m classdef Xorshift32 handle properties x uint32 end methods function obj Xorshift32(seed) if nargin 1 seed uint32(1); end obj.x uint32(seed); end function r next(obj) r obj.x; r bitxor(r, bitshift(r, 13, uint32)); r bitxor(r, bitshift(r, -17, uint32)); r bitxor(r, bitshift(r, 5, uint32)); obj.x r; end function u unit(obj) r obj.next(); u double(bitshift(r, -8)) / 16777216.0; end end end参数说明13, 17, 5这一组移位量是经过验证的三元组组合起来能保证2^32 - 1的满周期。不要随手改成16, 8, 3之类的数位运算发生器对移位参数极其敏感改错一个数周期可能掉到几千万测试时还不容易发现。xorshift 的另一个优势是状态极小且保存方便并行场景下给每个线程分一颗互不相同的种子内存开销几乎可以忽略。2.5 MT19937C 与 MATLAB 各自的标准随机数发生器梅森旋转Mersenne TwisterMT19937是目前主流语言默认随机数发生器的事实标准。名字里的 19937 指的是它的周期是2^19937 - 1并且保证 623 维均匀分布。C11 把它收进randomMATLAB 的rng默认用的也是它。C 侧最小用法#include random #include iostream int main() { std::mt19937 gen(42); // 指定 32 位种子 std::uniform_real_distributiondouble dist(0.0, 1.0); for (int i 0; i 5; i) { std::cout dist(gen) \n; } return 0; }MATLAB 侧对应的最小写法rng(42, twister); % 显式指定梅森旋转算法 u rand(1, 5);这里有一个容易踩的细节两边虽然都是 MT19937生成的内部整数状态序列可以做到一致但 C 的uniform_real_distribution和 MATLAB 的rand在把内部整数映射到[0,1)浮点数时用了不同的策略所以double输出不完全相等。跨语言对拍时不要直接比较浮点数先比较整数状态或者老老实实用自己实现的算法。2.6 五种随机数发生器的周期与适用场景对比算法状态大小周期适用场景注意点LCG4 字节最多2^32嵌入式、教学、跨语言对拍低位周期短取高 24 位Park-Miller8 字节2^31 - 2老系统兼容、低功耗设备种子不能为 0防溢出平方取中法不定极短且退化仅供课堂演示不要用于生产xorshift324 字节2^32 - 1游戏、并行分片、高性能场景移位参数必须经过验证MT199372.5 KB2^19937 - 1仿真、统计计算、默认选择状态大恢复现场成本高3. 种子、周期与参数选择五种随机数发生器用错参数的后果3.1 C 与 MATLAB 的种子机制直接给整数还是给状态向量C 的std::mt19937接受一个 32 位整数做种子内部再用初始化算法把这个整数展开成 624 个 32 位状态MATLAB 的rng(seed, twister)也做类似的事情但它的种子会附带一个内部 shuffle 步骤。结论是就算两边种子整数相同、调用次数相同得到的整数序列也不保证逐位一致。想让 C 与 MATLAB 的随机数发生器输出完全对齐最简单的方法是选择状态转移简单、不依赖复杂初始化的算法比如前面写的 LCG32 和 Xorshift32。给随机数发生器设置种子时还有一个经常被忽略的规则种子不能是 0。Park-Miller 和 xorshift 在实现里都做了seed seed ? seed : 1的保护原因很简单这两种算法的状态为 0 时整个序列会永久停留在 0。LCG 相对宽容一些但如果c和种子凑巧让状态落入某个短循环同样会造成周期缩短。3.2 溢出与取模LCG 参数一旦选错周期直接缩水先看一个反例m 2^32, a 65539, c 0这是老的 RANDU 算法它在三维空间中的随机点会全部落在一组平面上周期和均匀性都极差。另一个反面教材是把a写成偶数且c 0这时状态永远不会出现奇数随机数的信息量直接少了一位。完整的周期条件只有三个c与m互质a - 1能被m的所有质因子整除如果m是 4 的倍数a - 1也必须是 4 的倍数。工程上没时间验证这些条件时直接用 2.1 节给出的1664525和1013904223这两个值已经被反复验证过完整周期成立。溢出问题在 C 侧最隐蔽。看下面这行代码uint64_t bad_next(uint32_t x) { return (x * 16807) % 2147483647; // x 先被提升为 int可能溢出 }x是uint32_t16807是int两者求积时uint32_t会先把int提升为uint32_t但如果x来自一个有符号变量情况就会变成int乘法先溢出再转换。修复方式是给常量加后缀return (x * 16807u) % 2147483647u;类型后缀看起来是小事但在随机数发生器里一个中间溢出的整数乘法足以让整条序列从某个位置开始错位。对拍脚本通常要跑到几十万步才能发现错位那时候排查成本就很高了。3.3 状态恢复随机数发生器的序列如何保存与续跑蒙特卡洛模拟跑了一半宕机需要从断点继续这时不能只保存种子要保存生成器当前状态。MT19937 的状态长达 2.5 KB但标准库提供了序列化能力LCG 和 xorshift 只需要保存一个整数。C 里可以用std::ostream把std::mt19937的状态整体写出#include random #include fstream std::mt19937 gen(42); std::ofstream out(mt_state.bin, std::ios::binary); out gen; // 把 624 个状态字整体写入文件MATLAB 侧对应的做法是用rng输出结构体s rng(42, twister); save(mt_state.mat, s); % 恢复时load(mt_state.mat); rng(s);提示MATLAB 的rng结构体里包含Type,Seed,State三个字段跨 MATLAB 版本恢复时一般兼容但不要手工去改State数组那会让生成器进入未定义状态。续跑模拟的正确姿势是把状态快照和业务数据放在同一份日志里而不是重新跑一遍再赌一次随机序列。4. 怎么看五种随机数发生器是否够随机均匀性检验4.1 先跑均值与方差这步能筛掉一半弱发生器拿到一个随机数发生器我通常先让它生成 100 万个[0,1)样本看均值和方差。均匀分布的均值理论值是0.5方差理论值是1/12 ≈ 0.08333。C 侧统计代码写起来很快#include vector #include numeric double mean(const std::vectordouble v) { return std::accumulate(v.begin(), v.end(), 0.0) / v.size(); } double variance(const std::vectordouble v, double m) { double s 0.0; for (double x : v) s (x - m) * (x - m); return s / v.size(); }用这个函数分别喂给五种发生器的unit()如果均值偏差超过0.001说明低位截断或者状态转移有问题。平方取中法在这个测试里会立刻现形它可能一时均值正常但序列长度一旦进入循环方差就会异常。4.2 卡方检验的 C 与 MATLAB 实现均值方差只能看一阶矩分布形态要用卡方检验。做法是把[0,1)区间等分成 64 个箱统计每个箱里的样本数再比较观测频数和理论频数。C 侧计算卡方统计量#include vector #include cmath double chi2_stat(const std::vectorsize_t counts) { size_t total 0; for (size_t c : counts) total c; double expected static_castdouble(total) / counts.size(); double stat 0.0; for (size_t c : counts) { double d static_castdouble(c) - expected; stat d * d / expected; } return stat; }逻辑说明64 个箱自由度是 63卡方统计量应该在 63 附近波动。如果统计量超过 100基本可以判定分布不均匀。平方取中法和参数错误的 LCG 在这一步过不去。MATLAB 有现成的chi2gof用法更直接rng(42, twister); x rand(1, 1000000); [h, p] chi2gof(x, Edges, linspace(0, 1, 65));参数说明Edges定义分箱边界linspace(0, 1, 65)把区间切成 64 个等宽箱。h 0表示不能拒绝均匀性假设p 0.05是常见参考线。注意[0,1)区间的样本很少落在 1 上边界问题可以忽略但如果用[0,1]闭区间生成器要确认最后一箱是否被异常地填满。4.3 序列相关性与交叉验证分布均匀不等于序列独立。一个简单的相关性测试是计算相邻样本的相关系数理想值接近 0。C 可以用一阶自相关公式更省事的是在 MATLAB 里画 scatter 图rng(42, twister); u1 rand(1, 1000); u2 rand(1, 1000); scatter(u1(1:end-1), u1(2:end), 3, filled);视觉上LCG 会在散点图里出现规则网格xorshift 通常看不出结构MT19937 在 623 维以内都保持均匀。把这五种发生器的散点图并排放一起会非常直观地理解“格点结构”和“高维均匀性”的差别。4.4 跨语言对拍脚本验证两套实现是否真的等价当 C 和 MATLAB 各写了一个版本时最值得做的是逐位对拍。将对拍脚本放在 CI 里每次改动随机数发生器就能立刻发现不一致# 生成 C 侧前 10000 个整数序列 ./lcg_dump 7 cpp_seq.txt # 生成 MATLAB 侧序列 matlab -batch lLCG32(7); for i1:10000, fprintf(%u\n, l.next()); end matlab_seq.txt # 对比 diff cpp_seq.txt matlab_seq.txt没有输出就代表一致。这个做法能验证两个实现的算法逻辑是否真正同步比肉眼对比浮点数可靠得多。5. 让 C 与 MATLAB 跑出同一组随机数序列以及选型技巧5.1 为什么不要直接对拍 MT19937 的浮点输出很多人第一反应是让 C 和 MATLAB 都调用标准 MT19937用同一粒种子对拍。实际对上后会发现前几个 float 就不一样两边虽然核心算法都源自松本与西村的设计但标准库的uniform_real_distribution和 MATLAB 的rand在浮点映射时采取了不同的策略导致输出不可比。跨语言对拍请直接从内部整数状态下手。MT19937 的std::mt19937::operator()和 MATLAB 的RandStream底层状态可以导出但操作起来不如自写算法直观。5.2 同一算法、同一参数、同一状态最省事的跨语言对齐方案手动实现的 LCG 是最容易做到跨语言完全一致的。C 侧直接用无符号 32 位回绕MATLAB 侧用uint32bitshift两边只要种子相同、调用次数相同输出的整数就完全相同。实际中我会把 LCG 的next()结果前 8 个值打印出来作为回归测试的 golden 基线C 和 MATLAB 各跑一份CI 里自动 diff。// lcg_dump.cpp #include cstdio #include LCG32.h int main() { LCG32 lcg(7); for (int i 0; i 8; i) { std::printf(%u\n, lcg.next()); } return 0; }对应 MATLABl LCG32(7); for i 1:8 fprintf(%u\n, l.next()); end这两份代码的输出可以完全相同。原因在于它们做的是同一组无符号整数的乘法与截断不涉及浮点舍入也不涉及不同标准库的实现差异。这个技巧特别适合给跨语言数据管道做随机数对照实验时使用。5.3 生产项目里的选型倾向我一般在项目里这样选默认用 MT19937因为周期足够长、统计性质经过大规模验证并行模拟里改用 xorshift每个线程占用 4 字节状态分片代价极低嵌入式场景用 Park-Miller整型运算在低功耗芯片上表现稳定LCG 只用于跨语言对拍和教学平方取中法永远不碰。最后留一个可执行的小技巧下一次接手的项目如果还在裸调rand()先把全局调用改成显式std::mt19937或 MATLABrng(seed, twister)再跑一遍上面的对拍脚本没有异常再继续调业务逻辑。这两种语言虽然标准库不同但你手上已经有了五种可用的实现和一套验证脚本替换成本其实很低。本文还有配套的精品资源点击获取 
简介这套资源提供五种随机数发生器的C与MATLAB实现适合需要学习伪随机数算法、进行模拟或数据分析的开发者与学生。涵盖平方取中法、乘积取中法、Mersenne Twister、ISAAC和PCG并配有平方取中、线性同余、组合LCG等具体代码文件。压缩包共12个文件包含6个m脚本、1个cpp源文件、编译依赖文件、可执行程序及Code::Blocks工程文件整体约320KB便于直接查看算法实现与运行效果。已有228人学习下载。通过对比这两种语言下的实现读者不仅可以掌握不同随机数生成器的原理与参数选择还能了解C标准库 和MATLAB内置随机函数之外的自主构造方法为编写可靠模拟程序和进一步研究随机算法打下基础。示例中展示了多种生成器的核心结构与中间结果代码文件分类清晰适合入门到进阶循序渐进地实践。1. 五种随机数发生器为什么一套模拟里要塞五种算法很多人写模拟时拿rand()一把梭种子设好循环跑完结果复现不出来就去怀疑编译器版本。等到把 C 的rand()换成库里的梅森旋转再和 MATLAB 对拍时又会发现同一颗种子在两边的输出完全不同模拟结果又翻了一轮。其实随机数发生器是个有脾气的组件周期、内存、可复现性、并发下的独立性都不一样选错一个后面所有统计结论都要跟着遭殃。下面这套方案不是只给五种能跑的代码而是把 C 侧的模板类直接放进业务代码、MATLAB 侧的 classdef 用来做快速原型和数据复用两种语言各自实现同一套算法方便交叉验证。五种算法覆盖了最快的、最省内存的、最稳定的以及最不该在生产里用的教学案例。看完后你会愿意花十分钟把项目里的随机数发生器显式替换掉而不是继续依赖某个编译器自带的全局状态。2. 五种随机数发生器的 C 与 MATLAB 实现从线性同余到 MT199372.1 LCG线性同余发生器用 3 行代码建立随机数发生器线性同余发生器是所有伪随机数发生器里最朴素的一种递推式只有一行x_{n1} (a * x_n c) mod ma是乘数c是增量m是模数。当c 0时它退化为乘法发生器所以 LCG 也是后面 Park-Miller 的基础。它的特点是状态量只有一个整数保存现场、并行分片都非常方便缺点是周期最多只有m而且连续随机点的散布图里存在明显的格点结构。下面这份 C 代码直接使用 32 位无符号整数的自然回绕把m 2^32的取模运算省掉// LCG32.h #pragma once #include cstdint class LCG32 { uint32_t x; // 当前状态也是唯一的种子状态 public: explicit LCG32(uint32_t seed 1) : x(seed) {} uint32_t next() { x 1664525u * x 1013904223u; // 溢出回绕等价于 mod 2^32 return x; } double unit() { // 取高 24 位转 [0,1)避开低位的短周期问题 return (next() 8) * (1.0 / 16777216.0); } };逻辑说明无符号整数乘法溢出后自动截断mod 2^32由编译器完成不需要手写%。unit()刻意右移 8 位而不是直接除以 2^32原因在于 LCG 的低位比特周期较短比如最低位的周期是 2取高位才能保留均匀性。同样的行为写成 MATLAB 类需要在每个赋值处显式转uint32语法上比 C 啰嗦但思路一致% LCG32.m classdef LCG32 handle properties x uint32 end methods function obj LCG32(seed) if nargin 1 seed uint32(1); end obj.x uint32(seed); end function r next(obj) obj.x uint32(1664525) * obj.x uint32(1013904223); r obj.x; end function u unit(obj) r obj.next(); u double(bitshift(r, -8)) / 16777216.0; end end end参数说明a 1664525, c 1013904223是 Numerical Recipes 里推荐的参数组合对模数2^32满足完整周期条件即c与m互质、a-1能被 4 整除。如果自己随意换成别的数值周期会从2^32缩水到几万甚至几百均匀性也会跟着崩掉。这两组参数在 C 与 MATLAB 两侧是同一份后文做跨语言对拍时以它作为起点最稳。2.2 Park-Miller 乘式发生器避开大素数求模的溢出陷阱Park-Miller 是 LCG 的一个著名特例c 0, m 2147483647递推式为x_{n1} 16807 * x_n mod 2147483647。由于模数m是梅森素数2^31 - 1它天然拥有m - 1的满周期前提是种子不能落在0上。这个算法在上世纪八十年代是 Park 和 Miller 推荐的“最低标准”发生器很多老系统的rand()都沿用它的变体。直接写16807 * x % 2147483647看似简单但在 32 位平台上有隐患16807 * x的最大值接近16807 * 2147483646 ≈ 3.6 × 10^13远超 32 位整数范围。现代 PC 上用uint64_t可以一了百了// ParkMiller.h #include cstdint class ParkMiller { uint64_t x; public: explicit ParkMiller(uint64_t seed 1) : x(seed ? seed : 1) {} uint32_t next() { x (x * 16807ull) % 2147483647ull; return static_castuint32_t(x); } double unit() { return (next() 8) * (1.0 / 16777216.0); } };如果目标平台连 64 位乘法都嫌贵可以用 Schrage 分解把一次大数乘法拆成小整数运算避免中间结果溢出uint32_t schrage_next(uint32_t x) { const uint32_t a 16807, m 2147483647; const uint32_t q m / a; // 127773 const uint32_t r m % a; // 2836 int32_t t static_castint32_t(a * (x % q) - r * (x / q)); if (t 0) t m; return static_castuint32_t(t); }逻辑说明Schrage 方法利用m a * q r的恒等式把x拆成x q * (x / q) (x % q)从而让每一步乘法都落在 31 位整数范围内。它解决的正是 32 位平台上做伪随机数发生器最容易踩的溢出坑结果对了但中间乘积先溢出了。MATLAB 侧不需要这么小心因为uint64是内置类型可以照搬第一版 C 代码state uint64(1); state mod(state * uint64(16807), uint64(2147483647));2.3 平方取中法演示原理可以生产环境别碰平方取中法由冯·诺依曼在 1940 年代提出思路很直观取当前数的平方再截取平方结果中间的几位作为下一个状态。比如四位种子1234平方是1522756补零成八位01522756取中间四位得到5227再平方、再取中间如此循环。用 C 写一个十进制演示版本// middle_square_demo.cpp #include cstdint #include cstdio uint32_t middle_square(uint32_t x, int digits 4) { uint64_t sq static_castuint64_t(x) * x; uint64_t mask 1; for (int i 0; i digits; i) mask * 10; // 平方补零后取中间 digits 位等价于先除以 10 的幂再取模 return static_castuint32_t((sq / (mask / 10)) % mask); } int main() { uint32_t x 1234; for (int i 0; i 20; i) { x middle_square(x); std::printf(%04u\n, x); } return 0; }这里用了十进制取位方便肉眼观察序列退化。实际工程如果真要实现建议改用二进制位操作把 32 位状态平方后取中间的 16 位或 32 位。无论哪种写法它都有一个致命问题序列很容易进入循环甚至收敛到常数。种子换成0000后会立刻变成全零后面每一次平方取中间都是零。保留这个算法的意义在于教学它能很直观地解释“为什么随机数发生器需要精心设计状态转移”以及“为什么不能拿上一轮结果简单算一算就假装随机”。部署到线上代码前请先把这一节划掉。2.4 xorshift位运算型高性能随机数发生器xorshift 是 George Marsaglia 在 2003 年提出的一族位运算发生器。名字来源于它的核心操作异或xor加移位shift。没有乘法没有取模只靠位运算所以在嵌入式平台和游戏引擎里非常流行。32 位版本的状态转移写成 C 只有三行// Xorshift32.h #include cstdint class Xorshift32 { uint32_t x; public: explicit Xorshift32(uint32_t seed 1) : x(seed ? seed : 1) {} uint32_t next() { x ^ x 13; x ^ x 17; x ^ x 5; return x; } double unit() { return (next() 8) * (1.0 / 16777216.0); } };MATLAB 里没有运算符重载要用bitxor和bitshift显式表达% Xorshift32.m classdef Xorshift32 handle properties x uint32 end methods function obj Xorshift32(seed) if nargin 1 seed uint32(1); end obj.x uint32(seed); end function r next(obj) r obj.x; r bitxor(r, bitshift(r, 13, uint32)); r bitxor(r, bitshift(r, -17, uint32)); r bitxor(r, bitshift(r, 5, uint32)); obj.x r; end function u unit(obj) r obj.next(); u double(bitshift(r, -8)) / 16777216.0; end end end参数说明13, 17, 5这一组移位量是经过验证的三元组组合起来能保证2^32 - 1的满周期。不要随手改成16, 8, 3之类的数位运算发生器对移位参数极其敏感改错一个数周期可能掉到几千万测试时还不容易发现。xorshift 的另一个优势是状态极小且保存方便并行场景下给每个线程分一颗互不相同的种子内存开销几乎可以忽略。2.5 MT19937C 与 MATLAB 各自的标准随机数发生器梅森旋转Mersenne TwisterMT19937是目前主流语言默认随机数发生器的事实标准。名字里的 19937 指的是它的周期是2^19937 - 1并且保证 623 维均匀分布。C11 把它收进randomMATLAB 的rng默认用的也是它。C 侧最小用法#include random #include iostream int main() { std::mt19937 gen(42); // 指定 32 位种子 std::uniform_real_distributiondouble dist(0.0, 1.0); for (int i 0; i 5; i) { std::cout dist(gen) \n; } return 0; }MATLAB 侧对应的最小写法rng(42, twister); % 显式指定梅森旋转算法 u rand(1, 5);这里有一个容易踩的细节两边虽然都是 MT19937生成的内部整数状态序列可以做到一致但 C 的uniform_real_distribution和 MATLAB 的rand在把内部整数映射到[0,1)浮点数时用了不同的策略所以double输出不完全相等。跨语言对拍时不要直接比较浮点数先比较整数状态或者老老实实用自己实现的算法。2.6 五种随机数发生器的周期与适用场景对比算法状态大小周期适用场景注意点LCG4 字节最多2^32嵌入式、教学、跨语言对拍低位周期短取高 24 位Park-Miller8 字节2^31 - 2老系统兼容、低功耗设备种子不能为 0防溢出平方取中法不定极短且退化仅供课堂演示不要用于生产xorshift324 字节2^32 - 1游戏、并行分片、高性能场景移位参数必须经过验证MT199372.5 KB2^19937 - 1仿真、统计计算、默认选择状态大恢复现场成本高3. 种子、周期与参数选择五种随机数发生器用错参数的后果3.1 C 与 MATLAB 的种子机制直接给整数还是给状态向量C 的std::mt19937接受一个 32 位整数做种子内部再用初始化算法把这个整数展开成 624 个 32 位状态MATLAB 的rng(seed, twister)也做类似的事情但它的种子会附带一个内部 shuffle 步骤。结论是就算两边种子整数相同、调用次数相同得到的整数序列也不保证逐位一致。想让 C 与 MATLAB 的随机数发生器输出完全对齐最简单的方法是选择状态转移简单、不依赖复杂初始化的算法比如前面写的 LCG32 和 Xorshift32。给随机数发生器设置种子时还有一个经常被忽略的规则种子不能是 0。Park-Miller 和 xorshift 在实现里都做了seed seed ? seed : 1的保护原因很简单这两种算法的状态为 0 时整个序列会永久停留在 0。LCG 相对宽容一些但如果c和种子凑巧让状态落入某个短循环同样会造成周期缩短。3.2 溢出与取模LCG 参数一旦选错周期直接缩水先看一个反例m 2^32, a 65539, c 0这是老的 RANDU 算法它在三维空间中的随机点会全部落在一组平面上周期和均匀性都极差。另一个反面教材是把a写成偶数且c 0这时状态永远不会出现奇数随机数的信息量直接少了一位。完整的周期条件只有三个c与m互质a - 1能被m的所有质因子整除如果m是 4 的倍数a - 1也必须是 4 的倍数。工程上没时间验证这些条件时直接用 2.1 节给出的1664525和1013904223这两个值已经被反复验证过完整周期成立。溢出问题在 C 侧最隐蔽。看下面这行代码uint64_t bad_next(uint32_t x) { return (x * 16807) % 2147483647; // x 先被提升为 int可能溢出 }x是uint32_t16807是int两者求积时uint32_t会先把int提升为uint32_t但如果x来自一个有符号变量情况就会变成int乘法先溢出再转换。修复方式是给常量加后缀return (x * 16807u) % 2147483647u;类型后缀看起来是小事但在随机数发生器里一个中间溢出的整数乘法足以让整条序列从某个位置开始错位。对拍脚本通常要跑到几十万步才能发现错位那时候排查成本就很高了。3.3 状态恢复随机数发生器的序列如何保存与续跑蒙特卡洛模拟跑了一半宕机需要从断点继续这时不能只保存种子要保存生成器当前状态。MT19937 的状态长达 2.5 KB但标准库提供了序列化能力LCG 和 xorshift 只需要保存一个整数。C 里可以用std::ostream把std::mt19937的状态整体写出#include random #include fstream std::mt19937 gen(42); std::ofstream out(mt_state.bin, std::ios::binary); out gen; // 把 624 个状态字整体写入文件MATLAB 侧对应的做法是用rng输出结构体s rng(42, twister); save(mt_state.mat, s); % 恢复时load(mt_state.mat); rng(s);提示MATLAB 的rng结构体里包含Type,Seed,State三个字段跨 MATLAB 版本恢复时一般兼容但不要手工去改State数组那会让生成器进入未定义状态。续跑模拟的正确姿势是把状态快照和业务数据放在同一份日志里而不是重新跑一遍再赌一次随机序列。4. 怎么看五种随机数发生器是否够随机均匀性检验4.1 先跑均值与方差这步能筛掉一半弱发生器拿到一个随机数发生器我通常先让它生成 100 万个[0,1)样本看均值和方差。均匀分布的均值理论值是0.5方差理论值是1/12 ≈ 0.08333。C 侧统计代码写起来很快#include vector #include numeric double mean(const std::vectordouble v) { return std::accumulate(v.begin(), v.end(), 0.0) / v.size(); } double variance(const std::vectordouble v, double m) { double s 0.0; for (double x : v) s (x - m) * (x - m); return s / v.size(); }用这个函数分别喂给五种发生器的unit()如果均值偏差超过0.001说明低位截断或者状态转移有问题。平方取中法在这个测试里会立刻现形它可能一时均值正常但序列长度一旦进入循环方差就会异常。4.2 卡方检验的 C 与 MATLAB 实现均值方差只能看一阶矩分布形态要用卡方检验。做法是把[0,1)区间等分成 64 个箱统计每个箱里的样本数再比较观测频数和理论频数。C 侧计算卡方统计量#include vector #include cmath double chi2_stat(const std::vectorsize_t counts) { size_t total 0; for (size_t c : counts) total c; double expected static_castdouble(total) / counts.size(); double stat 0.0; for (size_t c : counts) { double d static_castdouble(c) - expected; stat d * d / expected; } return stat; }逻辑说明64 个箱自由度是 63卡方统计量应该在 63 附近波动。如果统计量超过 100基本可以判定分布不均匀。平方取中法和参数错误的 LCG 在这一步过不去。MATLAB 有现成的chi2gof用法更直接rng(42, twister); x rand(1, 1000000); [h, p] chi2gof(x, Edges, linspace(0, 1, 65));参数说明Edges定义分箱边界linspace(0, 1, 65)把区间切成 64 个等宽箱。h 0表示不能拒绝均匀性假设p 0.05是常见参考线。注意[0,1)区间的样本很少落在 1 上边界问题可以忽略但如果用[0,1]闭区间生成器要确认最后一箱是否被异常地填满。4.3 序列相关性与交叉验证分布均匀不等于序列独立。一个简单的相关性测试是计算相邻样本的相关系数理想值接近 0。C 可以用一阶自相关公式更省事的是在 MATLAB 里画 scatter 图rng(42, twister); u1 rand(1, 1000); u2 rand(1, 1000); scatter(u1(1:end-1), u1(2:end), 3, filled);视觉上LCG 会在散点图里出现规则网格xorshift 通常看不出结构MT19937 在 623 维以内都保持均匀。把这五种发生器的散点图并排放一起会非常直观地理解“格点结构”和“高维均匀性”的差别。4.4 跨语言对拍脚本验证两套实现是否真的等价当 C 和 MATLAB 各写了一个版本时最值得做的是逐位对拍。将对拍脚本放在 CI 里每次改动随机数发生器就能立刻发现不一致# 生成 C 侧前 10000 个整数序列 ./lcg_dump 7 cpp_seq.txt # 生成 MATLAB 侧序列 matlab -batch lLCG32(7); for i1:10000, fprintf(%u\n, l.next()); end matlab_seq.txt # 对比 diff cpp_seq.txt matlab_seq.txt没有输出就代表一致。这个做法能验证两个实现的算法逻辑是否真正同步比肉眼对比浮点数可靠得多。5. 让 C 与 MATLAB 跑出同一组随机数序列以及选型技巧5.1 为什么不要直接对拍 MT19937 的浮点输出很多人第一反应是让 C 和 MATLAB 都调用标准 MT19937用同一粒种子对拍。实际对上后会发现前几个 float 就不一样两边虽然核心算法都源自松本与西村的设计但标准库的uniform_real_distribution和 MATLAB 的rand在浮点映射时采取了不同的策略导致输出不可比。跨语言对拍请直接从内部整数状态下手。MT19937 的std::mt19937::operator()和 MATLAB 的RandStream底层状态可以导出但操作起来不如自写算法直观。5.2 同一算法、同一参数、同一状态最省事的跨语言对齐方案手动实现的 LCG 是最容易做到跨语言完全一致的。C 侧直接用无符号 32 位回绕MATLAB 侧用uint32bitshift两边只要种子相同、调用次数相同输出的整数就完全相同。实际中我会把 LCG 的next()结果前 8 个值打印出来作为回归测试的 golden 基线C 和 MATLAB 各跑一份CI 里自动 diff。// lcg_dump.cpp #include cstdio #include LCG32.h int main() { LCG32 lcg(7); for (int i 0; i 8; i) { std::printf(%u\n, lcg.next()); } return 0; }对应 MATLABl LCG32(7); for i 1:8 fprintf(%u\n, l.next()); end这两份代码的输出可以完全相同。原因在于它们做的是同一组无符号整数的乘法与截断不涉及浮点舍入也不涉及不同标准库的实现差异。这个技巧特别适合给跨语言数据管道做随机数对照实验时使用。5.3 生产项目里的选型倾向我一般在项目里这样选默认用 MT19937因为周期足够长、统计性质经过大规模验证并行模拟里改用 xorshift每个线程占用 4 字节状态分片代价极低嵌入式场景用 Park-Miller整型运算在低功耗芯片上表现稳定LCG 只用于跨语言对拍和教学平方取中法永远不碰。最后留一个可执行的小技巧下一次接手的项目如果还在裸调rand()先把全局调用改成显式std::mt19937或 MATLABrng(seed, twister)再跑一遍上面的对拍脚本没有异常再继续调业务逻辑。这两种语言虽然标准库不同但你手上已经有了五种可用的实现和一套验证脚本替换成本其实很低。本文还有配套的精品资源点击获取