面向32位处理器的53→32位模运算优化方案探讨

核心思路：避免除法，利用m=2^n -d的快速模性质 #

MRG32k3a第二个递推的模为m2=2^32 -209（属于m=2^n -d形式，n=32，d=209），针对ARM、RISC-V及GPU这类32位处理器，可通过以下两种无除法的低成本策略实现高效模运算：

策略1：直接对递推中间值做快速模（无需预调整状态） #

递推式为x = c0*s0 - c1*s1（其中c0=1403580，c1=810728，s0/s1为32位状态值且0≤s0,s1<m2），已知|x|<2^53，直接计算x mod m2的步骤：

1. 处理符号，转为等价正数运算

   • 若x≥0，直接以无符号64位格式处理；
   • 若x<0，计算y = -x（无符号64位），最终结果为m2 - (y mod m2)（若y mod m2≠0，否则结果为0）。

2. 对正数执行快速模运算
   以正数z（z<2^53）为例：

   • 拆分z为高32位hi和低32位lo：hi = z >> 32，lo = (uint32_t)z（因z<2^53，hi<2^21，仅21位有效）；
   • 计算临时值temp = hi * d + lo（利用2^32 ≡d mod m2的性质，z ≡ temp mod m2）；
   • 比较temp与m2：
     • 若temp ≥ m2，r = temp - m2；
     • 否则r = temp；
   • r即为z mod m2。

所有操作仅涉及64位加减乘、移位和比较，无任何除法操作，适配支持64位乘法的32位处理器。

策略2：预计算状态乘积的模，避免大中间值 #

利用模运算分配律(a - b) mod m = (a mod m - b mod m) mod m，拆分递推步骤：

1. 分别计算两个乘积的模m2
   对c0*s0和c1*s1分别执行快速模：

   • 以c*s为例（c为常量，s为32位状态值）：
     • 计算64位乘积prod = (uint64_t)c * s；
     • 拆分prod为hi = prod >>32，lo = (uint32_t)prod；
     • 计算temp = hi * d + lo；
     • 若temp ≥m2，mod_result = temp -m2；否则mod_result = temp。

2. 计算最终递推结果

   • 令a = c0*s0 mod m2，b = c1*s1 mod m2；
   • 若a ≥b，new_state = a -b；
   • 否则new_state = a -b + m2。

此策略将大中间值分解为两个小乘积的模运算，步骤更简洁，且完全避免了大数值的符号处理，进一步降低运算开销。

关键优势 #

两种策略均无需32位除法操作，仅依赖32位处理器普遍支持的64位乘法、加减和移位指令，相比现有带除法的偏移方案，性能提升显著，同时完全保证MRG32k3a的伪随机数质量。

内容的提问来源于stack exchange，提问作者njuffa