随机数生成与相关算法解析

### 随机数生成与相关算法解析 #### 1. 随机数生成基础代码逻辑 在随机数生成的过程中，会根据不同参数的变化来计算一些有用的量。以下是一段代码示例： ```c if (n != nold) { // If n has changed, then compute useful quantities. en=n; oldg=gammln(en+1.0); nold=n; } if (p != pold) { // If p has changed, then compute useful quantities. pc=1.0-p; plog=log(p); pclog=log(pc); pold=p; } sq=sqrt(2.0*am*pc); // 以下是拒绝法，使用洛伦兹比较函数 do { do { angle=PI*ran1(idum); y=tan(angle); em=sq*y+am; } while (em < 0.0 || em >= (en+1.0)); // Reject. em=floor(em); // Trick for integer-valued distribution. t=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0) -gammln(en-em+1.0)+em*plog+(en-em)*pclog); } while (ran1(idum) > t); // Reject. This happens about 1.5 times per deviate, on average. bnl=em; if (p != pp) bnl=n-bnl; // Remember to undo the symmetry transformation. return bnl; ``` 这段代码的逻辑是，当`n`或`p`发生变化时，会重新计算一些相关的量。然后使用拒绝法生成随机数，通过一系列的条件判断和计算，最终得到符合要求的随机数。 #### 2. 随机比特生成 在某些特殊应用场景，如实时信号处理中，需要快速生成随机比特。这里介绍两种基于“模 2 本原多项式”的方法。 ##### 2.1 模 2 本原多项式 模 2 本原多项式是系数为 0 或 1 的特殊多项式。例如： \(x^{18} + x^{5} + x^{2} + x^{1} + x^{0}\) 可以简写成\((18, 5, 2, 1, 0)\)。每个\(n\)阶的模 2 本原多项式都定义了一个从\(n\)个前导比特得到一个新随机比特的递推关系，能产生最大长度的序列。 ##### 2.2 方法一（Method I） 方法一在硬件实现上较为简单，只需要一个\(n\)位长的移位寄存器和几个异或（XOR）门。对于上述的本原多项式，递推公式为： \(a_0 = a_{18} \land a_{5} \land a_{2} \land a_{1}\) 以下是对应的 C 语言代码： ```c int irbit1(unsigned long *iseed) { unsigned long newbit; newbit = (*iseed >> 17) & 1 ^ (*iseed >> 4) & 1 ^ (*iseed >> 1) & 1 ^ (*iseed & 1); *iseed=(*iseed << 1) | newbit; return (int) newbit; } ``` 该方法通过一系列的位操作来获取随机比特，但在 C 语言中实现时相对繁琐，需要进行全字掩码操作。 ##### 2.3 方法二（Method II） 方法二不太适合直接的硬件实现，但非常适合在 C 语言中使用。对于上述本原多项式，其规则为： \(a_0 = a_{18}\) \(a_5 = a_5 \land a_0\) \(a_2 = a_2 \land a_0\) \(a_1 = a_1 \land a_0\) 以下是对应的 C 语言代码： ```c #define IB1 1 #define IB2 2 #define IB5 16 #define IB18 131072 #define MASK (IB1+IB2+IB5) int irbit2(unsigned long *iseed) { if (*iseed & IB18) { *iseed=((*iseed ^ MASK) << 1) | IB1; return 1; } else { *iseed <<= 1; return 0; } } ``` 方法二的优点是可以将异或操作通常作为一个全字异或操作来完成。 ##### 2.4 部分模 2 本原多项式列表 | 多项式 | 表示形式 | | --- | --- | | \(x^1 + x^0\) | \((1, 0)\) | | \(x^2 + x^1 + x^0\) | \((2, 1, 0)\) | | \(x^3 + x^1 + x^0\) | \((3, 1, 0)\) | | ... | ... | | \(x^{100} + x^{8} + x^{7} + x^{2}\) | \((100, 8, 7, 2)\) | 需要注意的是，不要将这些例程生成的连续比特用作大的随机整数的比特，或者作为随机浮点数尾数的比特，它们在这些用途上的随机性不佳。 #### 3. 基于数据加密的随机序列 传统的使用数据加密标准（DES）生成随机数的方法在软件实现时速度非常慢。这里介绍一种更快更简单的算法，虽然在加密意义上可能不安全，但能生成质量相当的随机数。 ##### 3.1 算法原理 DES 是一种“块积密码”，通过迭代应用一种高度非线性的比特混合函数来处理 64 位输入。为了快速生成随机数，我们需要一个能在高级计算机语言中快速计算的非线性函数\(g\)。 函数\(g\)的设计思路是计算 16 位输入半字的三种不同 32 位乘积，然后通过快速操作（如加法或异或）将这些乘积与一些固定常数组合成一个 32 位结果。 以下是相关代码实现： ```c #define NITER 4 void psdes(unsigned long *lword, unsigned long *irword) { unsigned long i,ia,ib,iswap,itmph=0,itmpl=0; static unsigned long c1[NITER]={ 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L}; ```