小名开开

在春天的光影里喘息

※ 本文有 AI 参与写作:Kimi M3.1 Flash 根据作者大纲构建了本文第一版,Deeepseek V4.1 Flash 和 GLM-5.3-Flash 帮助改写和清理逻辑。但三个模型都无法很好地处理话题转换时的过渡小段,且在文章过长时不再遵守 AGENTS.md 和 skill.md。

在 《葬送的芙莉莲》 里,精灵大法师芙莉莲有一手鉴定魔法,看一眼就知道眼前的箱子是不是伪装的宝箱怪。故事中,这手魔法的准确率高达 99%。

然后她一次又一次地开,一次又一次地被吞。

这个设计好笑的地方在于:99% 已经是一个高到离谱的数字,现实世界里的检测大多不敢吹到 99%。而动画里芙莉莲的失败率看着比设定的 1% 高得多,多到让人怀疑这个 99% 是不是抄错了行。

先别管动画抄没抄错。有意思的是「99%」这个数字本身,到底是什么。

一、99%,几乎什么都不是

无意义的 99%

假设漫长的时光里,大部分真宝箱早被别的冒险者探索过了。荒野里剩下的每 100 只宝箱中,只有 1 只是真宝箱。那么倘若这所谓「鉴定魔法」只是个假魔法,永远一律判定为宝箱怪,它是不是也能打出 99% 的准确率?

$$\frac{99}{100} = 99\%$$

一分不差。而这么一个假魔法,干了什么有用的事?一次都没有。它甚至比闭着眼睛乱猜更没价值:乱猜至少偶尔能撞上,而它把荒野里仅存的那一个真宝箱,精确地排除掉了。

这不是巧合,而是一整类情况。把这片荒野的不同时期排开来看:

时期 真宝箱 假定它是假的:准确率 它能找到的宝箱
无人踏足时 满地都是 $0\%$ $0$
探索中期 每 100 箱有 1 个 $99\%$ $0$
后来者 每 1000 箱有 1 个 $99.9\%$ $0$
终于被翻干净 一个不剩 $100\%$ $0$

四行里三行的准确率在 99% 以上,四行的战果全是零。「99% 的准确率」可以是对的,也可以同时是一句空话。这一点不需要任何公式,在荒野上走一圈就明白了。

大魔法师芙莉莲自然知道这手鉴定魔法是真的。可真魔法又能如何呢,99% 只是个虚妄的数字,这一点她同样是清楚的。更何况,宝物的珍贵足以让她情愿被多咬上几口。

那个假想的魔法压根没在做鉴定这件事,可它的分数偏偏很高。问题出在「准确率」的定义方式:它把所有判断一视同仁地数一遍,报对算一分。于是这个分数的高低,不仅取决于鉴定魔法有多灵,更取决于这批箱子里真宝箱和宝箱怪各占多少。而这个比例从来不是固定的:刚到一片野外,宝箱多;走了一天,宝箱少。同一手鉴定魔法,在不同地方可以有天差地别的准确率。

所以要评价这手鉴定魔法,第一步不是算它的准确率,是先问:想靠它知道什么?

她用它来干什么

芙莉莲站在一个箱子前,掏出鉴定魔法。她只问一件很具体的事:「这个箱子里有没有东西?」

「找到宝物」这件事对鉴定魔法提出两种要求,两件看似相似实则截然无关的事:

  • 有宝的时候,它能不能认出来,这是「灵敏度」。衡量它,等于把全部真宝箱摆到鉴定魔法面前,看它捞得出多少。
  • 没宝的时候,它会不会乱说,这是「特异度」。衡量它,等于把全部假箱子摆上去,看它有多少会瞎报。

这两个方向可以彻底背离。这手魔法完全可能把有货的全认出来,同时把没货的也全说成有货。「准确率」不够用,不是因为它算错了,是因为它把两件本来非分开不可的事情,揉成了一个数。

还有一个问题:这两种判断,她更关心哪一种?她在迷宫中走了一整天,遇到 100 只箱子,其中 1 只真宝箱。

  • 一种情况:她认出了那个宝箱,打开,带走一件东西。这一趟有收获,而那件东西只在这里有。
  • 另一种情况:她把宝箱怪当成了宝箱,被咬了一口。只有发型乱一点,明天还能接着找。

这两件事要记的,不是同一本账。前者记的是「拿到一件遗落的宝物,错过就不会再有第二次机会」,后者记的是「付掉一次无伤大雅的代价」。挨咬的代价可以反复支付,挨三口和挨一口差不太多;收获不能反复,那件东西只在那一个箱子里,也只在这一次。所以结论很干脆:

芙莉莲要的是宝物。被咬一口,只是去拿它要付的价钱。

「准确率」数的不是「有没有拿到」,而是「总共有几个判断对了」。而她要的那个量是:它说有货的时候,值不值得走过去。这两件事在空箱被错报成「有货」多到一定程度之后会分道扬镳,而「准确率」这个数字看不见这件事。

于是要把两边分开数

把真实情况和鉴定回答交叉起来,一共四个格子:

格子 名字 意思
有货 + 说有货 真阳性 真宝箱
有货 + 说没货 假阴性 漏掉了
没货 + 说有货 假阳性 被咬了
没货 + 说没货 真阴性 躲对了

四个格子各自对应一种性质不同的收益或损失:同样「说对了」,左上角让她心满意足,右下角让她毫发无损;同样「说错了」,一个挨咬,一个错失。合成一个数的那一刻,信息就没了。

灵敏度,在所有真有货的箱子里有多大比例被正确认出来:

$$Se = P(+\mid\text{有货}) = \frac{TP}{TP + FN}$$

特异度,在所有没货的箱子里有多大比例被正确说成没货:

$$Sp = P(-\mid\text{没货}) = \frac{TN}{TN + FP}$$

有了四格表,还差一个量:它说有货的时候,值不值得走过去。这个量叫阳性预测值,下文记作 PPV。

到这里,「99% 准确率」就显出了它的问题:它可以指 $Se$,可以指 $Sp$,可以指两者都 99%,还可以指某种把它们揉起来的东西。四个数字都叫「准确率」,都不算撒谎,但回答的是四个不同的问题。

揉在一起的数字为什么会塌

把四格揉成一个最常用的数:报对算一分,

$$\mathrm{Acc} = \frac{TP + TN}{N}$$

设这片地图上的宝物率为 $p$,这个数可以重写成一个加权平均:

$$\boxed{\mathrm{Acc} = p\cdot Se + (1-p)\cdot Sp}$$

看它的两个极端:宝物极其罕见时,$\lim_{p\to 0}\mathrm{Acc} = Sp$,准确率塌成特异度。原因很直白:罕见时绝大多数格子是「说没货」,默认答案几乎全对,分数自然漂亮。开头那个假想的撒谎鉴定,在这里现出了原形:它的准确率高,只是因为「说不是」在这个地方恰好是对的。

可芙莉莲要的不是「总共报对几个」,是「它说有货的时候,箱子里到底有没有」。这个量得单独算。

还是这片荒野,鉴定的都是同一手魔法。探索中期,10000 个箱子里有 100 个真货;快被翻完的时候,10000 个箱子里只剩 1 个真货。把两段时光并排:

探索中期 快被翻完
先验 $p$ $0.01$ $0.0001$
$Se = Sp$ $99\%$ $99\%$
$\mathrm{Acc}$ $99\%$ $99\%$
阳性预测值 $50\%$ $0.98\%$

Acc 一模一样,阳性预测值差了 51 倍。变化只来自一件事:它被用在了什么地图上。这个错误有个专门的名字,叫基率谬误:只看「有货的时候它报对的条件概率」,忘了先问「这里到底还有多少货」。

二、漏掉和误伤不是同一种错

挨咬的代价可以反复支付,错过的不能。既然两笔账轻重悬殊,判断就不能只数对错,得把损益直接摆开来记。设开到宝物得到收益 $b$,开错挨咬损失 $c_{FP}$。她可以去开箱,也可以放弃:不开的期望收益是零,开的期望收益是 $b$ 乘「说有货时真有货的概率」,减去 $c_{FP}$ 乘「说有货时其实没货的概率」。整理成比值,开箱的条件是:

$$\boxed{\frac{P(\text{有宝物}\mid\text{说有宝物})}{P(\text{没宝物}\mid\text{说有宝物})} > \frac{c_{FP}}{b}}$$

「她不怕被咬,她怕错过」这句话在数学上的全部内容,就是右边这个 $\dfrac{c_{FP}}{b}$。

如果开箱子的代价是零呢?$c_{FP} = 0$,判据右边是零,无论这片地图多空,都应该开。而这是对的。「反正开了不亏」,那当然一直开。

在一个「开箱完全免费」的世界里,「见到箱子就开」是可以论证的最优行为,而不是某个角色的性格缺陷。

这句话反转了批评的方向:要批评一个人乱开箱,必须先说清楚他的 $c_{FP}$ 不为零。芙莉莲的真实约束是,她的 $c_{FP}$ 不是零,但也差不了多少——她可以承受连续几次被吞,承受不了错过那件东西。

代价不会消失,只会换形态。代价明显不对称的场景至少有三种:

场景 $c_{FP}$ $c_{FN}$ 现实后果
市集地摊 小 中 收摊了白忙一天,但能再摆一天
传说难度区段 中 极大 错过的东西不会再刷,且会卡住进度
炼金熔炉 极大 小 投入的材料不可逆,失败就报废

同一套开箱判据,在这三张桌子上的答案完全不同:地摊前尽管随便挑,传说区里值得反复确认,而熔炉边恰恰要反过来——宁可错杀,不可错放。判据的形状没变,变的是代价的相对大小;每一个场景的阈值 $p^{*}$,画在哪里由代价决定,鉴定只通过 $LR_+$ 一个数进去。

三、宝物率不是常数

贝叶斯公式

魔法报「有宝物」的箱子里,真有宝物的占多少?设宝物率为 $p$:

$$PPV = P(\text{有宝物}\mid\text{报有}) = \frac{Se\cdot p}{Se\cdot p + (1-Sp)(1-p)}$$

它讲的事情一句话能念出来:分子是报「有」的箱子里真宝箱的那一份,分母是全部报「有」的箱子,真宝箱和被冤枉的空箱都算在里面。它是全文最重要的一个式子,一共只有三个成分——灵敏度、特异度、先验——缺一不可。来源就是贝叶斯公式:把分母 $P(\text{报有})$ 用全概率公式展开,再代入 $P(\text{报有}\mid\text{没货}) = 1-Sp$ 与 $P(\text{没货}) = 1-p$,即得。

它有个对称的量,阴性预测值,回答「说没宝物,敢不敢信」:

$$NPV = \frac{Sp(1-p)}{Sp(1-p) + (1-Se)p}$$

先验越低,它越贴近一。这是筛查里出名的不对称:先验低的地方,阴性很可信,阳性很可疑。

这两个量分别服务两种人。对输不起的冒险者,判「没宝物」就绕开走是他唯一的选择,漏掉的真宝箱不进开箱的账,NPV 就是他的全部世界。芙莉莲不在这本账里——判定对她从来不是开不开的开关,她用的是鉴定魔法,然后无视鉴定魔法。漏掉在她这里几乎不存在,存在的只是挨咬。

而野外的存货也只有这样,才能留给芙莉莲去尝试。判「有宝物」又真是宝箱的,大部分当场就被人开走了;能一路留到她面前的真宝箱,恰恰是一次一次被判「没宝物」的那些。她翻找的货架,是前面所有人按各自的判定筛剩的阴性堆。宝物率一年比一年低,不是箱子长了腿,是真货那一堆先见了底。

让 p 跑起来

为了让 PPV 达到某个值 $q$,先验至少要多少?令 $q$ 等于上式右边,两边乘分母、把含 $p$ 的项移到一边、提取公因式,解出:

$$\boxed{p_q = \frac{q(1-Sp)}{Se - q(Se + Sp - 1)}}$$

代 $Se = Sp = 0.99$,令 $q = 0.9$,得 $p_{0.9} \approx 0.083$。想让一个 99%/99% 的鉴定给出九成可信的结论,地图上的宝物率至少要到 8.3%,一个相当高的门槛。再取 $q = 0.99$,得 $p_{0.99} = 0.5$:先验得有一半是真宝箱,PPV 才追得上灵敏度。

还有一个退化情形:$Se + Sp = 1$ 时 $p_q = q$,魔法只是把先验原样搬过去,一个数字都没有改。$Se + Sp$ 这个和是衡量「信号含量」的粗略指标,它在 1 附近时这手魔法几乎没用。

而先验趋近零时,$PPV \approx \dfrac{Se}{1-Sp}p$,先验把 PPV 压向零,$Se$ 和 $Sp$ 决定的只是压下去的坡度。所以单独吹嘘其中任何一样都没有意义。

稀有度是指数

前面反复出现「宝物率很低」,这句话在数学上意味着什么?一个真宝箱留到今天,不是躲过了一关,是躲过了很多关。

芙莉莲不是第一个站到这些箱子面前的人。漫长的岁月里,每个箱子前面都来过别的冒险者。他们不像芙莉莲,被咬是真会没命的:只要判成宝箱怪,碰都不会碰。设每个路过的冒险者把真宝箱留下的概率是 $c$,一个箱子先后被 $m$ 个互不相干的冒险者路过还能留到今天:

$$P(\text{留到今天}) = c^{m}$$

代数字看它塌得多快:$c = 1/20$ 时,$m = 2$ 是 $2.5\times10^{-3}$,$m = 4$ 是 $6.25\times10^{-6}$,$m = 8$ 是 $3.91\times10^{-11}$。所以「稀有」不是一个形容词,是一个指数。

路过的冒险者 $m$ $c=1/20$ 时的留存概率
1 $0.05$
2 $2.5\times10^{-3}$
4 $6.25\times10^{-6}$
8 $3.91\times10^{-11}$

这带来三个后果:

  1. 稀有度不均匀。一个地段被多少人路过全看它的名气,难度相邻的两片地图,宝物率可以差一个数量级。
  2. 稀有度可以被制造。在一个箱子前面多摆一道独立的鉴定,稀有度立刻乘一个 $c$。
  3. 中段最危险。$10^{-4}$ 落在「说不准但不至于没人信」的区间里,会被大量尝试;$10^{-8}$ 落在「根本没人会去找」的区间里,不会被。于是人们系统性地把注意力集中在最不该集中的那一档上。

第三条是本文后半段所有故事的引擎。

还有一件事在现实里的用处比看起来大:四次放行分给四个陌生人,和四次放行落在同一个人头上,对「箱子能否留到今天」来说是同样罕见的两件事。决定稀有度的是独立判定的道数,不是做判定的人数。所以「多拉几个人一起看就更保险」在数学上站不住,除非那些判定本身相关。

把先验抬上去

先验低不是不可改变的。抬 $p$ 有四种正当手段,按成本从低到高排:

手段 机制 宝箱版本
打听 用别的信息源筛一次 问路过的商人这层有没有宝
看相 用与鉴定无关的独立特征 箱子的成色、铭牌、封蜡
换地图 只去宝物率更高的区段 只进传说难度区
加一道鉴定 换一个原理的鉴定魔法 找第二个鉴定魔法复核

代一组数字:普通区段先验 $0.005$,一个 $Se = Sp = 0.99$ 的鉴定魔法给出的 PPV 只有 33%;先用一次前置筛选把先验抬到 $0.026$,同一个魔法给出 72%。差距 2.2 倍,而鉴定魔法一个字没改。

提高先验是性价比最高的动作,因为它不需要动鉴定魔法,只需要动「谁来开箱」。

四、挖出来是杂物箱:第三种结局

前面所有结局都是二元的:有宝物 / 没宝物。现实里还有一类情况不在这个划分里。

荒野里还有这样一种箱子:鉴定魔法说「有宝物」,判定没有错。挖开之后,里面确实有东西——不是空的,也没有宝箱怪,而是满满一箱瓶瓶罐罐,件件寻常,拉到集市上也换不来几个铜板。

检测完全正确,而整个动作仍然是错的。它之所以需要单独讨论,是因为它的严重性被四格表遮蔽了:在四格表里这一箱被记为 $TP$,属于「成功」那一列。统计上它漂亮,账上它是个坑。

设鉴定后被挖开的箱子总数为 $D_{\text{det}}$,其中「不鉴定也会被别人挖走并用上」的为 $D_{\text{surf}}$,空挖率定义为

$$\mathrm{OF} = \frac{D_{\text{det}} - D_{\text{surf}}}{D_{\text{det}}}$$

麻烦在于 $D_{\text{surf}}$ 是反事实,「如果不挖会发生什么」没法数。但不需要它精确,需要的是一个能推翻它的证据。记录两个数:鉴定制度铺开后的年挖箱总数 $I_{\text{after}}$,和其中真正进入市场并成交的数量。一个不挖也会自然出现的箱子,也会自然地流入市场成交,所以 $D_{\text{surf}}$ 不超过铺开前的成交量 $V_{\text{before}}$,于是

$$\mathrm{OF} \ge \frac{I_{\text{after}} - V_{\text{before}}}{I_{\text{after}}}$$

代一组数字,$I_{\text{before}} = 35$,$I_{\text{after}} = 96$:至少 63.5% 是白挖。这个推导只用了两个数字和一个不等式,每一个前提都写在纸面上。它给的是下界,真实值通常更高。

这个不等式能证明「有白挖」,证明不了「有多少」和「是哪一类」。要加强它,需要更细的证据。一条是分布不变:如果被挖箱子的成色分布和处理成本分布在制度铺开前后都几乎一样,那么多出来的那批箱子在行为上无法与原有的区分,就不能声称它们「本来也会被挖」。另一条是分层:把箱子按「东西能不能出手」「要花多少工夫」分档,如果每一档都同比例上涨,那是整体普涨,也就是白挖;如果只有最低那一档在涨,那只是正常的品质下沉。

空挖的代价也不止「白挖一次」:

代价成分 宝箱版
直接 搬运费、处理费
机会 那一天本可以用来找传说箱
时间 从出发到空手而回的全过程
记录 这个人手上有一次「挖了个空」的经历
心理 未来的箱子会更难让人相信

后两项是会累积的:一次空挖不是一次性成本,它是一条会跟着这个人走的记录。所以空挖的真实代价不是空挖率本身,而是它乘上一个随时间增长的因子。

最后记一个幻想世界里数不出来的东西:怪物吃人,被挖开的箱子越多,遇到怪物的次数越多,死亡人数应该跟着涨。如果挖箱数涨而死亡数不动,多出来的挖箱就有问题了。这个「谁都不在意的、但确实存在的」外部终点,现在只是个类比。后面会遇到一个真实的领域,那里有真实可靠的外部终点,同样的推理一字不差地成立。

五、开一次不够:级联,和一个反例

多做几次

单道灵敏度 $Se$、各道相互独立时,$m$ 道全部漏掉的概率是 $(1-Se)^m$,累积灵敏度:

$$Se_m = 1 - (1-Se)^m$$

增长很快:$Se = 0.8$ 时,一道 80%,两道 96%,三道 99.2%。这就是「多道确认」的全部合法性来源:便宜的、灵敏度中等的鉴定,做三次比做一次强得多。反过来,单道灵敏度 20% 的鉴定,要做 21 道才能到 99%。

但这个公式有一个前提:各道独立。同一个箱子反复鉴定恰恰不独立——它的成色、封蜡状态决定了每一次鉴定的读数,第一道和第二道共享了几乎全部的物理信息。设相关性为 $\rho$,$m$ 道的有效道数不是 $m$,而是 $\dfrac{m}{1+(m-1)\rho}$。代 $Se = 0.8$、$m = 3$、$\rho = 0.5$:独立假设给出 99.2%,考虑相关性后只有 91.1%。真正安全的做法是每次用不同的手法,针对箱子的不同侧面——先听,后照,再拆封蜡。

级联:芙莉莲真正在做的事

$1-(1-Se)^m$ 描述的是「对全体箱子重复鉴定」,而现实里的多道确认是级联的:第一道说有宝物的,才进第二道。这个差别是决定性的,因为级联的每一道都在重新估计先验。

设定:$N = 10000$ 个箱子,$p = 10^{-4}$,1 个真宝箱。$Se = Sp = 0.99$,各道独立。第一道对全体做,阳性池 $0.99 + 99.99 = 100.98$ 个,PPV 约 0.98%。第二道只做这 100.98 个,阳性池缩到 1.98 个,PPV 49.5%。第三道只做 1.98 个,阳性池 0.98 个,PPV 98.98%。

0.98% → 49.5% → 98.98%。一个百分点的起点,三道之后是九成九。

用几率形式交叉验证一遍。先验几率是 $OD_0 = 10^{-4}/(1-10^{-4})$,每道鉴定的 $LR_+ = 0.99/0.01 = 99$,三道独立则几率连乘:$OD_3 = OD_0\times99^3 = 97.04$,转回概率得 0.98980——和级联链算出的结果一致到小数点后五位。两条路对上了:级联的有效性确实来自每道先验的重新估计,而「全体重复」的 $1-(1-Se)^m$ 和它不是一回事。

再做一个反例:同样三道,但都对全体做,任何一道报「有宝物」就算有。累积灵敏度还是 97%,可一个不是宝箱的箱子三道至少错报一次的概率是 2.97%,阳性池涨到 298 个,PPV 只有 0.33%,比一道的 0.98% 还低三倍。在低先验场景给同一个对象做三道独立鉴定,等于给了它三次被冤枉的机会。

方案 阳性池 PPV 累积灵敏度
一道 $100.98$ $0.98\%$ $99.0\%$
三道全体 $297.94$ $0.33\%$ $97.0\%$
三道级联 $0.98$ $98.98\%$ $97.0\%$

三道级联和三道全体的灵敏度几乎相同,阳性池差了 300 倍。这一行对比是本文最想让人记住的东西:

同样的灵敏度提升,可以花 300 倍的代价去换。区别只在于有没有把中间那 99 个空箱挡在门外。而挡住它们的不是更好的鉴定,是流程设计。

学徒和批量:把成本压下来

芙莉莲的鉴定魔法贵。便宜的替代品是让学徒看一眼,$Se \approx Sp \approx 0.95$。单独用不行:代进 PPV 公式,$p = 0.01$ 时只有

$$\frac{0.95\times0.01}{0.95\times0.01 + 0.05\times0.99} \approx 16\%$$

但学徒能一次看很多个箱子:把 $b$ 个箱子混成一批整批判一次,判阳的批再逐个细查。批处理的总成本在

$$b^{*} = \frac{1}{\sqrt{p}}$$

处最低。$p = 10^{-4}$ 时最优批是 100,总成本只有逐个细查的 2%。

批处理还顺手做了一件事:富集。判阳批里真宝箱的比例是

$$p_{\text{pool}} \approx \frac{1}{b}$$

在成本最优的 $b^{*} = 1/\sqrt p$ 上,细查池的先验恰好是 $\sqrt p$:$p = 10^{-4}$ 的地图,批处理把细查池的先验抬到 1%,抬高了 100 倍。批处理本身就是一个前置筛选,和「打听、看相」是同一个数学现象。

但批处理不是免费的。一个真宝箱拖着整批 $b - 1$ 个无辜箱子进细查,这叫污染;批越大,单个真宝箱的信号被摊得越薄,可能整批判阴,那一个真宝箱被漏掉,这叫稀释。池化的正确性依赖「混合之后信号仍然可检出」,这个前提的成立范围必须实测,不能假设。

六、把这台机器搬出幻想世界

在奇幻故事里挑一个开宝箱的小段子,展开这么一大篇计算,像是在认真地搞笑。

至少,认真是真的。

因为现实中也有一种东西,同样重要到错过就损失巨大,值得不惜代价把它检测出来。

它不是宝物。但一路读下来,应该也能猜到了:它是罕见病,而医学筛查就是鉴定魔法。

逐项对应

前几节攒下来的每一样东西,在现实里都有一个对应物:

  • 真宝箱,对应病人;宝箱怪、空箱、杂物箱,对应健康人。
  • 鉴定输出的「有宝物 / 没宝物」,对应检测结果的阳性 / 阴性。
  • 灵敏度与特异度,名字照旧;$PPV$ 与 $NPV$,对应阳性预测值与阴性预测值。
  • 先验 $p$(宝物率),对应患病率。
  • 宝物价值 $b$,对应早发现早治疗的收益。
  • 挨咬代价 $c_{FP}$ 与错过代价 $c_{FN}$,对应假阳性与假阴性各自的代价。
  • 被放行 $m$ 次才留到今天的箱子 $c^m$,对应罕见病。
  • 学徒目视加批量,对应廉价初筛加批量检测。
  • 外部终点「尸体数量」,对应死亡率。
  • 空挖,对应过度诊断。

把这份清单填进去之后,有一件听起来夸张但确实是事实的事:前几节的每一个公式,一个符号都不需要改。

现实里的检测没有一台是假的,搬过来的是数学,不是怀疑。这不是巧合。回头数一数,全部推导只用到六样东西:两个不交且并为全体的类,一个二元判定,两种错误,一个先验,一笔代价,一条判定规则。没有一条依赖「箱子里装的是什么」。只要一个领域同时满足这六条,它就自动获得本文的全部结论。

这份对应不是幻想

第三节埋过:一个真宝箱要在漫长的岁月里留到今天,它面前的每个冒险者都得放行,概率是 $c^m$。现实里有一件事和这个结构完全一致:一个人要发病,必须同时从父母双方各拿到一个罕见的致病基因。

以常染色体隐性遗传病为例,设致病等位基因频率为 $q$,携带者频率为 $c \approx 2q$,则发病率是 $q^2 \approx c^2/4$。这就是 $c^m$,其中 $m = 2$。代几个真实数字:脊髓性肌萎缩(SMA)发病率约 $1/10000$,反推携带者频率约 $1/50$;苯丙酮尿症(PKU)在欧洲人群中的携带者频率约 $1/30$ 到 $1/60$,报道发病率恰在 $1/10000$ 量级。罕见不是偶然的社会现象,它是 $c^m$ 这个函数的性质。

于是第一节开头那个循环闭合了:稀有 $\Rightarrow$ 先验低 $\Rightarrow$ PPV 低 $\Rightarrow$ 阳性结果不可信 $\Rightarrow$ 需要二轮、三轮确认试验。而四种抬先验的手段在这里全部有名字:打听是问诊和暴露史,看相是体格检查和影像,换地图是高危人群门诊,加一道鉴定是换一种原理的检测。批处理的两种新错误也有对应:污染,是一个阳性标本带崩整批复查;稀释,是标本混得太多、信号低于检出下限。

七、四格表进诊室

四个数字,关心的人不一样

数字 回答什么 谁关心它 典型误解
灵敏度 $Se$ 有病的时候仪器灵不灵 筛查项目负责人 以为它等于「阳性结果可信」
特异度 $Sp$ 没病的时候仪器会不会乱报 实验室与临床 以为它等于「阴性结果可信」
阳性预测值 $PPV$ 仪器说有病,真有病的概率 拿到阳性报告的人 以为它只由仪器决定
阴性预测值 $NPV$ 仪器说没病,真没病的概率 拿到阴性报告的人 以为高 $NPV$ 就能排除
总体准确率 $Acc$ 什么也回答不了 广告牌 以为它有任何诊断意义

最常见的一处错位是:$Se$ 回答的是 $P(\text{阳性}\mid\text{有病})$,$PPV$ 回答的是 $P(\text{有病}\mid\text{阳性})$。一对互逆的条件概率,在人群构成改变时会朝相反方向移动。

同一根验孕棒,从 92% 到 9%

家庭验孕棒的 $Se$ 和 $Sp$ 都在 99% 上下,广告写的「准确率 99%」说的就是这个。差别全在使用它的队列:

使用场景 大致先验 $p$ $PPV$
停经六周、自述有性生活 $10\%$ $92\%$
同房后几天、来查「到底有没有事」 $1\%$ $50\%$
育龄女性健康体检普筛 $0.1\%$ $9\%$

同一根验孕棒,同一批生产标准,从 92% 到 9%。而广告说的 99% 是哪一行,没人说明。

临床里同样的账天天在算。急诊胸痛患者做主动脉夹层的增强 CT,$Se \approx 98\%$,$Sp \approx 99\%$,够好了;可夹层在胸痛患者里的比例只有千分之二上下,算出来 PPV 约 16%。剩下六分之五的阳性是各种形态学上的临界情况,影像科给它们起了个诚实的名字:「不能排除」。它不是诊断,也不是排除,是一个中间态;临床真正负责的工作,是把这些假阳性从「疑似夹层」的队列里挪出去。

血库里的一场接力

1985 年,第一个 HIV 抗体检测(ELISA)进入血库筛查。它不差:$Se$ 和 $Sp$ 都在 99.5% 上下。可献血人群的感染率只有约万分之五,代进 PPV 公式:

$$PPV = \frac{0.995\times0.0005}{0.995\times0.0005 + 0.005\times0.9995} \approx 0.09$$

9%。每报出 11 个阳性,10 个是假的。这个数字远低于血库能接受的水平,但血库不能不筛——假阴性放出去的是一袋能传染 HIV 的血。出路就是级联:ELISA 阳性的标本全部送去做第二道原理不同的确证试验(Western blot),它把病毒蛋白逐条分开识别抗体,特异性高得多;后来又加入第三道核酸检测,直接找病毒的 RNA,顺带把「感染了但还测不出抗体」的窗口期也压掉了几天。

三道级联下来,阳性才敢用来拦截一袋血。这是级联在现实里最成熟的一个样板,而它的每一环,都对应前面算过的那一层:第一道负责不漏(灵敏度优先),后面的道数负责不冤枉(特异性优先),先验在每一道之间重新抬升。

窗口太窄的时候

级联的每一环都要花时间,而有一类疾病连时间都不给。婴儿型克拉伯病(Krabbe 病)是一种常染色体隐性的代谢病,有效治疗只有脐血造血干细胞移植一种,而且必须赶在症状出现之前做,通常认为出生后六个月以内是有效窗口。麻烦在于,这个病的可检出时点(酶活性测得极低)和症状出现时点几乎重合:等到症状出现再确诊,等于确诊即宣告失败。所以对它,新生儿筛查是唯一有意义的形态——必须在孩子还完全正常的时候就知道。

算一下账。患病率约 $1/80000$,酶活性检测取 $Se = 0.95$、$Sp = 0.99$。对已经有过患儿的家庭,先验约 $1/300$,代进 PPV 公式约 24%;对未治疗即死亡的病,$c_{FN}$ 近乎无穷大,代进开箱判据,实际先验高出阈值两个数量级,筛查压倒性地值得做。

但这个案例真正难的部分不在算术。一个孩子被确诊时六个月大,完全健康;移植可能救他,也可能在婴儿期把他杀死;不移植,他会在一年内开始不可逆的神经退化。筛查把一个原本由自然承担的痛苦,转成了一道必须由父母在两周之内做完的选择题。这是筛查伦理里最棘手的部分,公式帮不上忙。

年龄如何决定一个检测有多准

产前筛查是这套数学最好的展示窗口,因为先验 $p$ 是一个可以直接查表的已知量,而且它随母亲年龄变化得极其剧烈。

取无创产前筛查(NIPT)对 21 三体的性能:$Se = 0.995$,假阳性率 $5\times10^{-4}$,仪器这一侧进步巨大。把四个年龄代进 PPV 公式:

母亲年龄 先验 $p$ $PPV$
20 岁 $1/2000$ $49.9\%$
25 岁 $1/1400$ $58.7\%$
35 岁 $1/400$ $83.3\%$
40 岁 $1/100$ $95.3\%$

同一个检测,同一批试剂,同一个实验室,PPV 几乎翻倍,而一个物理量都没变,变的只是母亲的生日。

18 三体和 13 三体更糟:它们的灵敏度本身只有 0.90 左右,先验也更低,算出来 PPV 分别约 26% 和 13%。这直接产生了一条指南规则:NIPT 对 13 和 18 三体报阳性,必须做侵入性产前诊断才能做任何不可逆的决定。理由不是「这个结果不可靠」,而是「这个结果的 PPV 只有 13% 到 26%」。这两句话听起来像,但只有后一句能被引用、被质疑、被检验。

单次检测的假阳性看着温和,多年累积起来就是另一回事。Elmore 等人 2005 年发表的研究里,连续 10 年每年做一次钼靶筛查的 962 名女性中,49.2% 至少收到过一次假阳性结果,其中 145 人经历了活检。单次假阳性率约 5%,若 10 次之间完全独立,累积概率是 $1-0.95^{10} \approx 40\%$;实测的 49.2% 更高,因为同一个人的重复筛查之间存在正相关。同一批数据,收益是每 1000 人检出 11 到 13 例乳腺癌,成本是每 1000 人约 150 次活检。「检出 1 例,付出约 12 次活检」,这个比值比「99% 准确率」有信息得多。

韩国甲状腺:一次教科书级的现场实验

第四节那个空挖率公式,在现实里等着一组现成的数字。

2011 年起,韩国把甲状腺超声检查纳入国民癌症筛查项目。2012 年约 76.8 万人接受检查,发现 3104 例甲状腺癌;2013 年约 103.5 万人,发现 5865 例。发病率从筛查前大约每 10 万人 35 例,涨到 2014 年的约 96 例。

关键的一行在下面:同期甲状腺癌死亡率基本没有变化,维持在每 10 万人 4 例上下。

发病率涨了近三倍,死亡率一动不动。这正是空挖率不等式要抓的东西:筛查前「本来也会出现」的病例必然包含在历史发病率里,多出来的部分必须另找解释。代进去,$I_{\text{before}} = 35$,$I_{\text{after}} = 96$:

$$\mathrm{OF} \ge \frac{96-35}{96} \approx 0.635$$

至少 63.5% 是过度诊断。换一组口径更宽的统计,$6 \to 35$,算出至少 83%。后续研究给出的估计范围大致在 50% 到 90% 之间。

被多检出来的绝大多数是很小的甲状腺乳头状癌,影像上可疑,生物学上几乎不会进展。代价是数万台切除术、数万人终身服用左甲状腺素、数百万次随访超声——最后这一项还挤占着别人的影像资源,机会代价在这里是零和的。值得强调的是,这不是哪台机器坏了,检测本身没有问题;是把全人群倒进检测器这个动作出了问题。普查没有入口约束,它提高的不是先验,是假阳性。

神经母细胞瘤筛查是这个案例的姊妹篇。1990 年代美、日、欧几项大规模试验结果一致:筛查确实检出了更多病例,但死亡率不变,死亡年龄分布不变,各期生存率也不变。三条判据同时成立,说明新检出的病例和原有病例在临床行为上无法区分。这大概是医学史上过度诊断最干净的例子,代价是几万个孩子接受了手术和化疗,其中一部分本可以只是观察。

一个检测毫无问题、项目仍然失败的案例

先天梅毒在美国消失了四十年,然后回来了。2013 到 2014 年降到历史最低点,每 10 万活产不到 10 例;之后曲线掉头向上,到 2022 年前后重新回到每 10 万活产 100 例以上。

第一反应是「检测变差了」。但梅毒的血清学检测几十年没变,实验室的标准化程度比二十年前更高。变的是检测之外的东西:产前检查的覆盖率、无症状梅毒的识别率、无保险和无稳定住所人群的比例。

把机制写成公式:设有效筛查强度为 $Se_{\text{eff}} = c\cdot Se\cdot\alpha$,其中 $c$ 是筛查覆盖率,$\alpha$ 是查出后真正完成治疗的比例。取 $Se = 0.95$、$\alpha = 0.90$,覆盖率从 95% 掉到 60%,有效灵敏度就从 81% 掉到 51%,检测性能一个数字没改。病例数还正比于患病率 $p$ 和查出不治的传播概率 $\beta$——五个因子里只有 $Se$ 属于实验室。

筛查项目的有效性等于检测性能乘以可达性乘以依从性乘以治疗可得性。四个因子里,只有第一个属于技术。

原因很简单:$Se$ 从 0.90 提到 0.99 需要一次实验室级别的技术突破,而覆盖率从 0.90 掉到 0.70,只需要一次预算削减。

塔斯基吉:$c_{FN} = \infty$ 被卡住了

1932 年起,美国公共卫生服务部在阿拉巴马州塔斯基吉县观察梅毒的自然史,数百名黑人男性患病而不被告知、不被治疗。1947 年青霉素已成为标准治疗,项目仍决定不给。

这个案例常被当作伦理丑闻引用,但它同时是第二节那个开箱判据的一个标本。梅毒三期致死,$c_{FN}$ 是无穷大,而判据阈值

$$p^{*} = \frac{c_{FP}}{b\cdot LR_+ + c_{FP}} \to 0$$

塌到零:任何阳性都必须按患者处理,没有例外,不存在任何「值得再权衡」的余地。所以塔斯基吉不是「当初不知道有青霉素」,也不是「权衡之后选了不治」,它是在 $c_{FN} = \infty$ 的情形下做出了 $c_{FN} = 0$ 的决策。这是逻辑上的失败,不是知识上的失败。

疾病定义没变,治疗能力变了,而决策框架没变。前面说过,检测的定义、诊断阈值和治疗方案三者是耦合的,改掉任何一个,另外两个都必须重新审视。这套结构直到 1974 年《国家研究法案》要求所有联邦资助的人体研究通过伦理委员会审查,才被制度性堵住——而那已经是事件暴露之后的事。制度总是跟着丑闻走,而不是提前于丑闻。

八、筛查项目是怎么拼出来的

1968 年,世界卫生组织出版了 Wilson 和 Jungner 的《疾病筛查的原则与实践》,二十几页的小册子,地位相当于这个领域的 Koch 推定列表。十条原则里最值得注意的是:没有一条要求「筛查必须被证明有益」。其余几条倒是和前文的量一一对应:要求疾病有一段可识别的无症状期,那是筛查的唯一立足点;要求对已识别的患者有可接受的治疗,因为没有治疗,$c_{FN}$ 就是无穷大;要求查找病例的费用与整体医疗支出相平衡,那是成本效果分析的位置。

这个缺陷后来被补上,现代标准要求随机对照试验证据优先。门槛实际变高了:美国处于某种形式筛查推荐中的疾病,从 1995 年前后的约 110 种降到 2012 年的约 50 种,删掉的大部分不是被证明有害,而是被证明「证据不足」。

RCT 之所以必要,是因为筛查的收益必须用死亡率来量,而死亡率是最不容易被过度诊断污染的终点;之所以昂贵,是因为样本量由效应量的平方决定:死亡率相对下降 20% 的项目,两组各需约 1700 人;下降 10%,各需约 7000 人;再往下,样本量很快把项目判死。罕见病筛查做不了 RCT,不是因为效应量小,而是因为事件数根本凑不出来——发病率 $1/50000$ 的病,要凑 30 个死亡事件需要一百五十万人和一二十年随访。

只讲失败会让人以为筛查都不行,所以值得记几个被广泛认为成功的项目。先天梅毒的孕期筛查加青霉素预防:检测便宜稳定,治疗一次性且代价极小,婴儿的结局可以被观察,筛查效果能被独立验证。血库的 HIV 与肝炎筛查:单看 PPV 只有 9%,但假阳性的后果只是扔掉一袋血,假阴性的后果近乎无穷大,而且同一份血天然被多个独立的检测系统交叉覆盖——那是级联的实践形态。乙型肝炎表面抗原筛查加疫苗:筛查的输出能直接转化成一次便宜、有效的干预。归纳起来就三条:治疗便宜有效,窗口够宽,假阳性有便宜的出路。三者齐备,筛查才立得住。

如果只从这篇文章带走一份操作性清单,可以是这几条:

  1. 看到「准确率 X%」,先问哪个方向:$Se$、$Sp$、$PPV$ 还是 $Acc$。
  2. 看到任何 $PPV$,先问先验。没有先验的 $PPV$ 是一个没有单位的数字。
  3. 低先验下,$Sp$ 比 $Se$ 值钱:提高 $Se$ 的边际收益是 $p$ 量级,提高 $Sp$ 的边际代价是 $(1-p)$ 量级。
  4. 多轮筛查的正确形式是级联,不是叠加。
  5. 批量检测是富集操作:先验从 $p$ 抬到约 $1/b$,成本在 $b = 1/\sqrt p$ 处最低。
  6. 阈值是决策,不是性能。美国儿童血铅干预线从 40 降到 10,测量数据没变,变的是「多高的血铅值得干预」这个决策。
  7. 排序不产生信息:没有额外信息时,从 100 个箱子里挑「最可疑」的一个,命中率就是 $p$。
  8. 警惕金标准本身有误差,那会让所有误差估计都偏乐观。
  9. 过度诊断三判据:发病率升而死亡率不动;生存曲线不变;各期生存率不变。
  10. 筛查项目最脆弱的环节是可达性,不是检测。

九、回到那只箱子

把芙莉莲的期望收益算一遍

现在回到芙莉莲。设开到宝物的收益 $b = 100$,挨咬代价 $c_{FP} = 1$,$Se = Sp = 0.99$。开一次箱的期望收益是 $E = 101\cdot PPV - 1$,对几个先验算一遍:

先验 $p$ $PPV$ 期望收益 该不该开
$0.01$ $50.0\%$ $+49.5$ 疯狂开
$0.001$ $9.0\%$ $+8.1$ 开
$0.0001$ $0.98\%$ $-0.01$ 五五开
$0.00005$ $0.49\%$ $-0.50$ 别开

零点在 $p \approx 1.01\times10^{-4}$。把 $c_{FP}$ 换成 5 再算一遍:期望收益变成 $105\cdot PPV - 5$,零点上移到 $5\times10^{-4}$。这个敏感性值得摆出来:在极低先验的场景里,「挨一次咬的代价有多大」这个纯主观判断,比「鉴定魔法有多准」更能决定行为。精度从 99% 提到 99.9%,相当于把 $LR_+$ 从 99 提到 999,开箱阈值只降到原来的十分之一;而「一次失败能忍」和「一次失败就退出」之间,是 $10^{-4}$ 到 $5\times10^{-4}$ 的差距。修正信念的收益,不如修正偏好。

她做的正是级联

一个 $Se = Sp = 99\%$ 的鉴定魔法,在宝物率 $10^{-4}$ 的环境里,单次阳性只有 0.98% 可信。而级联那一节算过那串数字:0.98% → 49.5% → 98.98%。所以芙莉莲做的事情,在结构上等价于级联确证,她的确证方式是把箱子打开看。至于「一次又一次被吞」,从上面的表看,是 $c_{FP}$ 取值偏低时必然的均衡状态:期望收益为正,但单次成功率不到 5%。她在持续为正期望的赌局付费。

「99% 准确率」没有错,只是它回答的问题和观众以为的不是同一个。观众以为它说的是「我开箱有九成把握」,它实际说的是「我判断宝箱时,看错的概率是 1%」。从前者到后者,要走一遍条件概率,要走一遍先验,要走一遍代价比较。而绝大多数误读,都发生在以为可以跳过这三步的地方。

医院里没有宝箱怪,但有检验单。有检验单的地方,就有 $PPV$。

如果她能改一件事

假设给芙莉莲三个选项。起点 $p = 10^{-4}$,用几率形式算:每做一轮鉴定,后验几率乘一个 $LR_+ = \dfrac{Se}{1-Sp}$,$n$ 轮之后

$$OD = \frac{p}{1-p}\cdot LR_+^n, \qquad PPV = \frac{OD}{1+OD}$$

方案 每轮 $LR_+$ 轮数 $PPV$ 相对成本
原始 99 1 $0.98\%$ 1
换更准的魔法(99.9%) 999 1 $9.1\%$ 极高
加一道独立鉴定 $\times 99$ 2 $49.5\%$ 中
便宜鉴定重复三次($Se = 0.80$,$Sp = 0.98$) $\times 40$ 3 $20.4\%$ 低

加一道独立鉴定胜出得毫无悬念,但背后的结构规律比这个具体比较更通用。对上式取对数:对数几率是加法,每一次检测贡献一个 $\ln LR$,与先验相加。所以「信息」是可以相加的资源,它的单价是 $\ln LR$ 除以成本:

方案 $\ln LR$ 假设成本 单价
换更准的魔法 $6.9$ $100$ $0.069$
加一道独立鉴定 $4.6$ $50$ $0.092$
便宜鉴定重复 $3.7$ $2$ $1.84$

便宜鉴定的单价高出 20 倍以上。只要它允许重复足够多次,它一定能赢。

在先验极低的地方,信息必须便宜。精确度的提升在成本上是指数级的,而把同一份信息重复采多遍是线性叠加的。所以在「十次里九次都是坏消息」的检测环境里,正确的策略是「用便宜的手段多做几次」,而不是「花钱做一次准的」。芙莉莲真正做对的事不是「她的鉴定魔法很准」,而是「她反复确认」。

一些豆知识

完美使用与一般使用

1970 年代口服避孕药上市时,广告上的「超过 99% 有效」指完美使用条件下的有效率;一般使用条件下失败率约每年 7%,实际有效率约 93%,差距几乎全部来自忘记吃药。同一句话里的两个 99%,一个指条件概率,一个指无条件概率。

「罕见」的定义非常随意

欧盟的定义是患病人数少于 5000,日本和台湾用年发病率低于万分之一,美国用患病人数少于 20 万。同一疾病在三个定义下的罕见程度可以差几十倍,而这直接决定它能不能进筛查清单、能不能进医保目录、有没有药厂愿意做临床试验。

一个数字巧合

本文前面算出,$Se = Sp = 0.99$ 的鉴定魔法在宝物率低于 $10^{-4}$ 时不值得开箱;又算出,常染色体隐性遗传病的发病率是携带者频率的平方,携带率 $1/50$ 对应发病率 $1/10000$。两个 $10^{-4}$ 在这里碰上了。当然这是巧合,一个讲物品稀有度,一个讲基因频率。

但如果把「宝箱」理解成「致病等位基因」,把「鉴定」理解成「症状」,那么医学史上有很长一段时期,人类处理罕见病的方式确实就是这样:靠症状去找,等于拿着一个 99% 准确的检测器,在宝物率万分之一的地图上乱逛。直到新生儿筛查把「确认」这个环节从「症状」换成「实验室指标」,这个局面才改变。前面那个级联算例,就是这次替换的数学写照。

参考(节选)

  • Wilson JMG, Jungner G. Principles and Practice of Screening for Disease. WHO Public Health Papers 5. 1968.
  • Elmore JG, et al. Screening for Breast Cancer. N Engl J Med. 2005; 352(15): 1495-1504.
  • Welch HG. Overdiagnosis: What It Is and What It Means. Johns Hopkins University Press. 2011.
  • 关于韩国甲状腺癌筛查的前后对比分析,New England Journal of Medicine. 2016.
  • 关于 Krabbe 病纳入新生儿筛查的早期报道,含纽约州项目启动后首批病例的经过。
  • 国家卫生健康委员会,第一批罕见病目录(2018 年)与第二批罕见病目录(2023 年)。

今天遇到一个问题,Firefox 打开 B 站直播,没有画面和声音,而弹幕和聊天一切正常。

控制台有几百行黄色的「Referrer 策略」警告,但查证后确认是无害噪音:Firefox 对跨站请求强制使用更严格的默认 referrer 策略,B 站声明的老策略被忽略,每加载一个跨站资源就记一条日志,请求本身照常发出。

初步排查

  • 同一局域网的 Win10 和 MacBook 上,同版本 Firefox 158.0.1 一切正常;排除 IP 级风控问题。
  • 本机 Edge 正常;排除网络链路和 CDN 问题。
  • B 站点播正常;排除解码、渲染整条媒体管线的问题,故障只在直播播放器的初始化。
  • 本机全新安装、零配置的 Zen-browser(Gecko 内核)一样黑屏;锁定 Gecko 内核。
  • 清空 Cookie 缓存、新建配置文件、关闭跟踪保护,hosts、防火墙、代理逐一核对,均无影响。

至此故障的触发条件已经收得相当窄:要「本机 × Gecko 内核 × 直播」三个条件同时成立,缺一不可。问题从「为什么放不出视频」变成「初始化卡在哪一步」。

深入排查

用调试协议在页面加载前往直播间注入钩子,包装 fetch 和 XMLHttpRequest 记录调用。播放器拿到了播放地址:接口返回 200,avc、hevc 的流列表一项不缺,然后就没有然后了,页面里没有创建任何 video 元素,也没有向 CDN 发出一条拉流请求,初始化卡死。在页面里对照测了一圈媒体 API:

API 结果
requestMediaKeySystemAccess('com.widevine.alpha', …) 1 毫秒正常返回
mediaCapabilities.decodingInfo(…) 永不返回

最简单的复现代码,任意非 B 站网页的控制台粘贴即可,四秒后见分晓:

1
2
3
4
5
6
7
8
9
10
(async () => {
const t = Date.now();
let ok = 0, rej = 0, ms = 0;
navigator.mediaCapabilities.decodingInfo({
type: 'media-source',
video: { contentType: 'video/mp4; codecs="avc1.64001F"', width: 1920, height: 1080, bitrate: 8e6, framerate: 60 }
}).then(() => { ok = 1; ms = Date.now() - t; }, () => { rej = 1; ms = Date.now() - t; });
await new Promise(r => setTimeout(r, 4000));
console.log(ok ? `正常,${ms}ms 返回` : rej ? `${ms}ms 被拒绝` : '挂死,四秒未返回');
})()

钩住 decodingInfo 后发现,播放器初始化时调了它四次,探测 HEVC 和 AV1 的支持情况,四个 Promise 全部悬空。编码探测不返回,播放器就永远停在选编码这一步,弹幕和聊天走另一条路,所以一切如常。写了个 150 毫秒超时保底的垫片脚本注入,直播立刻能放,因果链闭合。

原因

在撰写向 Bugzilla 提交的反馈时,发现已经有相关问题汇报了:这是 Mozilla 自己的回归,见 bug 2082265,Twitch 和 Prime Video 的播放器是同款受害者:

  1. Firefox 158 修 bug 2051510(MediaCapabilities 在远程进程上报编码支持之前抢先回答,把不支持的 HEVC 报成支持)时,把每个视频能力查询都改成了等待所有远程媒体进程上报完毕再回答,等待对象包括 media engine CDM 工具进程。
  2. 在部分 Windows 机器上,这个进程在 LPAC 沙箱下永远启动不完,于是查询永远在等,Promise 永远悬空。
  3. 修复进了 159.0b4 和 160,158 标记 wontfix 不修。

顺带排除了所有本机嫌疑:rdd、utility 进程、硬件解码、WMF 的开关逐个关掉再测,全部照挂;两台机器六千二百多项生效配置全量 diff,零相关差异;进程里没有第三方注入 DLL;显卡驱动从 591.86 换到最新 621.01 也没用。

解决办法

按省事程度排:

  1. about:config 里把 media.wmf.media-engine.enabled 设为 0,实测 27 毫秒恢复。副作用是关闭 WMF Media Engine 播放路径,B 站用不到它,哪天 DRM 场景出问题改回 2 即可。
  2. 升级 Firefox 159.0b4 及以上,根治。
  3. 或者用油猴脚本给 decodingInfo 加超时垫片顶着。我已经写了个脚本:firefox B 站直播无视频声音修复脚本。

至于这台机器为什么落在「部分机器」里,沙箱里那个进程为什么偏偏启动不完,bug 下面也没人给出答案,就这样吧。

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

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

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

黄巢

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

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

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

据史载:

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

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

起义初期,黄巢与王仙芝合兵一处,转战山东河南,沿途饥氓争赴,队伍迅速壮大。王仙芝战死后,黄巢收编其众,成为义军首领。此后他避实击虚,挥师南下,在岭南休整补充。次年北伐,势如破竹,于 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。游戏内赛季征服要求之一是一小时内完成全部剧情,经过研究后使用猎魔人成功完成,记录过程如下:

配装:

头肩胸手腿脚:恐惧之地机轮甲 三件 + 暗影装束 三件
腰带护腕:克雷姆的强力腰带 + 沃兹克臂甲
饰品指环:地狱火护符 + 黄道黑曜石指环 + 瑞秋的行窃之戒
武器副手:杨的反曲弓 + 希瑞的九号箭袋
萃取:要塞弩机 + 皇家华戒 + 猎手之怒,替换: 黎明
宝石:屯宝者的恩惠(25级) + 闪电华冠(25级) + 迅捷勾玉(25级)
被动:战术优势 + 热力追击 + 完美主义者 + 暗夜追踪者 + 血腥复仇
药水:无尽恐惧药水
技能:蓄势待发-集中心智 + 烟雾弹-飘忽不定 + 暗影之力-任意 + 复仇-恨意迸发 + 追踪箭-吞噬箭 + 扫射-暗影游移

备注:

  • 配装包括两个方向,其一是直接增加速度,以及保护增速 buff(弩机)。另一个是通过减 CD 减耗维持烟雾弹和扫射可以持续不断地释放。因此部分装备和被动可以替换甚至不用,只要能满足能量需求即可。但增速方向的装备几乎是必须的,只有瑞秋的行窃之戒+无尽恐惧药水因为触发场合不多,可以省略。保留 猎手之怒 是因为在长达一小时的游戏过程中,必然要打很多次手动吞噬箭以保持 20 层动能,使用猎手之怒可以降低这部分的时间。
  • 由于过程较多涉及鼠标左键点击互动,因此把吞噬键放到 1-4 号上比较好,左键建议绑定复仇技能。

剧情流程:

节约时间要点:

  • 剧情进度大部分是人物对话触发的,但也有一些是地形位置触发的。对于后者,可以跳过对话直接跑到目的地节约时间。
  • 按 T 回城时,会在原地留下传送门。而打开地图传送时,不会留下传送门,也就不会顶掉之前的传送门。利用这点可以在『先 A 点再 B 点再回 A 点』这种流程中为第二次去 A 点节约时间。
  • 暗黑3 地图分为完全随机生成的地图,以及固定大地图内含随机刷新地块两种类型,哪怕是随机的地图也是有一定规律的,熟悉地图可以节约不少时间。
  • 在游戏设置里为跳过对话设置顺手的快捷键,过程中需要频繁使用。
  • 每章剧情要点如下:

A1

  1. 旅馆清完怪,等待与莉娅交谈,别急着出门。
  2. 尸母女王在一路向上向右的传送点旁,不用杀路上三个奖励目标。(后同,奖励目标都不考虑)
  3. 艾德莉亚的小屋 - 隐藏的地窖里,等待与莉娅交谈。
  4. 大教堂一层和四层,都是一个大方形主路上长出支路,目标出口在支路上。
  5. 救完凯恩,立刻使用 T 回城。
  6. 不要堵在地窖门口,要给铁匠开门让路,清完地窖与铁匠交谈。
  7. 墓地固定在哭泣山谷最上面,三个墓穴中,正确的一条路进门不远有个三岔口。
  8. 拿完王冠不要用 T 回城,必须打开地图,点新崔斯特姆来回城。
  9. 铁匠修好王冠后,用传送门回到救凯恩的地方(ABA 式节约时间)。
  10. 消灭柯麦克周围的人后,直接去开他的装备箱,然后直接跑到楼梯处,不必交谈(免谈式节约时间)。
  11. 大教堂四层和一层,都是一个大方形主路上长出支路,目标出口在支路上。
  12. 沿苦难旷野略靠边的位置找卡兹拉洞穴。卡兹拉洞穴是小型随机刷新地块,不会刷新在地图中间的大型随机块上。
  13. 拿到宝剑碎片后不需要听麦格坦跟你讲下一个碎片在哪,直接回城。你是奈非天,你啥都知道。
  14. 拿完烂木林两个圣物,地图传送回到沉没神庙。
  15. 回城对话-船夫交谈-沃桑小教堂门口清怪-对话-地窖-传送回城-对话-旅馆-对话-传送沃桑,这段连续操作要背一下。
  16. 蛛后洞窟地图是固定的,一个地图块,连一条极长直线路,连四个地图块2x2,其中一个地图波浪路连到第六个地图块,再极长直道进 BOSS 房。
  17. 和卡瑞拉对话 - 杀蛛 - 点地上毒液 - 对话。
  18. 救完卡瑞拉去南方高地出门不必交谈。卡瑞拉遗失马车是小型刷新点,沿地图边找。
  19. 痛苦大厅一层和三层的出口支路一定是大铡刀长直路。

A2

  1. 每次出营地都有和艾席拉的对话。
  2. 不用和下坡路上的魔女对话,直接跑到路尽头。
  3. 石头打开后,不需要救迎面那波铁狼,直接靠地图边走去找两个仪式房。A2 的地图都是固定形状大地图,只有刷新地块变化。
  4. 卡辛哨站进屋,打开铁狼笼子后,马上出去走到哨站入口准备打架。
  5. 达格尔绿洲为上下两片大地图,每片中间一个大随机地块,周围几个小随机地块。被遗忘的废墟可能在任意一个小随机地块里。
  6. 拿完精魄通过传送石回到绿洲,然后跑到地图顶部远古水道门口,按 T 回城。对话完后进传送门,直接进水道。
  7. 凄凉沙漠是固定大地图带随机地块,同样,两个目标地块都是地图较靠边的小随机地块。藏书馆则在最顶部。
  8. 打完库勒记得拿灵魂石再走。
  9. 救人时,所有地图黄圈都要靠近晃一晃,黄点的车子要点一下。

A3

  1. 开局直接按 T 回城,不要和泰瑞尔废话。
  2. 直接去点五个烽火和三个投石机。
  3. 架设投石机,先沿路把三个任务都开了,出现三个进度条。然后从第三个往回挨个推。
  4. 要塞深渊地图随机度比较大,尽量选择固定一个方向走到底。
  5. 军械库清完怪,等待与莉娅交谈。
  6. 进入战场,不需要交谈。恶魔弩炮在地图四角,抛石机正中。
  7. 亚瑞特巨坑是随机地图,尽量往一个方向走。天谴者之塔和诅咒者之塔是固定地图,转圈就行。
  8. 消灭辛迪娅,打碎罪恶之心后直接下去到门口等。到阿兹莫丹都是固定地图,打完阿兹莫丹记得拿灵魂石。

A4

  1. 第四章开场希望花园一二层是固定地图,腐化物位置也是固定的。
  2. 命运图书馆在希望花园一层最左边。
  3. A4 没有特别节约时间的操作,银色高塔跑图也不复杂。

A5

  1. A5 有较多的左键点击开门互动行为,合理分配技能对应按键可以节约一些时间,当然仍要以习惯为先。
  2. 泰瑞尔跑动的路上也可以对话,不必等他跑到位。
  3. 威斯特玛城区地图是随机的,尽量找长直道沿一个方向走。
  4. 吉迪恩大街清完炼魂阵后等一下,救左边杂物堆里的附魔师。
  5. 棘草荒园出口一定在地图边缘,靠边走。
  6. 贵族陵园清完大炼魂阵后往回走,和路上的附魔师对话。
  7. 杀完厄兹尔等和马萨伊尔聊天完再走。
  8. 杀完艾德利亚等洛拉斯出来聊完天再走。
  9. 在第三个门捡到第一个攻城符文后,不用等英普瑞斯帮开门,直接地图传送到左上角废弃的攻城营地。
  10. 永恒战场的两个攻城符文事件绑定了固定地块,其一在大地图边缘,另一在大地图较大完整地块,有一个圈。
  11. 混沌要塞一层不要和泰瑞尔对话。
  12. 混沌要塞二层入口消灭门口守门怪直接到圆弧尽头等,大地图是 3x3 的块,沿左边走。

这是对暗黑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%。

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

阅读全文 »

今天是 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
// 代码略

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

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

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

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

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

0%