32、随机序列生成算法与并行化研究

随机序列生成算法与并行化研究

1. 线性递推快速跳跃算法

1.1 伪随机数生成器输出问题

在伪随机数生成器中，使用 $(x_i, \cdots, x_{i + k - 1})$ 作为第 $i$ 个输出的 $k$ 位整数并不合适，因为连续输出存在重叠。不过，这种线性反馈移位寄存器（LFSR）在流密码（伪随机比特生成器）中会配合合适的非线性输出函数 $o : S \to F_2$ 使用。

1.2 梅森旋转算法（Mersenne Twister）

1.2.1 算法描述

设 $w$ 为机器的字长（如 $w = 32$），行向量空间 $W := F_2^w$ 可看作字的集合。对于固定整数 $n > m$，梅森旋转算法（MT）通过以下递推式生成 $W$ 中的元素序列 $x_0, x_1, \cdots$：
$x_{j + n} := x_{j + m} \oplus (x_{j}^{w - r} | x_{j + 1}^{r})A$，$j = 0, 1, \cdots$
其中，$(x_{j}^{w - r} | x_{j + 1}^{r})$ 表示 $x_j$ 的高 $(w - r)$ 位 $(x_{j}^{w - r})$ 与 $x_{j + 1}$ 的低 $r$ 位 $(x_{j + 1}^{r})$ 的拼接，$w \times w$ 矩阵 $A$ 定义如下：对于任意 $w$ 维行向量 $x$，
$xA =
\begin{cases}
\text{shiftright}(x) & \text{如果 } x \text{ 的最低有效位（LSB）为 } 0 \
\text{shiftright}(x) \oplus a & \text{如果 } x \text{ 的最低有效位（LSB）为 } 1
\end{cases}$
这里，$a$ 是一个合适选择的常数。该生成器的状态转移函数为 $f(x_0^{w - r}, x_1, \cdots, x_{n - 1}) = (x_1^{w - r}, x_2, \cdots, x_n)$，状态空间为 $S = F_2^{nw - r}$。

1.2.2 流行实例 MT19937

最流行的实例 MT19937 使用参数 $n = 624$，$w = 32$，$r = 31$，其序列的最大周期等于梅森素数 $2^{19937} - 1$。由于其高速和良好的分布特性，MT19937 被广泛用作标准伪随机数生成器（PRNG），但此前缺乏高效的跳跃方法。

1.2.3 跳跃算法

对于 MT 生成器，如果选择输出函数 $o : (x_0^{w - r}, x_1, \cdots, x_{n - 1}) \to x_1$ 的最低有效位，那么 $o(k)$ 的逆像可以在时间复杂度 $O(k)$ 内计算。因此，如果使用 Karatsuba 多项式乘法，跳跃可以在 $O(k^{1.59})$ 时间内完成；如果使用快速傅里叶变换，则可以在 $O(k \log k)$ 时间内完成。

1.3 跳跃方法速度比较

进行了一个实验，比较三种跳跃方法的速度：
1. 直接实现霍纳方法（Horner）。
2. 带参数 $q$ 的滑动窗口方法（SW）。
3. 基于 Karatsuba 算法的多项式乘法新方法（PM）。

将这些方法应用于具有梅森指数 $k = 19937, 21701, 23209, 44497, 86243, 110503, 132049$ 的 MT 生成器。对于每个 $k$ 值，在两台不同的计算机上重复以下操作 1000 次：生成一个在 $F_2[t]$ 中次数小于 $k$ 的随机多项式 $g(t)$ 和一个在 $S = F_2^k$ 中的随机状态 $s$，然后使用三种算法计算 $f^J(s)$ 并测量所需的 CPU 时间。结果如下表所示：

计算机	$k$	Horner CPU (sec)	SW $q$	SW 内存 (KB)	SW CPU (sec)	PM CPU (sec)
Intel Core Duo 32 - bit	19937	17.7	7	312	5.1	6.7
	21701	20.8	8	679	5.9	7.4
	23209	23.7	8	726	6.5	7.9
	44497	84.0	8	1391	20.4	19.8
	86243	309.3	8	2696	71.7	52.3
	110503	445.4	9	6908	113.6	71.1
	132049	648.1	9	8254	158.3	89.8
AMD Athlon 64 3800 + 64 - bit	19937	8.6	7	312	2.7	2.2
	21701	10.1	7	340	3.1	2.4
	23209	11.4	7	363	3.5	2.6
	44497	40.4	6	348	11.7	6.0
	86243	148.3	6	674	44.8	14.5
	110503	242.2	5	432	73.2	19.7
	132049	345.0	5	516	106.2	24.7

从结果可以看出，对于 32 位计算机，当 $k \geq 44497$ 时，PM 方法比 SW 方法快；对于 64 位计算机，当 $k \geq 19937$ 时，PM 方法更快。这与计算复杂度近似结果一致，SW 方法的复杂度为 $O(k^2)$，PM 方法的复杂度为 $O(k^{1.59})$。

1.4 结论

基于多项式乘法的跳跃算法在状态空间维度 $k$ 足够大时优于滑动窗口方法。新的 PM 方法时间复杂度为 $O(k^{1.59})$，而滑动窗口方法为 $O(k^2)$。实验结果证实了这一点，并且表明这个足够大的 $k$ 大致对应于 MT 最流行实现中使用的 $k$ 值。更重要的是，新的 PM 方法空间复杂度为 $O(k)$，远小于滑动窗口方法的 $O(k^2q)$。不过，新方法的主要限制是需要一个高效的算法来计算 $o(k)$ 的逆。

2. 伪随机序列的并行生成

2.1 研究背景

在现代通信方案中，高速生成伪随机序列是一个重要问题。通过抽取（decimation）可以对序列进行缩放以实现并行性，从而得到子序列生成器。子序列生成器和抽取在过去已广泛应用于线性反馈移位寄存器（LFSR），但对于具有非线性反馈的自动机，相关研究仍在进行中。

2.2 研究目标

研究如何将带进位反馈移位寄存器（FCSR）转换为子序列生成器，考察两种转换解决方案：一种基于 $\ell$ 序列（即具有最大周期的 FCSR 序列）的抽取特性，另一种基于多步实现。

2.3 抽取的作用

抽取是将 1 位输出生成器转换为子序列生成器的主要工具，它可以提高伪随机序列生成器（PRSG）的吞吐量或降低功耗。

2.3.1 提高吞吐量

PRSG 的吞吐量 $T$ 定义为 $T = n \times f$，通常 $n = 1$。抽取可以实现吞吐量的优化：$T_d = d \times \gamma f$，其中 $0 < \gamma \leq 1$ 是原始自动机频率的降解因子。当 $\gamma d > 1$ 时，抽取可以提高吞吐量。在理想情况下，$\gamma = 1$，$d$ - 抽取将使吞吐量乘以 $d$。

2.3.2 降低功耗

CMOS 设备的功耗可以通过 $P = C \times V_{dd}^2 \times f$ 估算。序列抽取可以通过交织子序列来降低设备频率。子序列生成器将以频率 $\frac{\gamma f}{d}$ 时钟，输出将与一个以频率 $\gamma f$ 时钟的 $d$ 输入多路复用器组合。原始功耗可以降低 $\frac{\gamma}{d}$ 倍，其中 $\gamma$ 必须接近 1 以保证最终序列以频率 $f$ 生成。

2.4 LFSR 和 FCSR 的基本概念

2.4.1 LFSR

LFSR 是一种生成线性递归序列的自动机。设 $s(x) = \sum_{i = 0}^{\infty} s_i x^i$ 是 LFSR 生成序列 $S = (s_0, s_1, s_2, \cdots)$ 的幂级数，则存在两个多项式 $p(x)$ 和 $q(x)$，使得 $s(x) = \frac{p(x)}{q(x)}$，其中 $q(x)$ 是由自动机反馈位置定义的连接多项式，$Q(x) = x^m q(\frac{1}{x})$ 是特征多项式。如果 $q(x)$ 是本原多项式，则 LFSR 的输出序列是具有最大周期 $2^m - 1$ 的 $m$ - 序列。LFSR 有两种表示形式：Galois 和 Fibonacci 结构，如下图所示：

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A(x7):::process --> B(x6):::process --> C(x5):::process --> D(x4):::process --> E(x3):::process --> F(x2):::process --> G(x1):::process --> H(x0):::process
    I(x7):::process --> J(x6):::process --> K(x5):::process
    K --> L(x3):::process
    K --> M(x4):::process
    L --> N(x2):::process --> O(x1):::process --> P(x0):::process
    A1(Galois setup):::startend
    A2(Fibonacci setup):::startend
    style A1 fill:#ffffff,stroke:#000000,stroke-width:2px
    style A2 fill:#ffffff,stroke:#000000,stroke-width:2px

两种结构都使用模 2 加法。Fibonacci 结构的特点是多输入单输出反馈函数，而 Galois 结构有多个具有公共输入的反馈函数。在两种结构中，$x_0$ 表示 LFSR 的输出单元。同一序列可以由具有相同特征多项式但不同初始化的 Fibonacci 或 Galois 结构的 LFSR 生成。序列 $S$ 的线性复杂度定义为能够生成该序列的最小 LFSR 的大小。

2.4.2 FCSR

FCSR 由 Klapper 和 Goresky 引入，与 LFSR 不同，FCSR 使用带进位加法，需要额外的内存来存储进位。其非线性更新函数使其在对线性敏感的领域（如流密码）中特别有用。二进制 FCSR 的输出对应于有理数 $\frac{h}{q} \leq 0$ 的 2 - 进展开。与 LFSR 类似，FCSR 也有 Fibonacci 和 Galois 结构，如下图所示：

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A(x2):::process --> B(x3):::process --> C(x4):::process --> D(x5):::process --> E(x6):::process --> F(x7):::process
    G(x4):::process --> H(x6):::process --> I(x3):::process --> J(x0):::process --> K(x1):::process --> L(x2):::process --> M(x5):::process --> N(x7):::process --> O(x1):::process --> P(x0):::process
    Q(c /2):::process
    R(mod2):::process
    S():::process
    A1(Galois setup):::startend
    A2(Fibonacci setup):::startend
    style A1 fill:#ffffff,stroke:#000000,stroke-width:2px
    style A2 fill:#ffffff,stroke:#000000,stroke-width:2px

在本文中，考虑 $\ell$ 序列，即具有最大周期 $\phi(q)$ 的序列（其中 $\phi$ 是欧拉函数），这等价于 $q$ 是素数幂且 2 模 $q$ 的乘法阶为 $\phi(q)$。为了简化，还限制生成的序列是严格周期的，这对于 Galois FCSR 总是成立，但对于 Fibonacci FCSR 不一定。序列的 2 - 进复杂度定义为能够生成该序列的最小 FCSR 的大小，在周期情况下，这等价于 $q$ 的比特长度。

2.5 LFSR 序列的子序列生成

2.5.1 基于 LFSR 合成的方法

该方法为每个子序列关联一个 LFSR，基于 LFSR 序列抽取的已知结果，可应用于 Fibonacci 和 Galois 表示。设 $S$ 是由特征多项式 $Q(x)$ 在 $F_2$ 上不可约且次数为 $m$ 的 LFSR 生成的序列，$\alpha$ 是 $Q(x)$ 的根，$T$ 是 $Q(x)$ 的周期。$S$ 的 $d$ - 抽取得到的子序列 $S_i^d$ 可以由具有以下性质的 LFSR 生成：
- $\alpha^d$ 在 $F_{2^m}$ 中的最小多项式是所得 LFSR 的特征多项式 $Q^{\star}(x)$。
- $Q^{\star}(x)$ 的周期 $T^{\star}$ 等于 $\frac{T}{\gcd(d, T)}$。
- $Q^{\star}(x)$ 的次数 $m^{\star}$ 等于 $Q(x)$ 在 $\mathbb{Z}_{T^{\star}}$ 中的乘法阶。

在实践中，可以使用 Berlekamp - Massey 算法确定特征多项式 $Q^{\star}(x)$。子序列由 $d$ 个由特征多项式 $Q^{\star}(x)$ 定义但初始值不同的 LFSR 生成。在 LFSR 的情况下，$m^{\star}$ 总是小于或等于 $m$。

2.5.2 基于多步 LFSR 的方法

该方法由 Lempel 和 Eastman 首次提出，通过改变存储单元之间的连接并复制反馈函数，在一个时钟周期内对 LFSR 进行 $d$ 次时钟操作，得到一个线性互连的移位寄存器网络。该方法对于 Galois 和 Fibonacci 结构有所不同。

对于 $m$ 位 Fibonacci LFSR 转换为每周期生成 $d$ 位的自动机，使用以下方程：
$next_d(x_i) = x_{i - d \bmod m}$
$(x_i) {t + d} =
\begin{cases}
f(X {t + i - m + d}) & \text{如果 } m - d \leq i < m \
(x_{i + d})_t & \text{如果 } i < m - d
\end{cases}$
其中，$next_d(x_i)$ 是连接到 $x_i$ 输出的单元，$f$ 是反馈函数。第一个方程对应于存储单元连接的转换，将原始 LFSR 中满足 $i \bmod d = k$ 的所有单元 $x_i$ 聚集形成一个子移位寄存器，第 $k$ 个子移位寄存器的最后一个单元的内容对应于第 $k$ 个子序列 $S_k^d$。第二个方程对应于反馈函数的转换，合成需要仅建立寄存器在时间 $t + d$ 和 $t$ 的状态之间的关系。

对于 Galois 结构，转换由以下方程描述：
$(x_i) {t + d} =
\begin{cases}
(x_0) {t + d - m + i} \oplus \sum_{k = 0}^{m - 2 - i} a_{i + k} (x_0) {t + d - k - 1} & \text{如果 } m - d \leq i < m \
(x {i + d}) t \oplus \sum {k = 0}^{d - 1} a_{i + d - 1 - k} (x_0) {t + k} & \text{如果 } i < m - d
\end{cases}$
其中，$q(x) = 1 + a_0 x + a_1 x^2 + \cdots + a {m - 2} x^{m - 1} + x^m$。虽然这个方程没有直接提供寄存器在时间 $t + d$ 和 $t$ 的状态之间的关系，但可以很容易地推导得到更实用的公式。

2.5.3 两种方法的比较

方法	存储单元	逻辑门
LFSR 合成	$d \times m^{\star}$	$d \times wt(Q^{\star})$
多步 LFSR	$m$	$d \times wt(Q)$

基于 LFSR 合成的方法证明了子序列生成器合成存在解决方案，该方案的内存成本和门数量都取决于抽取因子 $d$。而 Lempel 和 Eastman 提出的多步 LFSR 方法在子序列生成器合成中使用固定数量的存储单元。此外，使用 Berlekamp - Massey 算法定义的子序列生成器不适合降低 LFSR 的功耗，而多步 LFSR 可以实现更好的功耗降低效果。

2.6 FCSR 序列的子序列生成

2.6.1 基于 $\ell$ 序列抽取特性的方法

探索 $\ell$ 序列的抽取特性来设计子序列生成器，使用 FCSR 合成算法。然而，与 LFSR 的情况相比，基于抽取特性的解决方案导致的结果成本要高得多。

2.6.2 基于多步实现的方法

展示如何在 Fibonacci 和 Galois 配置中实现多步 FCSR。在实现过程中，需要考虑进位的传播对设计的影响。

2.7 总结与展望

在伪随机序列的并行生成研究中，对于 LFSR 序列的子序列生成已经有了较为成熟的方法，包括基于 LFSR 合成和多步 LFSR 的方法，并且对这两种方法在存储单元和逻辑门使用上进行了比较。而对于 FCSR 转换为子序列生成器的研究，虽然提出了基于 $\ell$ 序列抽取特性和多步实现的两种解决方案，但基于抽取特性的方案成本较高，多步实现方案需要考虑进位传播的影响。

未来的研究可以进一步优化 FCSR 子序列生成的方法，降低成本，提高性能。例如，深入研究进位传播的规律，找到更有效的处理方法，以减少其对设计的不利影响。同时，可以探索更多的序列转换和并行生成策略，以满足不同应用场景对伪随机序列高速生成和低功耗的需求。

3. 综合比较与应用建议

3.1 线性递推跳跃算法比较

算法	时间复杂度	空间复杂度	适用情况
滑动窗口方法（SW）	$O(k^2)$	$O(k^2q)$	状态空间维度 $k$ 较小时
基于多项式乘法的方法（PM）	$O(k^{1.59})$	$O(k)$	状态空间维度 $k$ 足够大时

从复杂度和实验结果来看，当 $k$ 较小时，滑动窗口方法可能更合适，因为其实现相对简单。但随着 $k$ 的增大，基于多项式乘法的方法在时间和空间复杂度上都具有明显优势，尤其是在需要快速跳跃且对空间有一定限制的场景中。

3.2 LFSR 和 FCSR 子序列生成方法比较

生成器类型	子序列生成方法	存储单元	逻辑门	功耗降低效果
LFSR	基于 LFSR 合成	$d \times m^{\star}$	$d \times wt(Q^{\star})$	较差
LFSR	基于多步 LFSR	$m$	$d \times wt(Q)$	较好
FCSR	基于 $\ell$ 序列抽取特性	成本高	-	-
FCSR	基于多步实现	需考虑进位传播	-	-

对于 LFSR，基于多步 LFSR 的方法在存储单元使用和功耗降低方面表现更优。而对于 FCSR，基于多步实现的方法有一定的潜力，但需要解决进位传播带来的设计问题。

3.3 应用建议

• 线性递推跳跃场景 ：如果需要在状态空间维度较小的情况下进行跳跃操作，可优先考虑滑动窗口方法。而在状态空间维度较大，且对跳跃速度和空间有要求时，基于多项式乘法的方法是更好的选择。
• 伪随机序列并行生成场景 ：
  • 对于 LFSR，如果主要关注功耗降低和使用固定数量的存储单元，建议使用基于多步 LFSR 的方法生成子序列。
  • 对于 FCSR，目前基于多步实现的方法具有研究价值，可进一步探索如何优化进位传播的处理，以实现高效的子序列生成。

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A(线性递推跳跃):::process --> B{状态空间维度 k 大小}:::decision
    B -->|k 较小| C(滑动窗口方法):::process
    B -->|k 足够大| D(基于多项式乘法的方法):::process
    E(伪随机序列并行生成):::process --> F{生成器类型}:::decision
    F -->|LFSR| G{关注重点}:::decision
    G -->|功耗降低和固定存储单元| H(基于多步 LFSR 的方法):::process
    G -->|其他| I(基于 LFSR 合成的方法):::process
    F -->|FCSR| J(基于多步实现的方法):::process

综上所述，在不同的应用场景中，需要根据具体的需求和条件选择合适的算法和方法，以实现高效的伪随机序列生成和操作。未来的研究可以围绕 FCSR 子序列生成方法的优化以及新的序列生成和处理策略展开，以推动伪随机序列技术在更多领域的应用。