在上一篇文章中，我提到了 V8 引擎存在的一个 bug，对于数组的长度计算有点问题，导致在部分算法执行过程中出现 Invalid array length 异常而退出。有趣的是，Object 也会出现同样的 Invalid array length，尽管 Object 本身并没有 length 属性。

我写了个最简单的复现代码：

1
2
3
4
5
6
7
8
9

const n = 2 ** 24 - 1; // 16777215
// const dict = Object.create(null) // same bug
const dict = [];
for (let i = 0; i < n; i++) {
    let key = i * Math.floor(Math.log(i));
    dict[key] = i + 1;
}
const maxkey = (n - 1) * Math.floor(Math.log(n - 1)); // 279097901
console.log(dict[maxkey]); // should be 16777215

在最新版本的 nodejs 和 chrome 控制台中，均出现相同的报错。但在 firefox 和 safari 的控制台中都可以正常执行，这两个浏览器都有自己的独立 JS 引擎。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

node object-dict-test.js
/Users/kaikai/Coding/object-dict-test.js:7
    dict[key] = i + 1
              ^

RangeError: Invalid array length
    at Object.<anonymous> (/Users/kaikai/Coding/hexo-blog/source/attach/2025/08/object-dict-test2.js:7:15)
    at Module._compile (node:internal/modules/cjs/loader:1738:14)
    at Object..js (node:internal/modules/cjs/loader:1871:10)
    at Module.load (node:internal/modules/cjs/loader:1470:32)
    at Module._load (node:internal/modules/cjs/loader:1290:12)
    at TracingChannel.traceSync (node:diagnostics_channel:322:14)
    at wrapModuleLoad (node:internal/modules/cjs/loader:238:24)
    at Module.executeUserEntryPoint [as runMain] (node:internal/modules/run_main:154:5)
    at node:internal/main/run_main_module:33:47

Node.js v24.6.0

查了一下原因：

1. V8 对 整数键（即符合数组索引规则的属性名）的 Object 有特殊处理。当普通对象使用整数作为 key 时，V8 通常会将其视作“类似数组”的元素：无论是 Array 还是普通对象，都有一个内部的 “elements” 存储空间来保存这些整数下标对应的值。这解释了 bug 为什么同时出现在了 Array 和 Object 上。

2. V8 在管理大量元素时区分快速 Fast 和字典 Dictionary 两种存储模式。当数组规模较小并密集时，使用快速模式；若遇到稀疏或非常大的索引，就可能切换到字典模式。在快速模式下，向新索引赋值会自动更新 length 并尽量扩容底层数组；超过阈值后会退化到字典模式，用哈希表存储元素。同时 V8 为了优化性能，设置了多个内部变量用于判断数组下不同的数据类型并进行优化。官方博客

3. 在隐式转换过程中，V8 忘了对内部变量进行判断，于是一个原本合法的 key 值因为大于内部变量限制，而导致了 Invalid array length 异常。这个属于 V8 长期存在的 bug，至今没有修正。

再深入的分析通过搜索已经找不到答案了，需要去查看庞大的 V8 源代码，超出了我的时间余裕和兴趣能力，就这样吧。

问题

在知乎上看见了个有意思的问题：

孤独素数：
一个素数 P 通常可以分为两个殆素数之和，其序数与对应的殆素数的素因子序数同时满足：
\begin{equation*}
\begin{cases}
P = \prod p_i + \prod q_i \\
\pi(P)=\prod \pi(p_i) + \prod \pi(q_i)
\end{cases}
\end{equation*}
例如 $11 = 2*3 + 5$，同时有 $\pi(11) = 5 = \pi(2)\pi(3) + \pi(5) = 1\cdot2 + 3。$
但也存在一些素数，任意分法均不能满足上述规则，称为孤独素数。

看见这题我是懵逼的，懵逼树上懵逼果，第一步直接撞死在黎曼猜想了，想不到在 $\pi(p)$ 还没有明确公式解的情况下怎么去判断孤独素数的规律。但至少这个看起来是个简单而有趣的编程练习题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

// 欧拉筛法获得 n 以内的素数表
function EulerSieve(n) {
    const primes = [];
    let eularboard = new Array(n + 1).fill(true);
    for (let i = 2; i <= n; i++) {
        if (eularboard[i]) {
            primes.push(i);
        }
        for (let j = 0; j < primes.length; j++) {
            if (i * primes[j] > n) break;
            eularboard[i * primes[j]] = false;
            if (i % primes[j] === 0) break;
        }
    }
    return primes;
}
const primes = EulerSieve(10 ** 6);

// 计算殆素数的素因子序数乘积，并使用使用 dp 避免重复的因式分解
const dp = [];
function factorOrdersProduct(n) {
    if (n === 1) return 0;
    if (dp[n]) return dp[n];

    const init = n;
    const sqrt_n = Math.floor(Math.sqrt(n));
    let indices = 1;
    for (let i = 0; i < primes.length; i++) {
        const p = primes[i];
        if (p > sqrt_n) break;
        while (n % p === 0) {
            indices = indices * (i + 1);
            n = n / p;
            if (dp[n]) {
                indices = indices * dp[n];
                n = 1;
                break;
            }
        }
        if (n === 1) break;
    }

    if (n > 1) {
        const idx = primes.indexOf(n);
        if (idx === -1) throw new Error(`Math crisis ${n}`);
        indices = indices * (idx + 1);
    }

    dp[init] = indices;
    return indices;
}

// 寻找孤立素数，将每个素数分解成两个殆素数，然后计算殆素数的素因子序数乘积
const lonelyPrimes = [];
for (let i = 0; i < primes.length; i++) {
    const p = primes[i];
    let lonely = true;
    for (let m = 2; m < p / 2; m++) {
        if (factorOrdersProduct(m) + factorOrdersProduct(p - m) === i + 1) {
            lonely = false;
            break;
        }
    }
    if (lonely) {
        lonelyPrimes.push(p);
    }
}

console.log(lonelyPrimes);
// 2, 3, 211, 277, 541, 631, 653, 1259, 2797, 2897, 4391, 5279, 5323, 5381, 5527,  6113, 6547, 7177, 7523, 7993, 8147, 8513, 9157, 9661, 10211, ...

实测下来 Apple M1 上计算 $10^6$ 以内孤独素数约 3.5 秒，$10^7$ 以内约 180 秒。

优化

整体来看上述代码是清晰的，先用欧拉筛获得素数数组，再根据孤独素数的定义逐个验证。由于殆素数分解是必然有解的，实际上就是验证第二条件序列数的运算是否成立。通过 dp 暂存已经计算过的结果，可以大大提高效率。但直观来看，还有一些可以优化的小点：

1. 在规模不大的时候，可以使用字典对象储存素数与其序数的键值对，替代 primes.indexOf(n)。一个对象可以容纳 2^24 个键值对，即最大素数为 $p_{(2^{24})}$，查询得知 $p_{16777216}=310248241$，足够大了。更大规模也可以使用若干个字典对象来应对。
2. 如果一个殆素数因子序数乘积已经超过了当前素数序数，可以提前结束。但需要考虑每次增加的判断带来的损耗是否小于节约的时间。考虑到：$\pi(p) \approx p / \ln p$，殆素数 m 大概率为合数， $\pi(m) <\prod_i(m_i/ln m_i), m = \sum_i m_i$，m-q 同理。而 $lnp$ 显然是大于 1 的，因此可以认为这种提前结束的情况几乎不存在。在极端情况下，m=2，p 与 p-2 为孪生素数，也只是等于序列数的和。因此这一条更大可能是负面的。
3. 使用 UInt32Array 等定长数组优化，这个我不是很喜欢。我知道它有效，但总觉得将语言特性也视为算法本身的优化的话，汇编就是唯一可用的语言了。

以及一个可能存在的大优化点：

在欧拉筛获得素数的过程中，我们已经把所有的合数都计算过一遍了

这意味着 factorOrdersProduct(n) 函数的计算过程，与筛法取素数的计算过程本身是有大量重合的。在题目的定义下，素数 p 的两个加数几乎就覆盖了从 1-n 的所有合数，也就是筛法求素数中被筛去的数，以及筛选这个动作本身也是有价值的。

值得一提的是，在因数分解相关算法中有个 SPF(Smallest Prime Factor) 算法。它与欧拉筛是完美配合的。简单地说，SPF 分为两步：

1. 在筛法取素数过程中，保留合数而不是直接筛去，并记录每个数最小质因数。
2. 在因数分解过程中，使用类似递归的方法，根据最小质因数快速计算各合数的全部因子。

在第 1. 步中将 SPF 与欧拉筛结合，并将第 2 步中的计算因子数组改为计算因子序数乘积，即可得到更高效的搜搜代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

const n = 10 ** 6;
const spf = new Uint32Array(n + 1).fill(0);
const primes = [];
const primeIndex = new Map();
const product = new Uint32Array(n + 1);

// 构建SPF和素数序号
for (let i = 2; i <= n; i++) {
    if (spf[i] === 0) {
        spf[i] = i;
        primes.push(i);
        primeIndex.set(i, primes.length); // 记录素数序号
    }

    for (let j = 0; j < primes.length; j++) {
        const p = primes[j];
        const ip = i * p;
        if (p > spf[i] || ip > n) break;
        spf[ip] = p;
        if (i % p === 0) break;
    }
}

// 直接计算因子序数乘积
product[1] = 1;
for (let i = 2; i <= n; i++) {
    const p = spf[i]; // 最小质因数
    const prev = i / p;
    product[i] = product[prev] * primeIndex.get(p);
}

// 寻找孤立素数
const lonelyPrimes = [];
for (let idx = 0; idx < primes.length; idx++) {
    const p = primes[idx];
    const target = idx + 1; // primeIndex[p]
    const half = Math.floor(p / 2);
    let lonely = true;
    for (let m = 2; m < half; m++) {
        if (product[m] + product[p - m] === target) {
            lonely = false;
            break;
        }
    }
    if (lonely) lonelyPrimes.push(p);
}

优化结果是 $10^6$ 以内约 0.8 秒，$10^7$ 以内约 52 秒，$10^8$ 约 8900 秒。

在 JS 实现中，对于相同的查询请求，性能排序大约是 Uint32Array > Array $\approx$ Object > Map。但 V8 引擎存在 bug，导致 Array 或 Object 在用作 primeIndex 时出现了预料之外的错误，不得不用 Map()。这就是为什么我不喜欢基于语言特有数据结构的优化。而相同的 python 代码虽然可以天然地支持更大的数字而无需考虑这些异常，也无需考虑默认 2GB 的内存上限，但代价是慢得多的计算速度以及单位计算量下更大的内存占用。在计算规模上升后还不如 nodejs 实用。

观察

在有了高效的孤独素数搜索算法以后，可以对其序列进行一些观察。

1. 分布。

   从计算结果看，孤独素数很稀疏，且存在较大的波动。但计算到 1 亿，以 10 万为一段统计其出现的频率，发现和素数本身的分布趋势是大致一致的。

   

   绿色为素数数量分布参考线 $\pi(p)=p / \ln p$

2. 解的分布。

   对于非孤独素数，第一个解往往在 m 较小时就出现了。一般地，对于 p = m + q 的情况，上述算法是从 m = 2 开始逐项验证的，等效于从 p 左边的所有数中取任意两个，使 $\pi(p)$ 有解。先来看一下所有数的因子序数积的分布情况。

   

   前 1000 个整数的因子序数积的分布

   

   前 30000 个整数的因子序数积的分布

   从图中我们发现，整数的因子序数积的分布出现了明显的原点出发放射状图样，最上面一条显然是质数和其对应的序列数。随着 x 轴放大点数密集，放射线变得更为清晰。因此我们进行局部放大并标注具体数据后如下：

   

   29000 到 30000 个整数的因子序数积的分布及数值标注

   经过简单计算，就可以发现素数线以下的第二条线是 $3\times素数$，再下面第三条线是 $5\times素数$。

   根据图像发现一个推论，若素数 p' 与 p $\frac{p'}{p} \approx \frac{1}{3}$，则 $\frac{1}{3} < \frac{\pi(p')}{\pi(p)} < \frac{1}{2}$，所以才有 $3\times素数$ 线稳定在素数线下。

   当 p 足够大时，$\pi(p)\approx p / ln(p)$，假设存在一个与 p 几乎相等的殆素数 3p'。则有：
   $$ \frac{\pi(3)\cdot\pi(p')}{\pi(p)}\approx \frac{2p}{3ln(p/3)} / \frac{p}{lnp} = \frac{2}{3}\cdot\frac{lnp}{lnp-ln3}$$
   令 $t=lnp$, $c=ln3$，则比值 $r(t)=\frac{2}{3}\cdot\frac{t}{t-c}$，其中$t>c$ ，即$p>3$。

   • 若 $r(t)>1$，则 $2t>3t-3c \Rightarrow t<3c \Rightarrow lnp<ln27 \Rightarrow p<27$
   • 若 $r(t)=1$，则 $p = 27$
   • 若 $r(t)<1$，则 $p > 27$

   也就是说 在 p 足够大时，$\pi(p) > \pi(3)\cdot\pi(p')$，因此殆素数 3p' 构成的因子序数积的线一定在 $\pi(p)$ 下方。

   以同相方法，容易比较 2p'、5p' 等不同殆素数构成的线条的相对位置。进一步地，n 阶殆素数也可以以类似的方法确定相对顺序。

   特殊地，因为$\pi(2)=1$，所以 $2^n$ 是最底部恒等于 1 的横线。

3. p = m + q

   根据上一条的结论，对于一个殆素数，每增加一个因子，它的因子序数积都是变小的，因为对任意 p 都有 $\pi(p)<p$。所以通常来说，其中一个值是较小的质数，会让等式成立的概率更大。

   如果两个都是殆素数，最优情况也至少是 $p=3p_1+5p_2$，如果 p 足够大，此时要求 $p/lnp = 2p_1/lnp_1+4p_2/lnp_2$，若 m 与 n 大小类似，容易导致 $p_1, p_2$ 过小。

   也就是说，$\pi(p)=\prod \pi(p_i) + \prod \pi(p_j)$ 等式要成立，可以调节的空间大小就是 $p/lnp - (p_i/lnp_i + p_j/lnp_j)$，这个当然没有稳定解，但从趋势上容易理解，m q 差异越大的，殆素数因子越少的，越容易得到解。

4. 特例与猜想

   • 孪生素数中较大的一个，一定不是孤独素数，这结论还挺文学的。
   • 有无穷多个孤独素数。
   • 指定 m 为任何 2 以上整数，都可以找到 p 与 q 使得公式成立。

得不出什么普遍性的结论，思而不学民科一下。

这是对暗黑3游戏过程中的一些计算结果的整理，也包括之前发布于论坛上的帖子的内容整合。由于与游戏内容高度相关，本文不会解释一些游戏名词与常识，默认读者已知晓。同时，本文提及的数值计算与运作方式已经过本人测试，由他人测试后的结论或历史结论均会单独注明。

伤害加成与独立计算

计算增幅时的基底

总的来说，计算装备更替技能调整的伤害变化，总是要 新属性 / 旧属性 两者比较的。因为 D3 装备众多，大多数时候你并不需要全身属性全算一遍，只直接比较受到影响的部分属性就行。下文关于增伤类的讨论就会详细讲述哪些互不影响。

这时候唯一要注意的就是，当你已经有部分属性时，在同属性上再加些数值，带来的实际增长比率要比看上去小。这些 已有的属性 作为分母，缩小了装备新增/改善带来的增伤，使实际增益不如装备上描述的那么高了。

你已经有 20000 敏捷，再加 5000 只增长了 25%；
你已经有 50 级困者之类(+30%)，再升到 150级(+60%)，实际只增长了 30% / (1+30%) = 23%。
你已经有 30% 面板的技能伤%，再加个 50 级太级石(+40%)，实际只增长了 40% / (1+30%) = 31%

你已经有的部分会降低新增部分带来的实际收益，我们把这些旧装备 已有的数值 称为基底值。

下文我们会看到，有些属性的基底值只要考虑这件物品属性本身就行，另一些属性的基底值可能还要考虑别的装备的影响。

伤害分类与倍率

得益于早年的科普帖：https://bbs.nga.cn/read.php?tid=8111152，『A 伤』的概念深入人心。但是这帖子里，作者只是语焉不详地提到：『除了其它伤以外都是 A 伤』，实际应用时仍有很多迷惑的地方。尽管该文明确提到了『相同类型是先相加，再和不同类型相乘』，但该文作者对分类的表述不甚严谨，分类有时用于乘法间隔，有时又用于整理性质的分组，甚至使用了大类小类这种二级分类，同时将各分类主观命名 ABCD。这些不严谨的表述可能造成理解上的歧义。

考虑到习惯因素，本文有时仍会沿袭 ABCD 类这个称呼方便对应解释，但会严格厘清并定义伤害分类的区别，并尝试溯源 A 类多且杂的原因，并给出列表。

从计算角度讲，我们把 “类” 作为乘法的分隔，存在多个条目以加法形式共同作于增伤的，称为一类。同类内部需要相加(additive)，不同类间需要相乘(multiplicative)。于是相比于原文，这里的分类的称呼稍有变化：

1. A 类仍然以 A 伤称呼，因为其多且杂，但会给出列表及判定方法。
2. B / C / D 类不再以字母称之，其有更直接明确的名字，“元素伤” / “精英伤” / “亡灵伤” 等等。
3. 将“伤害类别” 概念与 “同类中的条目/词缀数量” 明确分隔为两个概念。同一类中可以有多条增伤，也可以只有一条。 如果某类里目前只有一条增伤，那就是独自成类，也就是大家一般说的独立增伤了。
4. 放弃出于整理或文学性质的分组，严格将类作为计算用词汇，同类内部相加(additive)，不同类间相乘(multiplicative)。

玩家的每一次攻击，其伤害可以被精确计算如下：

\begin{align*}
攻击伤害 = DPH &\times 技能自身描述 \times (1 + 技能伤害加成 + 指定的技能加成) \times (1 + 暴击几率\times暴击增伤) \\
&\times (1 + 主属性增伤 ) \times (1 + 元素属性增伤) \times (1 + 精英目标增伤) \times (1 + 亡灵目标增伤) \times \cdots \\
&\times (1 + 武器特效) \times (1 + 副手特效) \times (1 + 指环特效) \times (1 + 腰带特效) \times \cdots\\
&\times (1 + 困者宝石特效) \times (1 + 贼神宝石特效) \times (1 + 至简宝石特效) \times \cdots\\
&\\
= DPH &\cdot 技能描述 \cdot (1 + 技能伤 + 指定技能伤)(1 + 暴率\times暴伤)\cdot\prod(1+各类增伤)
\end{align*}

DPH： 这是一件武器，上面有很多不同的条目，每个条目都为玩家的全局属性加一些数值。其中 DPH (Damage per Hit) 是武器的一项属性，表示一次攻击造成的伤害。在图中 2036-2364 确定了武器的 DPH，虽然存在上下限，但通常来说计算伤害时使用其平均值即中间值 $\frac{2364+2036}{2} = 2200$ 即可。

注意不要看大白字，那是『每秒伤害』，是综合了每次攻击造成的伤害，和每秒攻击次数得到的数值，但攻击频率和要计算的每次伤害无关。它会影响到最终的总输出，但不会影响到每次伤害的计算。

技能自身描述：

每次陨石会造成 740% 的伤害，并且在3秒内额外造成 235% 的伤害。也就是 $2200 * 740\% + 2200 * 235\%$。

1 + 技能伤害加成 + 指定的技能加成

在人物的详细信息数据面板里，我们可以看到第二行『技能伤害加成』20%，最后一行『陨石术伤害加成』52%。这两个数字是需要合并计算的，也就是说，如果你放的是陨石术，那么会在这两条一共享受到 72% 的增伤，如果你放的是其它技能，则只享受 20% 的增伤。

这也是旧科普所说的 『A 伤』。所有你选择的技能、装备、被动天赋，如果其增伤体现在这个人物面板的这两个条目内，就属于 A 伤。其伤害增加值需要合并到一个整体中，增量比例并没有条目描述的那么高。比如魔法师主动技能 魔法武器、魔星-烈焰火花，被动技能 玻璃大炮、坚定意志，猎魔人的箭术高超（弓），传奇宝石的太极石增伤，都是 A 伤。某些效果并不会展示在面板中，而是对怪物产生影响的，也是 A 伤，比如魔法师被动技能 冷血。对于后者属于哪一类，需要具体测试才能确认。

1 + 暴击几率✕暴击增伤

详细信息面板中的六七两行数据，61.5% 和 +423%。因为存在『不暴击』的情况，其实际公式是：

\begin{align*}
伤害加成 &= 暴击率 * 暴击伤害 + 不暴击 * 不暴击伤害 \\
&= CR \times (1+CE) + (1-CR) \times 1 \\
&= 1 + CR \times CE
\end{align*}

注意整个暗黑3里，所有的增伤描述都是以『额外增加 xx%』展示的，当增伤项值为 0 时，表示其没有额外增伤效果，而不是让整个乘法链条归 0。

1 + 主属性增伤

详细信息面板的第一行数据『智力伤害加成』33909%，增加了 339 倍伤害。不同人物的主属性不一，非主属性不带来增伤效果。

1+其它增伤

详细信息面板已经列举了精英伤和元素伤，它们内部都是相加(additive)的，也就是说，像技能增伤一样，有多个位置都提供了元素伤害，会被一并加和到这个条目中，然后再以乘法和其它条目相乘(multiplicative)进行最终增伤的计算。

传奇宝石（太极石除外）、单件装备增伤等都不在面板中显示，主要因为只有多条相加(additive)的数据，才需要额外显示方便一并展示，而只有单条相乘(multiplicative)的数据，不需要额外显示，也没有那么多空展示。

以目标的存活时间来计算真实增伤

乘数性增伤和加法增伤计算很直观，基数是什么，多少百分比。但其实这都是表像，归根结底，更高的伤害都是为了更快地干掉怪物冲更高的层，因此增伤的真实计算方式其实是：

$$增伤倍率 = 怪物原存活时间 / 怪物新存活时间 - 1$$

比如原本怪能活 300 秒(5 分钟)，现在怪只能活 60 秒，那你的实际火力无疑是原来的五倍，也就是增伤了 400%。

对于全额百分比增伤的伤害，增伤比例和怪存活时间的缩短是一致的，因此直接讨论增伤即可。但有些场景下则不能直接算，比如 伏击、脆弱 等与敌方血量或其它非全程状态相关的增益。

伏击： 原本怪有 $1000$ 血你每秒攻击 $10$ 血，怪能扛 $100$ 秒。伏击提供了前 $250$ 血时每次攻击 $14$ 滴，后 $750$ 血不变，则怪能存活 $250/14+750/10 = 92.857$ 秒。因此伏击的真实增伤为 $100/92.857-1 = 7.69\%$

有死者遗物箭袋时的伏击： 原本怪有 $1000$ 血你每秒攻击 $10$ 血，遗物箭袋下怪能存活 $400/10+600/20=70$ 秒。伏击提供了前 $250$ 血 $40%$ 增伤，则怪的存活时间为 $250/14+150/10+600/20=62.857$ 秒。伏击的真实增伤为 $70/62.857 = 11.36%$。多重套配伏击就是这么来的。

脆弱： 直接杀死血量低于 85% 的敌人。原本活 $100$ 秒的敌人现在活 $85$ 秒，等效增伤是 $100/85-1 = 17.647\%$。如果有受罚者，则这个技能的等效增伤会减少，具体计算放在后面。

受罚者之灾： 受罚者之灾的存活时长的缩短和存活时间/攻击次数以及多怪分担数量相关，攻击又分为像多重这种每下均匀的，和刀扇这种极端的。考虑到大部分配装都是不断重复主力伤害技能，因此我根据前一种大概跑了一下随机过程模拟，结论如下图：

在数值模拟中我设定了极端的上限 1000 次攻击，通常来说一趟大秘境大约需要花费 4-5 分钟左右击杀 boss。以每秒 1.5 次攻击的频率，总共需要进行 350-450 次攻击，则受罚者加成大约在 180% 左右。

有受罚者之灾时的脆弱：

由于受罚者之灾是逐渐增伤的，因此原本 脆弱 等效的最后 15% 处于高受罚叠层的状态，本身的存活时间并不占全部存活时间的 15%。具体计算如下：

在有受罚者之灾（以 150 级石头每次增伤 2.3% 为例）的情况下，总共 n 次攻击的总伤害为 $\frac{1 + 1 + (n-1)2.3\%}{2}\cdot n$，而需要完成其 85% 伤害的 $n'$ 次攻击的总伤害为 $\frac{1 + 1 + (n'-1)2.3\%}{2}\cdot n'$。即求解：

$$n'^2+\frac{1977}{23}n' = 0.85\cdot(n^2+\frac{1977}{23}n)$$

而要求的是 $k=n/n'$ 沿 n 的变化，这个式子虽然具体数值唬人，但无论是具体 n 对应的 k ，或者 k 的极限，都不算困难。在前文所述的 350-450 次攻击时，脆弱 的等效增伤约为 9% ~ 9.5%，而极限则可以直接消去一次项，得

$$\lim_{n',n \to +\infty} k =\frac{1}{\sqrt{0.85}}\approx 1.08465$$

即，在有受罚者之灾的情况下，脆弱 的等效增伤不低于 8.465%，不高于 15%，通常在 9% 多一些。

有受罚者之灾时的伏击：

我原以为和 死者遗物 类似，受罚者之灾 缩短了后半程的存活时间，会变相提升 伏击 的价值。计算方法完全类似，原本 n 次攻击的总伤害为 D，现在前 25% 血量需要 a 次攻击，，后 75% 血量需要 b 次攻击，则：

\begin{equation*}
\begin{cases}
D = \frac{1 + 1 + (n-1)2.3\%}{2}\cdot n \\
D \cdot 25\%= 1.4\cdot\frac{1 + 1 + (a-1)2.3\%}{2}\cdot a \\
D \cdot 75\%= \frac{1 +a\cdot2.3\%+ 1 + (a+b-1)2.3\%}{2}\cdot b \\
E = \frac{n}{a+b}-1
\end{cases}
\end{equation*}

但当实际计算时，才发现 $\lim_{n \to +\infty} E \approx 1.03775$，在通常的 350-450 次左右也只有 4%，甚至比原本的更差了。细究原因，主要是因为受罚者本身的伤害严重依赖于叠层，因此虽然 伏击 形式上是额外削减了前期 10% 的怪物生命，但相当于把整个叠层进度整体往后挪了，真正削减的实质上是最后 10% 血量对应的最高受罚叠层攻击时间，因此其最终折算增伤甚至不如平稳伤害的 7.69%。

当然，大密境有 2/3 以上的时间是在清理小怪。对小怪而言受罚的效果几乎没有，此时伏击仍然发挥了原有的作用，要根据自己实际的冲层情况选择。

受罚者之灾的叠层规则与计算

150 级宝石提供每次攻击 2.30% 的伤害叠加。这个伤害是加法叠加的，也就是第一次攻击 100% 伤害的话，第二次就是 102.3%，然后是 104.6% 106.9% 109.2% 111.5%……。

但是，每次攻击后有 CD，CD 期间再有其它攻击不再叠层 ，典型比如猎魔人 扫射。扫射一次技能施放有三次实际攻击，也就是伤害为 100% 一次, 102.3% 3次, 104.6% 3次，类推。

内置 CD 具体时间为 0.9/面板攻速，确保大多数技能每一次释放都能叠一层(不妨仔细想想)。某次改版后对 杨的反曲弓这类特定技能加速也有效，杨弓多重受罚的内置 CD = 0.9/面板攻速/2，确保杨弓每次攻击也都能叠一层。

内置 CD 具体多长，与有效叠层的那次攻击的攻速有关。而实际叠层还与子弹飞行速度有关，怪的距离够远，杨弓多重后的闪避射击有可能不叠加。扫射带出来的技能必然不叠加。火多重的导弹不叠加。但如果多重空了，则导弹叠加。

有玩家曾分析称，巫医的火墙之类高频小额伤害的，可以偷一些额外的叠层出来。比如你一个技能本来是1秒释放一次，虽然内置 CD 是 0.9秒，但你下一次放技能还是得 1 秒以后，实际只能一秒叠一层。但如果你这个技能是 1 秒放一次，实际伤害是分 10 次计算每次 0.1 秒。那么在过了 0.9 秒以后，你第一发技能的第 10 次结算就能叠层，第二发技能的第8次计算能叠层，依此类推，你的 10 发技能就可以叠出 11 层。类似的情况还有魔法师塔拉夏陨石：陨石术有伤害，砸地后续有地面伤害，又有风暴护甲电人，又有魔星补输出，各种伤害高频地打向boss，其叠层也比正常的攻速快一些。

一个可以击中多个怪的技能，例如多重，也只会对其中的一个怪叠一层(默认为第一个击中的怪)，其它怪仍然有伤害但不叠受罚。对于箭雨这类一个技能有多次独立范围伤害的，多波伤害也和巫医一样，需要在 0.9 倍攻击时间后的下一波才能再叠一层。

简单来说：在可触发叠层时，为第一个受到伤害的敌人叠加一层，然后就进入 CD 期。直到 CD 期结束进入下一次叠层计算，CD 时长为攻击间隔的 0.9 倍。

A 类的特殊性与原因考究

先看图：139 级太极石，面板直接增加 75.6% 技能伤害，是 A 类。

箭术高超(弓)，面板增加 8% 技能伤，是 A 类。

无独有偶……有仨……有n……魔法法爷的多个技能和被动，巫医的穿透迷障等被动也是 A 类。

稳固瞄准(猎魔人被动)，不在面板展示，独立增伤。伏击((猎魔人被动)，不在面板展示，独立增伤。两者也相互独立。

同样是被动，差在哪里？

在纷繁的内容中我最终还是找到了线索：伏击是 64 级被动。

我们知道，暗黑 3 原版并不算成功，而其资料片夺魂之镰则带来了大密境、天梯、橙装特效等一系列改进，指明了奈非天的飞升之路。大密境指数性的难度增长，也要求玩家伤害也有对应的指数增长，于是带来了大量新的增伤类目，且确定了增伤类目间的乘法运算规则。

那么原版的伤害加成呢？如你所料，都是加法。原版的计算逻辑就是把所有的增伤，都视为对全部技能的增伤，直接加到面板上。词缀提供的特定技能增伤，也单独加在面板的技能伤害上，而这两条之间，又仍是以加法叠加。

所以 64 级新增的伏击，作为全新的增伤被动以独立类目的形式加入了游戏，而箭术高超则仍是原版的加法逻辑，至今未作修改。

那么稳固瞄准，尽管是 16 级解锁被动，按上述推测应该也是夺魂之镰才加入的新被动。

为了验证这个观点，我查到原版暗黑3于2012年5月15发行，夺魂之镰于2014年3月25日发布。但确实稳固瞄准被动在2012年就有讨论帖了，并非 70 级添加的。猜想错了么？

经过一番搜索，我找到了这个视频(已从油管转载到 B 站)：https://www.bilibili.com/video/BV1qFVyzMEpf/?t=630

视频中非常明确地显示了，在原版 60 级的暗黑 3 中，稳固瞄准是 20% 以加法添加到所有技能伤害栏，而箭术高超则以加法 15% 添加，两者共同提升了 35% 的所有技能伤害。

猜想成立。

到 70 级后，稳固瞄准改为了独立增伤，箭术高超则不但没改独立增伤，还将数值从 15% 下调到了 8%，我猜这可能一半是有意为之，一半是忘了改。具体见下文A伤测试。

由于时间古老搜索不便，只做了稳固瞄准、力士护腕等少数的搜索验证，但可以断言暗黑3原版的增伤全是加法，暴击暴伤是唯一的乘法类目。“不同类别增伤相乘” 是从夺魂之镰资料片，随着大秘境的添加而添加的新计算逻辑。由于要新旧兼容，因此有大量的算法和文本要逐步修改，导致了一些虽然文字描述都是“所有伤害增加”，但实际效果却不一致的情况。

结论：A 伤多且繁杂，绝大多数都是原版暗黑3遗留问题，并且在版本更新中也在不断改少，但个别是有意为之，比如太极石和力士护腕。

A 类增伤的测试列表

其它直接体现在面板上技能伤和具体技能伤的装备就不测了，只提几个有说法的。

稳固瞄准： 这个被动为弓/双手弩/双持单手弩提供了三种不同加成。其中，双手弩 +50% 暴伤在人物有 62% 暴率 460% 暴伤的情况下，可以提供 8.1% 的增益。

如果人物双持单手弩，有双武器绿宝石提供暴伤的情况下，人物满暴击暴伤属性大约在 45% 550% 或更多，此时该被动 6% 的暴击能提供 9.5% 的增伤。同时还有回憎效果，以抵消武器词条不如箭袋多的劣势。

人物若带弓与箭袋，则这个被动为 A 伤，在三条技能增伤的情况下，只能提供 $8\% / 145\% = 5.52\%$ 的增伤，且不说还有太极等其它 A 伤增高基底，使其增益更少。

考虑到从 60 级到 70 级时，暴雪特意将该被动的弓增伤从 15% 下调到 8%，因此可以大胆推测，这个被动的 A 伤是暴雪忘了改类目，它应当也被改为独立增伤，这样三种不同的武器才有大致相当的收益。

另一个证据是，伏击 被动在无遗物箭袋的情况下，其收益也就是 $1-1/(0.25/1.4+0.75)=7.69\%$ 的收益。按暴雪一贯以来 “根据情况有不同选择” 的设计思路来看，稳固瞄准的弓类没理由不是独立增伤。

所以双枪机轮冰吞流其实是可以考虑一下这个稳固瞄准被动的。

太极石： 现有配装下对 DH 而言，其实除了三条技能词缀以外一般没有 A 伤了。150 级太极石能提供的增益是 80/145 = 55.17%。

相比之下，150 级困者是 60%，贼神是 16-80%。但是如果你技能伤小于 33.33%(比如箭袋没出+15技能伤)，那同等级太极是高于困者的，况且太极还有减伤。因此这是暴雪有意为之。

力士护腕： 其特效是对目标的一个持续 6 秒的减益效果，但属于 A 伤。因为是加在目标身上的减益效果，会对小队所有人生效，且多人的减益效果以加法叠加。这与一般的减益效果只取最高效果的不同。

如果玩家 1 对该怪施放了特效，那么所有人对怪的伤害都增加 30/145(假如所有人都有三条技能伤词缀)，如果四名玩家同时对该怪释放特效，则增伤 120/145。由于一人释放多人共享多人释放叠加共享的特性，它也只能分类到 A 伤，否则如果是乘法在多人游戏中就过于强大了(假设是乘法则增伤 1.3^4-1=185%)。因此这是暴雪有意为之。

深渊挖掘裤 与 杨先生的妖法裤：

前者是独立增伤，后者是 A 伤，但我找不到什么原因。我猜暴雪可能是认为妖法裤特效主要在于能量生成方面。在某些下能量生成多 25% 能否意味着伤害也能增加 25%？另一方面可能原因是，即使是 70 级引入了大量的乘法因子，但最初的数字也没有目前这样多次改版闭眼按 0 那么夸张，也没有克里森船长等铁匠套装填充六件套外的位置，妖法裤作为过渡装或者特定玩法是有一席之地的，哪怕其增伤属于 A 类。

某些职业被动： 如上所述，大多是 60 级遗留问题，且会在面板的全部技能伤位置显示，因此不再具体列举。

由 A 技能带出 B 技能时的增伤计算

跑马天拳、恐惧冰吞或者圣化扫射等构筑玩法中，反复释放的是 A 技能，但实际对敌人发动有效攻击的是 B 技能。这种情况下，你对 A 技能的增伤通常无法作用于 B 技能，你需要针对 B 技能配置增伤词条。典型误区之一是为恐惧冰吞添加扫射技能的增伤词条。但另一方面要注意的是，如果为 A 技能添加词条有助于更高效地触发 B 技能，那么该词条仍然是有价值的，例如在冲高层时为恐惧冰吞搭配维拉的遗赠。

减耗、减伤、减 CD

三减词条的独立计算规则

和伤害计算不一样，几乎所有的三减词条都是独立计算的。你的装备上有两条 +15% 技能伤词条，该技能的增伤是两个词条之和也就是 +30%。但如果你有两个减 8% CD 词条，你的总减 CD 并不是 16%，而是 $100\%-92\%\cdot92\%=15.36\%$。这个原因也很简单，你总不能凑两条 50% 减伤词条，就能拼出一个无敌效果吧。

$$总减量 =100\% - (1-三减词条1)\cdot(1-三减词条2)\cdot(1-三减词条3)\cdots$$

但其实，虽然三减词条的数值计算看起来相对复杂一些，但实际在装备选择上反而更简单，你完全不必考虑其它装备词条，不必考虑基底值。只要这个词条上写的减 10%，那它就是在你现有承受伤害的基础上减 10%，也就是增加了 11.11% 的生存时间。

另外，在详细信息面板上，你也可以看到已经计算完成的三减词条的最终计算结果，例如 CDR、减耗等，也不需要你真的逐一相乘。

护甲与抗性

两者独立地为人物提供减伤效果，其计算公式分别为：

$$护甲减伤 = 1-\frac{50 \times 敌人等级}{50 \times 敌人等级+护甲值}$$

$$对应属性抗性减伤 = 1-\frac{5 \times 敌人等级}{5 \times 敌人等级+对应属性抗性值}$$

• 通常来说，敌人是 70 级的。
• 护甲数值起到的效果正好是抗性的 10 倍。
• 这两者是完全独立的，也就是说，敌人对你的一次物理攻击，会被护甲减伤计算一次，然后会被物理抗性减伤再计算一次，如果还有其它减伤效果，再另行计算一次。
• 一次法术攻击，同样会被护甲减伤计算一次，然后被对应属性的抗性减伤计算一次，如果还有其它减伤效果，再另行计算一次。
• 要注意的是，『护甲和抗性提供了减伤』并不是『护甲和抗性就是减伤词条』，这两属性是先通过加法计算出一个总和值，再通过公式为人物计算一个减伤比例。这与三减词条逻辑不一样。
• 你可以通过智力主属性获得抗性，也可以直接通过抗性宝石、装备属性获得抗性值，护甲值也一样。同时还有其它一些对应的增加百分比的词条或效果。
• 对不同的职业，其计算公式是一样的，智力职业在护甲值上获得的减伤幅度和力敏职业完全相同，反之亦然，力敏职业在抗性上和智力职业也相同。

另外，在人物面板上也可以看到，每点力量敏捷主属性都会为人物提供 1 点护甲，每点智力主属性会为人物提供 0.1 点抗性。结合上述公式，可知每点主属性给到对应人物的减伤效果是完全一致的。

抗性减伤里的 A% + B%

有时我们会在抗性属性的鼠标提示里看到有如图所示的 A% + B% 两个减伤值。这个提示中间的 + 并不是把前后两个值加起来的意思（否则通常都会超过 100%）。而是说，这里有两条独立的减伤词条，前者 74.82% 是物理抗性对应的减伤，后者 61.08% 通常是 奥吉德、克里林船长 等套装提供的减伤。如果同时装备了两者，这里的减伤也已经是两者效果叠加计算后的值了。

仪式祭坛

一张网图：

三种喝药光圈与四种喝药圣坛效果

洋红色增伤光圈： 独立增伤 50%，持续 7 秒。增伤不如神目指环(70%-85%)，所以只有神目环位置不好，或者 boss 战无神目期间比较合适站红圈，但也存在神目与红环重叠的可能，此时能站就站，高达 277.5% 倍伤害。

草绿色降耗光圈： 降耗 50%，持续 7 秒。适合<span class="build>火多重、<span class="build>塔陨等对能耗要求高的构建。

天蓝色减 CD 光圈： 持续 7 秒。与普通 CDR 词缀不同，这个光环的运作方式是固定每 0.5 秒后减 1 秒所有技能 CD。 虽然长 CD 下基本等效于 -66% 的 CD 时间，但实际上会略小一些。例如对于一个 10 秒技能，在 6 个 0.5 秒后减去 9 秒 CD，再过 0.5 秒，理论上能再减1秒，实际也只有 0.5 秒可减。所以在光圈内需要 3.5 秒才能走完 CD。

光圈半径10码，并不一定刷在脚下，会在离人物 10-30 码的位置随机刷，甚至有可能刷在墙里面，但一定有边缘在墙内可供站立。

光圈的出现机率未实际测试，预期应为 1 : 1 : 1。

喝药只会出四种圣坛效果：飞驰 +25% 移速，狂暴 +25% 攻速，祝福 +25% 减伤，启迪 +25% 经验，但喝不出来强能圣坛与强盗圣坛的效果。四种情况出现机率预期应该是相同的，未实际测试。

仪式祭坛词条

仪式祭坛提供的各类效果词条和人物装备、巅峰等提供的效果都是一致的，伤害同样遵循同类相加异类相乘的规则，祭坛词条通常也会合并在人物面板详细信息一并显示。

传奇宝石

贼神的复仇之石

[ 图片来源 ] 由此可见主角至少有5米高，不愧是奈非天

容易判断的是屏幕左右边缘是 50 码。

我自己的测试结果（注意，图片中虽然在右边，但 0 线其实是在屏幕正中间的）：

结果显示， 贼神在近战范围就有最低增伤，而在 40 码外就满效果了。文字描述的『上限50码』，实际意思是『第50码的位置没有第六次伤害提档』，而不是『50码外才有满伤害』。

太极石

由于对 A 伤的嫌弃，现在大多数职业 build 中，一般都只有三条 15% 技能伤词缀是 A 伤了。所以尽管太极石也是 A 伤，但 150 级太极实际提供的增益是 80/145 = 55.17%。

相比之下，150 级困者是 60%，贼神是 16~80%。但是如果你技能伤小于 33.33%(比如箭袋没出多重)，那同等级太极是高于困者的，况且太极还有减伤。自行取舍。

受罚者之灾

上文已述，不再重复。

困者之灾 & 闪耀冰晶

很多冰系技能其实并不会为敌人添加 寒冷 效果，闪耀冰晶 就是这时候用的。但 寒冷 本身是带 60% 减速的（另外还减攻速），所以它是一种更强的减速效果，一般的减速效果通常只减移动速度。寒冷效果可以触发困者之灾，并不需要其它额外条件。困者之灾本身在 25 级时也会提供一个小范围的减速效果。

尽管如此，闪耀冰晶 并不会拿来和 困者之灾 一起用，因为有更多的办法达成减速而触发困者，而不是非得用 寒冷 效果。冰晶更多出现在多人组队时的辅助职业上。

转煞宝石

两条效果是独立的，第二条的减伤只在生命值 50% 以下时生效。要注意的是其抗性值增加 75% 并不是减伤直接增加 75%。按上面的公式，其有效减伤为：

$$\frac{350 + 你的抗性值\cdot 175\%}{350 + 你的抗性值}$$

如果你有 4000 抗，那其效果能接近 70%，如果你只有 1000 抗，那它效果为 55%。

闪电华冠

最强黄装搭配宝石。如果你刚到 70 级一身黄装开始打大密境，还没凑齐合适套装和对应宝石时，可以用它过渡一下。

职业：猎魔人

关于复仇全覆盖

前提：你已经有巅峰的 10% CDR(Cooldown Reduction)，头部钻石(无瑕的皇家钻石)12.5% CDR，以及萃取/满减黎明。

结论：你在装备上需要三条 CDR 词缀(所有技能的冷却时间缩短 x%)，可以从 肩、武器、箭袋、手套、指环、项链，魔女任意位置出，这三条只要加起来大于等于 21 就行。

无论是 10% 6% 5%(19.936秒)，或者 7% 7% 7%(19.953秒) 都可以满足复仇 CD 小于 20 秒。但 10% 5% 5%(20.149秒)或者 7% 7% 6%(20.168秒)就不行了。当然，更高的 CDR 仍然能带给你更余裕的操作时间，10% 8% 8% 能到 18.896 秒，8% 8% 8% 也能到 19.316 秒。

对，前文说了 CDR 每条词缀是独立计算不是直接加的。

但简化版结论更好记，从结果上也是正确的，尽管不同组合的实际结果不同，但都满足 <20 的要求。你可以认为只是巧合。严格来说是面板 CDR 数值 > 36.51%。

其实从计算结果来看，极限情况下，指环x2+肩出三条的“成本”反而是最低的，但条件苛刻掉率低。实际来看不如考虑船长套。

关于扫射攻速档位

背景： 攻速会加快无 CD 技能释放的频率从而提升伤害。但由于暗黑3本身的限制，一个攻击动作的动画必须在整数帧里播放完成，因此这些提升是分档位的。在同一档内，即使面板攻速数值上有差别，在战斗中实际技能释放频率却是一样的。因此攻速调整到适当档位，省出词缀提升其它伤害更合理。

对于猎魔人的扫射，还有一个额外的限制：由扫射带出来的技能，两次间隔之间最短不少于 9 帧。也就是说，如果你的扫射在第 1 帧和第 9 帧分别攻击了一次，不管是不是同一个敌人，第二次的攻击都不会带出另一个技能（比如机轮甲的冰吞箭，或者圣化赛季的多重、三刀），你需要等到第 17 帧的第三次攻击才会带出第二次伤害技。也就是说，随着攻速增加，你的有效伤害技能的施放速度从 60/9 帧掉到了 60/16，伤害反而下降了 43.75%。这意味着，理论上 9 帧间隔是最高频率，否则就应该考虑缩短间隔一轮的帧数。尽量把攻速堆到 6f 甚至 5f，而避开削减最大的 8f 7f 阶段。

前提： 巅峰有 10% 攻速。开了祭坛喝药随机圣坛 buff，随从带了圣坛手10分钟，认为喝药狂暴圣坛 25% 是接近常驻的 。低层割草药水随便喝，冲层本来就有赌图，进门喝药赌狂暴也很正常。

结论：

1. 对于弓(杨弓旅途终点)。装备+魔女攻速 >=19（通常为完美箭袋）。能 8f->7f，理论伤害提升+16.67%。如果无常驻药水狂暴，那么箭袋外，武器也要出攻速。

2. 对于单手弩(黎明要塞弩机)。方案一(建议)：武器出攻速，箭袋攻速 18 以上，不喝药 6f，喝药5f，喝药狂暴不要和加速塔叠加。方案二：武器不出攻速，箭袋攻速 18 以上。不点巅峰攻速，加速塔喝药狂暴可叠加。因为加速塔喝药狂暴叠加的情况毕竟相对少数，方案一还是高一档的，建议选方案一。

3. 对于家妮，武器完全不需要洗攻速，箭袋攻速也是有就行。不喝药7f，喝药狂暴6f，喝药狂暴+加速塔 5f。

4. 对于双手弩(暴雪重炮)。方案一，武器无攻速，箭袋+巅峰攻速<=26，药水狂暴不与加速塔叠加，但最好有一个。方案二：巅峰+10%。箭袋+武器攻速>=25，不带圣坛手，不喝药不开塔，必须喝药就赌不出，出了忍两分钟。理论上方案二是冰吞最高伤害。

职业：魔法师

关于奥陨能耗

陨石术 - 星辰契约 的伤害如其技能描述，是初始 40 点产生 740% 伤害 + 235% 地板伤害持续 3 秒，多余的每点奥能额外增加 20% 伤害。

由于地板伤害并不会叠加，在 1.4 攻速的情况下，连续扔陨石实际生效的只有 $740\%+235\%\cdot\frac{2}{9} = 795.95%$，很接近后续每点 20% 的伤害。所以对奥陨来说，如果没有其它因素，有能量就扔和满能量才扔的伤害是差不多的。但考虑到受罚者之灾和蕴火之心，尽量满攻速扔就变成了更好的选择。

进一步地，在 大维兹尔之伏 和 风暴护甲-风暴之力 以及其它减耗词条的帮助下，初始能量降到不足 17 奥能，而伤害是不变的，这使得初始奥能的伤害比更高，理论上最高伤害就变成了有能量就扔、能扔就扔，而不再是满能量才扔。

关于 A 伤

法师的技能 魔法武器 和 魔星-烈焰火花 以及被动 玻璃大炮 等都是 A 伤，会在面板上体现。冷血 不展现在面板上但也是 A 伤。前两个技能是 塔拉夏陨石 的构筑技能，套装要求必须四系法术，所以需要注意可以避开的也就是 玻璃大炮 和 冷血 了。

今天是 3.14，写写圆周率怎么算。

割圆术

一开始的思路很简单，既然是『圆周』的比值，那就寻找一个方法计算圆周的长度就完事了。

1. 首先准备一个圆内接正六边形，每边长是 r，周长是 6r，如果把六边形周长当作圆周长，圆周率就是 $\pi=\frac{6r}{2r} = 3$。别笑，历史上确实有以 3 当作圆周率的，『周三径一』曾在多处典籍出现。

   

2. 以正六边形的边作中垂线，构造一个正 12 边形，则正 12 边形的边长即 $x=L_{\overline{A_2A_8}}$。

   

   $\overline{A_2A_8}$是直角三角形 $\triangle A_2A_8C$ 的斜边，$\overline{A_2C} = \frac{r}{2}$，$\overline{A_8C} = r - \overline{OC}$。而$\overline{OC}$又是直角三角形$\triangle A_2OC$的直角边。因此应用两次勾股定理，容易得到：

   $$x = \sqrt{\left(\frac{r}{2}\right)^2+\left(r-\sqrt{r^2-(\frac{r}{2})^2}\right)^2}$$

3. 进一步在正 12 边形基础上继续作中垂线变为 24 边形、48 边形……类推。

   

   局部放大：

   

   在 n 次迭代操作后，可以得到圆的内接正 $6\cdot2^n$ 边形，其边长为可由$6\cdot2^{n-1}$ 边长递推而来：

   $$l_n = \sqrt{\left(\frac{l_{n-1}}{2}\right)^2+\left(1-\sqrt{1-(\frac{l_{n-1}}{2})^2}\right)^2}$$

   根号内的平方展开，公式化简后得到：$l_n=\sqrt{2-\sqrt{4-l^2_{n-1} } }$。所以：$\pi\approx3\cdot2^n\cdot l_n$

   尽管可以将 $\pi_n=6\cdot2^n\cdot l_n$，$\pi_{n-1}=6\cdot2^n\cdot\frac{l_{n-1}}{2}$ 代入得到 $\pi_n$ 的直接递推式，但结果反而会增加计算难度，没有意义。

   编程验证一下割圆法的递推公式：

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

   from decimal import Decimal, getcontext def print_pi_iterations_csv(n_iterations, precision=50): getcontext().prec = precision + 2 r,k,s = Decimal(1),Decimal(6),Decimal(1) pi_approx = (k * s) / (2 * r) print(f"0,{s:.{precision}f},{int(k)},{pi_approx:.{precision}f}") for n in range(1, n_iterations + 1): s = (Decimal(2) - (Decimal(4) - s**2).sqrt()).sqrt() k *= 2 pi_approx = (k * s) / (2 * r) print(f"{n},{s:.{precision}f},{int(k)},{pi_approx:.{precision}f}") print_pi_iterations_csv(n_iterations=100, precision=200)

   计算结果：

   n	$\pi$ （红色开始为误差项）	边数	边长
   0	3	6	1
   1	3.106	12	0.517638090205041524697797675248096656698137802639861027628006414630
   2	3.133	24	0.2610523844401031830968124557909780203874814096234645037812338966702656504596472
   3	3.139	48	0.1308062584602861336306311175503508828812646078901080356639753014926938148927779305724160345230399815
   4	3.14103	96	0.0654381656435522841273198526345762522299252545033923937561917675146374704265347854596049670897283379849162629786820274103417801265169
   5	3.14145	192	0.03272346325297356328594384696834610047132981567239244974814148723774665964804514057084743346984975274224223378392101609252042560186228338473297730652077981634239188734083635448788894004614630671707971
   6	3.141558	384	0.01636227920787425857039824658921801844533529416452341952057948786939579165064652866066255648348084605099544845686188429947694545964114075255218686278984440383734161800789305756187065766408033033714294
   7	3.141584	768	0.00818120805246957918924210834302383685031261696817992854723965236122383489144687206840513659491152519677943407736541665751720319595526819838604319161259091076464007027427841057546918033121536463776576
   8	3.1415905	1536	0.00409061258232819022882611784626426005671130508942161074227241896916750309821831955224474466566072674889796524473684900214896944666772227172928706453876775835266189348091823621540235396623950175141037
   9	3.14159211	3072	0.00204530736067660908238592229206387848510904769203534158371388824527333132200776714283048062007265192348301163404334423074239806419254304223189361723602246795475702449971333535955495781000908587365137
   10	3.14159252	6144	0.00102265381402739500247163869946902080272386682694276940575245689833248716463240037650358826374003409375522080937566141783691027096707879199337571242973168314629930390958384129075647765179770735929993
   11	3.141592619	12288	0.00051132692372483462812329904167729802119501834726478547125933748209551655293273596409563162021982192636730703620007115215205194042178721769005805718824184735509920188172391787369176696970438736673964
   12	3.141592645	24576	0.00025566346395130948052344901602532864686870304453467621730677964075264727555486053435606907311636775837244640506715770965493389536029465573114161194846506354894597613066543370721104363861111432890423
   13	3.1415926515	49152	0.00012783173223676626186947646649138260988152781056255781673371926384916560262546013047471577662181934711832084928074633138426333635453646621891256820530408464072964170436390615844992513167973181884419
   14	3.14159265306	98304	0.00006391586615102207116070807251821846928420566418226833293283142812113254271150192259083632846234511443992999462757190341371764220822115850525439907961017619817497180269701854339121002598869952981899
   15	3.14159265346	196608	0.00003195793307959090310938154255793333263503358446685194730153233192952836577298952906771466453949596139067707595400480370347309422927398029140060759615402672879750656729406381857083342148916810946482
   16	3.141592653556	393216	0.00001597896654030543499584362445350196006406104280692607674654698067095823491466166283932166777919474862554088504872494049081442239324689977187654126144557724034900570345338497245798640501861199132473
   17	3.1415926535814	786432	0.00000798948327021646542806668183879235631019681169477246348191864536472108682476590250406902878002067011109537606200591443120881102952774587693884967039607068647393920042977032508672304529703548793507
   18	3.1415926535877	1572864	0.00000399474163511620120530147346356461485004455676682492653080263097961031568786521018546178449577827905058028155793726100684082116620490208586760128562326052355254999492341178495600442520229714214498
   19	3.14159265358927	3145728	0.00000199737081755909666405925404488658903363571572547377195911130767842215943514172624435772225406061136635761994565450768012986355225672804758813211799003989963368820258474934638028323690684853391208
   20	3.14159265358966	6291456	0.00000099868540877967283970569170986518056233278254800139575240199441413157386950268773547426918651144098773169103279081732332903382750784965462242767295874351578603227570190990894498209480288483175579
   21	3.141592653589761	12582912	0.00000049934270438985198331235394158794637574452005163238561622949160536438778783833013780076841949385532244297830798431363949646238115306391824299463536681672767037675392104288262525145523660499757185
   22	3.141592653589785	25165824	0.00000024967135219492793708861548164863083528480050429414620649761964010850031010937136302482094374335613839348636992034690138588360189467559359010796462784072115098350334154845172067404838766858926926
   23	3.1415926535897912	50331648	0.00000012483567609746421172336255468185819030616389052416717182657674416467433744424882377430722467888793107943167913990376778051086967167245591527296295705271015718603906188572436811955215249723929612
   24	3.1415926535897927	100663296	0.00000006241783804873213625906312907314414694658977754689172890296684927941714902639547605794720789562158704287228431191826831833465509516927936598151213605097619366385055651234305254691020971269767748
   25	3.14159265358979311	201326592	0.00000003120891902436607192920429600309964886031266085579006387133906218733474982304042661535174860729942309322597890450240486561656308496161043473084985788277410293815655336611240206071283718092142279
   26	3.141592653589793207	402653184	0.00000001560445951218303643956123943486579303830946734590076137071151382758546729017465867897056956871062577565687205469967325499101603554656770693413534572719717163479597830054400212066111172965038686
   27	3.1415926535897932305	805306368	0.00000000780222975609151827915050614659739327282312314419213435000540384209982803362747958675796558062042982475357489548994098448825732511259496703453493157075790296236194140478050246191004140967384875
   28	3.1415926535897932365	1610612736	0.00000000390111487804575914699648887694425875179664923589163142535332766643652227643640373038294224713031271130104996335321398963832881871591336387327721771031647629833578335173087223319866660377209223
   29	3.1415926535897932380	3221225472	0.00000000195055743902287957442589891392782464098322741904183452748558916604225809692973822374008478147095680240506129083134376023950655446505210288531412793791591047020932205873679944497150982204577829
   30	3.14159265358979323834	6442450944	0.00000000097527871951143978732890626639587422864790677350546878071113342307743900198686359140290355881550128722580828503965413131817775505376202266065797748836863198053672189794374120960822078454085167
   31	3.141592653589793238432	12884901888	0.00000000048763935975571989367894773437693235284413627643353924138657991875227166863658134592478609629063235498679155112391032958494770696135492336364934405252910218154894519718632749329296292232873980
   32	3.141592653589793238455	25769803776	0.00000000024381967987785994684128569233584058123711119494820572430372864962831541128412896710085186107246109239553239175349774989214200004231103228126374017973510217647145205395451425934862543903141874
   33	3.1415926535897932384607	51539607552	0.00000000012190983993892997342086932431134209122043661067557410939165547834271901397817958184287679899746954872366931517378432756277219608864154084676071873439013814578810008899659458297578644652102792
   34	3.14159265358979323846216	103079215104	0.00000000006095491996946498671046297192359877068545345171015976479731683080387882390963230448827617887470480541069715613284937281391527142769874265021378174567842463712478657447014329047907537025078812
   35	3.14159265358979323846252	206158430208	0.00000000003047745998473249335523502468279035097713111976794487129248565177626429377675630894324491318143281681108283621034037740658494844040637391634904185368756479826080968237471660797772242533971580
   36	3.141592653589793238462613	412316860416	0.00000000001523872999236624667761795468151904619286610914234050926206938980743255828792757021971394828021509269854954204203807365348514101443729223110438242704353620510969922222720558589696914162834932
   37	3.141592653589793238462636	824633720832	0.00000000000761936499618312333880903263327500693447062322906813727064108287401977495862480364165191892572722821805378076342708109786024478510267128305685288129822216150408530005053668788474754596314193
   38	3.1415926535897932384626415	1649267441664	0.00000000000380968249809156166940452322820193894699000769679011251019721704206004157069165212410810544026360243081562173175012380943710784819833158848148454305274178869725472055344941764738359510503432
   39	3.1415926535897932384626429	3298534883328	0.00000000000190484124904578083470226247804652390846434085867764950648712663131010933239756397504390127793843615009932381929180666667653837973866688861742752372542199845768295329289064647703956013271496
   40	3.14159265358979323846264326	6597069766656	0.00000000000095242062452289041735113134701645625860333755562416727738678225630409133726241388412271419693280728455731432240919779384123153580275368398924018623402760598689881472309873103129157256079637
   41	3.141592653589793238462643354	13194139533312	0.00000000000047621031226144520867556568700737741734806466859775202820253268876017605259802738752786993037186448547714104069276714403847492681351960293481403671057123207625030082320476025045525187921347
   42	3.141592653589793238462643376	26388279066624	0.00000000000023810515613072260433778284519108236967983182064708458072713526423819809929748745762233445722522915255220218907083199050851571707542444274694226267953634788865545023196916578383333271638065
   43	3.1415926535897932384626433814	52776558133248	0.00000000000011905257806536130216889142280646539246564084611706836175233956662580323924238486502823822095273623016740503242823293473322752826943406712212347816633161124743887015182230112651127795772707
   44	3.1415926535897932384626433828	105553116266496	0.00000000000005952628903268065108444571142959822218603604003272493981728309761137840677254857193124035843137902790456916465278946093378200973737232257810799457850600559768648334606507173483923226798148
   45	3.14159265358979323846264338316	211106232532992	0.00000000000002976314451634032554222285571809480183716997213813631477682811379565938813807352206137928266371922361967513628494303350752813837220533755077556277267010508481067708258272280581697116368445
   46	3.141592653589793238462643383251	422212465065984	0.00000000000001488157225817016277111142785945936226160398008428988799695448379353411048669540675105529201156255393239241815105795960193162958466509541845372039556076088358306744001667409563608755116413
   47	3.1415926535897932384626433832723	844424930131968	0.00000000000000744078612908508138555571392978117629867936429404766032454532982972879927442015100333573102539703125512303180592336739443426803028897653497323464198804132243445293300830649524734564338552
   48	3.1415926535897932384626433832777	1688849860263936	0.00000000000000372039306454254069277785696489702504532435392851166970303119261182832676543366635866770056075491914594374169492690139868978298367458617935594850996354453930203928999116673294818087014542
   49	3.14159265358979323846264338327905	3377699720527872	0.00000000000000186019653227127034638892848244931713466026093694181479411041279007664518388039250219394917520581289134362666072356238423046393393267507079270873999818256870264565851827569558246606719889
   50	3.14159265358979323846264338327939	6755399441055744	0.00000000000000093009826613563517319446424122475914382989096505665488987955847187244503302737899350870943426865862471134189914127942916634993896055484923272330147659367044182796913164438902451505300070
   51	3.14159265358979323846264338327947	13510798882111488	0.00000000000000046504913306781758659723212061239214397741554460154588154282324605029445339776664665126602561797118126889216645671298658718824327096087054775868599403253919720208550301571372348300608040
   52	3.14159265358979323846264338327950	27021597764222976	0.00000000000000023252456653390879329861606030619764349651653005992524534679212430533767605152357040947956551476686121907768168532165192866105014965558860008759203488559386437922131721340234005270487190
   53	3.1415926535897932384626433832795011	54043195528445952	0.00000000000000011626228326695439664930803015309901818673435974985666074531862481319050223365839744480599681889753061151021039333721394536374942994115768269838484630595034485522996445808917007449174018
   54	3.14159265358979323846264338327950244	108086391056891904	0.00000000000000005813114163347719832465401507654953364817669171491508513164963273917601720652929341974206876005338536255726570022549748804281781239796629797911615377845161799814839230123476900245597338
   55	3.14159265358979323846264338327950277	216172782113783808	0.00000000000000002906557081673859916232700753827476989343953483745588691069860641116109055396819036476750547363087511217529594767943500465379157613992458109310350412218037800653716796508005520240576783
   56	3.141592653589793238462643383279502857	432345564227567616	0.00000000000000001453278540836929958116350376913738533038866604122773649845852696077719571424363288354984560163404176381957147290806540633568741458332329570497135696625855744102793092493269403397306521
   57	3.141592653589793238462643383279502877	864691128455135232	0.00000000000000000726639270418464979058175188456869271315294534842634237961791644978817963657368349018018547073553645174721809397062791589612692380420980430418549727383863300796950010128713119835412810
   58	3.1415926535897932384626433832795028825	1729382256910270464	0.00000000000000000363319635209232489529087594228434636257129921518973045610753984606903755555565090225823497728933861097425510012445275500660856570420314377959077737522021138499361073229058057856525305
   59	3.1415926535897932384626433832795028838	3458764513820540928	0.00000000000000000181659817604616244764543797114217318203500292521693513634109262568138724540009302627880665664054547671442602629245985300233979583099822160020500500483475659543003887476852580099590637
   60	3.14159265358979323846264338327950288409	6917529027641081856	0.00000000000000000090829908802308122382271898557108659111117062731122630670646165067155218116731953598635420518712225579467900141661189798523079951316477341015481355037440868192797508134947832305112068
   61	3.141592653589793238462643383279502884170	13835058055282163712	0.00000000000000000045414954401154061191135949278554329556729395924345799567022024256463341039252169509758470379152247370308692670281380402226420303760039400738209142733170889955014305024957720733537028
   62	3.1415926535897932384626433832795028841904	27670116110564327424	0.00000000000000000022707477200577030595567974639277164778511056032020960312473379843592387017238273841336013772018156620620139714636784838886552346958522729508544287095924157656559328774743942058403501
   63	3.1415926535897932384626433832795028841955	55340232221128654848	0.00000000000000000011353738600288515297783987319638582389273822774741487722356985886216283070820704775151719320297692305771229984585552789478670272471592185318112305923135816186725549723401189399465139
   64	3.1415926535897932384626433832795028841967	110680464442257309696	0.00000000000000000005676869300144257648891993659819291194639198232212119806943529938660652730685549751219967936668777898522258006789613584548518548947517399518980753460579735453010148010561666852518465
   65	3.14159265358979323846264338327950288419706	221360928884514619392	0.00000000000000000002838434650072128824445996829909645597319884971711231896692394593774390264752174589247798884844188386097084317662800257999539949531820854021372973920212589123057484043038413386480875
   66	3.141592653589793238462643383279502884197143	442721857769029238784	0.00000000000000000001419217325036064412222998414954822798659978217806262447498775999942703119802262260178073225171594371672767039910395837746754883514554387355534428542490502843047419903699298900650460
   67	3.1415926535897932384626433832795028841971628	885443715538058477568	0.00000000000000000000709608662518032206111499207477411399329993575396962036143460337853290058329403000825478551622790184696278207177431648889605252409000436176691963484486508085164556906087635528463600
   68	3.1415926535897932384626433832795028841971678	1770887431076116955136	0.00000000000000000000354804331259016103055749603738705699664997346010209869620989211161887341468235484256112337447060013347246770025475017772545318040100499719527221982164346924474197203014536962384233
   69	3.1415926535897932384626433832795028841971690	3541774862152233910272	0.00000000000000000000177402165629508051527874801869352849832498742794071041254151985860348959772059490108518983279739396677937494271519282335856564782606593038161921464824045057141901255243934791698939
   70	3.14159265358979323846264338327950288419716930	7083549724304467820544	0.00000000000000000000088701082814754025763937400934676424916249380120656283932533165465100141015772463551818630392263099429490062385425702670489272777099615392845908043683470262540481687443753631387209
   71	3.141592653589793238462643383279502884197169374	14167099448608935641088	0.00000000000000000000044350541407377012881968700467338212458124691150780737379448729299415778149104071588104247756832825705434634711181057142482354004448463241253905615654192091803661550138349834925883
   72	3.1415926535897932384626433832795028841971693929	28334198897217871282176	0.00000000000000000000022175270703688506440984350233669106229062345711696943116372132970566102529704265770576491705274450503252719017657177675877691511013902854190096207806219192633888747260021845335379
   73	3.1415926535897932384626433832795028841971693978	56668397794435742564352	0.00000000000000000000011087635351844253220492175116834553114531172872886793361517037525390327946746161632353791870268554275183884754782486331603778813068233615174199226884022215799223710907726838045600
   74	3.1415926535897932384626433832795028841971693990	113336795588871485128704	0.00000000000000000000005543817675922126610246087558417276557265586438573186906174890142708573558609834409560089188565508087950401784329091801122048811978107497632868024621903259103783240092289033277638
   75	3.14159265358979323846264338327950288419716939927	226673591177742970257408	0.00000000000000000000002771908837961063305123043779208638278632793219552817231264491493855962977459511403952943751000011500970438591506577337745170668333578535599951979023667877352456919294462872091966
   76	3.141592653589793238462643383279502884197169399350	453347182355485940514816	0.00000000000000000000001385954418980531652561521889604319139316396609809686587904376549740691013499079976873084270090861482240887454217130109700888814980605677101995723614993166563461344604800739186474
   77	3.141592653589793238462643383279502884197169399369	906694364710971881029632	0.00000000000000000000000692977209490265826280760944802159569658198304909003040486204625221934197345705522798618684369325162265879229605485469678366255030178352556066002693178577047232466989525115469101
   78	3.1415926535897932384626433832795028841971693993735	1813388729421943762059264	0.00000000000000000000000346488604745132913140380472401079784829099152455021488559854356404915684997373453194568910850150554234743333323887169027810355721265754992488889043086125818785113193573919393452
   79	3.1415926535897932384626433832795028841971693993747	3626777458843887524118528	0.00000000000000000000000173244302372566456570190236200539892414549576227575740319521183676701415789251813071691901508261310331929140249243806985456612550895329193795879589953572603154100824544871464065
   80	3.14159265358979323846264338327950288419716939937501	7253554917687775048237056	0.00000000000000000000000086622151186283228285095118100269946207274788113795994664709842522631154555946542345146881514528910460810254597164342807858892835620795414037008213423109321081194546364822645740
   81	3.14159265358979323846264338327950288419716939937509	14507109835375550096474112	0.00000000000000000000000043311075593141614142547559050134973103637394056899012895473577596850633110638350648736057102314237177980400366216536052889102776078611695252923253860326700438854434523571599969
   82	3.14159265358979323846264338327950288419716939937510	29014219670751100192948224	0.00000000000000000000000021655537796570807071273779525067486551818697028449633393126620840367198534402310258888355594288341333553345644955268599269768840268037991033640132662622651219876868571667511013
   83	3.1415926535897932384626433832795028841971693993751043	58028439341502200385896448	0.00000000000000000000000010827768898285403535636889762533743275909348514224832564737039425426334514586546996259218677535573509881948390470781402228825989395273182032651008058458132368741640505896617496
   84	3.14159265358979323846264338327950288419716939937510544	116056878683004400771792896	0.00000000000000000000000005413884449142701767818444881266871637954674257112418265890235838368509163216447481481489448816712110330223715773980089907117859221939072724239797636078782671271242814072375971
   85	3.14159265358979323846264338327950288419716939937510573	232113757366008801543585792	0.00000000000000000000000002706942224571350883909222440633435818977337128556209380885332434891172319848620488659729738164471724588802112783671406082911480437965024962904042577987603085039013447608945929
   86	3.14159265358979323846264338327950288419716939937510580	464227514732017603087171584	0.00000000000000000000000001353471112285675441954611220316717909488668564278104721435193031908950877204359837819737995801750320972363402779838300917945572908652501006153003542085374325991604412869361661
   87	3.141592653589793238462643383279502884197169399375105815	928455029464035206174343168	0.00000000000000000000000000676735556142837720977305610158358954744334282139052364591662367762396028262186118096103138740814467820927027954854394810262789410419235082714217186525080075111982148174778001
   88	3.141592653589793238462643383279502884197169399375105819	1856910058928070412348686336	0.00000000000000000000000000338367778071418860488652805079179477372167141069526182780089415357188087838593833946330836975399647327306680837620194187065419258155122989250888415010712734684797608033073625
   89	3.14159265358979323846264338327950288419716939937510582060	3713820117856140824697372672	0.00000000000000000000000000169183889035709430244326402539589738686083570534763091450576986613092803132734513835450326938323875340758736308820974542462674391428058499844637478820265381488838988104029032
   90	3.14159265358979323846264338327950288419716939937510582088	7427640235712281649394745344	0.00000000000000000000000000084591944517854715122163201269794869343041785267381545732855028173358746468046956525510777025489944130017542641677057978947214406421351026331982610650978815836832879531168973
   91	3.141592653589793238462643383279502884197169399375105820952	14855280471424563298789490688	0.00000000000000000000000000042295972258927357561081600634897434671520892633690772867373330945030916346733440713728590207285972872463543131778575672519330342530762911803363811511058158186027795432290378
   92	3.141592653589793238462643383279502884197169399375105820969	29710560942849126597578981376	0.00000000000000000000000000021147986129463678780540800317448717335760446316845386433804892579809401062455465663235945315460611537163618042257785088658544266429819636249671110253420808746597446520573378
   93	3.1415926535897932384626433832795028841971693993751058209735	59421121885698253195157962752	0.00000000000000000000000000010573993064731839390270400158724358667880223158422693216917224678316443392363825994914428934207508906198289830674985682660949635821384192403281631922844526642413278832413310
   94	3.14159265358979323846264338327950288419716939937510582097458	118842243771396506390315925504	0.00000000000000000000000000005286996532365919695135200079362179333940111579211346608460459637709689553823924642869271501663404845301205016530755451802616306227467630540663683104068369613595989640792411
   95	3.14159265358979323846264338327950288419716939937510582097485	237684487542793012780631851008	0.00000000000000000000000000002643498266182959847567600039681089666970055789605673304230460731173778259117213777111142880151658721675860020914535582465972145418673451986278984778163602362579726077072294
   96	3.141592653589793238462643383279502884197169399375105820974922	475368975085586025561263702016	0.00000000000000000000000000001321749133091479923783800019840544833485027894802836652115259229626755814834263320515134831240823898216087199538412524249058018818246139798994424707318178321345945256457982
   97	3.141592653589793238462643383279502884197169399375105820974939	950737950171172051122527404032	0.00000000000000000000000000000660874566545739961891900009920272416742513947401418326057633222818361243076588714252512839516036266280313248404349353781084532268698995350740374253514900417876053020966349
   98	3.1415926535897932384626433832795028841971693993751058209749432	1901475900342344102245054808064	0.00000000000000000000000000000330437283272869980945950004960136208371256973700709163028817062409803538495726488875624597744971172786690330281567563348535035541670308629447817987552336482887802186116224
   99	3.14159265358979323846264338327950288419716939937510582097494423	3802951800684688204490109616128	0.00000000000000000000000000000165218641636434990472975002480068104185628486850354581514408587579979633867542260906483321120854716349161878400707892481545467979558812567504614518139663047436054484489897
   100	3.141592653589793238462643383279502884197169399375105820974944503	7605903601369376408980219232256	0.00000000000000000000000000000082609320818217495236487501240034052092814243425177290757204300836874550011231007511825538341473499419058028357844460091683379454395130967959953704622673261060242184507762

4. 在圆的外切正六边形上不断向内收缩也可以用来计算 $\pi$，与内接多边形的逻辑基本一致，但值是从大往小推的，和内接多边形共同构成了对 $\pi$ 的双向逼近。

   • 外接六边形的边长由 ${L_0}^2 = {\frac{L_0}{2}}^2+r^2$ 得到，即 ${L_0} = \frac{2}{\sqrt{3}}$。
   • 边长递推公式由半角公式推得：
     $$t_n=\frac{2t_{n-1}}{1+\sqrt{1+t^2_n}}$$
     两者相夹：
     $$ 3\cdot2^n\cdot l_n \lt \pi \lt 3\cdot2^n\cdot t_n$$

无穷级数

割圆是个效率很低的算法。当然，最初除了割圆别无他法，自然无所谓效率问题。直到有一天，人们发现了无穷级数：

$$\frac{\pi}{4} = 1 -\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots = \sum_{n=1}^\infty\frac{(-1)^{n-1}}{2n-1}$$

这个公式称为莱布尼茨级数，最早由印度数学家马德哈瓦发现并应用。

※ 你没看错，真正的先驱往往成为无名的先烈。这个公式据称是莱布尼茨独立重新发现的。

严格来说，莱布尼茨级数的计算效率并不高，甚至比割圆法还低，每增加一位精度需要多计算 10 项。但公式本身启迪了人们，计算 $\pi$ 还有割圆术不断乘方开方以外的其它方法。同时，级数更方便并行计算，不受递推公式约束，每一项的计算也更简单，因此具有重大的历史意义。随后越来越多的 $\pi$ 级数被研究出来：

• 马德哈瓦级数，每项精度约 1.4 位：

  $$\pi = \sqrt{12}\sum_{k=0}^\infty\frac{(-1)^k}{3^k(2k+1)}$$

  最早实用于级数法计算 $\pi$ 的公式。

• 欧拉级数：

  $$\frac{\pi^2}{6} = \sum^{\infty}_{k=1}\frac{1}{k^2}$$

  任何一本高数书里都会有的公式。

• 牛顿级数：

  $$\pi=24\left(\frac{\sqrt{3}}{8}+\sum^{\infty}_{k=0}\frac{(2k)!}{(k!)^2(2k+1)16^{k+1}}\right)$$

  牛爵爷是大英帝国的雄狮。

• 反正弦函数泰勒展开

  $$\frac{\pi}{2} = 1 + \frac{1}{3} + \frac{1}{3}\frac{2}{5}+\frac{1}{3}\frac{2}{5}\frac{3}{7}+\cdots = \sum_{k=0}^{\infty}\frac{2^k(k!)^2}{(2k+1)!}$$

• Machin 公式：

  \begin{equation*}
  \begin{cases}
  \frac{\pi}{4}=4arctan\left(\frac{1}{5}\right)-arctan\left(\frac{1}{239}\right) \\
  arctan(x)=\sum^{\infty}_{k=0}\frac{(-1)^kx^{2k+1}}{2k+1}
  \end{cases}
  \end{equation*}

• 高斯-勒让德算法：

  \begin{equation*}
  \begin{cases}
  \pi \approx \frac{(a_n+b_n)^2}{4(1-\sum_{k=0}^{\infty}2^k(a_k-b_k)^2)} \\
  a_{n+1}=\frac{a_n+b_n}{2} \\
  b_{n+1}=\sqrt{a_nb_n}
  \end{cases}
  \end{equation*}

  计算机出现前夜的成果。

• 拉马努金级数：

  $$\frac{1}{\pi}=\frac{2\sqrt{2}}{99^2}\sum^{\infty}_{k=0}\frac{(4k)!(1103+26390k)}{(k!)^4\cdot396^{4k}}$$

  来自娜玛卡尔女神梦中的启示。

• Chudnovsky 级数：
  $$\frac{1}{\pi}=\frac{12}{640320^3}\sum^{\infty}_{k=0}\frac{(-1)^k(6k)!}{(k!)^3(3k)!}\cdot\frac{13591409k+545140134}{640320^{3k}}$$

  目前最高效的级数，每项可以推进 14 位精度。

• BBP 级数：

  $$\pi=\sum^{\infty}_{k=0}\frac{1}{16^k}\left(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}\right)$$

  这个级数具有革命性的意义：可以直接计算圆周率 16 进制（或 2 进制）下任意一位的值，而不需要从头开始逐项推进。因此尤其适合于分布式。

随着计算机的应用，单项是否易计算不出错已经不再是问题，计算效率变为了各类公式中最重要的考核标准。目前应用得最多的就是 Chudnovsky 级数和 BBP 级数。

一些豆知识

• 祖冲之与刘徽

  虽然课本里多提到祖冲之，但实际上割圆法是魏晋时期刘徽提出的。祖冲之据推测应该是在计算方法上有技术创新，同时明确提出了疏率 22/7 和密率 355/113。

• 东西方割圆

  前文演示的是通过内接圆和外切圆逼近 $\pi$ 的方法，但实际上刘徽的割圆术是使用内接多边形的面积来计算的，因为他先证明了 $圆面积 S = \frac{圆周 R}{2}\times半径 r$，面积计算有助于减少迭代时的计算量。纯靠线长计算是阿基米德干的。

• 割圆术/法的效率：

  $$k\approx2nlog_{10}2+log_{10}36\approx0.6n+1.56$$

• 椭圆周长公式为：

  $$C=4a\int_{0}^{\pi/2}\sqrt{1-{(1-\frac{b^2}{a^2})}sin^2\theta}\ d\theta$$

  $C$ 与其长轴 $a$ 短轴 $b$ 有关，但并不存在固定的比值。该积分没有简化形式，需要通过数值方法计算。但拉马努金提出了两个近似公式可以简化计算：

  $$(1), C \approx \pi\left[3(a+b)-\sqrt{(3a+b)(a+3b)}\right]$$

  $$(2), C \approx \pi(a+b)\left[1+\frac{3h}{10+\sqrt{4-3h}}\right], h = \frac{(a-b)^2}{(a+b)^2}$$

  前者最大误差 0.05%，后者最大误差 0.01%。

• 椭圆面积公式为：

  $$S=\pi ab$$

  面积与周长不同，面积是二维度量，在几何变换下的缩放是线性的，对形状的局部变形不敏感，因此公式就是如此简单直接。

  而周长是一维弧长的积分，依赖曲线每一点的曲率变化，因此需要专门的椭圆积分计算。

142857

近期刷到了一个视频，说 142857 这串数字如何的神秘，与古埃及有如何如何的关系。由于它分别乘以 1~6 都还是这串数，所以叫什么走马灯数，如何如何。抛开神秘学话题，这就是简单的 1/7 循环节：

不过『走马灯数』这个概念还是有点意思的。它最显著的特征是循环节长度 = n-1，同时附带的性质是当乘以 $1, 2, 3, \cdots, n-1$ 时，循环节顺序和长度不变，只是起始位置换了一下。因此写了段程序简单搜索一下还有没有其它的 n 也有这个特征。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

const loop = (n) => {
    let a = 1;
    let pool = new Set();

    while (true) {
        a = (a % n) * 10;
        if (a === 0) {
            return "";
        } // 除尽没有循环节
        let digi = (a / n) | 0;
        let b = a + digi;
        if (pool.has(b)) {
            // 当当前位的值和当前余数都出现过时，新的循环开始了
            let loop = [...pool];
            let index = loop.findIndex((i) => i === b);
            let str = loop
                .slice(index)
                .map((e) => e % 10)
                .join("");
            return str;
        } else {
            pool.add(b);
        }
    }
};

let start = 7;
while (start < 1000) {
    let str = loop(start);
    if (str.length == start - 1) {
        console.log(start, str);
    }
    start++;
}

结果出乎意料地多：

简化代码并扩大计算量搜索后，得到 100 万以内的走马灯数分布：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

const n_devide_len = (n) => {
    let a = 1,
        pool = new Set();
    while (true) {
        a = (a % n) * 10;
        if (a == 0) {
            return 0;
        }
        if (pool.has(a)) {
            return pool.size;
        } else {
            pool.add(a);
        }
    }
};

有点意思。

研究原因

对于任意由 1/n 形成的循环小数，$\frac{1}{n} = 0.\overline{******}$，设循环节部分为 $r$，其长度为 $l$，都可以写成一个无穷级数之和：

$$\frac{1}{n} = \frac{r}{10^l}+\frac{r}{10^{2l}}+\cdots+\frac{r}{10^{kl}} +\cdots = \sum_{k=1}^\infty\frac{r}{10^{kl}}$$

对这个无穷等比数列重新求和，$S=1/n$，$10^l * S-S = r$，得到 $S = \frac{r}{10^l-1} = \frac{1}{n}$，所以 $n=\frac{10^l - 1}{r}$ ，或者写成 $r=\frac{\overbrace{999\cdots999}^l}{n}$，这里的 r 一定是整数。这个方法同时也是将循环小数化为分数的标准方法。

• 当 $n = 7, l = 6$ 时，也就得到了 $142857 = \frac{\overbrace{999999}^6}{7}$，所以 $1/7 = 0.\overline{142857}$

• 当 $n = 3, l = 2$ 时，得到 $r = \frac{99}{3} = 33$，所以 $1/3 = 0.\overline{33}$，逻辑上其实也是对的，只是一般来说小数循环时只考虑最小周期，所以正确的写法应该是 $1/3 = 0.\overline{3}$。

• 当 $n = 11, l= 10$ 时，得到 $r = \frac{9999999999}{11} = 0909090909$，要注意这是循环节的计算因此首位的 0 不能去掉，所以 $1/11 = 0.\overline{0909090909}$，正写为 $1/11=0.\overline{09}$

对于 $1/6 = 0.1\overline{6}$ 这类存在非循环小数部分的，可以通过乘 10 等操作，把非循环部分移到整数上并抛弃，即化为 $10/6 \rightarrow 4/6 = 2/3 = 0.\overline{6}$，考察 $0.\overline{6}$ 即可。

同余

在上述 $n=3, n=11$ 两个例子中，虽然得到了逻辑上可行的循环节，但不符合通常的习惯，于在结果上进行了缩减，变为符合习惯的最小循环节。

• $10^k-1$ 本身就是 3 的倍数，所以循环节永远只有一位。

• 11 因其自身的特殊性质，11 = 10 + 1，天然地符合 10 为底的级数的结构，如果一个数字各位数值相同，那么奇数位数字对 11 都是同余的，偶数位也是。所以 $10^k-1 (mod\ 11)$ 只有两个结果 9 和 0。并且在两者间不断跳转。

可见，循环节长度即是最早出现同余时的 k 值，若要 1/n 的循环节长度为 n-1，需要 $10^k \ mod\ n$ 在 $k = 1, \cdots, n-1$ 时都不同余。

进制

从上可见，循环小数与进制密切相关。例如在 8 进制中，9（在 8 进制中也写为 11）的性质将和 10 进制中的 11 相似，也会出现奇位同余偶位同余的情况。所以 $1/11_{(8)} = 0.\overline{07}_{(8)}$，形式上完全一致。

而 $1/7_{(8)}$ 则变为了 $0.\overline{1}_{(8)}$，因为 $8^k-1_{(8)} = {\overbrace{77\dots7}^{k}}_{(8)}$ 各位恒为 7，也就不存在什么神秘数字循环了。

由于和进制本身高度相关，因此某个 1/n 是否有 n-1 位循环节，取决于它和进制 $R^k - 1$ 之间的相除结果，这个关系不知道数论中是否已有结论。

推论

1. n-1 位循环节天然存在，只是因为有可能内部重复导致长度缩减。而在 n-1 基础上的缩减不会破坏原有模式，因此最小循环节长度一定是 n-1 的因数。

2. 只有当 n 为素数时，才有可能实际出现 n-1 位循环节。由于余数需要遍历 1 到 n-1 的所有数。若 n = pq，余数会在 q-1 与 p-1 的共同作用下循环，但 p-1 与 q-1 至少含有公因素 2，不可能产生所有 {0, p-1} 与 {0, q-1} 之间的全部两两组合。

3. 一切循环节在若干次乘 10 排除不循环部分，整理为 $0.\overline{******}$ 形式后，都可以用 $a * 1/n$ 来表示。

4. 如果 n 为素数，而其循环节长度 $l$ 不为 n-1，则必有互相不能抵达的 $\frac{n-1}{l}$ 套循环，进入哪套循环取决于分子 $a$ 的值。因为素数没有因数，$a/n (a=1, 2,\dots, n-1)$ 必然有 n-1 个不同余数。

5. 如果 n 为素数，则 1/n 的循环节必然形如 $0.\overline{******}$，从小数点后就开始，即使起始是若干个 0 也是循环节的一部分。

6. 如果 n 是若干个不同素数乘积，则循环节长度是这些素数对应的循环节长度的最小公倍数。如果 n 是素数 p 的 k 次方时，循环节长度是 $l * p^{k-1},(p>3)$ 或 $l * p^{k-2},(p=3)$。

   $p=3$ 的特例在于 $p^2=9$ 仍然小于 10，没有触发进退位，即当 p^k 仍小于对应进制单位时，不往上叠加循环节长度。在 8 进制下就是正常的 $l*p^{k-1}$ 了。

7. 不存在比 n-1 更长的循环节。

8. 不同进制下的走马灯数具体是不一样的，例如 7 在 10 进制下是走马灯数，但在 8 进制下不是。2 在一切奇数进制下都是走马灯数，因为只要除不尽，循环节长度就至少为 1。

9. 除 2 以外，走马灯数的循环节长必度是偶数，将循环节的两段拆分后相加，一定是 $999\cdots$。例如 $1/7 = 0.\overline{142857}$，则 $142 + 857 = 999$。上述结论在其它进制下也成立，例如 8 进制下为 $777\cdots_{(8)}$。

   证明：设一个循环小数 $\frac{1}{n}$ 的循环节两段分别为 $\overline{r_1r_2}$，总长度为$2k$，则其既可以表达为

   \begin{align*} \frac{1}{n} = \frac{\overline{r_1r_2}}{10^{2k}} + \frac{\overline{r_1r_2}}{10^{4k}} + \cdots + \frac{\overline{r_1r_2}}{10^{pk}}+\cdots = \frac{\overline{r_1r_2}}{10^{2k}-1} &&(1) \end{align*}

   也可以表示为：

   \begin{align*} \frac{1}{n} =\frac{r_1}{10^k} + \frac{\overline{r_2r_1}}{10^{3k}} + \cdots + \frac{\overline{r_2r_1}}{10^{(p+1)k}}+\cdots &&(2) \end{align*}

   $(2)$ 式乘以 $10^n$ 再加 $(1)$ 式，得到：

   \begin{align*} \frac{(10^k+1)}{n} = r_1 + \frac{(r_1+r_2)(10^k+1)}{10^{2k}-1} = r_1 +\frac{r_1+r_2}{10^k-1} &&(3) \end{align*}

   再将 $(1)$ 式分子分母同除 $10^n+1$ 得到：

   \begin{align*}
   &&\frac{1}{n} = \frac{\overline{r_1r_2}}{10^{2k}-1} = \frac{\overline{r_1r_2}/(10^k+1)}{(10^{2k}-1)(10^k+1)} = \frac{\frac{\overline{r_1r_2}}{10^k+1}}{10^k-1} &&(4) \\
   \Rightarrow &&\frac{\overline{r_1r_2} * n}{10^k+1} = 10^k-1 &&(5)
   \end{align*}

   其中，$(4)$ 中分子 $\frac{\overline{r_1r_2}}{10^k+1}$ 不可能为整数，否则直接可以化简：

   \begin{align*} \frac{1}{n} = \frac{\frac{\overline{r_1r_2}}{10^k+1}}{10^k-1} = \frac{r_k}{10^k-1} &&(6) \end{align*}

   循环节长度直接由 2k 降为 k，这与之前的假设不符。因此在 $(5)$ 式中，$n$ 与 $(10^k+1)$ 必然存在公因数。而 n 为素数，则 $(10^k+1)$ 为 $n$ 的倍数。即 $(3)$ 式中 $\frac{(10^k+1)}{n}$ 为整数。

   因此，$r_1 + \frac{r_1+r_2}{10^k-1}$ 为整数。$r_1$ 本来就是提出来的前半循环节为整数，而 $\frac{r_1+r_2}{10^k-1}$ 是由两个长度为 k 的半循环节加和得到的分子，要求其为整数只可能是 1。所以 $r_1+r_2 = 10^k-1 = \overbrace{999\cdots}^k$

   得证。其它进制下的证明将底数 10 替换为对应进制即可。

计算程序

走马灯数搜索程序在得到以上结论后可以提升一些效率：

1. 使用 欧拉筛法，只在素数堆里进行验证。
2. 如果循环节小于 n-1，由于推论 1，则其至多 $\frac{n-1}{2}$，因此当计算到第 $\frac{n+1}{2}$ 项时，便可得出结论而不必计算后半段。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

function searchCarousel(n) {
    let primes = [];
    let eularboard = new Array(n + 1).fill(true);
    let carousels = [];

    const isCarousel = (p, radix = 10) => {
        let a = 1,
            pool = new Set(),
            end = (p + 1) / 2;
        while (pool.size < end) {
            a = (a % p) * radix;
            if (a == 0) {
                return false;
            }
            if (pool.has(a)) {
                return false;
            } else {
                pool.add(a);
            }
        }
        return true;
    };

    for (let i = 2; i <= n; i++) {
        if (eularboard[i]) {
            primes.push(i);
            if (isCarousel(i)) {
                carousels.push(i);
            }
        }
        for (let j = 0; j < primes.length; j++) {
            if (i * primes[j] > n) break; // 上限判断
            eularboard[i * primes[j]] = false;
            if (i % primes[j] === 0) break; // 跳过本轮后续筛选
        }
    }

    return carousels;
}

console.log(searchCarousel(1000));
// [7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193, 223, 229, 233,  257, 263,
// 269, 313, 337, 367, 379, 383, 389, 419,  433, 461, 487, 491, 499, 503, 509, 541, 571, 577,  593, 619,
// 647, 659, 701, 709, 727, 743, 811, 821, 823, 857, 863, 887, 937, 941, 953, 971, 977, 983]

其它进制

根据前文继续厘清思路，在其它进制下：

• 循环节长度的度量和进制数无关，表述形式不改变数值大小，在 10 进制下若为 16，在 8 进制下就写为 20，但都指的是同一个值。
• 但循环节的生成和进制有关，在当前位下不满足该进制条件下的整除，产生余数，就会退到下一位。
• 在其它进制下，$r = \frac{R^{n-1}-1}{n}$ 仍然成立，这是循环小数定义公式的变形。
• 仍然只需要在素数表中搜索走马灯数，因为 r 需要由 $n-1$ 个不同余的结果来确保循环长度上限。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

/* 计算各进制下的走马灯数 */
function searchCarousel(n) {
    let primes = [];
    let eularboard = new Array(n + 1).fill(true);
    for (let i = 2; i <= n; i++) {
        if (eularboard[i]) {
            primes.push(i);
        }
        for (let j = 0; j < primes.length; j++) {
            if (i * primes[j] > n) break; // 上限判断
            eularboard[i * primes[j]] = false;
            if (i % primes[j] === 0) break; // 跳过本轮后续筛选
        }
    }

    const isCarousel = (p, radix = 10) => {
        if (p === 2 && radix % 2 !== 0) return true; // 2 在奇数进制下必然是 1 位循环
        let a = 1,
            pool = new Set(),
            end = (p + 1) / 2;
        while (pool.size < end) {
            a = (a % p) * radix;
            if (a == 0) {
                return false;
            }
            if (pool.has(a)) {
                return false;
            } else {
                pool.add(a);
            }
        }
        return true;
    };

    for (thisRadix = 2; thisRadix < 100; thisRadix++) {
        let carousels = [];
        for (let i = 0; i < primes.length; i++) {
            if (isCarousel(primes[i], thisRadix)) {
                carousels.push(primes[i]);
            }
        }
        console.log(thisRadix, carousels.length);
    }
}
searchCarousel(10000);

搜索在各不同进制下， 10000 以内的走马灯数的数量，结果如下：

这里最有特点的是当进制为偶数的完全平方时，没有走马灯数，当进制为奇数的完全平方时，只有唯一一个走马灯数 2，循环节长度为 1。

特殊进制 $R=s^2$

无论在哪种进制下，循环节公式 $r = \frac{R^{n-1}-1}{n}$ 都仍然是成立的。其中 $R$ 为指定的进制数，$n$ 为测试数字， $r$ 为对应的循环节。根据前文所述，在一些情况下循环节会缩减。即当 r 的数字序列内部出现重复时，需要用更短的循环节代替。前文 $\frac{1}{11} = 0.\overline{0909090909} = 0.\overline{09}$ 便是例子。

而当 $R$ 为完全平方数 $s^2$ 时，将循环节公式展开为；

$$ r_{(R)} = \frac{(s^2)^{n-1}-1}{n} = \frac{(s^{n-1}+1)(s^{n-1}-1)}{n} = (s^{n-1}+1) \cdot \frac{s^{n-1}-1}{n}$$

因为 $s$ 是整数，注意到 $\frac{s^{n-1}-1}{n}$ 实际上是在 $s$ 进制下的 $\frac{1}{n}$ 的循环节表达式 $r_{s(s)}$，必定是有整数解的，并且在 $s$ 进制下的长度至多为 $n-1$ 。

现在将 $r_{s}$ 的数值从 $s$ 进制转到 $R$ 进制。因为$R=s^2$，$R$ 进制下一位数值可以容纳 $s$ 进制下两位。其长度在转换后正好一半 $\frac{n-1}{2}$，n 为素数。则：

$$r_{(R)} =(s^{n-1} + 1) * r_{s(s)} \overset{s\rightarrow R}{=\!=\!=} (R^{\frac{n-1}{2}} + 1) * r_{s(R)}$$

除 2 以外的 n 都是素数，所以 $\frac{n-1}{2}$ 是整数，$(R^{\frac{n-1}{2}} + 1)$ 是形如 $\overbrace{100\cdots01}^{(n+1)/2}$ 的整数，长度上正好使 $r_{s(R)}$ 重复两遍，而不会有叠加进位。

因此当进制 $R$ 为完全平方数时，任意 $n$ 的 $r_{(R)} = \frac{R^{n-1}-1}{n}$ 数字序列至少会内部重复一次，所以长度一定小于 $n-1$，也就不存在跑马灯数。

推论：当进制 $R=s^{(2^n)}$ 时，对应的最长循环节为 $\frac{n-1}{2^n}$ 向上取整。

费马小定律

在前文中数次使用 $r=(R^{n-1}-1)/n$，并且认定它是整数，是因为该式是从循环节的级数表示法计算出来的，而循环节本身一定是 n-1 位的整数，所以认为是整数。这个推论是单薄的，不能证明任意 $(R^{n-1}-1)/n$ 一定是整数。但其实关于这一点的正向证明早就在三百多年前就被书空白写不下学家完成了：

$(R^{n-1}-1)/n$ 为整数，其实本身就是 费马小定理 的一部分，即：

• 假如 $a$ 是整数，$p$ 是素数，那么 $a^p \equiv a (mod\ p)$
• 如果 $a$ 不是 $p$ 的倍数，则 $a^{p-1} \equiv 1 (mod\ p)$，这其实就是我们用的循环节公式。
• 如果 $a$ 是 $p$ 的倍数，则 $a^{p-1} \equiv 0 (mod\ p)$，这时就除尽不构成循环节了。

我们之前提到，进制的转换会影响循环节数字序列的变化，但不影响素数、整除等基于数本身的性质，因此这里的 R 在整除判定上反而和进制本身没什么关系，就是一个单纯的幂底数。

模运算与原根

当引入费马小定律后，我们进入了数论中模运算的领域。在前文中有一条推论 2：

只有当 n 为素数时，才有可能出现 n-1 位循环节。由于余数需要遍历 1 到 n-1 的所有数。若 n = pq，余数会在 q-1 与 p-1 的共同作用下循环，但 p-1 与 q-1 至少含有公因素 2，不可能产生所有 {0, p-1} 与 {0, q-1} 之间的全部两两组合。

其实是在数论的门外匆匆一瞥。思而不学则民科，在查询了一些资料后，发现其实这些走马灯数在数论中早已有深入的讨论。本文前半篇是从计算中发现现象，再从现象中思索原因。后半篇打算从原理出发，厘清前因后果全套链条。

欧拉函数 ϕ(n)

为了循序渐进，先来熟悉一下欧拉函数，其定义为：

$$\phi(n)=小于正整数 n 且与 n 互素（即最大公约数为1）的正整数的个数。$$

注意这些正整数不一定是素数，例如 8 与 9 是互素的，但两个都不是素数。

这个函数描述相比一般函数通俗得多，这也是数论领域的特色。尽管如此，这个函数实际上却并不用真的去逐个数个数，它确实是有公式通解的：

一般地，若 $n=p_1^{k_1}\cdot p_2^{k_2}\cdot \cdots \cdot p_m^{k_m}$，$p_m$为 n 的素因数，则：
$$\phi(n)=n\cdot\left(1-\frac{1}{p_1}\right)\left(1-\frac{1}{p_2}\right)\cdots\left(1-\frac{1}{p_m}\right)$$

• 当 n 为素数时，所有 $1\leq k < n$ 均有效，$\phi(7)=7-1=6.$
• 当 $n=p^k$ 为单一素数的幂时，, $\phi(8)=\phi(2^3)=2^3-2^2=4.$

欧拉公式便是上文 推论 6 的完整表述。略一思索也可以很快得出结论：『互素』其实就要求所有数都『避开』n 的所有素因数及其构成的合数集，这其实就是个变形的组合数与容斥问题，形式上类似于：

$$\phi(n)=n -\sum{\frac{n}{p_i}} + \sum{\frac{n}{p_ip_j}} - \sum{\frac{n}{p_ip_jp_k}} + \cdots + (-1)^{m-1}\frac{n}{p_1p_2\cdots p_m}$$

整理便可得到通项公式。

另外值得注意的是，欧拉公式中并没有显式地使用到幂数 $k_1 \cdots k_m$，实际上是因为 $n=p_1^{k_1}\cdot p_2^{k_2}\cdot \cdots \cdot p_m^{k_m}$ 中已经包含了所有的 $k$。

整数环与有限域

在模运算 $\ mod\ n$ 中，所有余数相同的整数都被视为等价的，即 $a+kn \equiv a \ mod\ n$，通常用最小非负余数 $[0],[1],[2],…,[n−1]$ 表示这些等价类。所有的余数类构成一个集合，称之为 n 的整数环，用奇怪符号 $\Bbb{Z}/n\Bbb{Z}$ 表示。

整数环也有类似一般逻辑的『加法』与『乘法』运算：

$$[a]+[b]=[a+b]\ (mod\ n)$$
$$[a]\cdot[b]=([a] \times [b])\ (mod\ n)$$

从进制的角度也很好理解，大于 n 的部分『进位』上去了，无法影响最后一位余数的运算。尽管模运算本身没有提到进制，但进制本身就是模的一种应用，通过两者的相似性有助于理解问题。

因为同余的是一类数而不是一个数，因此没法直接做除法运算，但可以定义一个 乘法逆元，当：

$$[a] \times [b] \equiv 1 \ mod\ n$$

时，$[a]$ 与 $[b]$ 就互为乘法逆元。

同样地，和乘法中 0 类似，模运算也会出现 $a \times b \equiv 0 \ mod\ n$ 的情况，零因子没有逆元。出现零因子的原因也很简单，当 n 为合数时，a b 两个非零元素相乘后结果能被 n 整除。

整数环仅仅是个称谓，任意整数 n 都有长度为 n 的整数环，并不能体现什么特殊性。但 n 不同会使得对应的整数环内部表现出不同的性质，则是深入探究的目标。

前面也看到，有些整数环内部存在零因子，有些则没有。判断方法很简单，当 n 为素数时，其整数环内部没有零因子，所有非零同余类都有对应的乘法逆元。这就和上面的 推论 2 一致了。此时整数环获得了一个新名字：有限域。

原根、欧拉定理与阶

当我们在关注循环小数的循环节时，本质上是对其不断循环取除的过程，将余数乘以进制，除以 n 以后再取余数，如此不断循环。如果将这个过程抽象一下，其实就是如下表述：

在模 n 的整数环中，对整数 $r$ 的幂次 $(r^1, r^2, \cdots, r^k)$ 不断进行模运算，则其或落入零因子，或者在一或多个非零同余间循环。

——实际上就是 r 进制下 1/n 的循环节计算过程，或者被整除，或者处于某个循环。

在数论课本中，对上述内容做了展开，有更严谨而完备的表述与证明过程：

欧拉定理

若 a 与 n 互素，则 $a^{\phi(n)} \equiv 1 \ mod\ n$

欧拉定理的证明过程大致分为四步：

1. 设模 n 的简化剩余系为 $S = {r_1, r_2, \cdots, r_{\phi(n)}}$ 满足 $gcd(r_i, n) = 1 $，且任意两项不同余，即 $r_i \not\equiv r_j$。

2. 将 $S$ 每个元素乘以 a，得到 $a\cdot S= {a\cdot r_1, a\cdot r_2, \cdots, a\cdot r_{\phi(n)}}$。此时若 $a\cdot r_i \equiv a\cdot r_j \ mod\ n$，则因 a 与 n 互素可消去，得到 $r_i \equiv r_j$，矛盾。 故 $a\cdot S$ 仍为简化剩余系，与 $S$ 包含相同的剩余类，仅顺序不同。

3. 由于两者同余，$S$ 和 $a\cdot S$ 的乘积相等：

   $$\prod_{i=1}^{\phi(n)} r_i \equiv P\ mod\ n \equiv a^{\phi(n)} \cdot P\ mod\ n \equiv \prod_{i=1}^{\phi(n)} (a \cdot r_i) $$

4. 由于 $P$ 是乘积，$r_i$ 与 n 互素，故 $P$ 也与 n 互素，两边同乘 $P$ 的逆元 $P^{-1}$，得到 $ a^{\phi(n)}\equiv 1 \ mod\ n$

阶与原根

由欧拉定理，可以定义 阶 与 原根：

1. $a$ 和 $n$ 为互素的正整数，则 $a$ 模 $n$ 的阶是满足 $a^k \equiv 1\ mod\ n$ 的最小正整数 $k$。
2. 若整数 a 与 n 互素，且其模 n 的乘法阶数等于 ϕ(n)，则称 a 是模 n 的一个原根。

或者使用另一个等效表述：

在模 n 的整数环中，若存在整数 $a$ 使得 $(a^1, a^2, \cdots, a^{\phi(n)})$ 能够遍历所有与 n 互素的剩余类，则 $a$ 称为模 n 的一个原根。

一次能遍历所有可能的余数，也就是循环节达到了最大值，这个最大值就是欧拉函数 $\phi(n)$。存在原根和最长循环是一体两面的。

应用

进一步地，若要 $\phi(n) = n-1$ 则 n 必须是素数，也就是说，走马灯数的分母 n 必须是素数。而当 $n=p $（素数）时，$\phi(n)=p-1$，欧拉定理简化为：$a^{p-1}\equiv 1 \ mod\ p $，这就是费马小定理！

到这里，就和前文的推论 2, 6, 7, 8 都对上了，甚至进一步地和 1, 4 也都能对上，一切全联系起来了！甚至我们还发现了之前没有注意到的 推论 10：n 可能有多个原根，则 1/n 可能在多个进制下有循环小数，并且所有可能的进制下，循环节的长度都是一样的。

接下来就是在欧拉定理的指导下，科学地判定走马灯数了：

1. 当 $n = 1, 2, 4, p^k, 2p^k$ 时存在原根，且阶为 $\phi(n)$。但只有 $\phi(n) = n - 1$ 时才能称为『走马灯数』，所以 n 必须是素数。
2. 当且仅当 $a^{n-1} \equiv 1\ mod\ n$ 时，a 即为原根，对应的 n-1 即为循环节长度。
3. 只需排查任意 $n-1$ 的真因子 $d$ 是否会使得 $a^d \equiv 1\ mod\ n$，如果有，则 $a$ 不是 $n$ 的原根。

程序优化

1. 添加 n-1 的真因子 d 的计算过程，这是个因式分解过程。考虑到我们之前已经计算出了素数表用于 n 的选择，这里可以复用。
2. 在实际计算中，$a^{\phi(n)}$ 仍然是个巨大的数，容易有溢出风险，但模运算支持乘法，因此可以把 $a^b\ mod\ n$ 的指数$b$拆分求模：
   $$((a\ mod\ n )^{b_1} \times (a\ mod\ n )^{b_2} \times \cdots \times (a\ mod\ n )^{b_k})\ mod\ n$$
3. 只验证所有 d 下 $a^d \not\equiv 1\ mod\ n$，有 false 直接退出。即排除推论 4 描述的情况。

据此我们写出了最终的走马灯数快速算法，直接应用欧拉定理，分解 n-1 的素因数，代替原来的试除法进行验证，并使用模运算的乘法分解式来计算 $a^b$ 以避免大数问题。这个算法在 n 较大时的效率比原来有了很大提升。旧算法搜索 100 万以内的走马灯数约需 20 分钟，新算法不到 1 秒，提升了 1000 倍。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

function searchCarousel(n, radix = 10) {
    let primes = [];
    let eularboard = new Array(n + 1).fill(true);
    let result = [];

    // 去重列举 n 的所有素因数
    function primeFactors(n) {
        const factors = new Set();
        let i = 0;
        while (n > 1 && primes[i] < n) {
            while (n % primes[i] === 0) {
                factors.add(primes[i]);
                n = n / primes[i];
            }
            i++;
        }
        if (n > 1) factors.add(n);
        return Array.from(factors);
    }

    // 计算 a^b mod n 的值
    function modExp(a, b, mod) {
        let result = 1n;
        a = BigInt(a) % BigInt(mod);
        b = BigInt(b);
        mod = BigInt(mod);
        while (b > 0n) {
            if (b % 2n === 1n) {
                result = (result * a) % mod;
            }
            a = (a * a) % mod;
            b = b / 2n;
        }
        return Number(result);
    }

    // 判断是否为全循环质数（原根验证法）
    function isCarousel(p, radix) {
        const phi = p - 1;
        if (modExp(radix, phi, p) !== 1) return false; // 费马小定理
        const factors = primeFactors(phi);
        for (const factor of factors) {
            const d = phi / factor; // 因为 factor 是独立最小素因数，使用对应的 d 可以最大化否决。
            if (modExp(radix, d, p) === 1) return false;
        }
        return true;
    }

    // 欧拉筛获得素数表
    for (let i = 2; i <= n; i++) {
        if (eularboard[i]) {
            primes.push(i);
            if (isCarousel(i, radix)) {
                result.push(i);
            }
        }
        for (let j = 0; j < primes.length; j++) {
            if (i * primes[j] > n) break; // 上限判断
            eularboard[i * primes[j]] = false;
            if (i % primes[j] === 0) break; // 跳过本轮后续筛选
        }
    }

    return result;
}

console.log(searchCarousel(1000000));
//[ 7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193, 223, 229, 233, 257, 263, 269, 313, 337, 367, 379, 383,
// 389, 419, 433, 461, 487, 491, 499, 503, 509, 541, 571, 577, 593, 619, 647, 659, 701, 709, 727, 743, 811, 821, 823, 857, 863, 887,
// 937, 941, 953, 971, 977, 983, 1019, 1021, 1033, 1051, 1063, 1069, 1087, 1091, 1097, 1103, 1109, 1153, 1171, 1181, 1193, 1217, 1223,
// 1229, 1259, 1291, 1297, 1301, 1303, 1327, 1367, 1381, 1429, 1433, 1447, 1487, 1531, 1543, 1549, 1553, 1567, 1571, 1579, 1583, 1607, 1619,
//...29400 more items]

进一步地，搜索非 10 原根即非 10 进制下的走马灯数时，除了 n 质数这一条件不变之外，其它也有一些应用定理可以优化的点：

1. 若 $g$ 的阶为 $\phi(n)$，则 $g^k$ 的阶为 $\frac{\phi(n)}{gcd(k,\phi(n))}$
2. 当且仅当 $gcd(k,\phi(n)) = 1$ 时，$g^k$ 的阶仍为 $\phi(n)$

故，假设 $g$ 是模 n 的一个原根，则所有其它原根可表示为 $g^k\ mod\ n$，其中 $gcd(k,\phi(n))=1$

所以算法上，应当从 $g=2$ 开始搜索第一个原根，先与 n 做互素判定，再应用上文的 isCarousel() 函数。找到 $g$ 后，再寻找 $k$ 令 $gcd(k,\phi(n))=1$，则其它原根即为 $g^k\ mod\ n$。gcd() 函数的计算比 modExp() 函数更高效。

• 以 n = 7 为例，$\phi(7) = 6$，
• 搜索到第一个 $g=3$，$3^6 \equiv 1\ mod\ 7$，且验证了 $3^2 \not\equiv 1\ mod 7$ 和 $3^3 \not\equiv 1\ mod 7$，所以 3 是原根。
• 搜索 $k$ 使 $gcd(k,6)=1$，且 $1<k<6$，有唯一解 $k=5$，计算 $3^{k=5}\ mod\ 7$，madExp(3,5,7) == 5，即 5 为第二个原根。注意 $g_2=5$ 与 $k=5$ 是两个独立的计算过程，虽然值一样。
• 对于其它进制，例如 7 应当在 10 进制下也是 6 位循环，是因为 $10\ mod\ 7 \equiv 3\ mod\ 7$。若要在程序的输出结果上体现出来更大的进制，只需要输出前加个循环即可。

1

// 代码略

江南水乡的婴儿呱呱坠地时，窗外的桃枝上正结着薄霜。接生婆用红绸裹住婴儿，在族谱写下"一岁"的瞬间，这个生命就完成了对时间的第一次突围。虚岁不是简单的数字游戏，而是一把丈量生命的游标卡尺，在文明的刻度上为成长预留出余量。

中国人对时间的认知始终保持着弹性。上古时期，先民们在龟甲上刻下"岁"字时，就赋予了它双重含义：既指周天星辰的运行周期，又指人类生命的生长韵律。《周礼》将三百六十五日划为"岁实"，又在每个朔望月设置"虚日"，这种虚实相生的智慧投射到年龄计算上，便形成了虚岁与实岁并行的独特系统。孩子尚未学会走路，却已在时间的账簿上拥有了完整的一年，如同春分前的惊蛰，用先声夺人的雷鸣预告生命的力量。

传统礼仪中的冠礼最能体现虚岁的深意。少年十五行冠礼，实则身体仅十四龄。当长辈将缁布冠缓缓加诸其首时，不仅是在完成形式上的成年仪式，更是在进行一场庄重的心理暗示。提前赋予的年龄如同试穿的礼服，让少年在宽出一寸的衣襟里，提前感知未来需要承担的重荷。就像苏州园林的借景艺术，虚岁借来未来的光阴，将成长的远景拉进现实的取景框。

这种时间预支的智慧，在当代依然焕发着生命力。上海弄堂里的阿婆教孙子打算盘，会故意把算珠拨快一个档位；岭南祠堂的族老主持祭祖，总让孩童站在比实际年龄高一级的台阶上。这些细微的"误差"，实则是为成长铺设的缓冲地带。当成都的茶馆里老人们谈论"虚岁三十该立事"时，他们深谙人生如戏，需要提前一年拉开帷幕，才能让主角从容登场。

生命的年轮从不是严丝合缝的同心圆。虚岁制在时间褶皱里藏着的这枚硬币，正面镌刻着"未雨绸缪"，背面铭刻着"留有余地"。它提醒我们，成长从来不是瞬间的蜕变，而是层层铺垫的渐进。就像黄山顶上的迎客松，在嶙峋山石间蜿蜒的根系，早在地面枝干成形前，就已在地下完成千百次生长的预演。

担心 TF 卡易损坏导致重新配置系统的麻烦，我准备定期使用 dd 对树莓派进行全卡备份：

1

sudo dd if=/dev/mmcblk0 of=/mnt/momoda/$(date +%Y-%m-%d)-sysbak.img

尽管使用的是 64GB 的 TF 卡，生成的镜像文件也是 64GB。然而，实际系统盘的数据占用并不多，大部分空间是空白。因此，我选择通过管道压缩备份数据的方式以节省存储空间，并测试了几种支持管道的压缩算法的性能。

测试环境：

系统存储被分为 /dev/mmcblk0p1 和 /dev/mmcblk0p2 两个分区，实际上是同一张 TF 卡，目前总空间使用率约为 10%。目标是定期对全盘进行备份以应对意外，同时通过压缩方式保存数据，减少备份所需的存储空间。由于备份操作在运行中的设备上进行，数据在此期间仍可能被读写，因此选择的压缩算法需要尽可能高效，以减少数据变动的影响。

全盘空间信息：

1
2
3
4
5
6
7
8
9

~ df ⮐
Filesystem          1K-blocks       Used  Available Use% Mounted on
udev                  3728364          0    3728364   0% /dev
tmpfs                  799744       3176     796568   1% /run
/dev/mmcblk0p2       59551920    5034336   51780988   9% /
tmpfs                 3998704          0    3998704   0% /dev/shm
tmpfs                    5120         16       5104   1% /run/lock
/dev/mmcblk0p1         522232      66728     455504  13% /boot/firmware
tmpfs                  799740          0     799740   0% /run/user/1000

首先生成一个未压缩的镜像文件：

1
2
3
4
5

~ sudo dd if=/dev/mmcblk0 of=/mnt/exDisk/sys.img bs=4M ⮐
14909+1 records in
14909+1 records out
62534975488 bytes (63 GB, 58 GiB) copied, 1459.15 s, 42.9 MB/s
# 如果压缩算法在树莓派的 CPU 上能达到 42.9 MB/s 的处理速度，就可以认为压缩完全不影响原备份速度。

由于大部分空间是空数据，如果以 63GB 为基数计算，压缩后的文件大小差异将显得微不足道。因此，下表中的百分比均以 zstd 默认压缩率的结果作为 100% 的基准进行比较。

同时，受限于树莓派的储存方式，该镜像放在了另一台 NAS 上，两台机器间的传输能力为千兆有线，纯传输时间最低约为 536 秒。

测试结果

• 总体来看，zstd 的性价比明显高于 xz 和 7z 等其它算法。
• lz4 速度极快。尽管没有多线程参数，单核性能仍快于其它算法的双核表现。然而，其较低的压缩率不适合作为备份用途，更适合处理持续性超大量数据的场景。
• 使用 lz4 -fast 时，压缩算法的数据吞吐量已经超过数据传输的上限。这是硬盘的瓶颈，不是算法的处理上限。
• 由于 TF 卡的读取速度更低，因此对于本文的需求而言，压缩时间低于 1459.15 秒时，是算法在等待数据读取，此时时间消耗的差异不再有区别。
• 综合考虑数据分布和性能表现，zstd -T2 -10 是当前测试中最合适的选择。可以进一步测试 -8、-12 等参数，以在 TF 卡读取时间内达到最佳压缩率。由于 TF 卡和 CPU 性能的差异，最佳方案会有所不同，需根据具体情况选择。
• 7z 和 xz 实际上用的都是 lzma 算法，两者的压缩率、性能和内存占用都没有什么实际差异。
• 在测试过程中，发现 xz -9 和 7z -mx=9 运行时的内存占用也高达 2GB 多，如果是小内存版本的树莓派也根本跑不起来。

简述

直接用浏览器打开 VMWare 的 CDS 地址：

https://softwareupdate.vmware.com/cds/vmw-desktop/

其中

fusion 目录是 MacOS 下的 VM Fusion，x86 和 M1 的 MacOS 都能用。
ws 目录是 VM Workstation。有 Windows 和 Linux 版。

一路点进去，选完版本号编译号系统版本后就可以下载 tar 包。用 7zip 等软件可以解压出 VMWare 的安装包。

比如发文时 ws 的最新版本号是 17.5.2，下载链接就是
https://softwareupdate.vmware.com/cds/vmw-desktop/ws/17.5.2/23775571/windows/core/VMware-workstation-17.5.2-23775571.exe.tar

记

自从被博通收购以后，VMWare 对个人用户免费了。虽然以前满地都是序列号，个人用户实质也是随便用的，但现在毕竟有官方直接授权了。但博通的帮助系统是真难用，找了半天找不到下载链接。然后想起 vmware 自带的更新是走一个 cds 的，尝试搜索了一下还在，这比在帮助系统里找方便多了。

记录一下。

记 II

遇到了 VMWare 17.5.1 跑 Windows XP 黑屏的 bug，搜索得知 VMWare 一直在不断地删减功能，未来甚至要把底层技术换成 KVM/HyperV，自己只做套皮。

又一个新版本不如老版本的软件。所以挖了一下把 CDS 上现有的老版本下载链接都提了出来。

1
2
3
4
5

wget -r -np -nd -e robots=off \
    --spider --wait=1 -l inf \
    https://softwareupdate.vmware.com/cds/vmw-desktop/ \
    2>&1 | grep -E '^--'  \
    | awk '{ print $3 }' >> links.txt

整理结果：

质数，也叫素数，是对 Prime Number 的不同翻译。素数的叫法大概率来自日语，素在日语里引申出了根本，源头以及不可再分割的意思，这一点古汉语里是没有的。因此诞生出了一大批和制词汇，比如元素，素材，要素等等。而 Prime Number 正好符合这个特性，叫素数再合适不过了。现代汉语由于引入了上面这些和制词汇，所以素也开始有了上述意思。至于港台地区保留了这个叫法也是因为延袭了民国时期的翻译习惯，没有做变更而已。

《为什么叫素数》 来源：Spenser Sheng
链接：https://www.zhihu.com/question/22456389/answer/967574484
至此，质数这个名词就在中国成了官方认可的科学术语，并一直没用至今。但是，素数这个说法并没有因此退出历史舞台。同时代的秀多著作，还是用了这个名词，比如华罗庚先生于 1940 年出版的《堆垒素数论》等。这直接影响了后来更多数学家在著书立说时没用素数这个名词（可能因为要考虑与他人表述的统一性），如陈景润研究“哥德巴赫猜想”也是讲素数。……在中小学教材上，目前是统一使用质数这个名词，但是标明了“质数，又叫素数”《梳理源流，叩问本质——对质数、合数名词的考证与思考》顾志能 https://www.doc88.com/p-1478574668277.html

如何计算 N 以内的素数，或者判断 p 为素数，在编程领域是孪生的两个问题，也是很多编程语言的入门题目。这个题目优点在于有相当多的解法和优化空间，尤其适合用来练习算法和程序优化。这里我们就来逐步优化一个素数计算程序。

基础实现

我们从最简单的开始，先写一个函数 isPrime(n) 判断 n 是否为素数，这个判断函数严格按素数的定义实现：

除了 1 和 n 之外，没有其他整数能整除 n

也就是说，将 2 ~ n-1 之间的每个整数都试一遍，看看是否能整除 n。再用 isPrime(n) 逐个找出 < N 的所有素数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

function isPrime(n) {
    // 判断 n 是否为素数
    for (let i = 2; i < n; i++) {
        if (n % i === 0) return false;
    }
    return true;
}

function calcPrimes(n) {
    // 找出 <N 的所有素数
    let primes = [];
    for (let i = 2; i <= n; i++) {
        if (isPrime(i)) {
            prime.push(i);
        }
    }
    return primes;
}

console.log(calcPrimes(100));
// [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

※ 本文演示代码为 JavaScript，方便在浏览器控制台内测试。但考虑到与正文内容保持一致，并未使用 JavaScript 内置的 .map(), .forEach(), .filter() 等遍历函数，而只用普通的 for loop if else 等代码。

简单优化

显然上面的代码是有很多无效的计算的，一些简单的数学知识就可以大幅地优化计算效率：

1. 所有的偶数都不是素数，除了 2。

1
2
3
4
5
6
7
8
9
10

function calcPrimes(n) {
    let primes = [2]; // for 循环从 3 开始，因此直接把特例 2 加入
    for (let i = 3; i <= n; i = i + 2) {
        // 从 i++ 优化为 i=i+2
        if (isPrime(i)) {
            prime.push(i);
        }
    }
    return primes;
}

2. 同上，所有 3 的倍数都不是素数，除了 3。因此，我们可以只检查 6k ± 1 的数。

1
2
3
4
5
6
7
8
9
10
11
12
13

function calcPrimes(n) {
    let primes = [2, 3];
    for (let i = 6; i <= n; i = i + 6) {
        // 计算 6k
        if (isPrime(i - 1)) {
            prime.push(i - 1);
        } // 判断 6k-1
        if (isPrime(i + 1)) {
            prime.push(i + 1);
        } // 判断 6k+1
    }
    return primes;
}

※ 理论上我们可以继续讨论 5 的倍数，但这会使 for 循环需要处理 30k ± a，且 a 有 1,7,13 等多个值。代码变得复杂，而优化效果一般。过于复杂的代码会让你不幸福。

3. 判断素数时，只需要检查到 sqrt(n) 即可。如果 n 有一个大于 sqrt(n) 的因子，必然对应一个小于 sqrt(n) 的因子，则之前的循环中已经判断过了。

1
2
3
4
5
6
7
8

function isPrime(n) {
    if (n < 2) return false;
    for (let i = 2; i * i <= n; i++) {
        // 从 i < n 优化为 i*i <= n
        if (n % i === 0) return false;
    }
    return true;
}

更多优化

进一步思考 isPrime 函数，我们发现，不必让所有小于 sqrt(n) 的数都试一遍，只需要用之前已经计算出来的素数尝试即可。但我们需要整体调整代码，使 isPrime() 函数可以调用 primes 数组。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

function calcPrimes(n) {
    let primes = [2, 3];

    function isPrime(n) {
        for (let i = 0; primes[i] ** 2 < n; i++) {
            // 在 isPrime 函数中可以直接使用 primes 数组
            if (n % primes[i] === 0) return false;
        }
        return true;
    }

    for (let i = 6; i <= n; i = i + 6) {
        if (isPrime(i - 1)) {
            prime.push(i - 1);
        }
        if (isPrime(i + 1)) {
            prime.push(i + 1);
        }
    }

    return primes;
}

console.log(calcPrimes(100));
// [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

筛法

以上的方法都是用较小的数字去尝试整除 N 以判断 N 是否为素数，这种方法称为试除法。从直观考虑，每个 N 都会被 2,3,5 ... 等数字都除一遍，比如 49 和 77 都会在 2, 3, 5 上各浪费三次计算。如果可以绕过无效除法，只用 7 去尝试 49，而 2,3,5 则不去浪费时间，耗时会大幅会降低。事实上这种方法是存在的，并且非常简单，称为『筛法』。

简单到解释起来只需要两句话：

• 从 2 开始，将所有 2 的倍数都划掉；从 3 开始，将所有 3 的倍数都划掉；从 5 开始，将所有 5 的倍数都划掉；
• 以此类推。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

function calcPrimes(n) {
    let primes = [];
    let numboard = new Array(n + 1).fill(true);

    for (let i = 2; i < n; i++) {
        if (numboard[i]) {
            primes.push(i); // 将素数加入 primes 数组
            for (let j = i * 2; j < n; j += i) {
                // 将 i 的倍数都划掉
                numboard[j] = false;
            }
        }
    }

    return primes;
}
console.log(calcPrimes(100));
// [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

这个筛法也被称为埃拉托斯特尼筛法，是为了纪念古希腊数学家埃拉托斯特尼而命名的。这个方法看似很简单，实际也确实很简单，是一种朴素的空间换时间方案。算法改除法为乘法，同时避免互质两数相遇，避免了很多无效运算。虽然当 N 非常大时，numboard 数组会占用很多内存，但时间优势更为明显。

※ 由于这个方法如此简单，普通人也能想到，不值得特地冠名。我猜测这个名字只是为了纪念历史名人，不能说明是他最早发明了这个方法。

同样，上述代码也有一些优化空间，比如：

• 划掉倍数时，可以从 j = i*i 开始，因为 2i, 3i ... (i-1)i 在对应的 2, 3 ... i-1(或其因数)那一轮时已经被划掉了。
• 当划到 i > sqrt(N) 时，则剩余未标记数也都是素数。假设有任意合数 m 满足 sqrt(N) < m < N，则 m 的素因数 p1,p2 ... 中最小的一个 p 也必然小于 sqrt(N)。而我们已经找到所有小时 sqrt(N) 的素数，并划掉其倍数了，所以找到 p 时必定已经划掉 m。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

function calcPrimes(n) {
    let primes = [];
    let numboard = new Array(n + 1).fill(true);
    let sqrt_n = Math.sqrt(n);

    for (let i = 2; i <= n; i++) {
        if (numboard[i]) {
            if (i > sqrt_n) break; // 根据上文证明，剩余未标记数也都是素数
            for (let j = i * i; j < n; j += i) {
                // j = i * 2 -> j = i * i
                numboard[j] = false;
            }
        }
    }

    for (let j = 2; j <= n; j++) {
        if (numboard[j]) {
            primes.push(j);
        } // 清点未标记数
    }

    return primes;
}

console.log(calcPrimes(100));
// [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

在算法进入筛法层次后，isPrime() 函数与 calcPrimes() 函数的功能区别也就不再明显。为了判断 N 是否为素数，需要至少计算到 sqrt(N) 的素数，而上文又证明剩余的未标记数也都是素数。差别只是最后判断一下 N 是否在 primes 数组中而已。

另外，在埃氏筛法中，6k ± 1 等手段不再有优化作用。构造一个『只有奇数的筛盘』是可以的，但后续的划除操作需要大量的脚标奇偶运算，反而会增加计算量。

位运算

在埃氏筛法中，内存占用是一个问题，尤其在筛法本身就是为了计算大 N 而生的。在实际应用中，有时也需要在有限的内存空间内尽量优化。注意到 numboard 数组只是用来标记是否为素数，只需要 0,1 两个状态，因此可以用位运算来减少内存占用。这使得内存的占用降为原来的 1/8。但此时算法的实现与编程语言也开始耦合，不是每种编程语言都支持直接的位运算，例如 JavaScript 的位运算就只是某种『模拟』，效率并不高。Javascript 也不支持直接修改某个 bit 位，而要经过多次字节操作来实现。尽管如此，我们仍可以编写一个 setBitFalse() 函数来模拟这个操作：

1
2
3

function setBitFalse(arr, i) {
    arr[i >> 5] &= ~(1 << (i & 31));
}

代码解释：

• JavaScript 里数值默认 32 位，通过 i >> 5 相当于 i/32 取整，可以得到对应的数组下标，即具体操作哪个元素。
• i & 31 将 i 与 31 (00001111) 进行按位与操作，相当于对 32 取余，得到具体操作该元素的第几个 bit 位。
• 1 << (i & 31) 将 1 左移 i & 31 位，得到一个只有第 i & 31 位为 1 的数，即只有指定 bit 位是 1 其它都是 0 的掩码。
• ~( ... ) 按位取反以后就得到了 1111....0...111 的反掩码，只有指定 bit 位为 0。
• &= 是按位与且赋值，将数组指定元素与反掩码进行按位与运算，其它位与 1 后值不变，指定 bit 位则变为 0。目标达成。

然后，我们需要写一个 isBitTrue() 函数来判断某个位是否为 1：

1
2
3

function isBitTrue(arr, i) {
    return (arr[i >> 5] & (1 << (i & 31))) !== 0;
}

类似地，函数计算数组下标和 bit 位置，生成一个只有指定 bit 位是 1 的掩码，将其与数组元素进行按位与操作。因为按位与操作遇 0 必 0，所以其它位都会被掩码数值置为 0。如果结果不为 0，则说明数组元素与掩码在该 bit 位上都是 1。

最后，我们还需要写一个 bitsToNums() 函数，将位运算后的结果数组转换为素数数组：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

function bitsToNums(arr,n) {
    let primes = [];
    for (let i = 0; i < arr.length; i++) {
        for (let j = 0; j < 32; j++) {
            if (arr[i] & (1 << j)) {
                primes.push(i * 32 + j);
            }
        }
    }

    let real_primes = [];
    for (let j = 2; j <= n; j++) {
        real_primes.push(primes[j])
    }
    // 因为是字节位长的整倍数，位运算可能会多一些数字，因此要二次筛选一下
    let real_primes = [];
    for (let j = 2; j <= n; j++) {
        if (primes[j] <= n) real_primes.push(primes[j]) else break;
    }

    return real_primes;
}

将几个函数代入原程序中的对应位置，我们可以得到一个位运算版本的埃拉托斯特尼筛选算法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

function calcPrimes(n) {
    let bitboard = new Array((n >> 5) + 1).fill(-1); // 位运算数组
    let sqrt_n = Math.sqrt(n);

    for (let i = 2; i <= n; i++) {
        if (isBitTrue(bitboard, i)) {
            if (i > sqrt_n) break;
            for (let j = i * i; j <= n; j += i) {
                setBitFalse(bitboard, j);
            }
        }
    }

    return bitsToNums(bitboard);
}

console.log(calcPrimes(100));
// [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

可以看到，位运算版本的代码更为复杂，原本只是根据下标取数组元素并判断 true/false 的操作，现在需要生成掩码，计算偏移量，转换位数等多次计算。只有在 N 极大，以至于必须减少内存占用的情况下，才需要使用位运算。在一般情况下并不值得作此妥协。

同时，对于 JavaScript / Python / java 等高级语言，位运算的效率并不高。即使实际情况确实拮据，也必须先考虑更换编程语言。比如在 C++ 中，std::bitset 这种数据结构，就特别适合于埃氏筛法结合位运算计算素数表。

分段筛法

说到内存不足，很容易就能想到另一种解决方案：分段筛法。在筛选过程中，我们把素数放进了单独的 primes 数组，而 numboard 数组则只是用来筛除合数。当 primes 找到大小为 p 的素数，就意味着 p*p 之前的 numboard 数据已经没有用了，如果我们能把这部分内存释放掉，就能大幅减少内存占用。

因此从逻辑上讲，当求 N 以内的素数时，我们只需要保留 sqrt(N) 以内的素数，并将整个 N 范围划为若干块，对每一块使用 primes 数组进行筛选。如果不考虑分块本身的内存分配回收时间，单只考虑数值与操作量的话，分段筛法的时间复杂度与原始筛法是相同的。对于每一个被筛选的元素，都会被其每一个独特的素因数筛选一次，无论是在全局 numboard 中，还是在局部 numboard_i 中，都一样。

基于前文的铺垫，似乎分块大小设置为 sqrt(N) 是直觉上合理的方案。但实践中，分块大小更多与内存大小有关。比如 nodejs 的默认内存上限为 2GB，尽管可以通过 --max-old-space-size 参数调整，但在个人电脑环境下再提升其实也并不会有什么质变。只要确保 primes 数组空间的前提下仍有足够的内存，就是合理的块大小。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

// JavaScript 数值默认 4Byte，1GB 内存可以存储 2^28 个数值。约 2.68 * 10^8 个。
const BOARD_SIZE = 10; // 分块大小，基于演示目的设置为 10

function calcPrimes(n) {
    let primes = [];

    function cleanBoard(s, e) {
        // 分块筛选区间 [s, e] 素数的子函数
        let board = new Array(e - s + 1).fill(true);
        let sqrt_e = Math.sqrt(e);

        for (let i = 0; i < primes.length; i++) {
            // 用先前得到的素数筛选本区间
            let p = primes[i];
            let start = Math.ceil(s / p) * p;
            for (let j = start; j <= e; j += p) {
                board[j - s] = false;
            }
        }

        for (let i = 0; i < board.length; i++) {
            // 本区间埃筛
            if (board[i]) {
                let p = i + s;
                if (p > sqrt_e) break;
                for (let j = p * p; j <= e; j += p) {
                    board[j - s] = false;
                }
            }
        }

        for (let i = 0; i < board.length; i++) {
            // 清点本区间的素数
            if (board[i]) {
                primes.push(i + s);
            }
        }
    }

    cleanBoard(2, Math.min(BOARD_SIZE - 1, n)); // 第一个区间特殊处理，去掉 0 和 1
    for (let i = 1; i < Math.ceil(n / BOARD_SIZE); i++) {
        cleanBoard(i * BOARD_SIZE, Math.min((i + 1) * BOARD_SIZE - 1, n));
    }

    return primes;
}

console.log(calcPrimes(100));
// [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

欧拉筛法

在埃氏筛法中，我们注意到，数字 6 会被 2, 3 划掉两次，30 会被 2, 3, 5 划掉三次。所有合数都会被其不同的素因数重复操作若干次，浪费还是有不少。但我们并不能去判断筛盘上某个具体的数是否已经划掉，因为『读取与判断』的时间本身可能比直接划掉时间还长。

进一步用力注意到，筛法的本质是从『验证素数』变成了『排除合数』，而且是有规律的排除合数。考察一个合数 6 = 2 x 3，我们希望在算法中，2 x 3 只发生一次，而不是在 2 与 3 的循环中各发生一次。对于 12 也一样，尽管它有 2 x 2 x 3, 4 x 3, 6 x 2 这三种分解方式，但我们希望只有其中一种发生，其它的被跳过。

仿佛一头雾水但哪里又有一些灵感，这好像是可能的，只要足够聪明合理地安排合数生成，好像可以做到？

事实确实如此。欧拉筛法就是为了解决这个问题而生的。具体操作是这样的：

1. 有一个自增的数字 i，从 2 开始。
2. 有一个数组 primes，存放已经找到的素数。
3. 如果筛板上 i 还没被划掉，则 i 是素数，将其加入 primes 数组。
4. 用每一个 i 乘以当前所有 primes 中的素数，得到的数都是合数，在筛板将它们标记为 false 划掉。
5. 如果 i 是 primes 中某个素数的倍数，那么 i * p 这个数会被标记 false，但后续直接跳过，进到 i+1 轮。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

function calcPrimes(n) {
    let primes = [];
    let eularboard = new Array(n + 1).fill(true);

    for (let i = 2; i <= n; i++) {
        if (eularboard[i]) {
            primes.push(i);
        }
        for (let j = 0; j < primes.length; j++) {
            if (i * primes[j] > n) break; // 上限判断
            eularboard[i * primes[j]] = false;
            if (i % primes[j] === 0) break; // 跳过本轮后续筛选
        }
    }

    return primes;
}

console.log(calcPrimes(100));
// [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ]

在欧拉筛中，我们通过 i * primes[j] 的方式来筛除合数，且当 i 本身为 primes[j] 的倍数时，跳过 i 与后续素数的计算环节。例如：

• 当 i = 4 时，本应划掉 4 * 2 与 4 * 3，但可以跳过 4 * 3，因为它会被后续的 6 * 2 划掉。
• 当 i = 35 时，会先划掉 35 * 5，而 35 * 7 会被跳过，因为它将在 i = 49 时被 5 划掉。

一般来说，如果 i = k * primes[j]，则接下来要划掉的 s = i * primes[j+1] = k * primes[j] * primes[j+1] 会在 i' = k * primes[j+1] 轮次时会被 primes[j] 划掉（因为 primes[j] 是 i' 的最小素因子），所以此时可以先跳过。

可以看到，在欧拉筛法中，由于 primes 是从小到大排列的，因此每个合数只在其最小素因数 m 对应的 i = n / m 轮次时，才会被划掉。

在这个算法中，合数的分解是唯一且有序的，因此欧拉筛与埃筛内存占用一样，但算法复杂度优化到了大约 O(n)。——严格来说，因为一些必须的额外判断，实际是 $O(nlog_2log_2n)$。

※ 数学先师的伟业至今仍被人传颂，有一个算法以他的名字命名。纵然身为一个 18 世纪数学家，那个时代还远未出现计算机。

素数的分布

有了高效的素数筛选算法，我们可以开始搜索大范围的素数了。计算 1 亿以内的素数，并每 10 万分为一组，统计每组的数量。可以看到，虽然有波动，但整体是逐渐减少的，并且这个图好像有点眼熟？

事实上，确实有一个 素数渐近分布定理，大致描述了素数的分布规律：

设不大于 n 的所有素数的个数为 π(n)，则 π(n) ≈ n / ln(n)

具体举例来说，当 n = 1 亿时，使用前文的程序计算出素数具体有 5761455 个，而 $n/ln(n)\approx 5428681$，两者比例为 1.061，较为近似。而当 n→∞ 时，值趋近于 1。

另外也可以得到一个推论，在 n 附近随机取一个数字，它是素数的概率是 $P(n)\approx 1/ln(n)$。函数下降不算快，就算 n = 1 亿，概率也有 1/18.42，5.43%。

这个定理已经被证明，但这个值本身只是一个渐近值。也就是说，虽然 $\lim\limits_{x \rightarrow 0}\frac{\pi(n)}{ n / ln(n)} = 1$，但 $\lim\limits_{x \rightarrow 0} [\pi(n) - n / ln(n)]$ 却是越来越大的。从函数图像上而言，两者只是越来越接近于平行，但距离却越来越远。

素数的分布规律是数论的一个重要研究方向。后来数学家又提出了几个更精确的估算公式，这一方向的顶点便是著名的黎曼猜想。

梅森素数与 Lucas-Lehmer 检验。

在了解素数的大致分布，并回顾欧拉筛法以后，我们发现即使是 O(n) 的算法也仍然有很大的局限性。内存限制是最初的现实问题，但也只是所有问题中最简单的一个，更大的限制仍然是来自于算法本身。我们要筛到至少 sqrt(N)，才能得到 N 是否为素数的结果。在 N 非常非常大时，O(n) 的算法仍然太慢了。

仔细思考发现，这是因为筛法天然地要求算法结果的素数是连续的，如果跳过某个数，那么后续的合数也可能被认为是素数，整体就会崩盘。那么是否存在『跳过一部分数，快速往更大的素数 N 前进』，以及在不知道小于 sqrt(N) 全部素数的情况下，判断 N 是否为素数的算法呢？

两者都是存在的，并且有多种方法。典型例子是梅森素数和配套的 Lucas-Lehmer 算法。

梅森素数是形如 $2^p - 1$ 的素数，并且简单可以证明，当 $2^p - 1$为素数时，p 也是素数。

当 p 是合数时，令 p = ab，有：
$2^p-1= (2^a)^b - 1 \overset{ x=2^a}{=\!=\!=} (x - 1)\cdot[x^{b-1}+x^{b-2}+\cdots+x+1]$ 为合数。

尽管反过来并不一定成立，当 p 是素数时，$2^p - 1$ 仍可能是合数，比如 2^11 - 1 = 2047 = 23 x 89。但即使如此，梅森素数依然有很多优点，比如：

• 比普通的序列增长得更快，可以跳过很多中间的数，尽快找到更大的素数。
• $2^p - 1$ 中包含素数的比例比全体自然中要高，因此更容易找到素数。
• 否定性验证很容易，只需要验证 p 是否为素数。正向验证难度也比一般的低。

Lucas-Lehmer 算法是验证梅森素数的一种方法。它构造了一个递推公式：

\begin{equation*}
s_i=\begin{cases} 4, & i=0 \\ (s_{i-1}^2-2) \mod (2^p-1), & i>0 \end{cases}
\end{equation*}

则当 $s_{p-2} \equiv 0 \pmod{2^p-1}$ 时，$2^p - 1$ 为素数。其算法复杂度为 $O((log\underset{2}{}n)^2)$。

※ 在大数环境下的计算需要编程语言和专门的数学包支持，代码就不放了

由于梅森素数拥有的诸多优点，目前已知的最大素数都是梅森素数，也有专门的 互联网协作项目 GIMPS 来众算更大的梅森素数。该项目前段时间正好确认了 M(57885161) 是一个新发现的素数。

质性检验

从梅森素数我们得到启发，在合理添加一些额外的限制条件后，素数的出现概率会提高，检验手段也可以针对性设计而变得简单，这对大素数的搜索非常有帮助。

在这个层次上，calcPrimes() 基本就不再是人们的目标，而 isPrime() 则不断地发展，成为专门的『质性测试』算法。这些算法的目标是尽可能快地判断一个数是否为素数，而不是找到所有素数。在这个领域，有很多经典的算法，比如：

Lucas-Lehmer 检验

上文提到的梅森素数检验算法。仅适用于 $2^p - 1$ 形式的素数。

Pépin 检验

验证费马数 $2^{2^n} + 1$ 是否为素数的算法。但由于费马数目前只发现 n=0,1,2,3,4 五个，后面都是合数，甚至有猜想认为 n > 4 时皆为合数。因此 Pépin 检验的应用范围有限。而欧拉和卢卡斯接力证明了费马数的素因数皆可表达为 $k^{2n+2} + 1$，为具体的因数分解也提供了不小的方便。

※ 感觉 Pépin 被费马坑了。

Miller-Rabin 素数测试算法

一个随机性的否定性算法。当测试结果为 false 则一定是合数，测试结果为 true 可能为素数。算法基于费马小定理和二次探测定理，在进行 n 次测试后，错误概率为 $1/4^n$。

费马小定理：如果 p 是一个素数，而 a 是任意不是 p 的倍数的整数，则 $a^{p-1} ≡ 1 (mod\ p)$。这意味着对于任意素数 p，选择一个不是 p 的倍数的整数 a，计算 $a^{p-1}\ \% \ p$，如果不等于 1，则 p 一定不是素数；等于 1，则 p 可能是素数 。
二次探测定理：假设 p 是一个素数，我们可以将 p-1 写为 $2^s * d$ 的形式，其中 d 是奇数。则对于任意整数 a，如果存在整数 x 满足：

1. $a^d ≡ 1 (mod\ p)$
2. 存在一个 i 满足 0 ≤ i < s，使得 $a^{2^i * d} ≡ -1 (mod\ p)$

以上两个条件之一，则 a 是一个模 p 的非平凡平方根，即 $a^2 ≡ x (mod\ p)$：

假设我们要测试数字 p = 17 是否为素数。

1. 选择不是 p 的倍数的整数 a，假设选择 a = 3。计算 $a^{p-1}\ \%\ p = 3^{16} \% 17 = 1$ 结果等于 1。根据费马小定理，17 可能是素数。

2. 根据二次探测定理，当 p = 17，有 $p-1=16=2^4*1$，因此 s = 4，d = 1。我们需验证任意整数 a，可以满足以下两个条件之一：

   • $a^d = a ≡ 1 (mod 17)$
   • 存在一个 i 满足 0 ≤ i < 4，使得 $a^{2^i} ≡ -1 (mod 17)$

   假若仍然选择 a=3，则 $3 \% 17 =3 \neq 1$，不满足条件一，继续检验条件二：

   • $i=0, 3^{2^0}\ \% 17 = 3$
   • $i=1, 3^{2^1}\ \% 17 = 9$
   • $i=2, 3^{2^2}\ \% 17 = 13$
   • $i=3, 3^{2^3}\ \% 17 = 16 = -1$，满足条件。因此 17 可能是素数。

继续类似地进行多轮测试，如果任意一轮验证失败，则 p 一定是合数，反之随着验证轮次增加，p 为素数的可能性便越来越大。

AKS 类质性测试

一个确定性的质性测试算法，复杂度为 $O(log^{6+\epsilon}n)$。AKS 是第一个被发表的一般的、多项式的、确定性的和无依赖的素数判定算法。算法基于一个简单定理：

当 $(x+a)^n\equiv(x^n+a)(\text{mod}\ n)$ 时，n 为素数。

后续 ASK 被不同科学家多次改进，因此有多个版本，统称为 AKS 类算法。AKS 类算法也是目前效率最高的确定性算法，只要算法给出 true 的结果就一定是素数。

※ ChatGPT 与 Copilot 对本文亦有贡献。

1. 不放回抽取

从一套充分洗匀的扑克中不放回地一直抽，直到抽到两张大小王都抽到，一共抽了多少张的期望值？

先说结论：

一般地，从 $n$ 张牌中不放回地抽取特定 $k$ 张，期望抽取次数是 $Exp = k \cdot (n+1)\ /\ (k+1)$

然后说证明：

首先，在充分洗匀的情况下，从牌堆顶部抽取和从牌堆中间任意位置抽取是等价的，所以可以假设牌堆是从左到右的展开一列，抽取动作等效于从左开始逐张提取卡片，直到所有目标牌都抽取到手里。想像 k 张目标牌已经排成一列，其顺序无关紧要因为最终都会被清空。这 k 张牌相隔包括左右外侧共产生 k+1 个空隙：

然后将 n-k 张非目标牌随机插入这 k+1 个空隙中，每个空隙插入的概率相等，所以期望张数是 $\frac{n-k}{k+1}$。

最后考察清点过程，显然最右侧外侧可以抛弃，所以需要抽取的期望张数是 $n - \frac{n-k}{k+1} = \frac{k \times (n+1)}{k+1}$。得证。

这类抽奖励在各类游戏中都有，其特点是：

• 奖池有限
• 不会抽到重复的奖品
• 通常玩家心仪的是奖池中的几个特定奖品，其它是无关紧要的

有些游戏会添加『升星级』的设置，某个抽卡人物 A 多次抽到后，会升级为 A+，A++ 等。这种情况完全可以视为若干个独立的 A 在奖池中，需要全部抽齐，计算时调整 n 与 k 的值即可。

—— 算下来『全部抽齐』确实是成本很高的目标。

2. 放回抽取

在数量为 n 的奖池中进行抽取，抽取结果是可重复的（抽完放回），抽完所有 n 件奖品的期望抽取次数是多少？

$$Exp= n \cdot (1 + 1/2 + 1/3 + \cdots + 1/n) = n \cdot \sum _{i=1}^n \frac{1}{i} = n \cdot H_n $$

$H_n$ 就是调和级数。

证明，从后往前考察抽取过程：

• 如果当前只剩最后一件奖品没抽到，则从 n 件奖品中抽取特定一件的几率是 1/n，期望抽取次数是 n。
• 如果剩最后两件奖品没抽到，只要从 n 件奖励中抽到这两件的任意一件，就会进入上一条的『只剩一件』状态。抽到这两件奖励之一的几率是 2/n，期望抽取次数是 n/2。
• 每抽到一件新的就会进入下一个状态。而在任意某个状态时，当前对应的期望抽取次数是总数 / 剩余新奖品数。
• 因此，总的期望抽取次数是所有状态的期望抽取次数之和。

调和级数是发散的，可以用 Excel 等工具计算期望值，也可以用近似公式 $H_n \approx \ln n + C + 1/2n$ 算。这里的 $C$ 就是欧拉常数，约为 0.5772156649。

进一步的，如果不需要抽完，而只是抽取其中的特定 k 件奖品，仍可以从最后一件倒数推算，其期望抽取次数是：

$$Exp = \frac{n}{1} + \frac{n}{2} + \cdots + \frac{n}{k-1} + \frac{n}{k} = n \cdot H_k$$

3. 放回抽取，不同概率

在各类游戏中，奖励物品往往有不同的『品质』等级，以及不同的出现概率。可以简单理解为

奖池分为 $m$ 个小奖池，每个小奖池进入概率分别为 $p_1, p_2, \cdots, p_m$，奖品数量分别为 $n_1, n_2, \cdots, n_m$。获得整个奖池全部奖励的期望抽取次数 $N$ 是多少？

对于每个小奖池，期望抽取次数是 $n_i \cdot H_{n_i}$，再除以小奖池本身的概率 $p_i$，即得 $N_i=n_i\cdot H_{n_i} \ /\ p_i$。但对于整个奖池，总期望抽取次数 $N$ 不是所有小奖池的次数之和，而是它们之中的最大值。因为当有 $n_i = p_i\cdot N$ 的抽取次数落到 $i$ 号小奖池时，也有对应的 $(1-p_i) N$ 溢出到其它小奖池。最大值保证溢出部分覆盖其它小奖池的期望抽取次数。随意举例如下：

总抽取次数 $N$ 会以 $0.01:0.09:0.3:0.6$ 的比例分配到四种品质中。观察表格，橙色传说所需的抽取次数最多，且溢出部分按比例仍可以覆盖其它所有奖池的有效抽取次数，也就是说抽到 2929 次时，期望可以让全部 10 件橙色传说抽取到手里。此时溢出的 2899 次已经完全覆盖了其它三种品质的有效期望抽取次数。所以 $$Exp = Max\left(n_i \cdot H_{n_i} \ /\ p_i\right)$$

进一步考虑发现，只要每个小奖池的期望有效抽取次数和设定的小奖池进入概率比例不一致，则一定会有一个小奖池的期望抽取次数最大。如果正好一致，也意味着每个都是最大值。这个结论是普适的。

4. 放回抽取，反还比例

在某些游戏中，抽取到重复奖品时，会通过『碎片』等给予一定的返还比例，若干个碎片可以合成需要的道具。先考虑单个奖池的情况：

假设奖池中有 n 件奖品，抽取概率相同。抽到重复奖品时，返还比例是 r(r<1)，那么抽取次数的期望值是多少？

没有傻逼会去合成已经存在的道具，所以返还可以视为抽取了 r 个新道具。同时由于 r <1，重复的价值仍是小于不重复的，因此在确保碎片可以合成奖池所有剩余奖品前，应当不合成任何奖励减少重复，使得整体收益最大化。再次从后往前考察抽取过程：

• 如果只差最后一件奖品没抽到，那么我仍需要从 n 件奖励中抽取，抽到这件奖励的几率是 1/n，期望抽取次数是 n。其中 1 次是有效抽取，n-1 次是重复抽取，返还 r(n-1) 碎片。
• 如果剩最后两件奖品没抽到，只要从 n 件奖励中抽到这两件的任意一件，就会进入『最后一件』状态，抽到这两件奖励的几率是 2/n，期望抽取次数是 n/2。其中 1 次是有效抽取，n/2-1 次是重复抽取，返还 r(n/2-1) 碎片。
• 以此类推，在 n 奖池中抽取到 t 个奖品时，对应的总抽取碎片为：

\begin{align*}
R(t) &= r\cdot(n-1) + r\cdot(\frac{n}{2}-1) + \cdots + r\cdot(\frac{n}{t-1}-1) +r\cdot(\frac{n}{t}-1) \\
&= r \cdot \sum_{i=1}^t \left( \frac{n}{i} - 1 \right) \\
&= r \cdot \left( nH_t - t \right) \\
\end{align*}

在某个节点 $t=k$ 时，获得的碎片就已经足够合成了所有剩余的奖品，即 $R(k) \geq n-k$。因为 $R$ 单调递增，$k$ 必然存在，此时对应的期望抽取次数即是 $n \cdot H_k$。$k$ 的通式求解超出了我的能力，但使用 Excel 或编程求解在有限的 n 范围内仍是方便的。

5. 放回抽取，反还比例，多个奖池

假设还是 $m$ 个分概率的小奖池，但返还的碎片是通用的，那么全部奖品的总抽取次数期望值是多少？

这里需要再添加一个参数，每个奖池的相对价格是不同的，设为 $c_i$。在 r 相同同，高价值奖池重复物品返回的碎片绝对数量更多。需要通过调整各奖池的 $c_i \cdot R(t)$，使各奖池的每抽价值归一化，后续就化为多个小奖池取最大值的已知问题了。即：

• 设总抽取次数为 $N$，则对于每个小奖池会分得 $k_i = N \cdot p_i$ 次。
• 各小奖池 $i$ 各抽取 $k_i$ 次，得到 $t_i$ 个奖品，满足 $n_i \cdot H_{t_i} \leq k_i$，虽然计算期望可以取等号，但调和级数无法依此简化计算。
• 此时产生的碎片为 $r\cdot(k_i-t_i)$ 即 $R(t_i)$，剩余未抽取的奖品数量为 $n_i - t_i$。
• 当总碎片数量大于等于剩余奖品所须数量时，N 得解。整理得到：

\begin{array}{l}
\left\{
\begin{aligned}
& k_i =N \times p_i, \\
&n_i \cdot H_{t_i} \leq k_i, \\
&\sum_{i=1}^m [R(t_i) \cdot c_i] \geq \sum_{i=1}^m [(n_i-t_i)\cdot c_i] \\
\end{aligned}
\right. \\
\\
\Rightarrow \large{ N = Max\left( k_i / p_i \right) }
\end{array}

• $c_i$ 为不同奖池间的价值权重，与碎片兑换各级别奖品的相对数量有关，与 $r$ 无关。
• $p_i c_i$ 为定义数据，调和级数只有离散解，每个小奖池计算的全局期望 $N_i$ 仍可能不一致。

6. 概率提升

在某些游戏中，当抽取未能获得奖励时，会提升下一次抽取的概率，直到抽取到奖励后重置概率。则真实概率是多少？

先从一个简单例子入手：

事件 A 的发生概率是 1/2，事件 B 发生概率为另 1/2，当 A 未发生时，下一次发生事件 A 的概率提到 100%。那么在多次重复中，A 事件的真实概率是多少？

直觉可能会认为是 75%。但进一步思考发现，事件 B 的下一轮必定跟随事件 A，而事件 A 的下一轮 A 与 B 各有 50% 几率发生。将多轮过程展开后会发现这是 $A:BA=1:1$ 的序列。因此单看事件 A 发生的真实概率为 2/3。

推广到 A 事件 1/3 的情况，可以得到：

在足够长的序列中，三者以 $A:BA:BBA=\frac{1}{3}:\frac{4}{9}: \frac{2}{9}$ 的比例出现，序列总长度为 $L = \frac{1}{3}\cdot n + \frac{4}{9}\cdot 2n+ \frac{2}{9}\cdot 3n$，其中 A 在每个组合中都出现一次，因此 A 的真实概率是 $\frac{n}{L} = \frac{9}{17}$。

由此观察到一般规律，对于初始概率为 $a$，提升步长为 $p$，达到 100% 时的步数为 $k$，即 $a+kp\geq 1 \Rightarrow k\geq(1-a)/p$ 时，真实概率 $P(A)$ 是以下式子的解：

\begin{align*}
L &= p(A)+2\cdot p(BA)+3\cdot p(BBA)+ \cdots + k \cdot p(\underbrace{BBB \ldots B}_{k-1}A) + (k+1) \cdot p(\underbrace{BBB \ldots B}_{k}A) \\
&= a + 2(1-a)(a+p) + 3(1-a)^2(a+2p) + \cdots + k(1-a)^{k-1}(a+(k-1)p) +(k+1)(1-a)^k\cdot1 \\
&= \sum_{i=1}^k \left[i\cdot(1-a)^{i-1}(a+(i-1)\cdot p) \right] + (k+1)(1-a)^k \cdot 1 \\
P(A) &= k/L
\end{align*}

这其实就是马尔克夫链与状态转移方程。

因为 $k$ 是有限的，所以式子不涉及无穷级数，但是计算仍然不方便，需要借助工具。

在 $a = p = 1/n$ 的特殊情况下，可以进一步简化为：

\begin{align*}
L &= \frac{1}{n} + \frac{n-1}{n}\cdot\frac{2}{n}\cdot2 + \frac{n-1}{n}\cdot\frac{n-2}{n}\cdot\frac{3}{n}\cdot3 + \cdots + \frac{n-1}{n}\cdot\frac{n-2}{n}\cdots\frac{1}{n}\cdot\frac{n}{n}\cdot n \\
& = \sum_{i=1}^n \left[\frac{i^2}{n^i} \cdot \prod_{j=1}^{i-1} (n-j) \right] \\
Q(n) &= 1/L \\
\end{align*}

手动计算几个较小值：
\begin{align*}
Q(2) &= \frac{1}{ \frac{1}{2} + \frac{2^2}{2^2}(2-1) } = \frac{2}{3} \\
Q(3) &= \frac{1}{ \frac{1}{3} + \frac{2^2}{3^2}(3-1) + \frac{3^2}{3^3}(3-1)(3-2) } = \frac{9}{17} \approx 0.5294 \\
Q(4) &= \frac{1}{ \frac{1}{4} + \frac{2^2}{4^2}(4-1) + \frac{3^2}{4^3}(4-1)(4-2) + \frac{4^2}{4^4}(4-1)(4-2)(4-3) } = 32/71\approx 0.4507 \\
Q(5) &= \frac{1}{ \frac{1}{5} + \frac{2^2}{5^2}\cdot4 + \frac{3^2}{5^3}\cdot4\cdot3 + \frac{4^2}{5^4}\cdot4\cdot3\cdot2 + \frac{5^2}{5^5}\cdot4\cdot3\cdot2\cdot1 } = 625/1569\approx 0.3983 \\
Q(6) &= \frac{1}{ \frac{1}{6} + \frac{2^2}{6^2}\cdot5 + \frac{3^2}{6^3}\cdot5\cdot4 + \frac{4^2}{6^4}\cdot5\cdot4\cdot3 + \frac{5^2}{6^5}\cdot5\cdot4\cdot3\cdot2 + \frac{6^2}{6^6}\cdot5\cdot4\cdot3\cdot2\cdot1 } = 324/899 \approx 0.3604 \\
Q(7) &= 117649/355081 \approx 0.3313, \\
Q(8) &= 131072/425331 \approx 0.3082, \\
Q(9) &= 4782969/16541017 \approx 0.2892, \\
Q(10)&= 1562500/5719087 \approx 0.2732,
\end{align*}

来个大的：

\begin{align*}
Q(100)&=\frac{\tiny{10587911840678754238354031258495524525642395019531250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000}}{\tiny{129277986730885202106151856642773382549665465071825090615034369359178917814747670780594211447306244001059786961886915926297133393516168977984471526892507817}} \\
&\approx 0.0819
\end{align*}

再来个图：

7. 保底

在某些游戏中，抽取到一定次数后，会有保底机制，即抽取次数超过一定值后，会直接给予某个奖励。

待续。