dzmingli_vs_floating 设计

设计目标

本包回答关于 DzmingLi/decimal@0.2.2 和 Luna-Flow/floating/decimal_gda@0.7.1 的两个问题:它们算出的精确十进制结果是否相同?当系数增长到数千位时它们算得有多快?只有在两个答案都正确时才报告速度,因此设计围绕一个不会因两个库彼此一致而被蒙蔽的正确性检查展开。

API 页列出各条目,教程运行它们。测量数据见性能章节。

数学背景

十进制值

有限十进制数是一个对 (c,s)∈Z×Z(c, s) \in \mathbb{Z} \times \mathbb{Z},它表示

v(c,s)=c⋅10−s.v(c, s) = c \cdot 10^{-s}.

映射 vv 不是单射:(c,s)(c, s) 与 (10c,s+1)(10c, s + 1) 表示同一个数。GDA 算术区分这样的对(指数是结果的一部分),而基准比较的是数。因此它需要每个等价类的一个规范代表;本包取 s≥0s \ge 0 且 cc 末尾零最少的那个对。

对 c≠0c \ne 0,记 d(c)d(c) 为 ∣c∣|c| 的十进制位数,于是

10d(c)−1≤∣c∣<10d(c).10^{d(c) - 1} \le |c| < 10^{d(c)} .

working_precision 把 dd 计算为 ∣c∣|c| 的十进制字符串长度,因此 d(0)=1d(0) = 1。

精度与向零舍入

精度为 pp 的 GDA 上下文恰好表示满足 ∣c∣<10p|c| < 10^{p} 且 ee 在指数范围内的数 c⋅10ec \cdot 10^{e}。运算先计算精确结果 xx,再把它舍入为这样的数。舍入模式为 Down(向零)时,结果为

round⁡p(x)=sgn⁡(x) ⌊∣x∣⋅10 p−E(x)⌋⋅10 E(x)−p,E(x)=⌊log⁡10∣x∣⌋+1,\operatorname{round}_p(x) = \operatorname{sgn}(x)\, \bigl\lfloor |x| \cdot 10^{\,p - E(x)} \bigr\rfloor \cdot 10^{\,E(x) - p}, \qquad E(x) = \lfloor \log_{10} |x| \rfloor + 1 ,

因此 ∣round⁡p(x)∣≤∣x∣|\operatorname{round}_p(x)| \le |x|,当且仅当 xx 至多有 pp 位有效数字时取等号。在固定的 kk 位小数处做同样的截断,记为

trunc⁡k(x)=sgn⁡(x) ⌊∣x∣⋅10k⌋⋅10−k.\operatorname{trunc}_k(x) = \operatorname{sgn}(x)\,\bigl\lfloor |x| \cdot 10^{k} \bigr\rfloor \cdot 10^{-k} .

带参考预言机的差分测试

设 I1,I2I_1, I_2 为被测实现,OO 为预言机,κ\kappa 为映到一个可判定相等的集合的规范化映射。对输入 xx,实现 kk 的判定为

valid⁡k(x)  ⟺  κ(Ik(x))=κ(O(x)).\operatorname{valid}_k(x) \iff \kappa(I_k(x)) = \kappa(O(x)) .

成对差分测试只检查 κ(I1(x))=κ(I2(x))\kappa(I_1(x)) = \kappa(I_2(x))。这种检查对共模错误视而不见——即两个实现给出同一个错误结果——而且失败时也不能说明是哪一方出错。有了预言机,每个实现都被单独评判,两者的分歧也可由各自的判定来解释。

设计决策

三个实现,测量其中两个

问题。 基准必须把“快”与“快但错”区分开。选项。 两个库之间的成对一致;以一个库作为另一个库的参考;独立的预言机。选择。 基于 BigInt 的独立精确预言机,在计时之外对每个数据集运行(ValidationCoverage::EveryDataset)。理由。 预言机小到可以逐行审查,与两个库都不共享代码,并且它发现了 DzmingLi 从 4,096 位起的 108 个失败,而成对一致只会把这些报告成“两个库不同”。只有当某规模下两个实现的所有校验都通过时,该规模才进入加速比图。

中立表示与零容差

问题。 两个库输出结果的方式不同(指数记法、末尾零、零的符号),而近似比较会掩盖真实错误。选择。 每个结果都解析为 DecimalValue,并按 canonical_string 是否相等来比较,即容差恰好为零。理由。 在下文的精度约定下,所有被测运算的结果都可精确表示,因此正确的实现必须逐位复现它;可靠性论证见正确性与不变量。代价是比较忽略了指数同类(cohort)和状态标志;官方 decTest 审计(tools/run_dzmingli_dectest_audit.sh)覆盖这些方面。

精确有限预言机

预言机只对系数使用整数运算。

加法与减法把两个操作数对齐到 s=max⁡(sℓ,sr)s = \max(s_\ell, s_r):

cℓ10−sℓ±cr10−sr=(cℓ10s−sℓ±cr10s−sr) 10−s.c_\ell 10^{-s_\ell} \pm c_r 10^{-s_r} = \bigl(c_\ell 10^{s - s_\ell} \pm c_r 10^{s - s_r}\bigr)\, 10^{-s} .

乘法将系数相乘、标度相加。除法把商写成整数分数并约分,

cℓ10−sℓcr10−sr=NM,N=cℓ10sr, M=cr10sℓ,N′=Ng, M′=∣M∣g, g=gcd⁡(N,M),\frac{c_\ell 10^{-s_\ell}}{c_r 10^{-s_r}} = \frac{N}{M},\quad N = c_\ell 10^{s_r},\ M = c_r 10^{s_\ell},\quad N' = \frac{N}{g},\ M' = \frac{|M|}{g},\ g = \gcd(N, M),

并且只有当 M′M' 不含 22 和 55 以外的素因子时才接受。然后它找出满足 M′∣N′10kM' \mid N' 10^{k} 的最小 kk,并返回 (N′10k/M′, k)(N' 10^{k} / M',\ k)。整数除法、余数、幂、完全平方数的平方根、FMA 和一元运算都由这些规则推出;API 页上的表格列出了它们。

截断类规则(DivideInteger、Quantize、Rescale、ToIntegralExact、ToIntegralValue)之所以向零截断,是因为 BigInt 除法如此。GDA 的 quantize、rescale 和 to_integral_* 按上下文的模式舍入,所以两个上下文都使用 Down;任何其他模式都会使预言机在这些运算上出错。在已发布的夹具中,这些运算本来就是精确的(见下一条决策)。

精度约定

问题。 每个运算都必须在能容纳精确结果的精度 pp 下执行,这样就永远不会发生舍入,零容差也才公平。更大的 pp 并非没有代价:GDA 除法的工作量会随精度增长,所以过大的 pp 可能会把结果不需要的工作也计入时间。选择。 按夹具根据操作数计算 pp,取精确结果位数的最小简单上界再加保护位:working_precision(op, ℓ, r) + 2,并对 Power 和 Fma 另行处理。理由。 下面的推导表明,该界对运行器生成的每类操作数都成立,并且也能精确容纳每个操作数,因此以 pp 解析从不舍入。

记 dℓ,dr,dtd_\ell, d_r, d_t 为操作数位数,m=max⁡(dℓ,dr)+∣sℓ−sr∣m = \max(d_\ell, d_r) + |s_\ell - s_r|。运行器使用如下操作数:Add、Subtract、Multiply、Fma、Quantize 和 Compare 取标度来自配置 (0,2)(0, 2)、(8,18)(8, 18)、(28,0)(28, 0) 的生成操作数;Divide、DivideInteger 和 Remainder 用标度为 00、88 或 2828 的生成操作数除以 22、88 或 2525;Power 求平方(指数 22);SquareRoot 取一个生成整数的平方;Quantize 使用量子 10−sℓ10^{-s_\ell},Rescale 使用指数 −sℓ-s_\ell,ScaleB 使用移位 00;其余一元运算取整数。

Add/Subtract:∣cℓ10s−sℓ±cr10s−sr∣<10m+10m≤10m+1⇒ ≤m+1 digits, p=m+4Multiply:∣cℓcr∣<10dℓ+dr⇒ ≤dℓ+dr, p=dℓ+dr+4Power (n=2):∣cℓ2∣<102dℓ⇒ ≤2dℓ, p=2dℓ+4Fma:∣cℓcr10σ−sℓ−sr+ct10σ−st∣<10m′+1, m′=max⁡(dℓ+dr,dt)+∣sℓ+sr−st∣⇒ p=dℓ+dr+dt+∣sℓ+sr−st∣+4≥m′+1Divide (b∈{2,8,25}):12=5⋅10−1, 18=125⋅10−3, 125=4⋅10−2⇒ ≤dℓ+3, p=dℓ+dr+4≥dℓ+5DivideInteger:∣trunc⁡(a/b)∣≤∣a∣<10dℓ−sℓ⇒ ≤dℓRemainder:∣r∣≤∣a∣, scale(r)=sℓ⇒ ≤dℓSquareRoot:cℓ=g2, d(g)≤⌈dℓ/2⌉⇒ p=dℓ+4Unary group:∣result coefficient∣≤∣cℓ∣⇒ ≤dℓ, p=dℓ+4Compare:result∈{−1,0,1}⇒ p=max⁡(dℓ,dr)+3\begin{aligned} \textbf{Add/Subtract:}\quad & |c_\ell 10^{s-s_\ell} \pm c_r 10^{s-s_r}| < 10^{m} + 10^{m} \le 10^{m+1} &&\Rightarrow\ \le m + 1 \text{ digits},\ p = m + 4 \\ \textbf{Multiply:}\quad & |c_\ell c_r| < 10^{d_\ell + d_r} &&\Rightarrow\ \le d_\ell + d_r,\ p = d_\ell + d_r + 4 \\ \textbf{Power}\ (n = 2):\quad & |c_\ell^2| < 10^{2 d_\ell} &&\Rightarrow\ \le 2 d_\ell,\ p = 2 d_\ell + 4 \\ \textbf{Fma:}\quad & |c_\ell c_r 10^{\sigma - s_\ell - s_r} + c_t 10^{\sigma - s_t}| < 10^{m'+1},\ m' = \max(d_\ell + d_r, d_t) + |s_\ell + s_r - s_t| &&\Rightarrow\ p = d_\ell + d_r + d_t + |s_\ell + s_r - s_t| + 4 \ge m' + 1 \\ \textbf{Divide}\ (b \in \{2, 8, 25\}):\quad & \tfrac{1}{2} = 5 \cdot 10^{-1},\ \tfrac{1}{8} = 125 \cdot 10^{-3},\ \tfrac{1}{25} = 4 \cdot 10^{-2} &&\Rightarrow\ \le d_\ell + 3,\ p = d_\ell + d_r + 4 \ge d_\ell + 5 \\ \textbf{DivideInteger:}\quad & |\operatorname{trunc}(a/b)| \le |a| < 10^{d_\ell - s_\ell} &&\Rightarrow\ \le d_\ell \\ \textbf{Remainder:}\quad & |r| \le |a|,\ \text{scale}(r) = s_\ell &&\Rightarrow\ \le d_\ell \\ \textbf{SquareRoot:}\quad & c_\ell = g^2,\ d(g) \le \lceil d_\ell / 2 \rceil &&\Rightarrow\ p = d_\ell + 4 \\ \textbf{Unary group:}\quad & |\text{result coefficient}| \le |c_\ell| &&\Rightarrow\ \le d_\ell,\ p = d_\ell + 4 \\ \textbf{Compare:}\quad & \text{result} \in \{-1, 0, 1\} &&\Rightarrow\ p = \max(d_\ell, d_r) + 3 \end{aligned}

表中每个 pp 都至少为 max⁡(dℓ,dr,dt)\max(d_\ell, d_r, d_t),因此操作数本身能被精确解析。之所以要对 Power 和 Fma 另行处理,是因为单靠 working_precision 得到 dℓ+dr+2d_\ell + d_r + 2,对指数 22 来说只有 dℓ+3d_\ell + 3,不足以容纳平方。

该约定只对这些操作数类别得到证明,而不是对任意输入。对一般的有限除数 q=2αq = 2^{\alpha},商的系数为 cℓ⋅5αc_\ell \cdot 5^{\alpha},约有 dℓ+0.699αd_\ell + 0.699\alpha 位,而界只增长 dr≈0.301α+1d_r \approx 0.301\alpha + 1。在 α=13\alpha = 13 时界失效:1/8192=0.00012207031251 / 8192 = 0.0001220703125 需要 1010 位,而 p=1+4+4=9p = 1 + 4 + 4 = 9。这样的夹具不会被悄悄接受;如下一节所示,它会被预言机比较拒绝。

计时范围

OperationOnly(arithmetic_only)对计时前已解析的操作数计时一次公开运算。FullPath(full_path)还会在准备好的上下文中解析规范操作数字符串,这正是接收文本的调用者需要付出的代价。两个范围运行相同的数据集、具有相同的指纹,因此它们的数字描述的是相同的输入。上下文构造、规范化、校验和报告都位于两个范围之外。

配对统计

对每个数据集 jj、重复 rr 和分块 bb,Mare Mark 为每个实现记录一个经校准的延迟。本包把共享 (j,r,b)(j, r, b) 的 DzmingLi 与 GDA 样本配对,并构造差值

Δi=tiGDA−tiDZ.\Delta_i = t^{\mathrm{GDA}}_i - t^{\mathrm{DZ}}_i .

若把样本建模为 t=μimpl+βb+εt = \mu_{\mathrm{impl}} + \beta_b + \varepsilon,其中 βb\beta_b 是该分块共享的漂移(频率变化、缓存状态),则 Δi=μGDA−μDZ+(ε−ε′)\Delta_i = \mu_{\mathrm{GDA}} - \mu_{\mathrm{DZ}} + (\varepsilon - \varepsilon'):分块效应被抵消。BalancedBlocks 轮换先运行的实现,因此顺序效应平均而言也会抵消。报告的量为

δ=100⋅median⁡iΔimedian⁡itiDZ %,speedup=median⁡itiGDAmedian⁡itiDZ,\delta = 100 \cdot \frac{\operatorname{median}_i \Delta_i}{\operatorname{median}_i t^{\mathrm{DZ}}_i}\ \%, \qquad \text{speedup} = \frac{\operatorname{median}_i t^{\mathrm{GDA}}_i}{\operatorname{median}_i t^{\mathrm{DZ}}_i},

当 δ≤−3\delta \le -3 时判定为 gda_faster,当 δ≥3\delta \ge 3 时为 dzmingli_faster,否则为 equivalent。中位数的崩溃点为 50 %,因此 Mare Mark 只报告而不剔除的离群值无法使其大幅偏移。3 %3\,\% 阈值是实际显著性的边际,而不是假设检验,报告也不含置信区间。11 Mare Mark 的 compare_paired 把 Δi\Delta_i 的四分位距存为其区间。它描述的是配对差的离散程度,而不是中位数的不确定性。

每个规模三个数据集、每个数据集 20 次确认性重复,一个有效规模就有 60 对。

特定运算的规模上限

DzmingLi 的 digit_count 在 32 位 Int 中计算 bit_length * 30103。乘积在下述条件下溢出

bit_length>231−130103≈71 337⟺digits≳71 337⋅log⁡102≈21 475.\text{bit\_length} > \frac{2^{31} - 1}{30103} \approx 71\,337 \quad\Longleftrightarrow\quad \text{digits} \gtrsim 71\,337 \cdot \log_{10} 2 \approx 21\,475 .

两个 nn 位操作数之积至多有 2n2n 位,因此乘法、FMA 和平方从 n≈10 738n \approx 10\,738 起就可能越过该界限。因此规模扩展运行器让所有运算止步于 10,000 位,只有 add、subtract、divide 和 compare 运行到 16,384 和 20,000 位。这些是由溢出条件得出的解析界;在 32,768 和 65,536 位上的运行复现了中止。

正确性与不变量

规范形式。 对 c≠0c \ne 0,normalize 返回满足 s′≥0s' \ge 0 且(s′=0s' = 0 或 10∤c′10 \nmid c')的 (c′,s′)(c', s'),并且 v(c′,s′)=v(c,s)v(c', s') = v(c, s),因为每一步都把 (10q,s)(10q, s) 换成 (q,s−1)(q, s - 1)。两个这样的对只有在相等时才表示同一个数:

c110−s1=c210−s2, s1<s2 ⇒ c2=c110 s2−s1 ⇒ 10∣c2 and s2>0,\begin{aligned} c_1 10^{-s_1} = c_2 10^{-s_2},\ s_1 < s_2 &\ \Rightarrow\ c_2 = c_1 10^{\,s_2 - s_1} \\ &\ \Rightarrow\ 10 \mid c_2 \text{ and } s_2 > 0 , \end{aligned}

这与 (c2,s2)(c_2, s_2) 的规范形式矛盾;而 s1=s2s_1 = s_2 迫使 c1=c2c_1 = c_2。canonical_string 写出符号、∣c′∣|c'| 的各位数字和小数点位置 s′s',它们都由数本身决定,因此字符串相等即数相等。

预言机除法。 N′/M′N'/M' 是既约分数。若 M′∣N′10kM' \mid N' 10^{k},则由于 gcd⁡(M′,N′)=1\gcd(M', N') = 1,有 M′∣10k=2k5kM' \mid 10^{k} = 2^k 5^k,所以 M′M' 不含 22 和 55 以外的素因子。反之,若 M′=2α5βM' = 2^{\alpha} 5^{\beta},则最小的这样的 kk 为 max⁡(α,β)\max(\alpha, \beta)。预言机在循环之前先检查因子分解,所以循环必然终止,且返回的标度是最短的精确标度。

整数平方根。 预言机迭代

xk+1=⌊xk+⌊n/xk⌋2⌋,x0=10⌈D/2⌉>n,x_{k+1} = \Bigl\lfloor \frac{x_k + \lfloor n / x_k \rfloor}{2} \Bigr\rfloor, \qquad x_0 = 10^{\lceil D/2 \rceil} > \sqrt{n},

其中 DD 是 nn 的位数。由算术-几何平均不等式 (x+n/x)/2≥n(x + n/x)/2 \ge \sqrt{n},每个迭代值都至少为 ⌊n⌋\lfloor \sqrt{n} \rfloor;当 xk>⌊n⌋x_k > \lfloor \sqrt{n} \rfloor 时有 n/xk<xkn / x_k < x_k,从而 xk+1<xkx_{k+1} < x_k。序列严格递减直到达到 ⌊n⌋\lfloor \sqrt{n} \rfloor,循环在此停止。随后预言机要求 x2=nx^2 = n,因此非平方输入会中止,而不是产生一个舍入后的根。

零容差的可靠性。 设 xx 为精确结果,pp 为夹具精度。

  1. 若 xx 至多有 pp 位有效数字,符合规范的 GDA 运算返回一个等于 xx 的数,因为在任何舍入模式下 xx 都是它自己的正确舍入值。于是 κ(I(x))=κ(O(x))\kappa(I(x)) = \kappa(O(x)):正确的实现永远不会被拒绝。
  2. 若 xx 多于 pp 位,向零舍入给出 ∣round⁡p(x)∣<∣x∣|\operatorname{round}_p(x)| < |x|,于是规范字符串不同,校验失败:精度不足永远不会被接受。
  3. 规范字符串相等即数相等,因此错误的结果永远不会被接受。

精度约定为每一类生成的操作数确立了 (1) 的前提。在这些类别之外 (2) 仍然成立,这使检查是失败即关闭(fail-closed)的。

确定性。 操作数由 Mare Mark 带种子的 derive_seed(seed, "<operation>:<digits>", profile) 导出,从不依赖时钟,因此相同的种子复现相同的语料和相同的指纹。generate_decimal 只使用数字 1..91..9,因此生成的系数恰好具有所要求的位数且没有末尾零。

预言机的代价。 对 nn 位系数,加法的主要开销是乘以 10∣sℓ−sr∣10^{|s_\ell - s_r|} 的对齐,乘法是一次 BigInt 乘积,除法是一次最大公约数加上 max⁡(α,β)\max(\alpha, \beta) 次乘以 1010,其中 M′=2α5βM' = 2^{\alpha} 5^{\beta}。所有这些都在计时区域之外运行。

被否决的方案

  • 只做成对一致性检查。 否决:它无法归因分歧,也会漏掉共同的错误。
  • 以一个库作为另一个库的预言机。 否决:GDA 库本身就是被测对象之一,而且基准发现的 DzmingLi 失败届时将无法与 GDA 失败区分。
  • 带容差的比较(ulp 或相对误差)。否决:每个被测结果在十进制下都是精确的,所以任何非零容差都只会掩盖错误。
  • 固定的大精度,例如 10510^5 位。否决:它可能使除法工作量超出结果所需,并把计时与一个任意常数绑定。
  • 每次运行使用随机操作数。 否决:结果必须可复现,指纹在不同计时范围之间也必须稳定。
  • 为 exp、ln 和 log10 做规模扩展基准。 否决:恒等夹具无法触及它们的算法,而一般结果需要一个独立的高精度超越函数预言机。它们只由 decTest 审计覆盖。

边界

  • 本包只检查数值。指数同类、结果的末尾零、状态标志、NaN、无穷大和带符号零都不在比较范围内;decTest 审计覆盖它们。
  • 除法只用有限除数 22、88 和 2525 测试。循环小数商和大分母不在基准范围内,预言机也会拒绝循环小数商。
  • Parse 与 Format 是恒等路径。它们不测试 GDA 的格式化;decTest 审计单独记录了 DzmingLi 的 329 个 toSci 失败。
  • OperandShape 只是标签;语料是每个规模三个确定性的标度配置,而不是随机分布。
  • 精度约定只对上面列出的生成操作数类别得到证明,而不是对任意输入。
  • 可执行程序只在 native 目标上运行;其他目标打印一条消息后退出。
  • DzmingLi/decimal@0.2.2 已被弃用,改为 moonbit-community/decimal;基准为了历史对比而固定使用它。

Footnotes

  1. Mare Mark 的 compare_paired 把 Δi\Delta_i 的四分位距存为其区间。它描述的是配对差的离散程度,而不是中位数的不确定性。 ↩