小名开开

在上一篇文章中，我提到了 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.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. 保底

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

待续。

PT 站大多使用了 nexusPHP 的方案，因此其 Seeding Bonus，通常称为魔力点数，基本也使用了相同的公式进行计算：

\begin{align*}
&A = \sum_i\left( 1-10^{ -\frac{T_i}{T_0} } \right) \cdot S_i \cdot \left( 1+\sqrt2 \cdot 10^{ -\frac{N_i-1}{N_0-1} } \right) \\
&B = B_0 \cdot \frac{2}{\pi}arctan\left( \frac{A}{L} \right) \\
\end{align*}

其中：

• 积分每小时计算并发放一次，后文不再赘述。
• $S_i$ 为第 i 个种子的大小，单位是 GB
• $T_i$ 为第 i 个种子发布起到现在经过的时间，单位为周，$T_0$ 为参数，$T_0 = 4$
• $N_i$ 为第 i 个种子当前的做种人数，$N_0$ 为参数，$N_0 = 7$
• $L$ 为参数，$L = 300$，$B_0$ 为参数，$B_0 = 50$，不同站点有时会在 $B_0$ 和 $L$ 上有所调整。
• $B$ 为每小时用户获得的点数，由于 $arctan$ 函数值有上限，因此 $B$ 恒小于 $B_0$

可以看到，对于 A，实质上是对每个种子的体积进行了加权求和，两个系数分别与时间和做种人数有关，实质上还是某种『有效做种体积』。因此先分析第二段公式，使用 Geogebra 绘制函数 $\frac{2}{\pi} \cdot arctan\left( \frac{A}{L} \right)$ 图像并取若干个样点：

函数图像是个有上限、增幅递减的增函数，在有效做种体积为 200GB 时，可以得到 $0.37B_0$ 的魔力值，在 1000GB 时为 $0.81B_0$。以下列举一些计算结果：

如果要达到特定的目标上限比，几个档位如下：

200 GB 以下大致可视为线性增长，2 TB 时可以达到 90% 上限分数，更高的分值需要极大的硬盘空间，逐渐不甚合理了。$L$ 参数的变化会等比例地影响上表中的加权体积数值。

下面来研究两个系数的影响，这两个系数实现了『实际做种体积』$\rightarrow$『有效做种体积』之间的换算。首先是时间系数，$f(T_i) = \left(1-10^{ -\frac{T_i}{T_0}}\right)$，由于 nexusPHP 本身对该分值是基于小数的，且『周』的展示不甚直观，这里直接改用『天』作为单位，换算 $T_0 = 28$，其图像为：

显然这也是和 arctan 一样的有上限、增幅递减的增函数，但增幅比 arctan 要大得多，靠近上限十分容易，也就是说到达 0.999 或更高值是个可行的目标。以下列举一些计算结果：

因为是 10 的幂运算，计算分段明显。每 28 天小数点后多一位 9，在寻找 free 种子攒积分时，选时间大于三个月 (>0.999) 的就行，『实际做种体积』与『有效做种体积』几乎一致。

最后是做种人数系数，$f(N_i) = \left(1+\sqrt2 \cdot 10^{ -\frac{N_i-1}{N_0-1}}\right)$，在其定义域内这是个有上下限的减函数，其图像为：

与前两个系数不同，$N_i$ 做种人数只能取正整数，因此其图像是离散的。以下列举一些计算结果：

可以看到，只有 8 人以下的种子才有可能达到 >1.1 的系数，而 20 人以上的种子，该系数的增益已经完全可以忽略了。由于做种人数是不可控的，因此可以简单划为两类，20 人以下有加成， 20 人以上无。又因为人数会变动，因此若要从这个系数上获得稳定收益，必须考虑 $N_i<8$ 的种子，以容许一些浮动范围。

通过以上两个系数的分析，可以得到一个结论，我称为『三八原则』：

三八原则：在 PT 站做种赚取积分时，选择时间大于三个月，做种人数小于 8 人的种子。

原则展开如下：

• 做种人数足够少，可以使实际 1GB 的硬盘空间达成 2.4GB 的做种效果，但条件苛刻。
• 但再多的做种人数也不过是系数为 1，没有奖励，但不会有折损。
• 三个月以前的种子，做种效果几乎 1:1 且容易达成。而一个月以内则折损可观，不是没得选尽量避免保新种。
• 积分有上限，有效做种体积 1-2 TB 已经可以达到 80%-90% 的上限了。再增长 10% 到 99% 则需要近 20TB 有效保种，即使都是长寿孤种也需要 8TB，意义不大。
• 形势所迫必须新种，则放置一周有 44% 的有效做种体积，准备 4TB 硬盘也就够了。
• 孤种才有奖励，新种只有折损。

这样可以最大化体积系数，从而用相同的硬盘空间获取更大的做种收益。当然，既然是赚取积分，肯定是要选择 free 种子。种子时间可以挑选，但做种人数加成八成是不必指望了。

其它一些 Bonus：

各站点会根据自己的实际调整参数，或添加其它的奖励机制，例如：

1. 每个种子 0.7 基本分，最多 14 个。

   该数值上限 9.8，是直接作为基础值的，因此该站的积分上限为 $9.8+B_0$。这个条件几乎不存在代价，保种 1-2 TB 通常来说都会超过 14 个种子。

2. 用户等级奖励，每升一级 1% 的比例奖励。

   这个奖励是在基础值之外的奖励，计算公式为 $(1+1\%\times Level) * (9.8+B_0)$

   用户等级奖励其实存在很大的隐性的代价。PT 站往往要求上传量远大于下载量才能升级，而 1% 的奖励却聊胜于无。这就像在银行存一百万拿 2% 的利息，银行要求你再存一千万升 VIP，而 VIP 的待遇仅仅是利息多 1%，并不是 $2\%\rightarrow 3\%$ 的 1%，而是 $2\%\rightarrow 2.02\%$ 的 1%。如果上传本有余裕，自然升级也就罢了，刻意去为 1% 奖励而花大力气冲没什么意义。

3. 为官方种子做种的额外奖励，上限为 25%。

   百分比类型奖励都是平行的，也就是 $(1+1\%\times Level + 25\%) * (9.8+B_0)$。
   注：该站点程序有 bug，多个官种奖励叠加反而会导致总 bonus 下降，因此只保一个 +25% 官种即可。

4. 每 100 个种子做种 24 小时以上开始 5% 的比例奖励，无上限，每个种子需要大于 5MB。

   $$(1+1\%\times Level + 25\% + \mathop{SeedsQuantity}\limits_{\substack{seedsize>5\text{MB} \\ seedinghour>24}} / 100 \times 5\%) * (9.8+B_0)$$

   该规则鼓励大量小体积种子做种，理论上来说，以 2TB 的目标计算，若全部靠 5MB 的种子达成，至多可以使用 40 万个种子，单这项 bonus 就能达到 $200B_0$，这是一个巨大的增益。客观上受限于种子总量有限及大小各异，不可能达到这个数量，但成为一个万种大佬还是有可能的。此时的 bonus 为 $$(1+25\%+\frac{10000}{100}*5\%) \times (9.8+0.9\times B_0) = 6.25 \times (9.8+0.9\times B_0)$$ 若 $B_0=50$，则值为 $342.5$。在该规则下，保种策略是尽量地将大种子更替为小种子，在总量有限而不变的情况下，提升种子数量。

   但是我们需要进一步讨论。更替小种虽有额外收益但通常也会产生额外的下载流量，需要靠魔力值兑换上传量来冲抵。一个种子的小时增益为 $5\% / 100 \times (9.8 + 0.9B_0) = 0.0274$，每天可以通过魔力值兑换的上传流量为 $0.0274 \times 24 \; / \; 425 \times 1024 = 1.5844$ MB，一个 5MB 的种子在 3.15 天回本，因为是独立计算，一百个 5MB 的种子也是 3.15 天回本。回本时间不受种子数量的影响，但与种子体积成正比，当种子体积到 50MB 时回本时间就超过一个月了。

   另一需要考虑的软性限制是 PT 与站点服务器的通讯。当做种数量足够多时，PT 客户端的 tracker update 操作变得十分频繁，甚至变相成为一种对服务器的 DOS 攻击。这也是热门大站往往 tracker update interval 很长的原因。大量种子对保种用户本身也是存在风险的，对同一地址的短时间频繁请求，在客户端到服务器的整条网络传输路径上的任何一个环节都有可能触发限制。服务器收不到 update 信息，用户就反而无法获得积分。

5. 使用 $250\times B_0$ 的价格购买 5% 比例的奖励徽章，可购买多个。

   同样是某种需要考虑投入产出的增益。该站点没有其它基础值奖励，以每小时 $0.9B_0$ 计算，购买后的增益为 $0.045B_0$，回本时间为 $250B_0 /0.045B_0 / 24 = 231.49$ 天。该站 50G 上传量为 $50B_0$，一个徽章相当于 250 GB 的上传量。从无穷长远来看，一切本金不会湮灭且无通胀环境下的投入都是值得的，算来第 232 天开始可就是净收益了。但是在合理的生命周期内，这个投入产出比是十分低的。因此这个徽章是一个『炫耀』的物品，而不是一个『赚钱』的物品。

   PT 站的生命周期是客观存在的，对于用户方，如果生活中发生了某些变化，帐号有时就冷落消失了。对于站点方，也有可能因为各种原因关闭跑路。长周期的收益更多只是在纸面上的。

   相聚离开都有时候，没有什么会永垂不朽。

6. 后宫加成，每个下线 5%。

   这个奖励的 ROI 计算更为复杂，暂不清楚是否有复利加成，即我后宫的后宫是不是我的后宫。单从实践经验来看，$50B_0$ 兑换邀请和赠与新人 $62.5B_0$ 启动资金，换得的是下线的 5% 增益。如果下线足够给力，则可以得到和徽章奖励一致的 $0.045B_0$。回本周期为 105 天。考虑到下线需要积累一段时间才能达到 $0.9B_0$，实际周期在 150 天左右。比起奖章仍合算一些。如果不提供启动资金，回本周期为 47 天，但新人启动时间会长不少，相当于额外再增加一个 0-50G 的启动时间。

   如果考虑风险因素，ROI 就不好说了。下线的行为是不可控的，且会影响到自己的账号。$112.5B_0$ 的本金浪费是一方面，其它权利影响是则另一方面。除非有站点相关数据作为支撑，否则这些风险很难量化。

7. 从奖池中抽道具，部分道具增加 10% 基础魔力值奖励，每次抽取动作消耗 $2.5B_0$。

   从 n 个道具中抽取特定的 k 个道具，假设无重复时，全部抽齐的期望次数为 $\frac{k(n+1)}{k+1}$，则成本为 $\frac{k(n+1)}{k+1} \times 2.5B_0$，收益为 $k \times 10\% \times 0.9B_0$，回本周期为 $1.1574 \cdot \frac{n+1}{k+1}$ 天，取决于奖池 n 和 k 的比例。事实上哪怕是 1/100 的比例，也比上一条后宫加成更合算。

   进一步考虑，第一个道具在抽出后开始生效，其成本为 $\frac{n+1}{k+1}\times 2.5B_0$，收益为 $10\% \times 0.9B_0$，回本周期同样为 $1.1574 \cdot \frac{n+1}{k+1}$ 天。也就是说，在整个过程中，每个道具的奖励比例一致，则回本周期也都是一致的。

   当有重复时，情况则有不同，第一个的期望成本为 $n/k\times 2.5B_0$，回本周期为 $1.1574 n/k$ 天。所有奖品抽到的期望次数为 $n\cdot H_k$，其中 $H_k=\sum_{i=1}^n\frac{1}{i}$ 为调和级数。越往后面，即使奖励相同，但抽取成本变大也会导致回本周期变长。最后一个的成本是 $n\times 2.5B_0$，回本周期对应为 $1.1574 n$ 天。

   如果除 k 外的其它道具也有一定的奖励，或者总数 n 难以计算，则 ROI 的计算只能通过少量抽样下的贝叶斯推断来进行了。

8. 上传流量的奖励，每 1GB 上传流量奖励 0.1 魔力值。

   PT 站本身上传量就是货币，这个其实属于某种重复奖励，且没有可优化或者可选择的余地。不讨论。

末尾一些啰嗦：

$f(A) = B_0 \cdot \frac{2}{\pi} \cdot arctan\left( \frac{A}{L} \right)$ 的导数为 $f^{'}(A) = \frac{2B_0}{\pi L} \cdot \frac{1}{1+\left( \frac{A}{L} \right)^2}$，为 $\frac{1}{1+x^2}$ 的缩放。当前每增加 1GB 的有效做种体积对应多少魔力值增量，可以通过导数来计算。

通过计算 $Total\_Bonus = N_i \cdot \left(\sqrt2 \cdot 10^{ -\frac{N_i-1}{N_0-1}}\right)$ 发现，在 3 人及以下做种时，PT 站额外奖励了更多的体积系数分到所有的人，而人数更多时 Bonus 总增益开始下降，奖励总量反而逐渐减少。算法并非是将固定的 Bonus 分发给所有人，而确实是针对孤种的单独加成：