Skip to content

线性同余伪随机数

线性同余法是一种基于确定性数学算法生成随机数序列的伪随机数生成方法。

算法概述

程序员通常对随机数比较熟悉,众所周知,计算机生成的随机数往往是伪随机的。那么,究竟什么是伪随机序列呢?

伪随机数是一列看似随机、但由确定性算法生成的数,也就是说,它们并非真正随机的。

既然随机性是通过算法模拟的,那么什么样的算法能接近真正随机的效果?它又是如何工作的呢?

其中一种较为简单的方法是线性同余生成器(LCG)。

LCG 基于如下递推公式:

Nj+1=(ANj+B)%MN_{j+1} = (A * N_{j} + B) \% M

LCG 中最重要的三个整数是:乘数 AA、增量 BB 与模数 MM

这些值 AABBMM 是由生成器设定的常数。LCG 的周期至多为 MM,不过在大多数情况下会小于 MM。为了获得最大周期,必须满足以下条件:

  • BBMM 互素;
  • MM 的所有素因子整除 A1A-1
  • M%4=0M \% 4 = 0,则 A1A-1 也必须能被 4 整除;
  • AABBN0N_{0} 都必须小于 MM,且 AABB 必须是正整数。

LCG 的优点是实现非常简单且生成速度快。然而,其缺点也很明显:对于 32 位整数,其最大周期仅限于 2322^{32},这对于需要高质量随机数的应用(如密码学)是不够的。

随机性分析

有了这些认识,人们可能会好奇:LCG 是如何模拟随机过程的?真正的随机性意味着任何数都可能出现,也就是说所有结果等概率出现,服从均匀分布。那么,LCG 是如何实现均匀分布的呢?

首先,我们需要理解 LCG 的周期。由上述定义可知,对于模数 MM,LCG 生成的序列最大周期为 MM。设递推公式为 Nj+1=(ANj+B)modMN_{j+1} = (A * N_{j} + B) \mod M。若周期为 TT,则 Nk+T=NkN_{k+T} = N_{k},又因为 Nk+1=(ANk+B)modMN_{k+1} = (A * N_{k} + B) \mod M,这意味着 NkN_{k} 的每个值唯一确定 Nk+1N_{k+1} 的值。因此,周期 TT 必须小于或等于 MM,因为对模 MM 取余后只可能有 MM 个不同的整数。第 (M+1)(M+1) 个数必然与前面某个数相同,而由于递推关系是一一对应的,后续序列将发生重复,从而周期 TMT \leq M

其次,乘数与模数互素很重要。这一点很容易理解:考虑 Nj+1=(ANj+B)modMN_{j+1} = (A * N_{j} + B) \mod M。若 AAMM 不互素,设 A=adA = a * dM=mdM = m * d,其中 d>1d > 1。展开递推公式得 Nj+1=ANj+B+kMN_{j+1} = A * N_{j} + B + k * M。由于 BB 是常数,可视为一个平移量,两边减去 BBNj+1B=(a(NjB)+aB+km)dN_{j+1} - B = (a * (N_{j} - B) + a * B + k * m) * d。这表明 Nj+1BN_{j+1} - B 的所有值都含有因子 dd。例如,若 d=2d = 2,则 MM 个可能的数只会产生偶数加上平移量 BB,即周期 TM/2T \leq M / 2。这会破坏序列的均匀分布。

几个简单例子可以说明这一点。设 A=3A = 3M=5M = 5N0=2N_0 = 2,则序列为 {2,1,3,4,2,1,}\{2, 1, 3, 4, 2, 1, \dots\},周期为 5。若设 A=6A = 6M=10M = 10N0=2N_0 = 2,则序列为 {2,2,2,}\{2, 2, 2, \dots\},周期为 1。

回到最初的问题:LCG 是如何实现均匀分布的?由上述讨论可见,只要精心选择参数,就可以构造周期为 MM 的序列,其中每个数在每个周期内恰好出现一次。这就保证了序列中每个数出现的频率相同,从而满足均匀分布的要求。

求解线性同余

在理解了随机分布与线性同余方法之间的关系后,让我们扩展知识,学习如何求解线性同余。这不仅是理论上的数学问题,还具有许多实际应用。例如:

在一条圆形跑道上,两只青蛙 AABB 分别从位置 xxyy 出发,AA 每次跳 mm 个单位,BB 每次跳 nn 个单位。跑道总周长为 LL。若两只青蛙同时出发,当它们落在同一点时称为相遇。问题是:两只青蛙最早在跳多少次后相遇?

根据这些条件,我们可以列出如下方程:

(x+mk)(y+nk)(modL)(x + m*k) ≡ (y + n*k) \pmod{L}

其中 表示对模 LL 同余。

展开得:

x+mk=y+nk+Lkx + m*k = y + n*k + L*k'

移项整理得:

xy=(nm)k+Lkx - y = (n - m)*k + L*k'

因此,

(xy)(nm)k(modL)(x - y) ≡ (n - m)*k \pmod{L}

nm=an - m = axy=bx - y = b
得到标准线性同余式:

akb(modL)a*k ≡ b \pmod{L}

定义:若 aabb 为整数,则形如 axb(modM)a*x ≡ b \pmod{M} 的表达式(其中 xx 为未知整数)称为线性同余。这里 ((modM))(\pmod{M}) 表示两边都对模 MM 取余。

求解线性同余,我们进行如下变换:

  1. 标准线性同余 axb(modM)a*x ≡ b \pmod{M} 等价于丢番图方程 ax+My=ba*x + M*y = b
  2. ddaaMM 的最大公因数,记 gcd(a,M)=dgcd(a,M) = d。若方程有解,则 dd 必须整除 bb。我们将方程改写为 a0x+M0y=b0a_0*x + M_0*y = b_0,其中 a=a0da = a_0*dM=M0dM = M_0*db=b0db = b_0*d。此时 a0a_0M0M_0 互素,即 gcd(a0,M0)=1gcd(a_0, M_0) = 1
  3. x=x0b0x = x_0*b_0y=y0b0y = y_0*b_0。方程变为 a0x0+M0y0=1a_0*x_0 + M_0*y_0 = 1,即等价于 a0x01(modM0)a_0*x_0 ≡ 1 \pmod{M_0}

最后,我们利用模逆元来计算 x0x_0

定理:若 aaMM 是互素的整数且 M>1M > 1,则 aa 关于模 MM 存在唯一的逆元。

该定理可证明如下:

由于 aaMM 互素,存在整数 sstt 使得 sa+tM=1s*a + t*M = 1。(若 aaMM 有公因数 d>1d > 1,则 sa+tM=dks*a + t*M = d*k,不可能等于 1)。对等式两边取模 MM,得 sa1(modM)s*a ≡ 1 \pmod{M}。因此,ssaa 关于模 MM 的模逆元。

aˉ\bar{a} 表示 aa 关于模 MM 的模逆元,则 x0aˉ(modM0)x_0 ≡ \bar{a} \pmod{M_0},即 x0aˉ+kM0x_0 ≡ \bar{a} + k*M_0

例如:

求解 5x7(mod9)5*x ≡ 7 \pmod{9}

  1. 由于 gcd(5,9)=1gcd(5,9) = 1,我们知道 55 关于模 99 存在逆元。用扩展欧几里得算法,可得 25+(1)9=12*5 + (-1)*9 = 1
  2. 因此,2255 关于模 99 的逆元。将原方程两边乘以 22,得 25x27(mod9)2*5*x ≡ 2*7 \pmod{9}
  3. 化简得 x14(mod9)5(mod9)x ≡ 14 \pmod{9} ≡ 5 \pmod{9}
  4. 因此,解为 x=5+9kx = 5 + 9k

总而言之,我们可以将一般的线性同余 axb(modL)a*x ≡ b \pmod{L} 转化为标准同余式 a0x01(modM0)a_0*x_0 ≡ 1 \pmod{M_0} 来求解。解出标准同余式后,再将解映射回去。若标准方程的解为 x0{x=p+qM0}x_0 \in \{x = p + q*M_0\},则原方程的解为 x{x=bdp+qM0}x \in \{x = \frac{b}{d}*p + q'*M_0\}

通过上述内容,我们概述了求解线性同余的过程,并加深了对线性同余方法的理解。

线性同余方法的不安全性

线性同余方法(LCG)是一种基于确定性数学算法的伪随机数生成技术。虽然它能满足一般伪随机数的需求,但不适用于高安全性应用,如密码学或需要强随机性的应用。以下是关于 LCG 不安全性的几点关键考虑:

  1. 确定性:LCG 是确定性的,即相同的种子与参数总是产生相同的随机数序列。这带来了安全风险,因为攻击者可以通过猜测种子与参数重现该序列,从而破坏随机性。
  2. 周期性:LCG 生成的随机数序列是周期性的。一旦序列达到周期,便开始重复。这种周期性会在需要长期随机性的应用中引发问题。
  3. 非均匀性:在某些情况下,LCG 生成的序列可能不服从均匀分布,尤其是当参数选择不当时,会导致统计偏差。
  4. 安全性:LCG 不适合用于密码学或高随机性的安全应用。对于此类应用,必须采用密码学安全的伪随机数生成器(CSPRNG)。

总之,尽管 LCG 可用于模拟、游戏或随机测试中的一般伪随机数生成,但不应在需要高安全性和强随机性的应用中使用。

代码实现

下面是用 TypeScript 实现线性同余生成器(LCG)的代码,以及一个测试其随机性效果的简单程序。

LCG 实现

class LinearCongruentialGenerator {
  private seed: number;
  private readonly a: number;
  private readonly c: number;
  private readonly m: number;

  constructor(
    seed: number,
    a: number = 1664525,
    c: number = 1013904223,
    m: number = 2 ** 32
  ) {
    this.seed = seed;
    this.a = a;
    this.c = c;
    this.m = m;
  }

  public next(): number {
    this.seed = (this.a * this.seed + this.c) % this.m;
    return this.seed;
  }

  public nextFloat(): number {
    return this.next() / this.m;
  }
}

// 使用示例
const lcg = new LinearCongruentialGenerator(12345);
console.log(lcg.next());
// 输出一个随机整数
console.log(lcg.nextFloat());
// 输出一个介于 0 和 1 之间的随机浮点数

随机性测试代码

我们可以通过生成一组随机数并绘制其分布来测试生成器的随机性。下面是一个简单的测试程序:

function testRandomness(
  generator: LinearCongruentialGenerator,
  iterations: number = 10000
): void {
  const results: number[] = new Array(10).fill(0);

  for (let i = 0; i < iterations; i++) {
    const randomValue = Math.floor(generator.nextFloat() * 10);
    results[randomValue]++;
  }

  console.log("随机性测试结果:");
  results.forEach((count, index) => {
    console.log(`${index}: ${"*".repeat(count / 100)}`);
  });
}

// 使用示例
const lcgTest = new LinearCongruentialGenerator(12345);
testRandomness(lcgTest);

该测试程序生成一组随机数并统计每个数出现的次数。通过观察输出的分布,我们可以对生成器的随机性作出初步评估。