decimal 设计

本页解释 Luna-Flow/floating/decimal 的算术模型：十进制浮点数是什么，同值类（cohort）与首选指数如何携带信息，交换编码如何把数字打包成比特，每个结果如何恰好舍入一次，由此得出哪些误差界，指数范围如何强制执行，以及初等函数如何认证。decimal API规定了每个函数；decimal 教程展示了它们的用法。

设计目标

decimal 对任意精度实现 IEEE 754-2019¹1 IEEE Std 754-2019, Standard for Floating-Point Arithmetic，第 3.3–3.5 条（十进制格式与编码）、第 4 条（属性与舍入）、第 5 条（运算）、第 7 条（异常）和第 9 条（推荐运算）。M. F. Cowlishaw 的 General Decimal Arithmetic 规范（1.70 版）以任意精度的形式给出了相同的模型；本文沿用其术语 coefficient（系数）、adjusted exponent（调整指数）、Etiny 和 clamp。的十进制算术，具有三项性质：

每个上下文运算都是正确舍入的。 结果是精确的数学结果，按所选方向一次舍入到上下文的精度和指数范围——基本运算与初等函数皆是如此。
不会静默丢失任何信息。 结果的指数（其量子（quantum））、零的符号、NaN 载荷以及每种异常条件，都是返回值或返回的 DecimalFlags 的一部分。
没有隐藏状态。 精度、舍入、指数范围和标志都是显式传入、显式返回的普通不可变值，与 Luna-Flow 其他部分一致。

数学背景

十进制浮点数

精度为 $p$ 、调整指数范围为 $[e_{\min}, e_{\max}]$ 的十进制浮点格式是如下数的集合

x = (-1)^s \cdot c \cdot 10^{q}, \qquad s \in \{0,1\},\quad c \in \mathbb{Z},\ 0 \le c < 10^{p},\quad E_{\text{tiny}} \le q \le E_{\text{top}},

再加上 $\pm\infty$ 和 NaN，其中 $c$ 是系数， $q$ 是指数或量子，且

E_{\text{tiny}} = e_{\min} - p + 1, \qquad E_{\text{top}} = e_{\max} - p + 1 .

非零 $x$ 的调整指数为 $\operatorname{adj}(x) = q + \operatorname{digits}(c) - 1 = \lfloor \log_{10} |x| \rfloor$ ：即 $x$ 写成科学记数法 $d_0.d_1d_2\ldots \times 10^{\operatorname{adj}(x)}$ 时的指数。非零 $x$ 在 $\operatorname{adj}(x) \ge e_{\min}$ 时为正规数，否则为次正规数；最小的正次正规数是 $10^{E_{\text{tiny}}}$ ，最大的有限值是

N_{\max} = (10^{p} - 1)\cdot 10^{E_{\text{top}}} = 10^{e_{\max}+1} - 10^{e_{\max}-p+1}.

DecimalContext 恰好存储 $p$ 、 $e_{\min}$ 、 $e_{\max}$ 、舍入模式、clamp（是否允许高于 $E_{\text{top}}$ 的指数；参见钳制）以及微小性（tininess）规则。交换格式如下：

格式	$p$	$e_{\max}$	$e_{\min}$	$E_{\text{tiny}}$	$E_{\text{top}}$	偏置 $=-E_{\text{tiny}}$	指数个数 $E_{\text{top}}-E_{\text{tiny}}+1$
decimal32	7	96	$-95$	$-101$	90	101	$192 = 3\cdot 2^{6}$
decimal64	16	384	$-383$	$-398$	369	398	$768 = 3\cdot 2^{8}$
decimal128	34	6144	$-6143$	$-6176$	6111	6176	$12288 = 3\cdot 2^{12}$

IEEE 754 固定 $e_{\min} = 1 - e_{\max}$ ，因此指数个数为 $E_{\text{top}} - E_{\text{tiny}} + 1 = e_{\max} - e_{\min} + 1 = 2e_{\max}$ ；标准选取 $e_{\max} = 3 \cdot 2^{w-1}$ 使该个数为 $3 \cdot 2^{w}$ ，这恰好是取值为 $0,1,2$ 的两个前导指数位再加 $w$ 个后续位所能编码的数量（参见编码）。把 $q$ 变为非负存储指数的偏置为 $-E_{\text{tiny}} = p - e_{\min} - 1$ ：对于 decimal64， $16 + 383 - 1 = 398$ 。

为什么用十进制

既约有理数 $n/m$ 在基数 $b$ 下有有限展开，当且仅当 $m$ 的每个素因子都整除 $b$ 。对 $b = 2$ ，唯一允许的分母是 2 的幂；对 $b = 10$ ，则为 $2^{i}5^{j}$ 。因此每个二进制浮点数都有有限的十进制展开，但 $0.1 = 1/(2\cdot 5)$ 在二进制中没有：与之最接近的 Double 是

0.1000000000000000055511151231257827021181583404541015625 = \frac{3602879701896397}{2^{55}} .

因此，以十进制定义的量——价格、费率、测量值、协议字段——可以由十进制浮点数精确表示，十进制舍入也发生在人或法规所指定的小数位上。代价是更大的摆动（wobble，见下文）以及更昂贵的数字运算。

同值类与量子

映射 $(s, c, q) \mapsto (-1)^s c\,10^{q}$ 不是单射。同一非零值的所有表示构成它的同值类（cohort）。若 $c$ 有 $d$ 位数字和 $k$ 个末尾零，则在不考虑指数范围时，其成员为 $(c\cdot 10^{j},\, q - j)$ ，其中 $-k \le j \le p - d$ ，因此该同值类有 $p - d + k + 1$ 个成员。例如，decimal32（ $c = 1$ 、 $q = 3$ 、 $d = 1$ 、 $k = 0$ ）中的 $1000$ 有七个成员 $1\text{E+}3, 10\text{E+}2, \ldots, 1000000\text{E-}3$ 。零对每个指数都有一个成员。

同值类成员携带着数值相等所不具备的信息：12.30 表示两位小数，1.2E+3 表示两位有效数字。因此 IEEE 754 为每个运算规定了一个首选指数，精确结果以指数最接近它的成员给出。首选指数源自精确结果自然所处的位置：

\begin{aligned} c_a 10^{q_a} \pm c_b 10^{q_b} &= \bigl(c_a 10^{q_a - m} \pm c_b 10^{q_b - m}\bigr)\,10^{m}, & m &= \min(q_a, q_b),\\ c_a 10^{q_a} \cdot c_b 10^{q_b} &= (c_a c_b)\,10^{q_a + q_b},\\ c_a 10^{q_a} / c_b 10^{q_b} &= (c_a / c_b)\,10^{q_a - q_b},\\ \sqrt{c\,10^{q}} &= \sqrt{c\,10^{q - 2\lfloor q/2\rfloor}}\;10^{\lfloor q/2 \rfloor},\\ x\cdot y + z &: \ \min(q_x + q_y,\ q_z). \end{aligned}

对于和与积，括号中的系数是整数，因此只要精确系数能放入 $p$ 位，就能达到首选指数：1.20 + 3.40 = 4.60 和 1.25 × 2.50 = 3.1250。对于商，只有当 $c_a / c_b$ 有有限十进制展开时精确结果才存在；此时在 $p$ 位允许的范围内尽量向 $q_a - q_b$ 移动，于是 2.400 / 1.2 = 2.00。不精确的结果总是用满 $p$ 位，即指数最小的成员。quantize 把指数作为显式参数，reduce_ctx/normalized 则选择指数最大的成员。

///|
test "design: preferred exponents" {
  let ctx = @decimal.DecimalContext::decimal64()
  let d = fn(s : String) { @decimal.Decimal::from_string(s).unwrap() }
  inspect(d("1.20").add_ctx(d("3.40"), ctx).0, content="4.60")
  inspect(d("1.25").mul_ctx(d("2.50"), ctx).0, content="3.1250")
  inspect(d("2.400").div_ctx(d("1.2"), ctx).0, content="2.00")
  inspect(d("0.0400").sqrt_ctx(ctx).0, content="0.20")
  inspect(d("1.5").fma_ctx(d("2.0"), d("0.25"), ctx).0, content="3.25")
}

舍入方向

设 $x > 0$ 为精确值，目标指数为 $t$ （即保留 $p$ 位数字的指数，对微小结果则为 $E_{\text{tiny}}$ ）。记

x \cdot 10^{-t} = c + f, \qquad c \in \mathbb{Z}_{\ge 0},\ 0 \le f < 1 .

每个舍入方向返回 $c$ 或 $c + 1$ （乘以 $10^{t}$ ）；选择取决于 $f$ 、 $c$ 的末位数字以及符号：

模式	IEEE 名称	返回 $c+1$ 的条件（ $f > 0$ ）
`Down`	roundTowardZero	从不
`Up`	—	总是
`Ceiling`	roundTowardPositive	$x > 0$
`Floor`	roundTowardNegative	$x < 0$
`HalfUp`	roundTiesToAway	$f \ge \tfrac12$
`HalfDown`	—	$f > \tfrac12$
`HalfEven`	roundTiesToEven	$f > \tfrac12$ ，或 $f = \tfrac12$ 且 $c$ 为奇数
`ZeroFiveUp`	—	$c \bmod 5 = 0$

对负的 $x$ ，对 $|x|$ 应用同一规则后恢复符号，因此 Ceiling 与 Floor 互换。每个模式 $\circ$ 都是单调的： $x \le y \Rightarrow \circ(x) \le \circ(y)$ 。初等函数的认证正依赖于单调性。

误差模型

设 $x$ 位于正规范围内， $10^{e} \le |x| < 10^{e+1}$ 。在该十进位段内，可表示数之间的间距为末位单位，即 $\operatorname{ulp}(x) = 10^{e - p + 1}$ 。就近舍入的误差至多为其一半，因此

\begin{aligned} \frac{|\operatorname{fl}(x) - x|}{|x|} \le \frac{\tfrac12\, 10^{e-p+1}}{10^{e}} = \tfrac12\, 10^{1-p} =: u , \end{aligned} \qquad\text{equivalently}\qquad \operatorname{fl}(x) = x(1 + \delta),\ |\delta| \le u .

该界在十进位段的底部附近达到。在顶部附近， $|x| \approx 10^{e+1}$ ，同样的绝对误差仅相当于 $\tfrac12 10^{-p}$ 的相对误差。因此，同一十进位段内最坏与最好相对误差之比为

\frac{\tfrac12\,10^{1-p}}{\tfrac12\,10^{-p}} = 10 = \beta ,

即基数 $\beta = 10$ 的摆动（wobble）。²2 Goldberg, “What every computer scientist should know about floating-point arithmetic”, ACM Computing Surveys 23(1), 1991, §1.2；Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.1–2.2。在二进制中摆动为 2，因此在相同存储下，十进制的最坏情形相对误差略差：decimal64 为 $u = \tfrac12 10^{-15} = 5\cdot 10^{-16}$ ，binary64 为 $u = 2^{-53} \approx 1.1 \cdot 10^{-16}$ 。定向舍入模式为 $|\delta| < 10^{1-p} = 2u$ 。机器 epsilon，即 1 与下一个更大数之间的距离，为 $10^{1-p}$ ——epsilon_contextual 返回该值，next_plus(1) 在 decimal64 中为 1.000000000000001。在 1 以下间距小十倍：next_minus(1) 为 0.9999999999999999。

对次正规结果，间距固定为 $10^{E_{\text{tiny}}}$ ，因此误差界变为绝对误差界 $|\operatorname{fl}(x) - x| \le \tfrac12 10^{E_{\text{tiny}}}$ 。由于每个上下文运算都是正确舍入的，标准模型 $\operatorname{fl}(a \circ b) = (a \circ b)(1+\delta)$ 、 $|\delta| \le u$ 对 $+,-,\times,/,\sqrt{\ }$ 、fma 以及每个结果为正规数的初等函数都成立，因此 Higham 的经典前向与后向误差分析可以在 $u = \tfrac12 10^{1-p}$ 下照搬。

设计决策

一种表示，多种上下文

问题。 应用需要 decimal32/64/128 和任意精度、IEEE 语义和 GDA 兼容性，而又不想在类型之间转换。

方案。 每种格式一个类型（类似硬件）；一个任意精度类型，格式由上下文携带；一个以格式为参数的类型。

选择。 一个 Decimal 类型，持有符号、任意长度的系数、指数、类别、一个 NaN 种类位和一个工作精度字段，再加上一个独立的不可变 DecimalContext。decimal64 计算就是在 DecimalContext::decimal64() 下的计算；交换编码器在编码前应用格式上下文。这样算术集中在一处，调用者可以用 50 位数字工作并在最后舍入到 decimal64，且与 General Decimal Arithmetic 模型一致。值中的精度字段只服务于无上下文的运算符和转换。

只在一处、只舍入一次

问题。 双重舍入——对已舍入的值再次舍入——可能改变正确舍入的结果。

选择。 每个有限的上下文结果都经过同一个最终化例程，它接收一个精确结果 $(s, C, Q)$ （ $C$ 可能远长于 $p$ ），并依次执行：舍入到 $p$ 位、上溢检查、舍入到次正规网格、次正规/下溢标志以及下折（fold-down）。运算内核计算精确整数，从不舍入。例外是精确结果为无限长的运算——除法、平方根、初等函数——下文将分别说明；它们每一个仍然依据精确信息决定最后一位。

舍入位与粘滞信息如何获得

把 $C$ （一个有 $D$ 位数字的非负整数）舍入到 $D - s$ 位，就是除以十的幂，

C = Q \cdot 10^{s} + R, \qquad 0 \le R < 10^{s},

并通过比较 $2R$ 与 $10^{s}$ 在 $Q$ 和 $Q+1$ 之间作出选择。这一次比较恰好携带了经典的舍入位与粘滞位的信息。记 $R = r\,10^{s-1} + R'$ ，其中舍入位为 $r \in \{0,\ldots,9\}$ ，余部为 $0 \le R' < 10^{s-1}$ 。则

2R - 10^{s} = 2(r - 5)\,10^{s-1} + 2R' ,

且由于 $0 \le 2R' < 2\cdot 10^{s-1}$ ：

\begin{aligned} r \ge 6 &\implies 2R - 10^{s} \ge 2\cdot 10^{s-1} > 0,\\ r = 5 &\implies \operatorname{sign}(2R - 10^{s}) = \operatorname{sign}(R'),\\ r \le 4 &\implies 2R - 10^{s} \le -2\cdot 10^{s-1} + 2R' < 0 . \end{aligned}

因此 $2R > 10^{s}$ 、 $2R = 10^{s}$ 和 $2R < 10^{s}$ 分别表示“高于、等于、低于中点”，这正是各 half 模式所需的全部信息； $R \ne 0$ 是定向模式所需的粘滞信息； $Q$ 的末位数字则是 ZeroFiveUp 和 HalfEven 所需的全部信息。当 $s \ge D$ 时整个系数都被丢弃， $Q = 0$ ，只有 $s = D$ 能达到中点：此时代码比较 $C$ 与 $5\cdot 10^{D-1}$ 。把 $10^{p}-1$ 变为 $10^{p}$ 的进位，通过一次精确除以 10 并使指数加一来消除。

只要 $s > 0$ 就引发 rounded，只要 $R \ne 0$ 就引发 inexact；末尾零在舍入前被去掉，使得丢弃零只引发 rounded。

除法

两个有限非零值的商 $a/b$ 依次采取以下三条精确路径之一。

除数是十的幂。 商为 $a$ ，指数作相应平移。
有限商。 用 $g = \gcd(c_a, c_b)$ 把 $c_a / c_b$ 约简为 $c_a'/c_b'$ 。商有有限十进制展开当且仅当 $c_b' = 2^{i}5^{j}$ ；取 $k = \max(i, j)$ ， $\frac{c_a'}{2^{i}5^{j}} = \frac{c_a'\,2^{k-i}\,5^{k-j}}{10^{k}},$ 因此精确结果是整数 $c_a' 2^{k-i}5^{k-j}$ ，指数为 $q_a - q_b - k$ ，与其他精确结果一样交给最终化例程。
无限商。 结果必然不精确。设 $a_c = \operatorname{adj}(c_a/c_b)$ ，它通过比较 $c_a$ 与乘以适当十的幂后的 $c_b$ 精确求得。将分子（或分母）乘以 $10^{p-1-a_c}$ ，使整数商恰有 $p$ 位： $N = Q D + R$ 。是否进位通过比较 $2R$ 与 $D$ 决定——与上文相同的中点测试，只是现在以精确余数作为粘滞信息——因此商只被舍入一次。

第三条路径用于扩展上下文中的正规结果。对次正规结果和子集上下文，代码计算 $p + \operatorname{digits}(c_b) + 2$ 位，按上下文模式舍入，再舍入到精度和 $E_{\text{tiny}}$ 。无上下文运算符 / 使用同样的带保护位方案，取 HalfEven。两次就近舍入并不总是正确的：若第一次舍入恰好落在第二次舍入的中点上，第二次舍入的平局规则就会对一个精确值本已确定的情形重新作出判决。因此该方案除这些近中点情形外都是精确的；例如运算符 / 在五位精度下对 $15/83294$ 会出现这种情况（得到 0.00018008 而非 0.00018009）。div_ctx 在正规结果上没有这一缺陷。

ZeroFiveUp 正是为了让这类两步方案安全而存在的。若 $x$ 先用 ZeroFiveUp 舍入到 $p + k$ 位（ $k \ge 1$ ），再用任意模式 $\circ$ 舍入到 $p$ 位，结果等于 $\circ(x)$ ：不精确的 ZeroFiveUp 结果末位既不是 0 也不是 5，因此它既不是 $p$ 位数，也不是 $p$ 位数的中点，并且相对于每个 $p$ 位数及其中点，它都与 $x$ 位于同一侧。证明见附录。³3 decimal 的舍入证明包含双重舍入引理、上溢表、下折界、认证引理和 NTT 界的完整证明。

平方根

设目标指数为 $t = \max(E_{\text{tiny}},\ \lfloor \operatorname{adj}(x)/2 \rfloor - p + 1)$ ，把操作数缩放为整数 $M = c\,10^{q - 2t}$ （在负分支中若 $q - 2t$ 为奇数，则先乘以 10）。整数平方根 $r = \lfloor \sqrt{M} \rfloor$ 用整数上的牛顿迭代计算，

a_{k+1} = \left\lfloor \frac{a_k + \lfloor M / a_k \rfloor}{2} \right\rfloor , \qquad a_0 = 10^{\lceil (\operatorname{digits}(M)+1)/2 \rceil} > \sqrt{M},

当 $a_k > \lfloor\sqrt M\rfloor$ 时迭代严格递减（由算术–几何平均不等式 $\tfrac12(a + M/a) \ge \sqrt{M}$ ，且 $a_{k+1} < a_k$ 当且仅当 $a_k^2 > M$ ），并在 $\lfloor \sqrt M \rfloor$ 处停止。余数 $M - r^2$ 即粘滞信息。中点测试为 $\sqrt{M} \gtrless r + \tfrac12 \iff 4M \gtrless (2r+1)^{2}$ ，而等号不可能成立，因为 $4M$ 为偶数而 $(2r+1)^2$ 为奇数。平方根永远不会恰好落在中点，因此 HalfEven、HalfUp 和 HalfDown 在其上结果一致。直接在次正规目标指数处舍入，避免了微小根的双重舍入。精确根会先被检测出来（把指数调为偶数后约简的系数是完全平方数），并以首选指数 $\lfloor q/2 \rfloor$ 返回。

指数范围

最终化例程作用于有 $p$ 位数字、指数为 $Q'$ 的已舍入系数 $C'$ ，即在无界指数范围下舍入后的值。

上溢

若 $\operatorname{adj}(C' 10^{Q'}) > e_{\max}$ ，结果上溢。IEEE 754 §7.4 规定，交付的结果是先在具有相同 $p$ 但指数无界的格式中对精确值舍入，然后饱和：从不增大量级的方向不会离开有限范围，可能增大量级的方向则得到无穷。因此，记 $N_{\max} = (10^{p}-1)10^{E_{\text{top}}}$ ：

模式	$x > 0$	$x < 0$
`HalfEven`, `HalfUp`, `HalfDown`, `Up`	$+\infty$	$-\infty$
`Down`	$+N_{\max}$	$-N_{\max}$
`Ceiling`	$+\infty$	$-N_{\max}$
`Floor`	$+N_{\max}$	$-\infty$
`ZeroFiveUp`	$+N_{\max}$	$-N_{\max}$

各 half 模式得到无穷，因为上溢的精确值至少为 $N_{\max} + \tfrac12 10^{E_{\text{top}}}$ （更小的值都会舍入为有限值而不上溢），它位于 $N_{\max}$ 与下一个十的幂之间的中点处或其外。ZeroFiveUp 会饱和，因为 $N_{\max}$ 的末位数字是 9。每次上溢都会引发 overflow、inexact 和 rounded。

///|
test "design: overflow depends on the rounding direction" {
  let ctx = @decimal.DecimalContext::decimal64()
  let big = @decimal.Decimal::from_string("9E+384").unwrap()
  let ten = @decimal.Decimal::from_int(10)
  inspect(big.mul_ctx(ten, ctx).0, content="inf")
  let down = ctx.with_rounding(@def.RoundingMode::TowardZero)
  let (sat, flags) = big.mul_ctx(ten, down)
  inspect(sat, content="9.999999999999999E+384")
  inspect(flags.overflow && flags.inexact, content="true")
}

次正规数、微小性与下溢

若精确结果需要低于 $E_{\text{tiny}}$ 的指数，则舍入到次正规网格：移位量变为 $s = E_{\text{tiny}} - Q$ ，并在保留少于 $p$ 位的情况下应用上述舍入规则。当结果的调整指数低于 $e_{\min}$ 时称其为微小的（tiny），判定可以基于精确值（BeforeRounding），也可以基于在无界指数下舍入到 $p$ 位后的值（AfterRounding，默认）。两种规则只在略低于 $10^{e_{\min}}$ 且向上舍入到它的值上有差异。微小结果引发 subnormal；微小且不精确的结果还会引发 underflow，这是 IEEE 754 §7.5 对默认异常处理的要求。舍入为零的结果获得指数 $E_{\text{tiny}}$ 和 clamped。

钳制

在 clamp 下（所有交换格式均如此），高于 $E_{\text{top}}$ 的指数不可表示，因为编码只容纳 $E_{\text{top}} - E_{\text{tiny}} + 1$ 个指数。 $Q' > E_{\text{top}}$ 且未上溢的结果会被下折（fold down）：系数乘以 $10^{Q' - E_{\text{top}}}$ ，指数设为 $E_{\text{top}}$ ，并引发 clamped。这永远不需要超过 $p$ 位：

\operatorname{digits}(C') + (Q' - E_{\text{top}}) = \bigl(\operatorname{adj} - Q' + 1\bigr) + Q' - (e_{\max} - p + 1) = \operatorname{adj} - e_{\max} + p \le p ,

这里用到了 $\operatorname{adj} \le e_{\max}$ 。值保持不变，只有同值类成员改变。在 decimal32 中，1E+96 存储为 1000000E+90。

///|
test "design: fold-down in decimal32" {
  let ctx = @decimal.DecimalContext::decimal32()
  let (x, flags) = @decimal.Decimal::from_string_ctx("1E+96", ctx)
  inspect(x.coefficient(), content="1000000")
  inspect(x.exponent10(), content="90")
  inspect(flags.clamped, content="true")
}

零没有可填充的数字，因此其指数只是被钳制到 $[E_{\text{tiny}}, E_{\text{top}}]$ （无 clamp 时为 $[E_{\text{tiny}}, e_{\max}]$ ），若有改变则引发 clamped。

标志作为返回值

问题。 在硬件中，IEEE 754 状态标志是粘滞的进程状态；GDA 又增加了陷阱。隐藏状态与 Luna-Flow“语义显式”的原则相冲突，也使并发或组合式代码变得脆弱。

选择。 每个上下文运算返回各自的 DecimalFlags。combine 是逐字段的 OR，因此标志集合构成以 DecimalFlags::new() 为单位元的交换幂等幺半群：在流水线上以任意结合方式累积它们，都得到同一集合，这恰好是去掉状态之后的粘滞标志语义。decimal_checked 封装了这种累积；decimal_gda 则作为独立模型实现粘滞状态与陷阱。五种 IEEE 异常对应 invalid_operation、division_by_zero、overflow、underflow 和 inexact；GDA 条件（rounded、subnormal、clamped、lost_digits、conversion_syntax、division_impossible、division_undefined、invalid_context）对它们作了细化。

quantize 与 same-quantum

x.quantize(y) 返回 $x$ 的值，指数恰为 $t = q_y$ ：

\operatorname{quantize}(c\,10^{q},\, t) = \begin{cases} c\,10^{q-t} \cdot 10^{t} & q \ge t \text{ (exact padding)},\\ \circ\!\left(c\,10^{q - t}\right)\cdot 10^{t} & q < t \text{ (rounding)} . \end{cases}

结果必须在该指数下可表示：新系数至多 $p$ 位， $t$ 必须位于 $[E_{\text{tiny}}, e_{\max}]$ 内，且结果的调整指数不得超过 $e_{\max}$ 。否则该运算无效。它从不替换为别的指数，因为指数本身就是契约（量化到分的金额必须有两位小数）。向上舍入系数可能多出一位（两位精度下的 $9.99 \to 10.0$ ），所以位数检查在舍入之后进行。same_quantum 是谓词 $q_x = q_y$ （对两个无穷或两个 NaN 为真）；在合并指数必须一致的值之前，应使用它进行检验。

交换编码

三种格式共享同一布局：一个符号位、一个 5 位的组合字段 $G$ 、 $w$ 个指数延续位以及 $10J$ 位的尾随有效数，其中 $p = 3J + 1$ ：

格式	$w$	$J$	$1 + 5 + w + 10J$
decimal32	6	2	32
decimal64	8	5	64
decimal128	12	11	128

带偏置指数 $E = q + \text{bias}$ 有 $w + 2$ 位，其最高两位只取 00、01、10；这就是上文推导出的 $3\cdot 2^{w}$ 个数。

DPD。 组合字段存放指数的最高两位和首位数字 $d_0$ ：若 $G = ab\,cde$ 且 $ab \ne 11$ ，则指数位为 $ab$ ，并且 $d_0 = cde \in [0,7]$ ；若 $G = 11\,cd\,e$ 且 $cd \ne 11$ ，则指数位为 $cd$ ，并且 $d_0 = 8 + e$ 。 $G = 11110$ 表示无穷， $G = 11111$ 表示 NaN，下一位区分信号 NaN 与静默 NaN。其余 $3J$ 位数字以密集打包十进制（densely packed decimal）存放在 $J$ 个 declet 中，每 10 位存 3 位数字。⁴4 M. F. Cowlishaw, “Densely packed decimal encoding”, IEE Proceedings — Computers and Digital Techniques 149(3), 2002。IEEE 754-2019 §3.5.2 给出了编码表；代码以布尔公式实现它们，符合性测试语料检查了全部 1024 个 declet。三位数字有 1000 个取值，10 位有 1024 个码字，效率为 $\log_2 1000 / 10 = 99.66\%$ 。称 0–7 的数字为小数字（3 位），8 或 9 为大数字（1 位）。declet $pqr\,stu\,v\,wxy$ 用 $v = 0$ 表示三个小数字，此时它们原样存放在 $pqr$ 、 $stu$ 、 $wxy$ 中； $v = 1$ 标记至少有一个大数字， $wx$ 以及随后的 $st$ 说明是哪几个。按大数字的个数计数，

\underbrace{8^3}_{v=0} = 512,\quad \underbrace{3\cdot 2\cdot 8^2}_{v=1,\ wx\ne 11} = 384,\quad \underbrace{3\cdot 2^2\cdot 8}_{wx=11,\ st\ne 11} = 96,\quad \underbrace{2^3}_{wx=11,\ st=11} = 8,

且 $512 + 384 + 96 + 8 = 1000$ 。前三种情形恰好使用 512、384 和 96 个码字。最后一种情形用 32 个码字表示 8 个值： $r$ 、 $u$ 、 $y$ 携带三位数字的低位， $p$ 、 $q$ 被忽略，因此有 24 个冗余码字：每个由三个大数字构成的值有四种编码，其中 $pq = 00$ 的那一种是规范编码。解码接受全部四种，canonical() 会将其改写。例如，125 是 declet 0010100101 = 0x0A5（三个小数字），999 是 0011111111 = 0x0FF，也可写作 0x1FF、0x2FF、0x3FF。

BID。 系数以二进制整数存储。若符号位之后的两位不是 11，则接下来的 $w+2$ 位是带偏置指数，其余 $10J + 3$ 位是系数。否则指数紧随 11 之后，系数为 $2^{10J+3}$ 加上其余 $10J+1$ 位（隐含“100”）。由于 $10^{7} - 1 < 2^{24}$ 、 $10^{16}-1 < 2^{54}$ 和 $10^{34}-1 < 2^{114}$ ，每个系数都能放下；系数 $\ge 10^{p}$ 的编码是非规范的，解码为零。

///|
test "design: redundant DPD declets decode and canonicalize" {
  let fmt = @decimal.DecimalInterchangeFormat::Decimal64
  let canonical = @decimal.DecimalInterchange::from_hex("#22300000000004FF", fmt).unwrap()
  let redundant = @decimal.DecimalInterchange::from_hex("#22300000000007FF", fmt).unwrap()
  inspect(canonical.to_decimal(), content="19.99")
  inspect(redundant.to_decimal(), content="19.99")
  inspect(redundant.is_canonical(), content="false")
  inspect(redundant.canonical().to_hex(), content="#22300000000004FF")
}

DecimalInterchange 保留原始比特，使非规范输入一直保留到调用者决定规范化为止；算术总是作用于 Decimal。

经认证的初等函数

问题。 对超越函数 $f$ ，对几乎每个十进制 $x$ ， $f(x)$ 都是无理数，因此只能近似；正确舍入的结果需要足以判定舍入的近似（即制表者困境）。

方案。 带先验误差界的固定精度求值（速度快，但正确性取决于该界对每个函数和参数都成立）；或采用带严格误差界的 Ziv 自适应策略⁵5 A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991；J.-M. Muller et al., Handbook of Floating-Point Arithmetic, 2nd ed., Birkhäuser 2018, ch. 12；包络模型参见 J. van der Hoeven, “Ball arithmetic”, 2009。。

选择。 在严格包络上运行 Ziv 循环。对运算 $f$ 和输入 $x$ ：

把 $x$ 精确转换为二进制区间： $\underline{x} = \nabla_w(x)$ 和 $\overline{x} = \Delta_w(x)$ ，to_bin_float，其中 TowardNegative 和 TowardPositive 取 $w$ 位，于是 $x \in [\underline{x}, \overline{x}]$ 。
用 ball_float 在该区间上计算 $f$ ，它返回一个对输入区间中每一点都成立的区间 $[L, U] \ni f(x)$ 。
把 $L$ 和 $U$ （二进分数，因而是精确的十进制数）精确转换为 Decimal，并用目标上下文分别舍入。若两者给出相同的表示（compare_total 相等）和相同的标志，则返回它。
否则增大 $w \leftarrow w + \max(32, \lfloor w/2 \rfloor)$ 并重复，至多 12 次；之后报告认证失败。

接受测试是可靠的，因为舍入是单调的： $L \le f(x) \le U$ 蕴含 $\circ(L) \le \circ(f(x)) \le \circ(U)$ ，若外侧两者是同一表示，中间的也是。标志同样可以传递，前提是 $f(x)$ 本身不可表示：此时有一个端点与公共结果不同，因此两者都不精确，而 overflow、subnormal 和 underflow 在不含零的区间上关于 $|f(x)|$ 单调。因此可表示的结果在循环之前就已处理（见下文）。

初始工作精度为 $w_0 = \max(128,\ 4\max(D, p) + 64)$ 位，其中 $D$ 为输入的位数：一位十进制数字需要 $\log_2 10 \approx 3.32 < 4$ 位，因此 $4\max(D,p)$ 位足以带余量地表示输入和目标，再多 64 位用于覆盖包络的损失。精度安排大约每步增长 $3/2$ 倍；12 步之后预算约为 $w_0 \cdot (3/2)^{12} \approx 130\,w_0$ 位。

有两类输入永远通不过一致性测试，会在循环之前判定：

精确结果。 若 $f(x)$ 是可表示的十进制数， $L < f(x) < U$ 在定向模式下会永远舍入到不同的相邻值。代码会检测精确情形（exp(0)、ln(1)、 $\log_{10} 10^{k}$ 、sinpi/cospi/tanpi 的整数与半整数参数、exp2/exp10 的整数参数、 $x^{1/2}$ 、整数次幂、…）；ball_float 对 $\log_2 8 = 3$ 、 $\sqrt[3]{8} = 2$ 和 $\operatorname{hypot}(3, 4) = 5$ 等精确的二进制结果返回点区间。两者都未检测到的精确结果——例如经由 power_ctx 计算的 $4^{1.5} = 8$ ——不会被识别，会耗尽整个细化预算。
超出范围的结果。 像 $[0, \text{tiny}]$ 这样的包络，其端点在舍入时会产生不同的标志。任何 $\ge 10^{e_{\max}+2}$ 的值都以相同方式上溢，任何非零且 $\le 10^{E_{\text{tiny}}-2} < \tfrac12 10^{E_{\text{tiny}}}$ 的值都以相同方式舍入（舍入为零或最小次正规数，仅取决于模式和符号），因此远处的端点会被替换为这样的代表值；二进制测试使用 $3.322 > \log_2 10$ ，以保证替换是保守的。对于 exp， $x \ge 3(e_{\max}+1)$ 会上溢，因为 $x \log_{10} e \ge 3 \cdot 0.434\,(e_{\max}+1) > e_{\max}+1$ ；由同样的估计， $x \le -3\,|E_{\text{tiny}} - 1|$ 会下溢到 $10^{E_{\text{tiny}}-1}$ 以下。出于同样的目的，power 用定向的 128 位算术界定 $y \log_2 x$ 。

初等函数拒绝 $p$ 或 $|e|$ 超过 999,999 的上下文（invalid_context），以此限制端点精确转换的规模。

系数内核

问题。 在十进制位置上舍入需要快速除以十的幂和快速计数位数；大精度需要次二次复杂度的乘法和除法。

选择。 包私有的 DecCoeff 要么是低于 $10^{9}$ 的内联 UInt，要么是以 $10^{9}$ 为基数的小端 limb 数组，没有前导零 limb，并缓存位数。基数 $10^9$ 是小于 $2^{30}$ 的最大的十的幂，因此 limb 乘积小于 $10^{18} < 2^{60}$ ；按九位的整数倍移位就是移动 limb，digits10 等于 limb 个数加上最高 limb 的位数。BigInt 只出现在公共边界上。

乘法分派：

形状	算法	开销
各一个 limb	内联	$O(1)$
大量零 limb（ $n_a' n_b' \cdot 4 < n_a n_b$ 个非零 limb）	稀疏乘积	$O(n_a' n_b')$
$n_{\text{long}} > 2 n_{\text{short}}$	按较短长度划分的平衡分块	$\frac{n_{\text{long}}}{n_{\text{short}}} M(n_{\text{short}})$
小规模	Comba 按列 / 教科书算法	$O(n^2)$
达到 Karatsuba 阈值起	Karatsuba	$O(n^{\log_2 3}) = O(n^{1.585})$
达到 Toom-3 阈值起	Toom-3	$O(n^{\log_3 5}) = O(n^{1.465})$
达到 NTT 阈值起	双素数 NTT	$O(n \log n)$

Comba 内核只在 $\min(n_a,n_b) \le 18$ 时才把一列 $n$ 个 limb 乘积累加到一个 UInt64 中：一列之和加上传入的进位至多为 $18(10^{9}-1)^{2} + 2\cdot 10^{10} < 1.8\cdot 10^{19} < 2^{64}$ ，而 19 个乘积可能超过 $2^{64} \approx 1.845\cdot 10^{19}$ 。

NTT 把系数拆分为以 $10^4$ 为基数的数字，并在模素数 $p_1 = 998\,244\,353 = 119\cdot 2^{23} + 1$ 和 $p_2 = 754\,974\,721 = 45 \cdot 2^{24} + 1$ 下进行卷积，这两个素数都有 $2^{23}$ 次单位根，因此长度至多为 $2^{23}$ 的变换都存在。卷积系数是至多 $m = \min(n_a, n_b)$ 个小于 $10^{4}$ 的数字乘积之和，因而小于 $m \cdot 9999^{2}$ ；它可以由中国剩余定理从其余数精确恢复，

z = r_1 + p_1 \bigl((r_2 - r_1)\, p_1^{-1} \bmod p_2\bigr),

只要 $m \cdot 9999^{2} < p_1 p_2 \approx 7.54 \cdot 10^{17}$ ，即 $m < 7.5 \cdot 10^{9}$ ——在变换长度上限以内总是成立。若用基数 $10^{9}$ 的数字，则需要 $m\cdot 10^{18} < p_1p_2$ ，这对任何 $m \ge 1$ 都不可能，所以 NTT 使用更小的数字。当这些界不成立时，内核回退到 Toom-3。

除法使用单 limb 除法、Knuth 的 Algorithm D、⁶6 D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., §4.3.1 (Algorithm D) and §4.3.3；R. Brent and P. Zimmermann, Modern Computer Arithmetic, Cambridge 2010, §1.3–1.4 and §2.4 (NTT)；C. Burnikel and J. Ziegler, “Fast recursive division”, MPI-I-98-1-022, 1998。 Burnikel–Ziegler 递归除法或 Newton 倒数除法。Newton 方法通过 $r_{k+1} = \lfloor r_k (2S - d r_k)/S \rfloor$ 为 $S = B^{n+1}$ 计算 $r \approx S/d$ ：由 $r_k = S/d - \varepsilon_k$ 可得 $S/d - r_{k+1} \approx d\,\varepsilon_k^{2}/S$ ，因此每一步正确的 limb 数翻倍。迭代必须从下方单调递增；若不是这样，或最终的商修正需要超过 $2n+8$ 步，该例程就回退到 Burnikel–Ziegler，而后者在不支持的形状上又回退到 Algorithm D。

各交叉点按目标平台测量，并存放在目标平台专属的文件中：

目标平台	Karatsuba mul/square	Toom-3	首个 NTT mul/square	Burnikel–Ziegler	Newton
native	96 / 48	1,152	1,728 / 640	从 2,816 起	已禁用
LLVM	96 / 96	2,048	4,096 / 2,048	2,048	4,096
Wasm / Wasm-GC / JS	96 / 96	4,096	8,192 / 4,096	2,048	4,096

（limb 为九位数字。）在 native 上，NTT 阈值还取决于变换长度——乘法依次为 1,728、2,816、4,608、7,680，然后是 8,192 个 limb；平方依次为 640、1,040、1,824、3,648、7,296，然后是 8,192——而对于更大的分块长度，Burnikel–Ziegler 的起点移到 5,120 和 10,240 个 limb。这些是分派边界：它们只改变开销，从不改变结果。native 的 Newton 路径已实现并经过测试，但处于禁用状态，因为 native 上的测量没有显示出交叉点。

正确性 / 不变式

表示。 有限的 Decimal 满足 $c \ge 0$ ；DecCoeff 的 limb 是规范的（无前导零 limb，位数精确）。零的符号保存在符号位中；coefficient() 从不携带符号。
单次舍入。 +、-、×、fma、sqrt、quantize、各种转换以及（对正规结果的）/ 的每个上下文结果，都是精确结果舍入一次所得；初等函数的结果只要返回就是正确舍入的，失败时会报告失败，而绝不给出近似值。
误差界。 对这些运算的正规结果， $\operatorname{fl}(x) = x(1+\delta)$ ，其中在 half 模式下 $|\delta| \le \tfrac12 10^{1-p}$ ，在定向模式下 $|\delta| < 10^{1-p}$ 。
精确性可见。 当且仅当返回值与精确结果不同时引发 inexact；只要有数字被丢弃就引发 rounded。
同值类保持。 能放下的精确结果以首选指数返回；quantize 要么返回指数 $q_y$ ，要么失败。
标志。 combine 满足结合律、交换律和幂等律，单位元为 new()。
序。 compare 是全预序（NaN 彼此相等且大于所有数， $-0 = +0$ ）；compare_total 是表示上的全序，并在非 NaN 值上细化 compare。
编码。 对规范比特先解码再编码是恒等映射；对适合该格式的值先编码再解码也是恒等映射，包括同值类、零的符号和 NaN 载荷（DPD）。
复杂度。 比较、加法、移位和单 limb 除法相对 limb 数为 $O(n)$ ；乘法和除法遵循分派表。

中点测试、ZeroFiveUp 双重舍入引理、上溢表、下折界、认证引理以及各内核界的证明汇集在附录中：

decimal 的舍入证明

符合性页面记录了有限的证据（固定的 IEEE 测试语料、穷举的 declet 检查、经 MPFR 认证的初等函数用例、四个目标平台）。

被否决的替代方案

二进制的 BigInt 系数。 在十进制位置上舍入，每次运算都需要除以 $10^{s}$ 并计数位数；用二进制系数时两者都很昂贵，而用基数 $10^{9}$ 的 limb 时，它们只是 limb 移动加一次小除法。
环境上下文与粘滞标志。 出于显式性和可组合性的考虑而被否决；标志幺半群提供了同样的信息。
遇到 NaN 就中止的 compare。 排序和泛型 Compare 代码会在含 NaN 的数据上中止。当前的 compare 是全预序，IEEE 语义可通过 compare_checked、compare_ctx、compare_signal_ctx 和 compare_total 获得。
带解析误差界的固定精度超越函数内核。 在小精度下更快，但需要针对每个函数、每种精度给出正确性证明；包络循环在构造上即正确，且其失败方式是显式的。
IEEE 与 GDA 共用一个包。 粘滞状态、陷阱及陷阱优先级会改变每个运算的类型；它们放在 decimal_gda 中。
规范化每个结果。 丢失同值类会使 12.30 与 12.3 无法区分，并破坏对量子敏感的协议。

边界

decimal 有意不做以下事情：

保持粘滞状态或陷阱（请使用 decimal_checked 或 decimal_gda）；
把无上下文运算符舍入到某个上下文：* 是精确的，+ 和 / 只舍入到操作数精度，且都不施加指数范围；
保证初等函数的认证一定成功：在细化预算耗尽后（预算的最后几步要处理非常宽的数，可能耗时很长），它会报告 CertificationFailure（try_*_ctx），或返回带 invalid_operation 的 NaN（*_ctx）；
在精度或指数超过 999,999 的上下文中计算初等函数；
在二进制转换中保留 NaN 载荷，或保留值精度与格式精度不同的 BID NaN 载荷；
暴露其系数表示、内核选择或阈值；
作出超出符合性页面上有限证据的符合性声明。

IEEE Std 754-2019, Standard for Floating-Point Arithmetic，第 3.3–3.5 条（十进制格式与编码）、第 4 条（属性与舍入）、第 5 条（运算）、第 7 条（异常）和第 9 条（推荐运算）。M. F. Cowlishaw 的 General Decimal Arithmetic 规范（1.70 版）以任意精度的形式给出了相同的模型；本文沿用其术语 coefficient（系数）、adjusted exponent（调整指数）、Etiny 和 clamp。 ↩
Goldberg, “What every computer scientist should know about floating-point arithmetic”, ACM Computing Surveys 23(1), 1991, §1.2；Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.1–2.2。 ↩
decimal 的舍入证明包含双重舍入引理、上溢表、下折界、认证引理和 NTT 界的完整证明。 ↩
M. F. Cowlishaw, “Densely packed decimal encoding”, IEE Proceedings — Computers and Digital Techniques 149(3), 2002。IEEE 754-2019 §3.5.2 给出了编码表；代码以布尔公式实现它们，符合性测试语料检查了全部 1024 个 declet。 ↩
A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991；J.-M. Muller et al., Handbook of Floating-Point Arithmetic, 2nd ed., Birkhäuser 2018, ch. 12；包络模型参见 J. van der Hoeven, “Ball arithmetic”, 2009。 ↩
D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., §4.3.1 (Algorithm D) and §4.3.3；R. Brent and P. Zimmermann, Modern Computer Arithmetic, Cambridge 2010, §1.3–1.4 and §2.4 (NTT)；C. Burnikel and J. Ziegler, “Fast recursive division”, MPI-I-98-1-022, 1998。 ↩