小名开开

在春天的光影里喘息

当我们摊开中国历史的长卷,目光总是不自觉地被那些醒目的坐标所吸引:王朝更迭、礼乐兴衰、苍生聚散、英雄沉浮。这些被墨迹重重勾勒的事件,构成了熟悉的历史叙事。

但如果俯下身来,把耳朵贴近那些泛黄的纸页之间,我们或许能听见另一种声音。那是运河上漕船划开晨雾的水声,是商帮会馆里算盘的轻响,是盐运码头纤夫们沉重的喘息,是无数名不登史的商人在暗夜中传递消息时的低语。这些声音太过细微,往往被时间的喧嚣所淹没,但从未真正消失。它们沿着某条隐秘的河道流淌,在某座商业城市的园林里盘旋,在某座宗祠的牌位间低徊,最终汇成一条看不见的暗流,与那条明面上的王朝更替线并行、交织、碰撞。

这条暗流,或许就是我们理解中国历史另一面的钥匙。

黄巢

待到秋来九月八,我花开后百花杀。冲天香阵透长安,满城尽带黄金甲。

在唐末乱世的历史图景中,黄巢是浓墨重彩的一笔。他起于盐商私贩,累举不第后,聚起数十万大军,攻入长安,自立为帝,建立了短暂的大齐政权。他对世家大族、门阀勋贵进行了毁灭性打击,是唐末走向五代十国割据局面的起点,也是世家落幕、士大夫阶层兴起的转折点。

但你有没有疑惑过,为什么他明明是私盐贩子,其起义却被定性为农民起义?

据史载:

唐末关东大旱,岁大饥,人相食。而地方官府催征愈急,迫使流民四起。私盐王仙芝率先起兵濮州,举『天补平均大将军』旗号,流民纷纷投奔,一时声势大振。

黄巢与王仙芝同乡,同为私盐贩首,本就是官府缉捕要犯。王起兵后,官府清剿日甚,黄巢进退维谷:若不奋起反抗,迟早家破人亡。眼见天下大乱、灾黎云集,他最终决意起兵响应。

起义初期,黄巢与王仙芝合兵一处,转战山东河南,沿途饥氓争赴,队伍迅速壮大。王仙芝战死后,黄巢收编其众,成为义军首领。此后他避实击虚,挥师南下,在岭南休整补充。次年北伐,势如破竹,于 880 年冬攻入长安,登基称帝,国号大齐,改元『金统』。

然黄巢军始终四方转战,未曾有稳固据地。后唐军纠合沙陀骑兵,合围进剿,大齐势渐蹙。中和三年,巢弃长安东撤,与唐军相持于陈州,终不支,败走狼虎谷。次年六月,自刎殉难,起义遂告覆亡。

从史料上可以得知:

其一,黄巢的起义军主体是数十万破产农民、流民与饥民,这些底层百姓构成了起义的核心战斗力,是运动得以燎原的根基。

其二,起义的核心矛盾直指唐末腐朽统治阶级。打击的是横征暴敛的朝廷、贪婪的贪官污吏与垄断权力的门阀贵族,依靠的是被压迫的底层百姓,深刻契合了农民阶级反抗封建压迫、追求生存权利的核心诉求。

其三,起义脉络始终贯穿「生存」与「复仇」的主题。黄巢最初的个人诉求虽源于摆脱身份困境,但在百万流民意志的裹挟与推动下,这些个人意志被升华为整个阶级的反抗意志,使他最终成为了农民阶级追求生存权利的集中代表。

所以,从性质上分析,这是一场农民起义的结论,确无争议。但如果只关注结论,我们会忽略关于黄巢本人的一些细节:

黄巢自幼读书习武、屡次赴长安应考进士,即使是在起事前半年仍然再次赴试。

黄巢「世鬻盐,富于赀」,拥有稳定的宗族亲信、贩盐徒党与护卫力量。

黄巢大军攻入长安后,虽然烧杀有之,但对「百工」往往采取收编而非单纯杀戮的态度。

黄巢的起义是无路可走后的绝地反击,是条件不成熟时的无奈之举,而非预设好的革命蓝图。

盐与武装

在唐代,盐铁专营是国家重要经济来源之一,从事私盐活动会遭到官府围剿。贩私盐一石以上即处死刑,持械拒捕者一律诛杀。因此,从事私盐活动必然与产盐区的基层社会组织紧密结合,并组织武装力量对抗巡盐捕盗与同行竞争。黄巢极有可能是通过控制某些盐业生产点,比如私自开辟的煎盐灶,或者收买官府监管下的盐亭来获取货源。如果只是一个二道贩子,他很难养活一支武装队伍。

而为了维持低成本和高产出,他必须介入生产环节。唐代的盐主要分池盐和海盐两种。黄巢老家在山东菏泽,靠近山东沿海和河北沿海的盐产区。当时的海盐生产主要是「煎炼」。这需要掌握卤水的浓度配比、熬制火候以及蒸发技术,并非简单的体力活,而是需要经验积累的手艺。私盐贩要获得高质量的盐,少时可以从盐户手中收购,但规模扩大后必须自己组织生产,必须有一批懂技术的灶户。这也是他注重百工的原因之一。

从这些看来,黄巢是一个合格的资产阶级生产组织者。他通过私人资本介入组织、生产、销售全链条。他能认知到工匠工艺的价值,同时建立武装力量保护商业利益,这些都是典型的资本主义商业活动特征。

求生路径

黄巢自幼读书习武、屡次赴长安应考进士。他能力有裕,但出身却被体制排斥。作为一个私盐商人,科举是他「洗白上岸」的唯一途径。通过科举进入官场获得政治地位,为家族的商业活动提供合法保护,是中国古代商人阶层典型的求生路径。

他所处的时代尚未形成完整的资产阶级群体,对自身所处阶级和利益诉求没有明确认知,没有任何资产阶级革命的路线与经验可供参考。黄巢及其家族的诉求,始终是依附皇权、洗白上岸,而非建立限制皇权、保护商人利益的新制度。在进入长安后,黄巢也曾下令惩办贪官污吏、救济贫民,但并未提出任何新的制度蓝图。没有保护商业贸易,没有建立新的社会秩序,只是延续了古代改朝换代的老路。

在早期盐铁专营的历史背景下,私盐的超额利润催熟了黄巢这类商人阶层。他们在经济上达到了资本主义的早期阶段,但因为是催熟阶层人数不足,又天然处于士农工商的底层,在政治上无法形成独立的阶级意识与政治诉求,也无法跳出皇权专制的循环。

暗线之始

然而要注意的是,黄巢及其后的一系列中国早期资产阶级在历史上进行过的尝试,与近代欧美资产阶级的早期发展逻辑高度相似,但具体的行为和后果则因为历史环境的不同而有所差异。任何在旧体制中崛起的新兴富裕阶层,最初的诉求都不是推翻现有秩序,而是试图通过合法渠道进入权力体系,以政治身份保护经济利益,在既有社会框架内获得承认与安全。

黄巢效仿耕读传家正统路径,企图通过科举跻身官场,从而庇护整个私盐集团的行为。与欧美资产阶级早期试图在封建王权中谋求议会席位、法律保障的诉求完全相同:他们首先想做秩序的参与者,而非颠覆者。黄巢最终走向起义,是新兴阶层的上升通道被彻底堵死后的被迫反抗,这同样与西方资产阶级的抗争逻辑一脉相承。

一个新兴阶层从产生到形成独立的阶级意识,再到发展出独立的政治诉求,开始与传统秩序进行斗争,从而塑造新的社会形态,这是一条相对固定的历史路径。它需要在不断的尝试和挫折中逐渐成熟,最终形成对现有秩序的根本挑战。这是一个漫长的过程,充满了曲折和反复,埋伏在历史的暗线上细长而悠远。在这条暗线上,我们还要继续追寻那些被历史尘埃掩盖的身影。我们将看到更多的他们,逐渐成长,形成自己的诉求,做出自己的选择,并对历史进程产生更深远的影响。

欧美成熟资产阶级的成功在于,他们在体制内晋升受阻后,能依托城市经济与阶级力量,走向议会斗争、制度革命,通过「换国王、定规则」掌握国家权力;而黄巢只能依附农民的力量,沿用农民起义的老路。他试图亲自下场夺权,但商业组织的逻辑与治国所需的权力平衡之间差距太远,最终只能重蹈改朝换代的老路。

但这次失败带来了一个深刻的教训。此后的新兴商业阶层不再试图取代皇权。他们学会了更隐蔽的方式:渗透、收买和包围皇权。

阅读全文 »

在上一篇文章中,我提到了 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 引擎。

chrome 控制台复现

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个整数的因子序数积的分布

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

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

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

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

    29000到30000整数的因子序数积的分布

    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

steadyaim60
archery60
details60

视频中非常明确地显示了,在原版 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 循环节:

1/7 Circle

不过『走马灯数』这个概念还是有点意思的。它最显著的特征是循环节长度 = 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++;
}

结果出乎意料地多:

1/n 循环节
1/7 142857
1/17 0588235294117647
1/19 052631578947368421
1/23 0434782608695652173913
1/29 0344827586206896551724137931
1/47 0212765957446808510638297872340425531914893617
1/59 0169491525423728813559322033898305084745762711864406779661
1/61 016393442622950819672131147540983606557377049180327868852459
1/97 010309278350515463917525773195876288659793814432989690721649484536082474226804123711340206185567
1/109 009174311926605504587155963302752293577981651376146788990825688073394495412844036697247706422018348623853211
1/113 0088495575221238938053097345132743362831858407079646017699115044247787610619469026548672566371681415929203539823
1/131 0076335877862595419847328244274809160305343511450381679389312977099236641221374045801526717557251908396946564885496183206106870229
1/149 0067114093959731543624161073825503355704697986577181208053691275167785234899328859060402684563758389261744966442953020134228187919463087248322147651
1/167 0059880239520958083832335329341317365269461077844311377245508982035928143712574850299401197604790419161676646706586826347305389221556886227544910179640718562874251497
1/179 0055865921787709497206703910614525139664804469273743016759776536312849162011173184357541899441340782122905027932960893854748603351955307262569832402234636871508379888268156424581
1/181 005524861878453038674033149171270718232044198895027624309392265193370165745856353591160220994475138121546961325966850828729281767955801104972375690607734806629834254143646408839779
1/193 005181347150259067357512953367875647668393782383419689119170984455958549222797927461139896373056994818652849740932642487046632124352331606217616580310880829015544041450777202072538860103626943
1/223 004484304932735426008968609865470852017937219730941704035874439461883408071748878923766816143497757847533632286995515695067264573991031390134529147982062780269058295964125560538116591928251121076233183856502242152466367713
1/229 004366812227074235807860262008733624454148471615720524017467248908296943231441048034934497816593886462882096069868995633187772925764192139737991266375545851528384279475982532751091703056768558951965065502183406113537117903930131
1/233 0042918454935622317596566523605150214592274678111587982832618025751072961373390557939914163090128755364806866952789699570815450643776824034334763948497854077253218884120171673819742489270386266094420600858369098712446351931330472103
1/257 0038910505836575875486381322957198443579766536964980544747081712062256809338521400778210116731517509727626459143968871595330739299610894941634241245136186770428015564202334630350194552529182879377431906614785992217898832684824902723735408560311284046692607
1/263 0038022813688212927756653992395437262357414448669201520912547528517110266159695817490494296577946768060836501901140684410646387832699619771863117870722433460076045627376425855513307984790874524714828897338403041825095057034220532319391634980988593155893536121673
1/269 0037174721189591078066914498141263940520446096654275092936802973977695167286245353159851301115241635687732342007434944237918215613382899628252788104089219330855018587360594795539033457249070631970260223048327137546468401486988847583643122676579925650557620817843866171
1/313 003194888178913738019169329073482428115015974440894568690095846645367412140575079872204472843450479233226837060702875399361022364217252396166134185303514376996805111821086261980830670926517571884984025559105431309904153354632587859424920127795527156549520766773162939297124600638977635782747603833865814696485623
1/337 002967359050445103857566765578635014836795252225519287833827893175074183976261127596439169139465875370919881305637982195845697329376854599406528189910979228486646884272997032640949554896142433234421364985163204747774480712166172106824925816023738872403560830860534124629080118694362017804154302670623145400593471810089020771513353115727
1/367 002724795640326975476839237057220708446866485013623978201634877384196185286103542234332425068119891008174386920980926430517711171662125340599455040871934604904632152588555858310626702997275204359673024523160762942779291553133514986376021798365122615803814713896457765667574931880108991825613079019073569482288828337874659400544959128065395095367847411444141689373297
1/379 002638522427440633245382585751978891820580474934036939313984168865435356200527704485488126649076517150395778364116094986807387862796833773087071240105540897097625329815303430079155672823218997361477572559366754617414248021108179419525065963060686015831134564643799472295514511873350923482849604221635883905013192612137203166226912928759894459102902374670184696569920844327176781
1/383 0026109660574412532637075718015665796344647519582245430809399477806788511749347258485639686684073107049608355091383812010443864229765013054830287206266318537859007832898172323759791122715404699738903394255874673629242819843342036553524804177545691906005221932114882506527415143603133159268929503916449086161879895561357702349869451697127937336814621409921671018276762402088772845953
1/389 0025706940874035989717223650385604113110539845758354755784061696658097686375321336760925449871465295629820051413881748071979434447300771208226221079691516709511568123393316195372750642673521850899742930591259640102827763496143958868894601542416452442159383033419023136246786632390745501285347043701799485861182519280205655526992287917737789203084832904884318766066838046272493573264781491
1/419 0023866348448687350835322195704057279236276849642004773269689737470167064439140811455847255369928400954653937947494033412887828162291169451073985680190930787589498806682577565632458233890214797136038186157517899761336515513126491646778042959427207637231503579952267303102625298329355608591885441527446300715990453460620525059665871121718377088305489260143198090692124105011933174224343675417661097852028639618138424821
1/433 002309468822170900692840646651270207852193995381062355658198614318706697459584295612009237875288683602771362586605080831408775981524249422632794457274826789838337182448036951501154734411085450346420323325635103926096997690531177829099307159353348729792147806004618937644341801385681293302540415704387990762124711316397228637413394919168591224018475750577367205542725173210161662817551963048498845265588914549653579676674364896073903
1/461 0021691973969631236442516268980477223427331887201735357917570498915401301518438177874186550976138828633405639913232104121475054229934924078091106290672451193058568329718004338394793926247288503253796095444685466377440347071583514099783080260303687635574837310195227765726681127982646420824295010845986984815618221258134490238611713665943600867678958785249457700650759219088937093275488069414316702819956616052060737527114967462039045553145336225596529284164859
1/487 002053388090349075975359342915811088295687885010266940451745379876796714579055441478439425051334702258726899383983572895277207392197125256673511293634496919917864476386036960985626283367556468172484599589322381930184804928131416837782340862422997946611909650924024640657084188911704312114989733059548254620123203285420944558521560574948665297741273100616016427104722792607802874743326488706365503080082135523613963039014373716632443531827515400410677618069815195071868583162217659137577
1/491 0020366598778004073319755600814663951120162932790224032586558044806517311608961303462321792260692464358452138492871690427698574338085539714867617107942973523421588594704684317718940936863543788187372708757637474541751527494908350305498981670061099796334012219959266802443991853360488798370672097759674134419551934826883910386965376782077393075356415478615071283095723014256619144602851323828920570264765784114052953156822810590631364562118126272912423625254582484725050916496945010183299389
1/499 002004008016032064128256513026052104208416833667334669338677354709418837675350701402805611222444889779559118236472945891783567134268537074148296593186372745490981963927855711422845691382765531062124248496993987975951903807615230460921843687374749498997995991983967935871743486973947895791583166332665330661322645290581162324649298597194388777555110220440881763527054108216432865731462925851703406813627254509018036072144288577154308617234468937875751503006012024048096192384769539078156312625250501

简化代码并扩大计算量搜索后,得到 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);
}
}
};

100 万以内的走马灯数分布

有点意思。


研究原因

对于任意由 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 以内的走马灯数的数量,结果如下:

进制 走马灯数数量 进制 走马灯数数量 进制 走马灯数数量 进制 走马灯数数量 进制 走马灯数数量
21 447 41 460 61 457 81 1
2 470 22 472 42 476 62 482 82 455
3 477 23 459 43 490 63 457 83 453
4 0 24 461 44 460 64 0 84 466
5 493 25 1 45 483 65 477 85 468
6 470 26 455 46 445 66 469 86 464
7 466 27 278 47 477 67 471 87 478
8 277 28 451 48 468 68 459 88 454
9 1 29 465 49 1 69 461 89 470
10 467 30 488 50 471 70 473 90 478
11 444 31 469 51 471 71 471 91 453
12 459 32 372 52 460 72 470 92 464
13 458 33 474 53 470 73 465 93 446
14 459 34 473 54 471 74 447 94 485
15 454 35 461 55 486 75 471 95 462
16 0 36 0 56 461 76 467 96 464
17 451 37 466 57 451 77 449 97 471
18 472 38 466 58 478 78 451 98 459
19 473 39 489 59 465 79 466 99 476
20 487 40 458 60 455 80 490 100 0

这里最有特点的是当进制为偶数的完全平方时,没有走马灯数,当进制为奇数的完全平方时,只有唯一一个走马灯数 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 -T2 -1 771.511 2574731933 89.66% 108.84%
zstd -T2 860.439 2365713277 100.00% 100.00% 默认压缩率档位 3
zstd -T2 -10 1274.784 2190541240 148.16% 92.60%
zstd -T2 -19 4416.055 2009856914 513.23% 84.96%
xz -T2 -0 1798.904 2293838248 209.07% 96.96%
xz -T2 4475.497 1924084292 520.14% 81.33% 默认压缩率档位 6
xz -T2 -9 5041.978 1803774144 585.98% 76.25%
xz -T2 -9e 8380.731 1801292000 974.01% 76.14% e: 消耗额外时间进一步提高压缩率
7z -mmt=2 -mx=1 1896.414 2330130495 220.40% 98.50%
7z -mmt=2 7938.620 1870578044 922.62% 79.07% 默认压缩率档位 5
7z -mmt=2 -mx=9 13600.655 1788204019 1580.66% 75.59%
lz4 -fast 606.257 3529871893 70.46% 149.21% lz4 没有多线程参数
lz4 640.828 3458840786 74.48% 146.21%
lz4 -best 2613.739 3038936614 303.77% 128.46%

测试结果

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

整理结果:

Version Workstation Player Version Fusion VMRC
e.x.p x86
3.0.2 x86
12.0.0 Windows Linux Windows Linux 8.0.0 x86
12.0.1 Windows Linux Windows Linux 8.0.1 x86
8.0.2 x86
12.1.0 Windows Linux Windows Linux 8.1.0 x86
12.1.1 Windows Linux Windows Linux 8.1.1 x86
12.5.0 Windows Linux Windows Linux 8.5.0 x86
12.5.1 Windows Linux Windows Linux 8.5.1 x86
12.5.2 Windows Linux Windows Linux 8.5.2 x86
12.5.3 Windows Linux Windows Linux 8.5.3 x86
12.5.4 Windows Linux Windows Linux 8.5.4 x86
12.5.5 Windows Linux Windows Linux 8.5.5 x86
12.5.6 Windows Linux Windows Linux 8.5.6 x86
12.5.7 Windows Linux Windows Linux 8.5.7 x86
12.5.8 Windows Linux Windows Linux 8.5.8 x86
12.5.9 Windows Linux Windows Linux 8.5.9 x86
8.5.10 x86
14.0.0 Windows Linux Windows Linux 10.0.0 x86
10.0.1 x86
14.1.0 Windows Linux Windows Linux 10.1.0 x86
14.1.1 Windows Linux Windows Linux 10.1.1 x86
14.1.2 Windows Linux Windows Linux 10.1.2 x86
14.1.3 Windows Linux Windows Linux 10.1.3 x86
14.1.4 Windows Linux Windows Linux 10.1.4 x86
14.1.5 Windows Linux Windows Linux 10.1.5 x86
14.1.6 Windows Linux Windows Linux 10.1.6 x86
14.1.7 Windows Linux Windows Linux
14.1.8 Windows Windows
15.0.0 Windows Linux Windows Linux 11.0.0 x86
15.0.1 Windows Linux Windows Linux 11.0.1 x86
15.0.2 Windows Linux Windows Linux 11.0.2 x86
15.0.3 Windows Linux Windows Linux 11.0.3 x86
15.0.4 Windows Linux Windows Linux
15.1.0 Windows Linux Windows Linux 11.1.0 x86
15.5.0 Windows Linux Windows Linux 11.5.0 x86
15.5.1 Windows Linux Windows Linux 11.5.1 x86
15.5.2 Windows Linux Windows Linux 11.5.2 x86
11.5.3 x86
15.5.5 Windows Linux Windows Linux 11.5.5 x86
15.5.6 Windows Linux Windows Linux 11.5.6 x86
15.5.7 Windows Linux Windows Linux 11.5.7 x86
16.0.0 Windows Linux Windows Linux
12.0.1 Windows Linux
12.0.2 Windows Linux MacOS
12.0.3 Windows Linux MacOS
12.0.4 Windows Linux MacOS
12.0.5 Windows Linux MacOS
16.1.0 Windows Linux Windows Linux 12.1.0 x86
16.1.1 Windows Linux Windows Linux 12.1.1 x86
16.1.2 Windows Linux Windows Linux 12.1.2 x86
16.2.0 Windows Linux Windows Linux 12.2.0 x86 arm64
16.2.1 Windows Linux Windows Linux 12.2.1 x86 arm64
16.2.2 Windows Windows
16.2.3 Windows Linux Windows Linux 12.2.3 x86 arm64
16.2.4 Windows Linux Windows Linux 12.2.4 x86 arm64
16.2.5 Windows Linux Windows Linux 12.2.5 x86 arm64
17.0.0 Windows Linux Windows Linux 13.0.0 x86 arm64 universal
17.0.1 Windows Linux Windows Linux 13.0.1 universal
17.0.2 Windows Linux Windows Linux 13.0.2 universal
17.5.0 Windows Linux Windows Linux 13.5.0 universal
17.5.1 Windows Linux Windows Linux 13.5.1 universal
17.5.2 Windows Linux Windows Linux 13.5.2 universal
17.6.0 Windows Linux Windows Linux 13.6.0 universal
17.6.1 Windows Linux Windows Linux 13.6.1 universal

质数,也叫素数,是对 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;
}
  1. 同上,所有 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 等多个值。代码变得复杂,而优化效果一般。过于复杂的代码会让你不幸福。

  1. 判断素数时,只需要检查到 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 对本文亦有贡献。

0%