19、随机数生成：从无理数到元胞自动机

随机数生成：从无理数到元胞自动机

在随机数生成的领域中，有多种方法可以探索。本文将深入探讨利用无理数和元胞自动机来生成随机序列的方法，分析其原理、实现步骤以及测试结果。

1. 利用无理数生成随机序列

通过对一些被认为是正态分布的常数（如 e、log(10)、φ 和 √3）的十六进制展开生成的二进制文件中的随机值，能够生成蕨类分形。从图中可以看出，每个常数生成的序列在这种情况下都很有用，并且每次都能得到相同的图形。这表明，如果能够快速生成许多（甚至大多数）无理数展开的任意位数，我们将拥有一种强大的方法来生成优质的随机序列。甚至可以将序列组合起来，或者像使用现有伪随机生成器的新种子一样使用新的无理数。

2. 利用阶乘生成随机序列

我们来探讨大阶乘的数字（忽略末尾累积的零）是否可以作为随机序列。首先，我们定义整数 n ≥ 0 的阶乘为 n! = n(n - 1)(n - 2)(n - 3) … (2)(1)，并规定 0! = 1。阶乘增长得非常快，例如：
| n | n! |
| — | — |
| 2 | 2 |
| 5 | 120 |
| 7 | 5040 |
| 11 | 39916800 |
| 13 | 6227020800 |
| 15 | 1307674368000 |
| 17 | 355687428096000 |
| 21 | 51090942171709440000 |
| 33 | 8683317618811886495518194401280000000 |

在 Python 中生成大阶乘很容易，因为 Python 支持任意精度的整数运算。经典的阶乘函数递归定义在 Python 中可能会导致栈深度超限，因此我们使用迭代定义：

def fact(n):
    if (n<2):
        return 1
    s = 1
    for i in xrange(1,n+1):
        s *= i
    return s

在研究大阶乘的随机性之前，需要处理阶乘末尾的零。可以将整数转换为字符串，从后往前数找到零开始的位置，但这种方法不够优雅。我们可以使用以下算法来计算 n! 末尾零的数量：

def trailing(n):
    s = 0
    while (n > 1):
        s += int(n/5.0)
        n = n/5.0
    return s

接下来，我们需要将大十进制数转换为二进制文件，以便用于测试程序。以下是完整的代码：

import sys
import os
def fact(n):
    if (n<2):
        return 1
    s = 1
    for i in xrange(1,n+1):
        s *= i
    return s

def trailing(n):
    s = 0
    while (n > 1):
        s += int(n/5.0)
        n = n/5.0
    return s

def main():
    n = int(sys.argv[1])
    oname = sys.argv[2]
    q = str(fact(n))[:-(trailing(n))]
    if (len(q)%8 != 0):
        q = q[:-(len(q)%8)]

    with open(oname,"w") as f:
        k=0
        while (k < len(q)):
            b = 0
            for i in range(8):
                v = 1 if int(q[k])>4 else 0
                b |= v<<(7-i)
                k += 1
            f.write("%s" % chr(b))

以 500,000 的阶乘为例，它有 2,632,342 位，其中末尾 124,999 位是零，转换后的输出文件为 313,417 字节。对该文件运行 ent 测试，结果如下：
- 熵：7.999373 比特/字节
- 最优压缩：文件大小无压缩
- 卡方分布：272.63，随机情况下超过此值的概率为 25.00%
- 算术平均值：127.6336（随机值为 127.5）
- 蒙特卡罗法计算的 π 值：3.149475458（误差 0.25%）
- 序列相关系数：-0.001209（完全不相关为 0.0）

为了得到更好的结果，我们可以合并多个阶乘的输出。由于 (n + 1)! = (n + 1)n!，我们可以对代码进行修改，更高效地计算一系列阶乘并合并它们：

def main():
    n = int(sys.argv[1])
    m = int(sys.argv[2])
    oname = sys.argv[3]

    for t in xrange(n,n+m):
        if (t == n):
            z = fact(t)
        else:
            z = t*z

        q = str(z)[:-(trailing(t))]
        if (len(q)%8 != 0):
            q = q[:-(len(q)%8)]

        with open(oname+"_%d" % t,"w") as f:
            k=0
            while (k < len(q)):
                b = 0
                for i in range(8):
                    v = 1 if int(q[k])>4 else 0
                    b |= v<<(7-i)
                    k += 1
                f.write("%s" % chr(b))

对 499,991! 到 500,110! 的阶乘输出进行合并，得到一个 37,609,440 字节的文件。再次运行 ent 测试，结果有了显著改善：
- 熵：7.999996 比特/字节
- 最优压缩：文件大小无压缩
- 卡方分布：229.62，随机情况下超过此值的概率为 75.00%
- 算术平均值：127.4932（随机值为 127.5）
- 蒙特卡罗法计算的 π 值：3.141752071（误差 0.01%）
- 序列相关系数：0.000359（完全不相关为 0.0）

然而，将该文件通过 dieharder 测试和 SmallCrush TestU01 测试集时，由于文件较小，部分测试结果不太理想。进一步对文件进行基本随机性测试，结果显示，除了一些由于文件大小导致不可靠的测试（如 gap 测试）外，合并后的阶乘值趋向于随机，但构建小数据文件所需的大量大阶乘使得这种方法在实际应用中存在困难。

3. 利用元胞自动机生成随机序列

元胞自动机是由离散的“单元”组成的数组，每个单元处于有限状态之一。通过对单元的初始配置应用规则，可以将其移动到新的配置，这个过程会不断重复。在这里，我们主要关注一维元胞自动机，每个单元有两种状态：存活（1）或死亡（0），Wolfram 将其称为基本元胞自动机。

对于每个单元，其状态由当前单元及其左右相邻单元的状态决定。由于单元状态是二进制的，对于任何中心单元，有八种可能的输入组合：
| 输入组合 | 含义 |
| — | — |
| 111 | 中心单元及其左右相邻单元都存活 |
| 110 | 中心单元和左相邻单元存活，右相邻单元死亡 |
| 101 | 中心单元和右相邻单元存活，左相邻单元死亡 |
| 100 | 中心单元存活，左右相邻单元死亡 |
| 011 | 中心单元死亡，左右相邻单元存活 |
| 010 | 中心单元死亡，左相邻单元存活，右相邻单元死亡 |
| 001 | 中心单元死亡，右相邻单元存活，左相邻单元死亡 |
| 000 | 中心单元及其左右相邻单元都死亡 |

我们可以用二进制表示规则，通过指定对应于当前世界三个比特状态的输出比特来确定规则。例如，一个规则可能是：
| 输入 | 111 | 110 | 101 | 100 | 011 | 010 | 001 | 000 |
| — | — | — | — | — | — | — | — | — |
| 输出 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 |

下面是一个简单的 Python 脚本，用于应用任何规则到初始的 16 位无符号状态，并观察其随时间的演化：

import sys

def pp(j,state):
    print "%04d:" % j,
    for i in range(16):
        print "1" if ((state>>(15-i))&1) else "0",
    print

def main():
    rule = int(sys.argv[1],16)
    state= int(sys.argv[2],16) & 0xFFFF
    steps= int(sys.argv[3])
    pp(0,state)
    for i in xrange(steps):
        nstate = 0
        for j in range(16):
            if (j == 0):
                if (rule & (1<<((state&3)<<1))):
                    nstate |= 1
            else:
                if (rule & (1<<((state>>(j-1))&7))):
                    nstate |= (1<<j)
        state = nstate
        pp(i+1,state)

运行这个脚本，对于规则 66 和初始状态 0x75ab，经过 16 次迭代后，系统会灭绝。如果我们放松边界条件，让单元环绕，代码修改如下：

for i in xrange(steps):
    nstate = 0
    for j in range(16):
        if (j == 0):
            if (rule & (1<<(((state&3)<<1)|((state>>15)&1)))):
                nstate |= 1
        elif (j == 15):
            if (rule & (1<<(((state>>14)&3)|((state&1)<<2)))):
                nstate |= (1<<15)
        else:
            if (rule & (1<<((state>>(j-1))&7))):
                nstate |= (1<<j)
    state = nstate
    pp(i+1,state)

运行修改后的代码，规则 66 和相同的初始状态下，系统不会灭绝，而是比特会无限向左循环。

Wolfram 声称规则 30 是混沌的，可以用于生成高质量的随机数。我们使用单个初始单元和环绕边界条件运行规则 30，发现中心单元（最初存活的单元）的比特序列被认为是随机的。早期版本的 Mathematica 甚至使用规则 30 作为随机数生成器。

为了生成更多的比特，我们将代码用 C 语言重写：

#define N 64
uint64_t apply_rule(uint8_t rule, uint64_t state) {
    uint64_t nstate=0, j;
    for(j=0; j<N; j++) {
        switch (j) {
            case 0:
                if (rule & (1<<(((state&3)<<1)|((state>>(N-1))&1))))
                    nstate |= 1;
                break;
            case (N-1):
                if (rule & (1<<(((state>>(N-2))&3)|((state&1)<<2))))
                    nstate |= (1ULL<<(N-1));
                break;
            default:
                if (rule & (1<<((state>>(j-1))&7)))
                    nstate |= (1ULL<<j);
                break;
        }
    }
    return nstate;
}

void main(int argc, char *argv[]) {
    FILE *f;
    uint8_t b, rule=30;
    uint32_t nsamp,i,j;
    uint64_t state=0x0000000100000000;

    nsamp = (uint64_t)atol(argv[1]);
    f = fopen(argv[2],"w");
    for(i=0; i<nsamp; i++) {
        b = 0;
        for(j=0; j<8; j++) {
            b |= ((state>>(N/2))&1)<<(7-j);
            state = apply_rule(rule, state);
        }
        fwrite((void*)&b, sizeof(uint8_t), 1, f);
    }
    fclose(f);
}

运行这个程序并将输出通过 ent 测试，结果非常好：
- 熵：8.000000 比特/字节
- 最优压缩：文件大小无压缩
- 卡方分布：223.39，随机情况下超过此值的概率为 90.00%
- 算术平均值：127.5002（随机值为 127.5）
- 蒙特卡罗法计算的 π 值：3.141616935（误差 0.00%）
- 序列相关系数：-0.000014（完全不相关为 0.0）

通过 dieharder 测试，得分 225（111 个测试通过，3 个弱通过，0 个失败）。在 TestU01 的 SmallCrush 测试集中，所有测试都通过，但运行时间较长，因为生成单个输出字节需要更新整个 64 位状态向量八次，且实现较为简单，每个比特位置本可以并行更新。

对于其他 255 个规则，我们对 rule30_rng.c 进行轻微修改，在命令行指定规则值。运行所有规则（[0, 255]），每次生成 4000 万字节的输出文件，使用 ent 计算得分。部分规则的得分如下：
| 规则 | 得分 |
| — | — |
| 89 | 0.0001523691 |
| 75 | 0.0001523691 |
| 225 | 0.0002438709 |
| 169 | 0.0002438709 |
| 86 | 0.0003192051 |
| 30 | 0.0003192051 |
| 101 | 0.0003796057 |
| 45 | 0.0003796057 |
| 135 | 0.0004629253 |
| 149 | 0.0004629253 |
| 137 | 1.1949211151 |
| … | … |
| 255 | 2.2360679775 |

我们选取得分最低的前 10 个规则，使用 20 亿字节的输出文件进行 dieharder 测试，结果如下：
| 规则 | 通过数 | 弱通过数 | 失败数 | 得分 |
| — | — | — | — | — |
| 30 | 111 | 3 | 0 | 225 |
| 86 | 111 | 3 | 0 | 225 |
| 135 | 109 | 5 | 0 | 223 |
| 149 | 109 | 5 | 0 | 223 |
| 45 | 108 | 6 | 0 | 222 |
| 101 | 108 | 6 | 0 | 222 |
| 75 | 107 | 7 | 0 | 221 |
| 89 | 107 | 7 | 0 | 221 |
| 169 | 80 | 5 | 29 | 107 |
| 225 | 80 | 5 | 29 | 107 |

实际上，规则 30 及其同构规则（86、135 和 149），以及规则 45 及其同构规则（75、89 和 101）本质上是相同的基本规则。以规则 30 为例，通过对索引比特进行镜像翻转、取补等操作，可以从规则 86、135 和 149 恢复到规则 30。

综上所述，利用无理数、阶乘和元胞自动机都可以生成随机序列，但每种方法都有其优缺点。在实际应用中，需要根据具体需求和场景选择合适的方法。同时，对于元胞自动机规则 30 及其同构规则，虽然在测试中表现良好，但由于实现效率问题，在实际应用中可能需要进一步优化。

总结

本文介绍了三种生成随机序列的方法，包括利用无理数、阶乘和元胞自动机。通过对这些方法的原理、实现步骤和测试结果的分析，我们可以看到每种方法都有其独特的特点和适用场景。在实际应用中，我们可以根据具体需求选择合适的方法，同时也可以进一步探索和优化这些方法，以提高随机序列的质量和生成效率。

展望

未来的研究可以集中在以下几个方面：
1. 优化算法实现 ：对于阶乘和元胞自动机方法，进一步优化算法实现，提高生成随机序列的效率。
2. 探索更多规则 ：除了规则 30 及其同构规则，深入研究其他元胞自动机规则，寻找更多能够生成高质量随机序列的规则。
3. 结合多种方法 ：尝试将不同的随机数生成方法结合起来，以获得更好的随机性能。
4. 应用拓展 ：将这些随机序列应用到更多的领域，如密码学、模拟实验等，验证其在实际应用中的有效性。

通过不断的研究和探索，我们有望找到更加高效、可靠的随机数生成方法，满足不同领域对随机数的需求。

随机数生成：从无理数到元胞自动机（续）

4. 不同随机数生成方法对比分析

为了更清晰地了解利用无理数、阶乘和元胞自动机生成随机序列的特点，我们对这三种方法进行详细对比。以下是一个综合对比表格：
| 生成方法 | 优点 | 缺点 | 适用场景 |
| — | — | — | — |
| 无理数 | - 理论上若能快速生成任意位数，可产生优质随机序列
- 可组合序列或使用新无理数作为种子 | - 实际中难以快速生成大量无理数展开的位数 | - 对随机序列理论研究有一定价值 |
| 阶乘 | - 可通过合并多个阶乘输出改善随机性
- Python 易于实现大阶乘计算 | - 构建小数据文件需大量大阶乘，实际应用困难
- 小文件测试结果受文件大小影响 | - 对随机性要求不高且数据量需求不大的场景 |
| 元胞自动机（规则 30 等） | - 部分规则（如规则 30）生成的随机序列测试结果良好
- 基本位操作，有并行硬件实现可能 | - 实现效率低，运行时间长
- 除规则 30 外，其他规则研究较少 | - 对随机序列质量要求高且有足够计算资源的场景 |

从这个对比表格可以看出，每种方法都有其独特的优势和局限性。在实际应用中，我们需要根据具体的需求和场景来权衡选择。

5. 元胞自动机规则可视化分析

为了更直观地观察元胞自动机规则的运行情况，我们可以对规则进行可视化。前面提到了规则 30 家族和规则 45 的相关图形，这里我们进一步分析其意义。

对于规则 30 家族（规则 30、86、135 和 149），图 7.2 展示了前 64 个状态的图形表示。从图中可以清晰地看到，规则 30 和 86、规则 135 和 149 分别是镜像关系，这也解释了为什么它们在 dieharder 测试中的得分相同。这种可视化有助于我们理解不同规则之间的内在联系，以及它们在生成随机序列过程中的相似性。

而对于规则 45，图 7.3 展示了 329 步的状态图（图像旋转 90° 至左，状态从左到右）。从图中可以看出似乎存在某种秩序，但随着状态的迭代，边界效应会向中心传播，可能会消除这种结构。这提示我们在研究规则 45 生成的随机序列时，需要考虑到状态迭代过程中结构的变化情况。

下面是一个简单的 mermaid 流程图，展示了元胞自动机规则可视化的基本流程：

graph TD;
    A[选择元胞自动机规则] --> B[设置初始状态];
    B --> C[迭代更新状态];
    C --> D[记录每个状态];
    D --> E[生成可视化图形];

6. 随机数生成方法的实际应用思考

随机数在许多领域都有广泛的应用，如密码学、模拟实验、游戏开发等。不同的应用场景对随机数的质量和生成效率有不同的要求。

在密码学领域，随机数的质量至关重要。密码系统需要高质量的随机数来生成密钥，以确保系统的安全性。对于这种场景，元胞自动机规则 30 及其同构规则可能是一个不错的选择，因为它们在测试中表现出了较好的随机性。但由于其实现效率低的问题，需要进行进一步的优化，例如采用并行硬件实现，以满足密码学对随机数生成速度的要求。

在模拟实验中，随机数用于模拟各种随机现象。如果实验对随机数的质量要求不是特别高，且数据量需求不大，阶乘方法可能是一个可行的选择。它相对容易实现，并且可以通过合并多个阶乘输出来改善随机性。

在游戏开发中，随机数用于增加游戏的趣味性和不确定性。对于这种场景，可以根据游戏的具体需求选择合适的随机数生成方法。如果游戏对随机数的质量和生成效率都有一定要求，可以考虑结合多种方法，例如将元胞自动机生成的随机数与其他伪随机数生成器相结合。

7. 未来研究方向的深入探讨

前面提到了未来研究的几个方向，这里我们进一步深入探讨每个方向的具体内容和挑战。

• 优化算法实现 ：对于阶乘方法，优化算法实现可以从减少计算量和提高存储效率入手。例如，可以研究更高效的大整数乘法算法，以减少计算大阶乘的时间。对于元胞自动机方法，并行硬件实现是一个重要的研究方向。可以利用 GPU 或 FPGA 等硬件平台，实现元胞自动机规则的并行计算，从而提高生成随机序列的效率。但并行硬件实现需要解决硬件资源分配、数据传输等问题，具有一定的技术挑战。

• 探索更多规则 ：除了规则 30 及其同构规则，元胞自动机还有 255 个其他规则。深入研究这些规则，需要建立一套有效的评估体系，以筛选出能够生成高质量随机序列的规则。可以结合数学分析和实验测试，对每个规则进行全面评估。同时，还需要研究规则之间的关系，探索是否存在其他同构规则家族。

• 结合多种方法 ：将不同的随机数生成方法结合起来，可以充分发挥每种方法的优势，获得更好的随机性能。例如，可以将无理数生成的随机序列作为种子，输入到元胞自动机中，以增加元胞自动机生成随机序列的随机性。但结合多种方法需要解决不同方法之间的兼容性问题，以及如何合理地融合不同方法生成的随机序列。

• 应用拓展 ：将随机序列应用到更多的领域，需要考虑不同领域的具体需求和特点。在密码学领域，需要确保随机数的安全性和不可预测性；在模拟实验中，需要保证随机数能够准确模拟实际的随机现象；在游戏开发中，需要考虑随机数对游戏平衡性和趣味性的影响。因此，在应用拓展过程中，需要进行大量的实验和验证，以确保随机序列在实际应用中的有效性。

总结与再次展望

本文全面介绍了利用无理数、阶乘和元胞自动机生成随机序列的方法，详细分析了每种方法的原理、实现步骤、测试结果以及优缺点。通过对比分析和可视化研究，我们对这些方法有了更深入的理解。

在未来的研究中，我们需要针对每种方法的局限性进行优化和改进，探索更多的可能性。同时，结合多种方法和拓展应用领域也是非常有意义的研究方向。相信通过不断的努力，我们能够找到更加高效、可靠的随机数生成方法，满足不同领域对随机数的多样化需求。

希望本文能够为对随机数生成感兴趣的读者提供有价值的参考，激发更多的研究和创新。让我们共同期待随机数生成领域在未来取得更大的突破！