decimal_gda 设计

本页解释 decimal_gda 所实现的数学,以及该包为何如此构建。教程展示了如何使用它,API 页面列出了所有名称。

设计目标

decimal_gda 以纯 MoonBit 值实现了 M. F. Cowlishaw 的 General Decimal Arithmetic Specification(1.70 版)11 M. F. Cowlishaw, General Decimal Arithmetic Specification, version 1.70 (2009), https://speleotrove.com/decimal/decarith.html。测试套件是同一作者的 dectest 集合,2.62 版。 。按规范,一个 GDA 运算产生的不只是一个数:它产生一个结果、一组条件、对上下文粘滞状态的更新,以及可能的陷阱——陷阱转移控制流,但有定义的结果仍然可用。目标是精确且可观察地建模所有这些,从而

  • 每个运算都是从(操作数,上下文)到(结果,引发的条件,下一个上下文,陷阱判定)的全函数,没有隐藏状态;
  • 结果与规范逐位一致,并以固定版本的官方测试套件检验(参见符合性);
  • 本包独立于 IEEE 754 包 decimal,因此一方契约的变化不会泄漏到另一方。

数学背景

数、同值类与调整指数

有限的 GDA 数是一个三元组 (s,c,e)(s, c, e),其中符号为 s∈{0,1}s \in \{0, 1\},系数为 c∈Nc \in \mathbb{N},指数为 e∈Ze \in \mathbb{Z},表示

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

映射 vv 不是单射:(0,250,−2)(0, 250, -2) 和 (0,25,−1)(0, 25, -1) 都表示 2.52.5。值相同的三元组构成一个同值类(cohort),GDA 保留三元组而不只是值,因为指数是有意义的(量子 10e10^e:“2.50”是精确到分的测量值)。零带符号,因此 (1,0,e)(1, 0, e) 是 −0-0。对 c>0c > 0,记 d(c)d(c) 为 cc 的十进制位数,并记

e^=e+d(c)−1\hat e = e + d(c) - 1

为调整指数,即首位数字的指数:10e^≤∣v∣<10e^+110^{\hat e} \le |v| < 10^{\hat e + 1}。除有限数外,还有 ±∞\pm\infty,以及携带符号和载荷 cc 的静默 NaN 与信号 NaN。

上下文与可表示集合

上下文固定一个精度 p≥1p \ge 1、一种舍入模式、一个指数范围 emin⁡≤emax⁡e_{\min} \le e_{\max} 以及一个钳制位。有限数在该上下文中可表示,当

d(c)≤p,e^≤emax⁡,e≥Etiny   where   Etiny=emin⁡−p+1,d(c) \le p, \qquad \hat e \le e_{\max}, \qquad e \ge E_{\mathrm{tiny}} \;\text{ where }\; E_{\mathrm{tiny}} = e_{\min} - p + 1,

且若开启钳制,还要求 e≤Etop=emax⁡−p+1e \le E_{\mathrm{top}} = e_{\max} - p + 1。值 EtinyE_{\mathrm{tiny}} 是最小单位的指数:e^=emin⁡\hat e = e_{\min} 且有 pp 位数字的数满足 e=emin⁡−p+1e = e_{\min} - p + 1,而低于 10emin⁡10^{e_{\min}} 的数(次正规范围)保持该单位,同时丢失前导数字。EtopE_{\mathrm{top}} 是满足 e^=emax⁡\hat e = e_{\max} 的最大 pp 位数的指数;钳制使指数集合恰好等于 IEEE 交换格式所能编码的集合,这就是 decimal32/64/128 预设开启钳制的原因。

算术规则:先求精确结果,再舍入一次

GDA 中的每个算术运算都由同样的两步规则定义:

  1. 计算精确的数学结果 xx,并在表示 xx 的三元组中选择指数最接近该运算理想指数 e∗e^\ast 的那一个;
  2. 若该三元组不可表示,则将其舍入到上下文,并引发描述所发生情况的条件。

由于第 1 步是精确的,第 2 步就是单次舍入,下文的每个界都是一次舍入的界。本包严格遵循该规则:每个运算从其操作数构造精确的系数和指数,然后调用同一个共享的最终化例程(normalize_decimal_parts_ctx_repr 之后接 finalize_finite_ctx),它是唯一对有限结果进行舍入、或引发与精度或范围相关条件的代码。

理想指数

理想指数是仅用操作数中已有的信息就能表示精确结果的最粗量子。推导很简短:

x1±x2=(c110e1−m±c210e2−m)10m,m=min⁡(e1,e2),x1×x2=(c1c2) 10e1+e2,x1÷x2=c1c2 10e1−e2,x=c 10e−2⌊e/2⌋  10⌊e/2⌋.\begin{aligned} x_1 \pm x_2 &= \bigl(c_1 10^{e_1 - m} \pm c_2 10^{e_2 - m}\bigr) 10^{m}, & m &= \min(e_1, e_2),\\ x_1 \times x_2 &= (c_1 c_2)\, 10^{e_1 + e_2}, & &\\ x_1 \div x_2 &= \frac{c_1}{c_2}\, 10^{e_1 - e_2}, & &\\ \sqrt{x} &= \sqrt{c\,10^{e - 2\lfloor e/2 \rfloor}}\;10^{\lfloor e/2 \rfloor}. & & \end{aligned}

对于加法,mm 是使两个操作数都为整数的最大指数,因此也是能保证和为量子整数倍的最大指数;括号内就是一个精确的整数系数。对于乘法,系数 c1c2c_1 c_2 在指数 e1+e2e_1 + e_2 下为整数,一般而言在更大的指数下则不是。对于除法,商 c1/c2c_1 / c_2 不一定是整数,因此只有在可达时才使用理想指数 e1−e2e_1 - e_2:精确的商以最接近 e1−e2e_1 - e_2、且保持整数系数不超过 pp 位的指数写出,不精确的商则用满 pp 位。对于平方根,⌊e/2⌋\lfloor e/2 \rfloor 是这样一个指数:它的平方是不超过 ee 的最大偶数指数。下表总结了本包所用的规则:

操作理想指数 e∗e^\ast
add、subtract、fma(求和部分)min⁡(e1,e2)\min(e_1, e_2)
multiply、fma(乘积部分)e1+e2e_1 + e_2
divide精确时为 e1−e2e_1 - e_2,否则为 pp 位
sqrt精确时为 ⌊e/2⌋\lfloor e/2 \rfloor,否则为 pp 位
remainder, remainder_nearmin⁡(e1,e2)\min(e_1, e_2)
divide_integer, to_integral_*整数值为 max⁡(e,0)\max(e, 0),商为 0
quantize, rescale所请求的指数
scalebe+ne + n
整数 nn 的 power,精确nen e
exp、ln、log10、非整数 powerpp 位(总是不精确,exp(0)、ln(1)、log10(10^k) 除外)
plus, minus, abs, applyee

零结果没有首位数字,因此其指数就是理想指数,钳制到 [Etiny,emax⁡][E_{\mathrm{tiny}}, e_{\max}](或 EtopE_{\mathrm{top}})之内。精确零和的符号为 ++,除非两个操作数都为负,或者模式为 Floor 且至少一个为负;这是 IEEE 规则 x−x=+0x - x = +0(向 −∞-\infty 舍入时除外)的十进制形式。

///|
test "ideal exponents" {
  let ctx = @decimal_gda.GdaContext::decimal64()
  let d = (s : String) => @decimal_gda.Decimal::from_string(s).unwrap()
  inspect(@decimal_gda.add(d("1.30"), d("1.2"), ctx).value(), content="2.50")
  inspect(@decimal_gda.multiply(d("1.30"), d("1.2"), ctx).value(), content="1.560")
  inspect(@decimal_gda.divide(d("2.40"), d("2"), ctx).value(), content="1.20")
  inspect(@decimal_gda.divide(d("1"), d("4"), ctx).value(), content="0.25")
  inspect(@decimal_gda.sqrt(d("1.00"), ctx).value(), content="1.0")
  inspect(@decimal_gda.subtract(d("1.0"), d("1.00"), ctx).value(), content="0.00")
}

舍入到精度

设精确结果为 (−1)sc 10e(-1)^s c\,10^{e},其中 d(c)=p+kd(c) = p + k、k≥1k \ge 1。除去低 kk 位,得 c=q 10k+rc = q\,10^{k} + r,其中 0≤r<10k0 \le r < 10^k。每种舍入模式都返回带增量 δ∈{0,1}\delta \in \{0, 1\} 的 (−1)s(q+δ) 10e+k(-1)^s (q + \delta)\,10^{e+k},八种模式的区别仅在于 δ\delta(should_increment_decimal_repr):

δDown=0,δUp=[r>0],δCeiling=[r>0∧s=0],δFloor=[r>0∧s=1],δHalfUp=[2r≥10k],δHalfDown=[2r>10k],δHalfEven=[2r>10k∨(2r=10k∧q odd)],δZeroFiveUp=[r>0∧q≡0(mod5)].\begin{aligned} \delta_{\mathrm{Down}} &= 0, & \delta_{\mathrm{Up}} &= [r > 0],\\ \delta_{\mathrm{Ceiling}} &= [r > 0 \wedge s = 0], & \delta_{\mathrm{Floor}} &= [r > 0 \wedge s = 1],\\ \delta_{\mathrm{HalfUp}} &= [2r \ge 10^k], & \delta_{\mathrm{HalfDown}} &= [2r > 10^k],\\ \delta_{\mathrm{HalfEven}} &= [2r > 10^k \vee (2r = 10^k \wedge q \text{ odd})], & \delta_{\mathrm{ZeroFiveUp}} &= [r > 0 \wedge q \equiv 0 \pmod 5]. \end{aligned}

ZeroFiveUp 的含义是“向零舍入,除非保留的末位数字是 0 或 5”;其目的是使之后用任意模式舍入到更少位数时不受第一次舍入的影响。2r2r 与 10k10^k 的比较在十进制 limb 上精确进行(gda_coeff_div_pow10_round_info_repr 返回商、是否 r>0r > 0 以及 2r−10k2r - 10^k 的符号)。若 q+δ=10pq + \delta = 10^p,结果重新规范化为 10p−1⋅10e+k+110^{p-1} \cdot 10^{e+k+1}。在舍入之前,最终化例程会先尝试去掉 cc 的末尾零;若精确系数过长仅仅是因为零,则结果是精确的,只引发 Rounded。

设 u=10e+ku = 10^{e+k} 为结果 x^\hat x 的末位单位。由于在 half 模式下 ∣δ 10k−r∣≤10k/2|\delta\,10^k - r| \le 10^k/2,其他模式下 <10k< 10^k,

∣x^−x∣≤12u    (half modes),∣x^−x∣<u    (others),|\hat x - x| \le \tfrac12 u \;\;(\text{half modes}), \qquad |\hat x - x| < u \;\;(\text{others}),

又由于 c≥10p+k−1c \ge 10^{p+k-1} 给出 ∣x∣≥10p−1u|x| \ge 10^{p-1} u,

∣x^−x∣∣x∣≤12 101−p    (half modes),∣x^−x∣∣x∣<101−p    (others).\frac{|\hat x - x|}{|x|} \le \tfrac12\,10^{1-p} \;\;(\text{half modes}), \qquad \frac{|\hat x - x|}{|x|} < 10^{1-p} \;\;(\text{others}).

这就是标准模型 fl(x∘y)=(x∘y)(1+ε)\mathrm{fl}(x \circ y) = (x \circ y)(1 + \varepsilon)、∣ε∣≤u|\varepsilon| \le \mathbf u 中的十进制单位舍入误差 u=12101−p\mathbf u = \frac12 10^{1-p}。22 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.2。十进制的摆动比二进制大:在一个十进位段内相对间距的变化因子为 10,而二进制为 2(Goldberg 1991, §1.2)。 只要移除了 k≥1k \ge 1 位数字,舍入就引发 Rounded;只要移除的数字中有非零者(r>0r > 0),就引发 Inexact。

该界以及本页其他舍入事实的完整证明见附录:

decimal_gda 的舍入证明

次正规结果与 EtinyE_{\mathrm{tiny}}

若精确结果满足 e^<emin⁡\hat e < e_{\min},则称其为微小的。微小结果可能仍至多需要单位 10Etiny10^{E_{\mathrm{tiny}}},因此舍入位置不是“保留 pp 位”,而是“保留不低于 10Etiny10^{E_{\mathrm{tiny}}} 的数字”:对精确指数 e<Etinye < E_{\mathrm{tiny}},最终化例程移除 k=Etiny−ek = E_{\mathrm{tiny}} - e 位数字,可能全部移除。此时绝对误差以次正规单位为界,

∣x^−x∣≤12 10Etiny    (half modes),|\hat x - x| \le \tfrac12\,10^{E_{\mathrm{tiny}}} \;\;(\text{half modes}),

但相对误差不再以 u\mathbf u 为界:随着 ∣x∣|x| 从 10emin⁡10^{e_{\min}} 降到 10Etiny10^{E_{\mathrm{tiny}}},精度从 pp 位逐渐降到 1 位。条件记录了这一点:每个微小结果都引发 Subnormal(GDA 以及本包的 GDA 函数在舍入之前依据精确的 e^\hat e 检测微小性),微小且不精确时引发 Underflow,舍入为零时引发 Clamped(此时零取指数 EtinyE_{\mathrm{tiny}})。无状态层还额外提供舍入后检测(DecimalTininessDetection::AfterRounding),供 IEEE 风格使用;它只改变哪些结果算作微小,从不改变值。

钳制

在 clamp 下,e>Etope > E_{\mathrm{top}} 但 e^≤emax⁡\hat e \le e_{\max} 的结果作为值可表示,但作为三元组不可表示。由于 e^≤emax⁡\hat e \le e_{\max} 意味着 d(c)+e−1≤emax⁡d(c) + e - 1 \le e_{\max},用 e−Etope - E_{\mathrm{top}} 个零填充系数可得

d(c 10e−Etop)=d(c)+e−Etop≤emax⁡+1−Etop=p,d\bigl(c\,10^{e - E_{\mathrm{top}}}\bigr) = d(c) + e - E_{\mathrm{top}} \le e_{\max} + 1 - E_{\mathrm{top}} = p,

因此填充后的三元组至多有 pp 位,且值不变。本包正是这样做的,并引发 Clamped;值不变,只有同值类成员不同。零以同样的方式钳制,即把其指数移入 [Etiny,Etop][E_{\mathrm{tiny}}, E_{\mathrm{top}}]。

上溢

若舍入后的结果满足 e^>emax⁡\hat e > e_{\max},则运算上溢:引发 Overflow、Inexact 和 Rounded,结果为 ±∞\pm\infty 或同号的最大有限数 Nmax⁡=(10p−1)⋅10EtopN_{\max} = (10^p - 1) \cdot 10^{E_{\mathrm{top}}}。取哪一个,可以把 ∞\infty 视为 Nmax⁡N_{\max} 之外的可表示数,再按该模式自身的方向得出:各 half 模式和 Up 远离零,因此到达 ∞\infty;Down 和 ZeroFiveUp 朝向零(ZeroFiveUp 只在末位为 0 或 5 时才远离零舍入,而 Nmax⁡N_{\max} 的末位是 9),因此停在 Nmax⁡N_{\max};Ceiling 对正结果给出 +∞+\infty,对负结果给出 −Nmax⁡-N_{\max},Floor 则与之镜像(overflow_to_infinity):

模式正向上溢负向上溢
HalfEven, HalfUp, HalfDown, Up+∞+\infty−∞-\infty
Down, ZeroFiveUp+Nmax⁡+N_{\max}−Nmax⁡-N_{\max}
Ceiling+∞+\infty−Nmax⁡-N_{\max}
Floor+Nmax⁡+N_{\max}−∞-\infty

条件、信号与陷阱

GDA 区分条件(发生了什么:除以零、结果被舍入、转换格式错误、…)、信号(条件所触发的具名事件)和陷阱(用户要求中断计算的信号)。本包为每个条件保留一个标志,并按同一条规则把条件映射到 GDA 信号:四个细分的无效条件 ConversionSyntax、DivisionImpossible、DivisionUndefined 和 InvalidContext 都发出 InvalidOperation 信号。形式化地,设 S\mathcal S 为十三个标志,ι⊂S\iota \subset \mathcal S 为五个无效族标志。对标志集合 RR 和信号 σ\sigma,

σ∈∗R  ⟺  {R∩ι≠∅σ=InvalidOperation,σ∈Rotherwise,\sigma \in^\ast R \iff \begin{cases} R \cap \iota \neq \emptyset & \sigma = \mathrm{InvalidOperation},\\ \sigma \in R & \text{otherwise,} \end{cases}

即 GdaFlags::contains。于是一个运算就是如下状态机

(x⃗,C)  ⟼  {Completed(v,C,∅)R=∅,Completed(v,C′,R)R≠∅, τ=⊥,Trapped(τ,v,C′,R)τ≠⊥,(\vec x, C) \;\longmapsto\; \begin{cases} \mathrm{Completed}(v, C, \emptyset) & R = \emptyset,\\ \mathrm{Completed}(v, C', R) & R \neq \emptyset,\ \tau = \bot,\\ \mathrm{Trapped}(\tau, v, C', R) & \tau \neq \bot, \end{cases}

其中 (v,R)=F(x⃗,π(C))(v, R) = F(\vec x, \pi(C)) 是仅由操作数和策略 π(C)\pi(C)(精度、舍入、指数范围、钳制、extended)计算出的结果与引发的标志,

C′.status=C.status∪R∪{InvalidOperation∣R∩ι≠∅},C′.π=C.π,C′.traps=C.traps,C'.\mathrm{status} = C.\mathrm{status} \cup R \cup \bigl\{\mathrm{InvalidOperation} \mid R \cap \iota \neq \emptyset\bigr\}, \qquad C'.\pi = C.\pi,\quad C'.\mathrm{traps} = C.\mathrm{traps},

而 τ\tau 是优先级列表 InvalidOperation、DivisionByZero、DivisionUndefined、DivisionImpossible、InvalidContext、ConversionSyntax、Overflow、Underflow、Subnormal、Inexact、Rounded、Clamped、LostDigits 中第一个满足 τ∈∗R\tau \in^\ast R 且 τ∈C.traps\tau \in C.\mathrm{traps} 的信号,若不存在则为 ⊥\bot(complete_gda 和 trapped_signal)。

由定义可直接得出三条性质,它们正是用户所依赖的:

  1. 值不依赖于状态或陷阱。 vv 和 RR 是 (x⃗,π(C))(\vec x, \pi(C)) 的函数;状态和陷阱只进入 C′C' 和 τ\tau。因此启用陷阱永远不会改变结果,并且 Trapped 可以携带有定义的结果。
  2. 状态单调且幂等。 C.status⊆C′.statusC.\mathrm{status} \subseteq C'.\mathrm{status},又由于 ∪\cup 满足结合律、交换律和幂等律,串联执行一系列运算之后的状态就是所有引发集合的并,与该序列如何分组无关。
  3. 陷阱的选择是确定性的。 优先级列表是信号上的全序,因此一个运算至多选择一个陷阱,无论其条件以何种顺序被检测到。
///|
test "status is the union of the raised sets" {
  let ctx = @decimal_gda.context(precision=3)
  let d = (s : String) => @decimal_gda.Decimal::from_string(s).unwrap()
  let a = @decimal_gda.divide(d("1"), d("3"), ctx) // Inexact, Rounded
  let b = @decimal_gda.divide(d("1"), d("0"), a.next_context()) // DivisionByZero
  let expected = a.raised().combine(b.raised())
  inspect(b.next_context().status() == expected, content="true")
  // Trapping changes the variant, never the value.
  let trapping = ctx.trap(Inexact)
  let t = @decimal_gda.divide(d("1"), d("3"), trapping)
  inspect(t.value() == a.value(), content="true")
  inspect(t is @decimal_gda.GdaOutcome::Trapped(Inexact, _, _, _), content="true")
}

设计决策

通过 GdaOutcome 传递状态,而非全局状态

问题。 GDA 规范把上下文描述为一个可变对象,运算会设置其状态标志。大多数实现(decNumber、Python 的 decimal)为每个线程保留一个当前上下文。

方案。 (a) 按引用传递的可变上下文;(b) 线程局部或全局的当前上下文;(c) 随每个结果一起返回的不可变上下文。

选择:(c)。 每个 GDA 函数都返回 GdaOutcome,调用者把 next_context() 传给下一个运算。理由就是上述性质。由于 FF 只读取策略,状态机分解为纯数值部分和纯簿记部分,两者都是引用透明的:表达式可以重新求值、记忆化或在另一个线程上运行,而不改变其结果或标志。测试运行器可以对上下文做快照,并在其下多次运行同一运算;frontend/gda_expr 中的 .decTest 运行器正是这样做的。MoonBit 也没有线程局部存储,所以 (b) 意味着进程全局状态,而这是 Luna-Flow 所排除的。代价是需要显式传递;decimal_gda_checked 包为线性流水线消除了这一负担。当没有引发任何条件时,输入上下文被原样返回,因此精确运算不会分配新的上下文。

保留陷阱的有定义结果

在 GDA 中,陷阱会转移控制流,但规范仍然定义了该运算本应交付的结果。把陷阱表示为错误(Result::Err、raise)会丢弃这个值。Trapped 保留了值、下一个上下文以及引发集合,调用者可以检查、记录它们或从中恢复;性质 1 保证它与没有陷阱时的值相同。

保留细分的无效条件

规范通过 InvalidOperation 报告 ConversionSyntax、DivisionImpossible、DivisionUndefined 和 InvalidContext。本包把它们保留为独立的标志,同时使 contains(InvalidOperation) 对每一个都为真,并在其中任何一个被引发时在状态中设置 invalid_operation。只了解八种 GDA 信号的程序看到的正是 GDA 的行为,而测试框架或诊断仍能区分格式错误的字面量与 0/00/0。把 InvalidOperation 放在优先级列表首位,使 InvalidOperation 陷阱能捕获全部四种,正如规范所要求的。

单一最终化例程,先处理特殊值

每个运算先判定特殊情形(NaN 传播、无效运算、无穷、精确零),然后构造精确的有限系数和指数,最后才调用最终化例程。没有任何系数算法会设置标志。这使条件始终是精确结果与策略的函数,与所选的乘法或除法内核无关,也使内核可以调优而不触及面向标准的行为。

独立于 decimal 的 GDA 包

IEEE 754-2008 十进制算术源于 GDA,因此两者在大多数有限结果上一致,但它们的契约不同:

方面decimal_gda(GDA 1.70)decimal(IEEE 754-2019)
精度每个上下文可取任意 p≥1p \ge 1格式的 pp(或自选)
舍入模式八种,包括 HalfUp、HalfDown、ZeroFiveUpIEEE 属性
条件上下文中的粘滞状态,带有定义结果的陷阱每个运算的标志随值一同返回
细分的无效条件、LostDigits、子集算术是否
微小性舍入前可选
初等函数sqrt, exp, ln, log10, powerIEEE 推荐的集合
交换格式DPDDPD 和 BID

共享的核心必须承载两种状态模型的并集,为一个标准所作的修改可能会悄然改变另一个标准的结果。因此本包拥有自己的值类型、系数内核、上下文、最终化例程和 DPD 编解码器,生产依赖扫描会检查它从不导入 decimal。目前系数阈值与 decimal 相同;这是测量得出的巧合,而非共享的依赖。

无状态层(DecimalContext、DecimalFlags、*_ctx 方法)是本包自身的引擎对外公开的部分。它像 IEEE 那样按运算返回标志,便于适配器使用(Luna-Flow/arithmetic trait 实现就使用它),但它实现的是 GDA 算术;IEEE 754 契约位于 decimal 中。

已证明等价的快速路径

有两类捷径绕开了通用机制,且不改变任何可观察的结果。

小的精确整数。 parse、add、subtract、multiply 和 fma 先检查:所有操作数是否都是指数为 0、系数小于 101810^{18} 的整数,精确结果是否为小于 101810^{18}、至多 pp 位的非零整数,其调整指数是否位于 [emin⁡,emax⁡][e_{\min}, e_{\max}] 内,钳制是否允许指数 0,以及上下文是否为 extended。在这些谓词成立时,通用路径不做舍入、不引发条件,并返回指数 0(这些运算的理想指数均为 0),因此捷径返回 Completed(v, C, none) 及相同的 vv。零结果被排除在外,因为其符号取决于舍入模式。

被吸收的加数。 当 x1x_1 与远小于它的 x2x_2 相加且 e^2<e^1−p+1\hat e_2 < \hat e_1 - p + 1 时,无需构造精确和。x2x_2 的每一位都低于 x1x_1 保留的最后一位,因此 x2x_2 只通过其符号及其与半个单位的比较来影响舍入:

x^=round(x1+x2)=(q+δ(sign⁡x2,  ∣x2∣≶12u))u,\hat x = \mathrm{round}\bigl(x_1 + x_2\bigr) = \bigl(q + \delta(\operatorname{sign} x_2,\; |x_2| \lessgtr \tfrac12 u)\bigr) u,

这就是上文的舍入表,只是把 rr 换成了属于同一比较类的粘滞代表值。exact_base_small_addend_result 正是这样求值的,其中 compare_magnitude_to_half_ulp 精确比较 x2x_2 与 u/2u/2。x1x_1 为十的幂且 x2x_2 符号相反的特殊情形(x1x_1 之下的单位小十倍)被排除,走通用路径。

除法

divide 单独处理 00、∞\infty 以及除以十的幂和小的精确除数的精确除法。否则它把被除数乘以 10t10^{t},其中 t=p+d(c2)+2t = p + d(c_2) + 2,使 Q=c110t/c2≥10p+2Q = c_1 10^{t} / c_2 \ge 10^{p+2} 至少有 p+3p + 3 位整数,按上下文模式把 QQ 舍入为整数 qq,再把 qq 舍入到 pp 位。商是否有限是精确判定的(has_finite_decimal_expansion_repr:c1/c2c_1/c_2 有限当且仅当 c2/gcd⁡(c1,c2)c_2 / \gcd(c_1, c_2) 没有 2 和 5 以外的素因子),据此为精确的商选择理想指数的同值类成员,对其余的商则强制取 Inexact 及 pp 位系数。

对定向模式,两次相继舍入等价于一次,因为截断可以复合:⌊⌊y/10j⌋/10m⌋=⌊y/10j+m⌋\lfloor \lfloor y/10^j \rfloor / 10^m \rfloor = \lfloor y/10^{j+m} \rfloor。对 half 模式,两者等价,除非第一次舍入制造出平局:qq 被丢弃的数字恰为 50⋯050\cdots0,而 Q≠qQ \ne q。

平方根

sqrt 先尝试精确根:去掉末尾零,使指数为偶数,用牛顿迭代 a←⌊(a+⌊c/a⌋)/2⌋a \leftarrow \lfloor (a + \lfloor c/a \rfloor)/2 \rfloor 求整数平方根 (s,ρ)(s, \rho)(s2+ρ=cs^2 + \rho = c),并在 ρ=0\rho = 0 时接受,然后在 pp 位以内向理想指数 ⌊e/2⌋\lfloor e/2 \rfloor 填充。否则直接在最终位置舍入。取根的指数为 f=max⁡(Etiny,⌊e^/2⌋−p+1)f = \max(E_{\mathrm{tiny}}, \lfloor \hat e/2 \rfloor - p + 1),它计算 s=⌊c 10e−2f⌋s = \lfloor \sqrt{c\,10^{e - 2f}} \rfloor,并通过比较被开方数与中点的平方来决定是否进位,这在整数中是精确的:

N≷s+12  ⟺  4N≷(2s+1)2.\sqrt{N} \gtrless s + \tfrac12 \iff 4N \gtrless (2s + 1)^2 .

在 ff 处只舍入一次(而不是先舍入到 pp 位、再在 EtinyE_{\mathrm{tiny}} 处舍入)避免了微小根的双重舍入。由于精确根已先行处理,比较只会用于判定无理根,而无理根不可能等于中点。该 GDA 函数总是采用 half-even 舍入。

整数次幂

对整数指数 nn,结果按 GDA 规范的规定计算:若精确幂能放入 pp 位,则以指数 nen e 精确返回;否则以工作精度 w=p+d(∣n∣)+2w = p + d(|n|) + 2(子集上下文中少一位)进行二进制幂运算,每次乘积后作 half-even 舍入,nn 为负时从 1/x1/x 开始,最后的乘积再舍入到上下文。∣n∣|n| 个因子中的每一个至多经过 ∣n∣−1|n| - 1 次舍入,因此工作结果为 xn(1+θ)x^n(1 + \theta),其中22 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.2。十进制的摆动比二进制大:在一个十进位段内相对间距的变化因子为 10,而二进制为 2(Goldberg 1991, §1.2)。

∣θ∣≤(1+uw)∣n∣−1−1≈∣n∣ uw<10d(∣n∣)⋅12 101−p−d(∣n∣)−2=12 10−1−p,|\theta| \le (1 + \mathbf u_w)^{|n|-1} - 1 \approx |n|\,\mathbf u_w < 10^{d(|n|)} \cdot \tfrac12\,10^{1 - p - d(|n|) - 2} = \tfrac12\,10^{-1-p},

它至多是结果末位的 0.050.05 个单位。加上最后一次舍入,在 half 模式下不精确的整数次幂与真值相差在末位 0.550.55 个单位以内。这接近但并非正确舍入,而规范对整数次幂也没有更高的要求。指数大到结果必然上溢或下溢的情形,会预先通过结果调整指数的界 n(e^+1)−1n(\hat e + 1) - 1 和 ne^n \hat e 检测出来。

正确舍入的 exp、ln、log10 与非整数 power

这些函数通过经认证的区间算术和 Ziv 式细化循环求值。33 A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991。关于球算术,参见 J. van der Hoeven, “Ball arithmetic”, 2009,以及 F. Johansson, “Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”, IEEE Trans. Computers 66(8), 2017。 对 f(x)f(x):

  1. 先处理精确情形。 exp⁡0=1\exp 0 = 1、ln⁡1=0\ln 1 = 0、log⁡1010k=k\log_{10} 10^k = k、特殊操作数、定义域错误,以及 power 能精确判定的情形(经由 sqrt 的 x1/2x^{1/2}、十的幂、1y1^y、必然上溢或下溢)永远不会进入循环。
  2. 包络输入。 xx 以定向舍入(to_bin_float 朝 −∞-\infty 和 +∞+\infty)转换为 ww 位的二进制球 [x−,x+][x^-, x^+],因此精确地有 x∈[x−,x+]x \in [x^-, x^+]。
  3. 包络输出。 ball_float 在球上计算 ff 并返回 [L,U]∋f(x)[L, U] \ni f(x)(对 log10,为 ln⁡\ln 除以缓存的 ln⁡10\ln 10 包络;对 power,exp⁡(yln⁡x)\exp(y \ln x) 经由 ln⁡x\ln x 的包络计算)。
  4. 认证候选值。 十进制近似 y~\tilde y(球心,或十进制级数求值的结果)被舍入为 rr。设 (a,b)(a, b) 为 rr 周围舍入到 rr 的实数构成的开区间:对 half 模式是与两侧相邻值之间的两个中点,对定向模式是 (r,r+)(r, r^+) 或 (r−,r)(r^-, r)。其端点是精确的十进制数;每个端点都被包络在二进制中,当满足下式时接受该结果 a≤a+<L≤f(x)≤U<b−≤b,a \le a^+ < L \le f(x) \le U < b^- \le b , 这就证明了 round(f(x))=r\mathrm{round}(f(x)) = r。第二种测试把两个精确的二进分数端点 LL 和 UU 都舍入到上下文;由于每种舍入模式都是单调的,round(L)=round(U)=r\mathrm{round}(L) = \mathrm{round}(U) = r 同样证明了 round(f(x))=r\mathrm{round}(f(x)) = r。
  5. 细化。 否则增大工作精度,w←w+max⁡(32,⌊w/2⌋)w \leftarrow w + \max(32, \lfloor w/2 \rfloor),至多十二次。

起始精度为 w0=max⁡(128,4D+64)w_0 = \max(128, 4D + 64) 位,其中 D=max⁡(d(c),p)D = \max(d(c), p):每位十进制数字四个比特超过了 log⁡210≈3.32\log_2 10 \approx 3.32,因此输入在传递时不会丢失十进制信息,并且还剩 64 个保护位。对于宽且未钳制、满足 p≤64p \le 64 的上下文中 e^=0\hat e = 0 的参数,会先用约 103D+12\frac{10}{3} D + 12 位作一次更廉价的尝试,若无法认证则回退到 w0w_0。

只要 f(x)f(x) 本身不是可表示数或中点,循环就会终止,因为包络会收缩到一点。对 exp、ln 和 log10,在第 1 步之后这总是成立:由 Lindemann–Weierstrass 定理,对有理数 x≠0x \ne 0,exe^x 是超越数,因此对有理数 x≠1x \ne 1,ln⁡x\ln x 是无理数,而 log⁡10x\log_{10} x 只在十的整数次幂处为有理数;可表示数和中点都是有理数。若预算仍然耗尽(例如某个精确值恰为中点且未在第 1 步捕获的 power),try_* 方法返回认证失败,GDA 函数返回带 InvalidOperation 的 NaN,绝不返回未经证明的结果。

GDA 函数 exp、ln、log10 和 sqrt 传入一份采用 HalfEven 舍入的上下文副本,因为规范把这些函数定义为无论上下文模式如何都按 round-half-even 正确舍入;power 则使用上下文的模式,这是规范及其测试向量所要求的,并且即使幂恰好是精确的,也把非整数次幂报告为带 pp 位数字的 Inexact(decimal_power_noninteger_gda_result)。另有两条 GDA 规则被逐字实现:这些函数只对 pp、emax⁡e_{\max} 和 −emin⁡-e_{\min} 都不超过 999,999 的上下文有定义(math_context_is_restricted,否则引发 InvalidContext);在子集算术中,ln 复现经典参考算法的结果,它可能比正确舍入的结果大末位一个单位(decimal_ln_subset_result),因为遗留的子集测试向量编码的正是那个结果。

系数表示与内核

系数要么是用于小于 101810^{18} 的值的 Small(UInt64),要么是以 10910^9 为基数、缓存了位数的 limb 数组,两者都是持久的:运算从不修改操作数的 limb。十进制 limb 使位数计数、去除末尾零、半单位比较以及 ZeroFiveUp 的末位测试都成为常数时间或线性时间,无需二进制到十进制的转换。乘法和除法按 limb 数在教科书算法、Karatsuba、Toom-3、双模数 NTT、Knuth 的 algorithm D、Burnikel–Ziegler 和 Newton 倒数除法之间选择,阈值按目标平台设定:

目标平台Karatsuba mul/squareToom-3首个 NTT mul/squareBurnikel–ZieglerNewton
native96 / 481,1521,728 / 640从 2,816 起关闭
LLVM96 / 962,0484,096 / 2,0482,0484,096
Wasm, Wasm-GC, JS96 / 964,0968,192 / 4,0962,0484,096

阈值是经测量的性能策略,而非语义:每个内核都返回精确的积或商,因此这种选择不可能改变结果或标志。测量数据参见性能。

正确性 / 不变式

  • 单次舍入。 除整数 power 和上文的除法情形外,每个有限结果都是精确结果舍入一次所得,因此在 half 模式下 ∣x^−x∣≤12u|\hat x - x| \le \frac12 u,其他模式下 <u< u,并且在次正规范围之外 ∣x^−x∣/∣x∣≤12101−p|\hat x - x|/|x| \le \frac12 10^{1-p}。
  • 理想同值类成员。 精确结果以最接近其理想指数的可表示指数返回,因此对精确数据的运算会保持其量子(2.50 × 3 = 7.50)。
  • 与状态无关。 F(x⃗,C)F(\vec x, C) 仅通过策略依赖于 CC:结果从不依赖于粘滞状态或陷阱。
  • 状态代数。 状态只会增长,串联执行一系列运算之后的状态就是各引发集合的并。
  • 陷阱确定性。 每个运算至多触发一个陷阱,按固定的全序选出;Trapped 携带的值与 Completed 所给出的相同。
  • 经认证的初等函数。 exp、ln、log10 或非整数 power 的有限不精确结果,只会连同其为正确舍入值的证明(附录中的引理 3 或 4)一起返回;若无法证明,则得到 NaN 和 InvalidOperation。
  • 全序。 compare_total 是三元组上的全序(依次比较符号、类别、值、指数、载荷),Decimal::compare 是全预序,其中所有 NaN 构成一个高于所有数的类。
  • 证据。 固定版本的 official 测试套件通过了 64,986/64,986 个合法的可执行标量用例,遗留的 official0 套件通过了 16,124/16,124 个(符合性)。这些都是有限的声明:上文的除法缺陷以及 API 页面所指出的取整差异,均未被任何固定用例覆盖。

被否决的替代方案

  • 可变或全局的当前上下文。 被否决,因为它使结果和标志依赖于求值顺序和隐藏状态(参见第一条设计决策)。
  • 把陷阱作为错误。 被否决,因为 GDA 定义的结果会丢失。
  • 与 decimal 共享引擎。 被否决,因为一个标准的改动可能改变另一个标准的结果;经测量的重复比意外的语义耦合更划算。
  • 用固定数量保护位的二进制浮点计算初等函数。 被否决,因为没有任何固定数量的保护位能保证正确舍入(制表者困境);带细化的认证则能保证,并在无法做到时显式失败。
  • 规范化每个结果。 被否决,因为量子是 GDA 数的一部分;规范化可通过 reduce 显式进行。

边界

  • 本包只实现 GDA 运算:没有三角函数、双曲函数或规范以外的其他函数,除无状态层的少数附加功能外也没有 IEEE 754 运算。
  • 它不解析 .decTest 文件、不运行测试套件,也不读取文件;这些由 frontend/gda_expr 和仓库工具负责。
  • 它不提供 BID 交换编码;只有 DPD。
  • 对整数次幂,它不提供超出 GDA 要求的正确舍入;在当前分支上,上文所述的 half 模式除法情形也不是正确舍入的。
  • 它不提供可变或全局的上下文,上下文也没有身份:字段相同的两个上下文可以互换。

Footnotes

  1. M. F. Cowlishaw, General Decimal Arithmetic Specification, version 1.70 (2009), https://speleotrove.com/decimal/decarith.html。测试套件是同一作者的 dectest 集合,2.62 版。 ↩

  2. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.2。十进制的摆动比二进制大:在一个十进位段内相对间距的变化因子为 10,而二进制为 2(Goldberg 1991, §1.2)。 ↩ ↩2

  3. A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991。关于球算术,参见 J. van der Hoeven, “Ball arithmetic”, 2009,以及 F. Johansson, “Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”, IEEE Trans. Computers 66(8), 2017。 ↩