尧图网站建设 
简介五种随机数发生器的C与MATLAB实现集中于一个压缩包内覆盖平方取中法、乘积取中法、梅森旋转算法、ISAAC和PCG既有入门级方案也有高性能算法。包内共有十二个文件以MATLAB的m脚本和C源文件为主另含可执行程序、工程配置与依赖文件压缩包约320KB轻便易用。目前已有228人学习使用。代码保留了每种算法的核心逻辑适合正在学习随机数原理的同学也适合在模拟、统计或机器学习任务中需要选用随机数工具的研究者。通过对照不同实现可以深入理解线性同余、位运算、置换等设计思路并能根据各算法的周期、统计特性与速度合理迁移到自己的项目中快速掌握随机数生成技术背后的数学构造。1. 解压后先看文件清单一份 C 与 MATLAB 并排的随机发生器源码解压五种随机数发生器-C与MATLAB代码(1).rar根目录里并排放着两套东西一堆.m脚本和一套 Code::Blocks 工程。脚本按算法命名SquareMidRand.m是平方取中MultiplyMidRand.m是乘积取中LinearModRandom.m是线性同余CombineLCG.m是组合线性同余FibRand.m是滞后斐波那契ConstMidRand.m是常数取中C 侧入口是main.cppSquareMid.cbp是工程描述文件obj/Debug与bin/Debug是构建产物目录。这种双语言并排的包不多见它的价值不在代码行数而在于能把同一条递推式在两种语言里写成完全一致的整数序列再用逐项比对验证自己对取位、取模、整数溢出的理解是否到位。做蒙特卡洛仿真、写信道模型、或者被rand()的短周期坑过一次的人适合从这份包开始动手。需要先说明的是包里的文件名和常见资料里说的「五种」并不完全对得上实际能数出六条递推式下面按文件为准来拆。2. 平方取中法与常数取中法取位对齐、种子退化与向量化写法2.1 平方取中法的递推式与十进制对齐平方取中法middle-square method的递推很直白把当前种子平方取结果的中间若干位作为下一个种子。写成公式设种子是 n 位十进制数、n 为偶数则X(n1) floor(X(n)² / 10^(n/2)) mod 10^n取位对齐是这类方法最容易翻车的地方。种子必须严格按 n 位补零比如 n4 时种子 1234 平方得 1522756要补成 8 位01522756再取第 3 到第 6 位得到 5227。少补一位后面所有结果全部错位。手算一次能把规则固定下来1234→5227→3215→3362→3030→1809→2724这串数验证了SquareMidRand.m的前几步。2.2 SquareMidRand.m 的循环实现与预分配MATLAB 里最直接的写法是循环加预分配。伪随机数发生器本质是串行递推第 n1 项依赖第 n 项不存在把循环向量化的空间能优化的只有数组预分配和减少临时变量。function [seq, seed] square_mid_rand(seed, n, len) % SQUARE_MID_RAND 平方取中法 % seed : 初始种子十进制 n 位数 % n : 种子位数必须为偶数 % len : 需要生成的个数 seq zeros(1, len); half n / 2; low 10^half; % 右移位数丢掉低位 mask 10^n; % 截断高位保证结果不超过 n 位 for k 1:len sq seed * seed; % 平方后最多 2n 位 seed mod(floor(sq / low), mask); % 取中间 n 位 seq(k) seed; end end参数说明low决定丢掉多少位低位mask负责把高位截掉。n必须取偶数否则「中间位」没有唯一取法这是实现层面的硬约束。调用示例rng(1); [s, ~] square_mid_rand(1234, 4, 10); disp(s);2.3 种子退化不动点与短周期把seed0100代进去100² 10000补足 8 位是00010000取中间四位回到0100序列原地锁死周期为 1。seed0000同理0 是不动点。这不是实现 bug而是平方取中法本身的缺陷每步都丢弃高位和低位熵持续流失一旦中间段进入全零或固定点就再也出不来。种子n4首次退化步数退化值周期010010100100001000011234未在 100 步内退化—约 10² 量级0101未在 100 步内退化—明显短于 1234上表的结论不是精确周期值而是量级判断4 位种子的状态空间只有 10⁴平方取中法实际能遍历的状态远小于这个数用于教学演示可以用于需要长序列的仿真不行。2.4 ConstMidRand.m 与取中类方法的通病常数取中法把「平方」换成「乘一个固定常数 C」再取中间位X(n1) middle(X(n) × C)。它的状态转移比平方取中多一个可调参数但退化机制完全一样取决于 C 的十进制尾部特征。取中类方法共有的问题是随机性依赖十进制位对齐低位规律性强且没有严格的可证明周期下界。我在实际项目里只在需要「可复现但不必严格均匀」的场合用它比如给排课、抽奖做演示数据。2.5 main.cpp 里的位运算版本C 侧不做十进制取位而是直接把「取中间 16 位」映射到位运算速度快且无十进制转换开销。#include cstdint #include cstdio // 32 位平方取中取 x*x 的中间 16 位 static uint32_t squareMid(uint32_t x) { uint64_t sq (uint64_t)x * x; // 必须升到 64 位否则 32 位乘法先回绕 return (uint32_t)((sq 8) 0xFFFFu); } int main() { uint32_t x 0x12345678u; for (int i 0; i 5; i) { x squareMid(x); printf(%u\n, x); } return 0; }参数说明 8丢掉低 8 位 0xFFFF丢掉高 40 位正好留下中间 16 位。最关键的坑在(uint64_t)x * x这个强制转换上如果写成uint32_t sq x * x乘法会在 32 位宽度内先回绕再赋给 64 位变量结果全错而且错得很安静。这几乎是 C 与 MATLAB 双端对拍时第一个分歧来源。3. 乘积取中、线性同余与组合 LCG满周期条件与 double 精度陷阱3.1 线性同余的递推式与 Hull-Dobell 定理线性同余发生器LCG是全篇最实用的一条X(n1) (a·X(n) c) mod m。它能否走满周期 m有明确的判定条件即 Hull-Dobell 定理c 与 m 互素a−1 能被 m 的所有素因子整除当 m 是 4 的倍数时a−1 也必须能被 4 整除。三个条件同时满足序列周期就等于 m这是取中类方法根本给不了的东西。参数组mac周期老式 C 库风格2³¹1103515245123452³¹Numerical Recipes2³²166452510139042232³²最小标准Park-Miller2³¹−11680702³¹−2组合 LCG 分支撑2147483563400140见 3.4选参数时的顺序建议是先定 m决定周期上限再按 Hull-Dobell 挑 a 和 c最后用谱检验看二维格点是否可接受。直接抄一组数就用往往拿到的是周期不满足或低位相关极强的劣质组合。3.2 乘积取中与纯乘同余的差别MultiplyMidRand.m属于乘积取中/乘同余这一族。纯乘同余就是 LCG 取 c0递推退化成 X(n1) a·X(n) mod m。它天生的毛病是永远取不到 0且低位比特的周期远短于整体周期比如 m2³² 时最低位只有周期 2。乘积取中在乘完之后多一步取中位操作把高位和低位一起牺牲掉换均匀性代价是周期无法用 Hull-Dobell 那套定理保证。包里的MultiplyMidRand.m与LinearModRandom.m放在一起正好方便对比 c0 与非零 c 的分布差异。3.3 MATLAB 的 double 精度陷阱与 Schrage 取模这是整份代码里最容易产生「看起来对、其实全错」的地方。MATLAB 默认数值类型是 double只有 53 位有效精度。当 m2³¹、a≈1.1×10⁹ 时a·X 的量级达到 2⁶¹超过 2⁵³mod运算的结果会被舍入误差污染而且不报任何警告。第一种解法是用uint64承载中间结果function [seq, x] lcg_rand(x0, a, c, m, len) % 用 uint64 做中间运算避开 double 的 53 位精度上限 seq zeros(1, len); x uint64(x0); A uint64(a); C uint64(c); M uint64(m); for k 1:len x mod(A * x C, M); % A*x 最大约 2^62仍在 uint64 范围内 seq(k) double(x) / double(M); end end参数说明A*x的上界是 a·(m−1)用 2³¹ 量级的参数时约 2⁶²uint64 装得下但如果把 m 提到 2⁶⁴这个写法立刻失效需要换成下面的 Schrage 分解。第二种解法是 Schrage 算法全程停在 double 内function y schrage_mod(a, x, m) % 精确计算 a*x mod m要求 0 x m0 a m q floor(m / a); % 商 r mod(m, a); % 余 y a * mod(x, q) - r * floor(x / q); if y 0 y y m; % 结果为负说明需要补一个模 end end参数说明a * mod(x,q)的上界约为 a·q ≈ mr * floor(x/q)也约为 m两项都在 double 能精确表示的小量级内所以结果无舍入误差。x必须严格小于m、a必须小于m否则补模一次的修正量不够。3.4 CombineLCG.m用组合消掉长程相关单条 LCG 再怎么调参二维格点上的平行线族是躲不掉的这是谱检验的经典结论。CombineLCG.m的思路是把两条周期相近但参数不同的 LCG 组合起来常见做法是取差值再模function u combine_lcg(seed1, seed2, len) % 两条独立 LCG 做差模 m1周期约为两周期之积 m1 2147483563; a1 40014; % 分支一 m2 2147483399; a2 40692; % 分支二 x1 uint64(seed1); x2 uint64(seed2); u zeros(1, len); for k 1:len x1 mod(uint64(a1) * x1, uint64(m1)); x2 mod(uint64(a2) * x2, uint64(m2)); z double(mod(int64(x1) - int64(x2), int64(m1))); u(k) z / m1; end end参数说明两个模数 m1、m2 互素且都是素数组合后的周期是 (m1−1)(m2−1)/2 量级远超单条 2³¹。int64转换是为了让差值可能出现的负数走「加一个模」的语义MATLAB 的mod对负数返回非负结果rem不会这里选mod是有意的。验收方法是把 (u(k), u(k1)) 散点画出来单 LCG 会呈现清晰平行线组合之后这类结构基本被打散。4. FibRand.m 与双端对拍Code::Blocks 工程配置到 mt19937 基准4.1 滞后斐波那契发生器的环形缓冲FibRand.m对应的是滞后斐波那契发生器X(n) (X(n−r) X(n−s)) mod m常见滞后对是 (17,5) 和 (24,55)。加法型不用乘法在只有移位和加法的硬件上很友好低位相关性也比乘同余弱但它需要先用另一条发生器把整个缓冲区填满冷启动阶段的质量取决于填充源。function [seq, buf] fib_rand(buf, lag1, lag2, m, len) % buf : 已初始化的环形缓冲长度 max(lag1,lag2) % lag1 : 较大滞后如 55 或 17 % lag2 : 较小滞后如 24 或 5 seq zeros(1, len); n numel(buf); idx 0; for k 1:len v mod(buf(mod(idx - lag1, n) 1) buf(mod(idx - lag2, n) 1), m); buf(mod(idx, n) 1) v; % 覆盖最旧的一个 seq(k) v; idx idx 1; end end参数说明mod(idx - lag1, n) 1是 MATLAB 下标从 1 开始导致的偏移修正漏掉这个1会越界或错位。lag1、lag2必须互素且都小于缓冲长度否则序列会退化成短周期。缓冲区填充建议用 3.3 节的lcg_rand不要用 MATLAB 自带的rand否则后续和 C 对拍时无法复现。4.2 SquareMid.cbp 的构成与命令行编译.cbp是 Code::Blocks 的 XML 工程描述记录编译目标、包含路径、链接库和构建配置SquareMid.depend与SquareMid.layout是 IDE 生成的依赖缓存和界面布局删掉后 IDE 会重建不影响编译结果obj/Debug、bin/Debug分别放中间文件和可执行文件。如果不想开 IDE直接命令行即可# 在工程根目录执行输出与 IDE 的 Debug 目标保持同路径 g -stdc17 -O2 -Wall -mfpmathsse -o bin/Debug/SquareMid main.cpp ./bin/Debug/SquareMid cpp_out.txt参数说明-O2不会改变整数算术结果-mfpmathsse是为了避免 32 位 x86 上默认走 x87 的 80 位扩展精度那会让最后的归一化除法在末位产生差异对拍时能把你折腾一晚上。-Wall打开后如果看到integer overflow in expression之类的提示基本就是 2.5 节那个没升位的问题。4.3 双端对拍的比较脚本对拍的要点是比原始整数而不是比归一化后的浮点。浮点除法在两侧的舍入路径不同会产生大量假分歧。% 读入 C 输出的整数序列与 MATLAB 生成的整数序列逐项比对 cpp load(cpp_out.txt); % C 侧用 printf(%u\n, x) 输出 mat load(mat_out.txt); % MATLAB 侧同样输出整数 n min(numel(cpp), numel(mat)); d find(cpp(1:n) ~ mat(1:n), 1); if isempty(d) fprintf(前 %d 项逐项一致\n, n); else fprintf(首个分歧在第 %d 项cpp%u mat%u\n, d, cpp(d), mat(d)); end参数说明find(..., 1)只取第一个分歧位置定位效率最高。如果分歧出现在第 1 项问题几乎一定在初始种子或取位方式出现在固定偏移处去看缓冲区填充或下标修正随机位置零星分歧优先怀疑精度和类型宽度。4.4 与 std::mt19937、PCG 的基准对照C 侧修完自己实现的那几条之后建议拿random里的标准件做基准这样对「好」的标准心里有数。#include random #include cstdio int main() { std::mt19937 gen(5489u); // 默认种子 std::uniform_real_distributiondouble d(0.0, 1.0); for (int i 0; i 5; i) printf(%.17g\n, d(gen)); return 0; }参数说明mt19937等价于 32 位版本的 Mersenne Twister周期 2¹⁹⁹³⁷−1归一化必须走uniform_real_distribution直接写gen() % N会引入取模偏差gen() / 4294967296.0虽然可用但要注意用 2³² 而不是 2³¹ 做除数。PCG 不在标准库中需要自备实现它的定位是在小状态量下拿到比 LCG 好的谱性质作为对照项加进去即可。发生器状态量周期量级C 标准库支持MATLAB 内置平方取中4 位十进制10² 量级无无经典 LCG32 位2³¹无无组合 LCG两条 31 位约 2⁶⁰无无滞后斐波那契55 个整数很大但难证明无无Mersenne Twister624 个整数2¹⁹⁹³⁷−1std::mt19937rand默认算法R2007a 起PCG64 位2⁶⁴无无MATLAB 侧的rand从 R2007a 起默认改用 Mersenne Twisterrng(seed)可以固定种子但它的序列和 C 的mt19937并不相同双端对拍时不能拿rand当参照系只能拿自己写的同一条递推式互比。5. 用卡方检验和二维格点给五种发生器打分5.1 卡方均匀性检验生成出来的序列不能光看着像随机的就收工。最省事的定量手段是卡方检验把 [0,1) 等分成若干箱统计落箱频数看和期望值的偏离是否超出抽样波动。MATLAB R2014b 之后用histcounts更老的版本用hist。function [h, p] chi2_uniform(seq, nbins) % seq : [0,1) 区间上的归一化序列 % nbins: 分箱数建议 10~50 counts histcounts(seq, 0:1/nbins:1); e numel(seq) / nbins; % 每箱期望频数 chi2 sum((counts - e).^2 / e); p 1 - chi2cdf(chi2, nbins - 1); h p 0.05; % 1 表示拒绝均匀假设 end参数说明样本量至少取分箱数的 5 倍以上否则卡方近似不成立nbins太小检不出局部偏差太大则每箱样本太少导致检验效力反而下降。p值落在 0.05 以下不代表发生器一定坏只说明这一段的分布与均匀假设不符需要换种子重测。5.2 二维格点图LCG 藏不住的平行线把相邻两项组成点对 (u(k), u(k1)) 画散点是判断 LCG 类发生器最直观的一招。u combine_lcg(12345, 67890, 5000); scatter(u(1:end-1), u(2:end), 2, filled); axis([0 1 0 1]); axis square; title(二维格点分布);参数说明点数取 5000 左右就够看出结构太少时线间距显示不出来。单条 LCG 会呈现明显的平行线族线的条数由谱检验的 ν 值决定组合 LCG 和 Mersenne Twister 在同样点数下视觉上接近云团平方取中和常数取中则会看到成片重复点和空洞这正是 2.3 节熵流失的图形化表现。5.3 一条实用的排查顺序遇到序列质量不达标时我一般按这个顺序过一遍先确认种子补位和取位偏移是否与设计一致再确认中间运算有没有掉进 double 的 53 位陷阱然后跑卡方和游程两项检验最后画二维格点。前两步能解决绝大多数「两端结果不一致」的问题后两步决定这条发生器能不能上生产。走完这四步包里那六条递推式各自的边界也就清楚了取中类适合做教学和可复现的演示数据线性同余族适合对速度敏感且周期要求明确的场景滞后斐波那契和组合 LCG 用于需要长序列的仿真真要对随机性做严格要求的场合直接上std::mt19937或自备 PCG 实现别在自制的取中法上硬扛。本文还有配套的精品资源点击获取 五种随机数发生器C++与MATLAB实现详解：从平方取中到组合LCG
简介五种随机数发生器的C与MATLAB实现集中于一个压缩包内覆盖平方取中法、乘积取中法、梅森旋转算法、ISAAC和PCG既有入门级方案也有高性能算法。包内共有十二个文件以MATLAB的m脚本和C源文件为主另含可执行程序、工程配置与依赖文件压缩包约320KB轻便易用。目前已有228人学习使用。代码保留了每种算法的核心逻辑适合正在学习随机数原理的同学也适合在模拟、统计或机器学习任务中需要选用随机数工具的研究者。通过对照不同实现可以深入理解线性同余、位运算、置换等设计思路并能根据各算法的周期、统计特性与速度合理迁移到自己的项目中快速掌握随机数生成技术背后的数学构造。1. 解压后先看文件清单一份 C 与 MATLAB 并排的随机发生器源码解压五种随机数发生器-C与MATLAB代码(1).rar根目录里并排放着两套东西一堆.m脚本和一套 Code::Blocks 工程。脚本按算法命名SquareMidRand.m是平方取中MultiplyMidRand.m是乘积取中LinearModRandom.m是线性同余CombineLCG.m是组合线性同余FibRand.m是滞后斐波那契ConstMidRand.m是常数取中C 侧入口是main.cppSquareMid.cbp是工程描述文件obj/Debug与bin/Debug是构建产物目录。这种双语言并排的包不多见它的价值不在代码行数而在于能把同一条递推式在两种语言里写成完全一致的整数序列再用逐项比对验证自己对取位、取模、整数溢出的理解是否到位。做蒙特卡洛仿真、写信道模型、或者被rand()的短周期坑过一次的人适合从这份包开始动手。需要先说明的是包里的文件名和常见资料里说的「五种」并不完全对得上实际能数出六条递推式下面按文件为准来拆。2. 平方取中法与常数取中法取位对齐、种子退化与向量化写法2.1 平方取中法的递推式与十进制对齐平方取中法middle-square method的递推很直白把当前种子平方取结果的中间若干位作为下一个种子。写成公式设种子是 n 位十进制数、n 为偶数则X(n1) floor(X(n)² / 10^(n/2)) mod 10^n取位对齐是这类方法最容易翻车的地方。种子必须严格按 n 位补零比如 n4 时种子 1234 平方得 1522756要补成 8 位01522756再取第 3 到第 6 位得到 5227。少补一位后面所有结果全部错位。手算一次能把规则固定下来1234→5227→3215→3362→3030→1809→2724这串数验证了SquareMidRand.m的前几步。2.2 SquareMidRand.m 的循环实现与预分配MATLAB 里最直接的写法是循环加预分配。伪随机数发生器本质是串行递推第 n1 项依赖第 n 项不存在把循环向量化的空间能优化的只有数组预分配和减少临时变量。function [seq, seed] square_mid_rand(seed, n, len) % SQUARE_MID_RAND 平方取中法 % seed : 初始种子十进制 n 位数 % n : 种子位数必须为偶数 % len : 需要生成的个数 seq zeros(1, len); half n / 2; low 10^half; % 右移位数丢掉低位 mask 10^n; % 截断高位保证结果不超过 n 位 for k 1:len sq seed * seed; % 平方后最多 2n 位 seed mod(floor(sq / low), mask); % 取中间 n 位 seq(k) seed; end end参数说明low决定丢掉多少位低位mask负责把高位截掉。n必须取偶数否则「中间位」没有唯一取法这是实现层面的硬约束。调用示例rng(1); [s, ~] square_mid_rand(1234, 4, 10); disp(s);2.3 种子退化不动点与短周期把seed0100代进去100² 10000补足 8 位是00010000取中间四位回到0100序列原地锁死周期为 1。seed0000同理0 是不动点。这不是实现 bug而是平方取中法本身的缺陷每步都丢弃高位和低位熵持续流失一旦中间段进入全零或固定点就再也出不来。种子n4首次退化步数退化值周期010010100100001000011234未在 100 步内退化—约 10² 量级0101未在 100 步内退化—明显短于 1234上表的结论不是精确周期值而是量级判断4 位种子的状态空间只有 10⁴平方取中法实际能遍历的状态远小于这个数用于教学演示可以用于需要长序列的仿真不行。2.4 ConstMidRand.m 与取中类方法的通病常数取中法把「平方」换成「乘一个固定常数 C」再取中间位X(n1) middle(X(n) × C)。它的状态转移比平方取中多一个可调参数但退化机制完全一样取决于 C 的十进制尾部特征。取中类方法共有的问题是随机性依赖十进制位对齐低位规律性强且没有严格的可证明周期下界。我在实际项目里只在需要「可复现但不必严格均匀」的场合用它比如给排课、抽奖做演示数据。2.5 main.cpp 里的位运算版本C 侧不做十进制取位而是直接把「取中间 16 位」映射到位运算速度快且无十进制转换开销。#include cstdint #include cstdio // 32 位平方取中取 x*x 的中间 16 位 static uint32_t squareMid(uint32_t x) { uint64_t sq (uint64_t)x * x; // 必须升到 64 位否则 32 位乘法先回绕 return (uint32_t)((sq 8) 0xFFFFu); } int main() { uint32_t x 0x12345678u; for (int i 0; i 5; i) { x squareMid(x); printf(%u\n, x); } return 0; }参数说明 8丢掉低 8 位 0xFFFF丢掉高 40 位正好留下中间 16 位。最关键的坑在(uint64_t)x * x这个强制转换上如果写成uint32_t sq x * x乘法会在 32 位宽度内先回绕再赋给 64 位变量结果全错而且错得很安静。这几乎是 C 与 MATLAB 双端对拍时第一个分歧来源。3. 乘积取中、线性同余与组合 LCG满周期条件与 double 精度陷阱3.1 线性同余的递推式与 Hull-Dobell 定理线性同余发生器LCG是全篇最实用的一条X(n1) (a·X(n) c) mod m。它能否走满周期 m有明确的判定条件即 Hull-Dobell 定理c 与 m 互素a−1 能被 m 的所有素因子整除当 m 是 4 的倍数时a−1 也必须能被 4 整除。三个条件同时满足序列周期就等于 m这是取中类方法根本给不了的东西。参数组mac周期老式 C 库风格2³¹1103515245123452³¹Numerical Recipes2³²166452510139042232³²最小标准Park-Miller2³¹−11680702³¹−2组合 LCG 分支撑2147483563400140见 3.4选参数时的顺序建议是先定 m决定周期上限再按 Hull-Dobell 挑 a 和 c最后用谱检验看二维格点是否可接受。直接抄一组数就用往往拿到的是周期不满足或低位相关极强的劣质组合。3.2 乘积取中与纯乘同余的差别MultiplyMidRand.m属于乘积取中/乘同余这一族。纯乘同余就是 LCG 取 c0递推退化成 X(n1) a·X(n) mod m。它天生的毛病是永远取不到 0且低位比特的周期远短于整体周期比如 m2³² 时最低位只有周期 2。乘积取中在乘完之后多一步取中位操作把高位和低位一起牺牲掉换均匀性代价是周期无法用 Hull-Dobell 那套定理保证。包里的MultiplyMidRand.m与LinearModRandom.m放在一起正好方便对比 c0 与非零 c 的分布差异。3.3 MATLAB 的 double 精度陷阱与 Schrage 取模这是整份代码里最容易产生「看起来对、其实全错」的地方。MATLAB 默认数值类型是 double只有 53 位有效精度。当 m2³¹、a≈1.1×10⁹ 时a·X 的量级达到 2⁶¹超过 2⁵³mod运算的结果会被舍入误差污染而且不报任何警告。第一种解法是用uint64承载中间结果function [seq, x] lcg_rand(x0, a, c, m, len) % 用 uint64 做中间运算避开 double 的 53 位精度上限 seq zeros(1, len); x uint64(x0); A uint64(a); C uint64(c); M uint64(m); for k 1:len x mod(A * x C, M); % A*x 最大约 2^62仍在 uint64 范围内 seq(k) double(x) / double(M); end end参数说明A*x的上界是 a·(m−1)用 2³¹ 量级的参数时约 2⁶²uint64 装得下但如果把 m 提到 2⁶⁴这个写法立刻失效需要换成下面的 Schrage 分解。第二种解法是 Schrage 算法全程停在 double 内function y schrage_mod(a, x, m) % 精确计算 a*x mod m要求 0 x m0 a m q floor(m / a); % 商 r mod(m, a); % 余 y a * mod(x, q) - r * floor(x / q); if y 0 y y m; % 结果为负说明需要补一个模 end end参数说明a * mod(x,q)的上界约为 a·q ≈ mr * floor(x/q)也约为 m两项都在 double 能精确表示的小量级内所以结果无舍入误差。x必须严格小于m、a必须小于m否则补模一次的修正量不够。3.4 CombineLCG.m用组合消掉长程相关单条 LCG 再怎么调参二维格点上的平行线族是躲不掉的这是谱检验的经典结论。CombineLCG.m的思路是把两条周期相近但参数不同的 LCG 组合起来常见做法是取差值再模function u combine_lcg(seed1, seed2, len) % 两条独立 LCG 做差模 m1周期约为两周期之积 m1 2147483563; a1 40014; % 分支一 m2 2147483399; a2 40692; % 分支二 x1 uint64(seed1); x2 uint64(seed2); u zeros(1, len); for k 1:len x1 mod(uint64(a1) * x1, uint64(m1)); x2 mod(uint64(a2) * x2, uint64(m2)); z double(mod(int64(x1) - int64(x2), int64(m1))); u(k) z / m1; end end参数说明两个模数 m1、m2 互素且都是素数组合后的周期是 (m1−1)(m2−1)/2 量级远超单条 2³¹。int64转换是为了让差值可能出现的负数走「加一个模」的语义MATLAB 的mod对负数返回非负结果rem不会这里选mod是有意的。验收方法是把 (u(k), u(k1)) 散点画出来单 LCG 会呈现清晰平行线组合之后这类结构基本被打散。4. FibRand.m 与双端对拍Code::Blocks 工程配置到 mt19937 基准4.1 滞后斐波那契发生器的环形缓冲FibRand.m对应的是滞后斐波那契发生器X(n) (X(n−r) X(n−s)) mod m常见滞后对是 (17,5) 和 (24,55)。加法型不用乘法在只有移位和加法的硬件上很友好低位相关性也比乘同余弱但它需要先用另一条发生器把整个缓冲区填满冷启动阶段的质量取决于填充源。function [seq, buf] fib_rand(buf, lag1, lag2, m, len) % buf : 已初始化的环形缓冲长度 max(lag1,lag2) % lag1 : 较大滞后如 55 或 17 % lag2 : 较小滞后如 24 或 5 seq zeros(1, len); n numel(buf); idx 0; for k 1:len v mod(buf(mod(idx - lag1, n) 1) buf(mod(idx - lag2, n) 1), m); buf(mod(idx, n) 1) v; % 覆盖最旧的一个 seq(k) v; idx idx 1; end end参数说明mod(idx - lag1, n) 1是 MATLAB 下标从 1 开始导致的偏移修正漏掉这个1会越界或错位。lag1、lag2必须互素且都小于缓冲长度否则序列会退化成短周期。缓冲区填充建议用 3.3 节的lcg_rand不要用 MATLAB 自带的rand否则后续和 C 对拍时无法复现。4.2 SquareMid.cbp 的构成与命令行编译.cbp是 Code::Blocks 的 XML 工程描述记录编译目标、包含路径、链接库和构建配置SquareMid.depend与SquareMid.layout是 IDE 生成的依赖缓存和界面布局删掉后 IDE 会重建不影响编译结果obj/Debug、bin/Debug分别放中间文件和可执行文件。如果不想开 IDE直接命令行即可# 在工程根目录执行输出与 IDE 的 Debug 目标保持同路径 g -stdc17 -O2 -Wall -mfpmathsse -o bin/Debug/SquareMid main.cpp ./bin/Debug/SquareMid cpp_out.txt参数说明-O2不会改变整数算术结果-mfpmathsse是为了避免 32 位 x86 上默认走 x87 的 80 位扩展精度那会让最后的归一化除法在末位产生差异对拍时能把你折腾一晚上。-Wall打开后如果看到integer overflow in expression之类的提示基本就是 2.5 节那个没升位的问题。4.3 双端对拍的比较脚本对拍的要点是比原始整数而不是比归一化后的浮点。浮点除法在两侧的舍入路径不同会产生大量假分歧。% 读入 C 输出的整数序列与 MATLAB 生成的整数序列逐项比对 cpp load(cpp_out.txt); % C 侧用 printf(%u\n, x) 输出 mat load(mat_out.txt); % MATLAB 侧同样输出整数 n min(numel(cpp), numel(mat)); d find(cpp(1:n) ~ mat(1:n), 1); if isempty(d) fprintf(前 %d 项逐项一致\n, n); else fprintf(首个分歧在第 %d 项cpp%u mat%u\n, d, cpp(d), mat(d)); end参数说明find(..., 1)只取第一个分歧位置定位效率最高。如果分歧出现在第 1 项问题几乎一定在初始种子或取位方式出现在固定偏移处去看缓冲区填充或下标修正随机位置零星分歧优先怀疑精度和类型宽度。4.4 与 std::mt19937、PCG 的基准对照C 侧修完自己实现的那几条之后建议拿random里的标准件做基准这样对「好」的标准心里有数。#include random #include cstdio int main() { std::mt19937 gen(5489u); // 默认种子 std::uniform_real_distributiondouble d(0.0, 1.0); for (int i 0; i 5; i) printf(%.17g\n, d(gen)); return 0; }参数说明mt19937等价于 32 位版本的 Mersenne Twister周期 2¹⁹⁹³⁷−1归一化必须走uniform_real_distribution直接写gen() % N会引入取模偏差gen() / 4294967296.0虽然可用但要注意用 2³² 而不是 2³¹ 做除数。PCG 不在标准库中需要自备实现它的定位是在小状态量下拿到比 LCG 好的谱性质作为对照项加进去即可。发生器状态量周期量级C 标准库支持MATLAB 内置平方取中4 位十进制10² 量级无无经典 LCG32 位2³¹无无组合 LCG两条 31 位约 2⁶⁰无无滞后斐波那契55 个整数很大但难证明无无Mersenne Twister624 个整数2¹⁹⁹³⁷−1std::mt19937rand默认算法R2007a 起PCG64 位2⁶⁴无无MATLAB 侧的rand从 R2007a 起默认改用 Mersenne Twisterrng(seed)可以固定种子但它的序列和 C 的mt19937并不相同双端对拍时不能拿rand当参照系只能拿自己写的同一条递推式互比。5. 用卡方检验和二维格点给五种发生器打分5.1 卡方均匀性检验生成出来的序列不能光看着像随机的就收工。最省事的定量手段是卡方检验把 [0,1) 等分成若干箱统计落箱频数看和期望值的偏离是否超出抽样波动。MATLAB R2014b 之后用histcounts更老的版本用hist。function [h, p] chi2_uniform(seq, nbins) % seq : [0,1) 区间上的归一化序列 % nbins: 分箱数建议 10~50 counts histcounts(seq, 0:1/nbins:1); e numel(seq) / nbins; % 每箱期望频数 chi2 sum((counts - e).^2 / e); p 1 - chi2cdf(chi2, nbins - 1); h p 0.05; % 1 表示拒绝均匀假设 end参数说明样本量至少取分箱数的 5 倍以上否则卡方近似不成立nbins太小检不出局部偏差太大则每箱样本太少导致检验效力反而下降。p值落在 0.05 以下不代表发生器一定坏只说明这一段的分布与均匀假设不符需要换种子重测。5.2 二维格点图LCG 藏不住的平行线把相邻两项组成点对 (u(k), u(k1)) 画散点是判断 LCG 类发生器最直观的一招。u combine_lcg(12345, 67890, 5000); scatter(u(1:end-1), u(2:end), 2, filled); axis([0 1 0 1]); axis square; title(二维格点分布);参数说明点数取 5000 左右就够看出结构太少时线间距显示不出来。单条 LCG 会呈现明显的平行线族线的条数由谱检验的 ν 值决定组合 LCG 和 Mersenne Twister 在同样点数下视觉上接近云团平方取中和常数取中则会看到成片重复点和空洞这正是 2.3 节熵流失的图形化表现。5.3 一条实用的排查顺序遇到序列质量不达标时我一般按这个顺序过一遍先确认种子补位和取位偏移是否与设计一致再确认中间运算有没有掉进 double 的 53 位陷阱然后跑卡方和游程两项检验最后画二维格点。前两步能解决绝大多数「两端结果不一致」的问题后两步决定这条发生器能不能上生产。走完这四步包里那六条递推式各自的边界也就清楚了取中类适合做教学和可复现的演示数据线性同余族适合对速度敏感且周期要求明确的场景滞后斐波那契和组合 LCG 用于需要长序列的仿真真要对随机性做严格要求的场合直接上std::mt19937或自备 PCG 实现别在自制的取中法上硬扛。本文还有配套的精品资源点击获取