bin_float 设计

bin_float 实现任意精度的 IEEE 754 二进制浮点算术。本页解释其背后的数学：值集、舍入函数及其满足的误差模型，各运算如何由精确的整数数据确定正确舍入的结果，如何处理指数范围、微小性和状态标志，为何 IEEE 余数是精确的，十进制转换与初等函数如何认证，以及为何快速整数内核不会改变结果。API 参考 列出了可调用接口，教程 展示了其用法。

设计目标

对 $p$ 从 1 到 $2^{28}$ 位的每种精度以及至多 $\pm(2^{30}-1)$ 的每种指数范围，bin_float 的每个运算都返回 IEEE 754-2019 对正确舍入运算所要求的值——即对精确实数 $f(x)$ 取 $\circ(f(x))$——并恰好给出相应的 IEEE 状态标志。同一套代码服务于两类用户：需要远超 Double 的二进有理值的任意精度数值计算，以及对 binary16、binary32、binary64 和 binary128 的逐位精确模拟（包括次正规数、五种舍入方向和两种微小性规则）。没有隐藏状态：精度、舍入和范围存放在不可变的 BinaryContext 中传入，标志则作为值传回。

数学背景

二进有理值与所存储的三元组

有限的 BinFloat 表示二进有理数

$$x = (-1)^s \cdot c \cdot 2^{e}, \qquad s \in \{0, 1\},\ c \in \mathbb{N},\ e \in \mathbb{Z},$$

其中 $c$ 是 BinCoeff，$e$ 是 exponent2()。该表示是规范的：若 $c \ne 0$ 则 $c$ 为奇数，若 $c = 0$ 则 $e = 0$。每个二进有理数恰有一种这样的形式（提出 $2^{\nu_2(c)}$ 即可），因此两个有限值数值相等，当且仅当它们的符号（对非零值而言）、系数和指数都一致。最高位的指数为

$$\operatorname{top}(x) = \lfloor \log_2 |x| \rfloor = e + \operatorname{bits}(c) - 1,$$

下文的每次比较、范围检查和舍入判定都以 $\operatorname{top}$ 和 $\operatorname{bits}$ 表述，而不使用浮点对数。每个值还携带满足 $\operatorname{bits}(c) \le p$ 的精度 $p$；它记录该值所属的格式，也是对该值进行普通运算时的默认精度。

IEEE 754 二进制格式

具有 $k$ 位的二进制交换格式包含一个符号位、一个 $w$ 位的带偏置指数字段 $E$ 和一个 $(p-1)$ 位的尾随有效数字段 $T$（IEEE 754-2019 第 3.4 条）。¹1 IEEE Std 754-2019，IEEE Standard for Floating-Point Arithmetic：第 3 条（格式）、4.3（舍入方向属性）、5（运算）、6（无穷、NaN、带符号零）、7（默认异常处理）。 设 $e_{\max} = 2^{w-1} - 1$，偏置 $e_{\max}$ 且 $e_{\min} = 1 - e_{\max}$，则一个编码表示

$$v = \begin{cases} (-1)^s\, 2^{E - e_{\max}} \bigl(1 + T\, 2^{1-p}\bigr) & 1 \le E \le 2^w - 2 \quad\text{(normal)},\\ (-1)^s\, 2^{e_{\min}} \bigl(0 + T\, 2^{1-p}\bigr) & E = 0 \quad\text{(subnormal or zero)},\\ (-1)^s\, \infty & E = 2^w - 1,\ T = 0,\\ \mathrm{NaN} & E = 2^w - 1,\ T \ne 0. \end{cases}$$

格式	$k$	$w$	$p$	$e_{\max}$	$e_{\min}$	最大 $\Omega$	最小正规数	最小次正规数
binary16	16	5	11	15	−14	$65504$	$2^{-14}$	$2^{-24}$
binary32	32	8	24	127	−126	$(2-2^{-23})2^{127}$	$2^{-126}$	$2^{-149}$
binary64	64	11	53	1023	−1022	$(2-2^{-52})2^{1023}$	$2^{-1022}$	$2^{-1074}$
binary128	128	15	113	16383	−16382	$(2-2^{-112})2^{16383}$	$2^{-16382}$	$2^{-16494}$

撇开编码不谈，该格式的有限值为

$$F(p, e_{\min}, e_{\max}) = \{\, M \cdot 2^{q} : M \in \mathbb{Z},\ |M| < 2^{p},\ q \ge e_{\min} - p + 1,\ |M| 2^{q} < 2^{e_{\max}+1} \,\}.$$

BinaryContext 恰好就是这个三元组，再加上舍入方向和微小性规则。其 $e_{\min}$ 和 $e_{\max}$ 是最高位指数，因此正规数满足 $e_{\min} \le \operatorname{top}(x) \le e_{\max}$，而 $2^{e_{\min}}$ 以下的网格具有固定的量子

$$\eta = 2^{e_{\min} - p + 1},$$

即最小的正次正规数。缺失的界由实现范围 $\pm(2^{30}-1)$ 代替，该范围足够大，使每一步指数运算都能放入 64 位中间量，每个存储的指数都能放入 Int；binary_precision_max $= 2^{28}$ 也使 $e_{\min} - p + 1$ 保持在范围内。

舍入函数

对于 $x \in \mathbb{R}$，令 $x^- = \max\{y \in F : y \le x\}$ 和 $x^+ = \min\{y \in F : y \ge x\}$，暂且将 $F$ 在 $\pm\Omega$ 之外以 $\pm\infty$ 延拓。BinaryRoundingMode 的六种舍入方向是映射 $\mathbb{R} \to F \cup \{\pm\infty\}$

$$\begin{aligned} \operatorname{RD}(x) &= x^-, \qquad \operatorname{RU}(x) = x^+, \\ \operatorname{RZ}(x) &= \operatorname{sign}(x)\,|x|^-, \qquad \operatorname{RA}(x) = \operatorname{sign}(x)\,|x|^+, \\ \operatorname{RNE}(x) &= \text{the nearer of } x^-, x^+, \text{ the one with even } M \text{ on a tie}, \\ \operatorname{RNA}(x) &= \text{the nearer of } x^-, x^+, \text{ the one of larger magnitude on a tie}, \end{aligned}$$

其中 RA（RoundAwayFromZero）不是 IEEE 属性，而是 @lf_arith.RoundingMode 与十进制核心共用的 GDA “round-up” 模式。上述每个 $\circ$ 的两条性质支撑了本页的大部分证明：

$$\text{(R1)}\ \ x \in F \implies \circ(x) = x, \qquad\qquad \text{(R2)}\ \ x \le y \implies \circ(x) \le \circ(y).$$

(R1) 成立，因为在 $F$ 上 $x^- = x^+ = x$。(R2) 成立，因为每个 $\circ(x)$ 都是 $x$ 的两个相邻值之一，其选择规则在 $x$ 穿过单元 $[x^-, x^+]$ 增大时只会从 $x^-$ 移到 $x^+$。

标准误差模型

设 $u = 2^{-p}$ 为单位舍入误差。取满足 $2^{t} \le |x| < 2^{t+1}$ 和 $t \ge e_{\min}$（正规范围）的 $x$。$F$ 在该二进制数量级区间内的点间距为 $2^{t-p+1}$，因此

$$\begin{aligned} |\operatorname{RN}(x) - x| &\le \tfrac12\, 2^{t-p+1} = 2^{t-p} \le 2^{-p} |x| = u|x|, \\ |\operatorname{RD}(x) - x|,\ |\operatorname{RU}(x) - x| &< 2^{t-p+1} \le 2u|x|, \end{aligned}$$

其中 RN 为 RNE 或 RNA。记 $\circ(x) = x(1 + \delta)$，便得到标准模型

$$\operatorname{fl}(a \circ b) = (a \circ b)(1 + \delta), \qquad |\delta| \le u \ \text{(nearest)}, \quad |\delta| < 2u \ \text{(directed)}.$$

将除以 $|x|$ 改为除以 $|\circ(x)|$，就近舍入的界可加强为 $|\delta| \le u/(1+u)$。²2 N. J. Higham，Accuracy and Stability of Numerical Algorithms，第 2 版，SIAM 2002，§2.2；D. Goldberg，“What every computer scientist should know about floating-point arithmetic”，ACM Computing Surveys 23(1)，1991。 在 $2^{e_{\min}}$ 以下，间距为常数 $\eta$，因此误差是绝对的：$|\operatorname{RN}(x) - x| \le \eta/2$ 且 $|\operatorname{RD}(x) - x| < \eta$。两种区段合在一起，得到带下溢项的模型，

$$\operatorname{fl}(a \circ b) = (a \circ b)(1 + \delta) + \epsilon, \qquad |\delta| \le u,\quad |\epsilon| \le \tfrac{\eta}{2},\quad \delta\epsilon = 0,$$

适用于就近舍入（对定向舍入为 $2u$ 和 $\eta$）。加法和减法从不需要 $\epsilon$：若 $a, b \in F$，则二者都是 $\eta$ 的整数倍，$a \pm b$ 亦然，而 $2^{e_{\min}}$ 以下的 $\eta$ 的倍数位于 $F$ 中。因此次正规的和是精确的，这正是渐进下溢能保持 $a - b = 0 \iff a = b$ 的原因。³3 J.-M. Muller 等，Handbook of Floating-Point Arithmetic，第 2 版，Birkhäuser 2018，§2.1 与 §4.3。

正确舍入是比该模型更强的性质：结果是唯一的那个点 $\circ(f(x))$，而不仅仅是 $u$ 范围内的某个点。bin_float 的全部实现都旨在返回这个点，因此上述模型对每个运算都成立，包括初等函数。

设计决策

由精确数据一次舍入

问题。 结果必须等于精确实数 $r$ 的 $\circ(r)$，但 $r$ 可能需要远多于 $p$ 的位数（两个 $p$ 位数的乘积有 $2p$ 位），甚至无穷多位（商、平方根、$e^x$）。

可选方案。 像不带 FMA 的硬件那样在更宽的格式中计算后再舍入一次；或保留固定数量的保护位；或依据精确信息决定舍入。

选择。 每个运算都计算 $r$ 的精确描述，并调用一个终结器。对于二进有理结果（和、差、积、fma、scaleb、remainder、转换），该描述是精确的整数绝对值 $m$ 和指数 $e$，满足 $r = \pm m 2^{e}$。终结器选取移位量

$$\sigma = \max\bigl(\operatorname{bits}(m) - p,\ (e_{\min} - p + 1) - e,\ 0\bigr),$$

即精度移位与移到次正规网格的移位二者中较大者，并以 $q = \lfloor m 2^{-\sigma} \rfloor$ 和 $0 \le f < 1$ 将 $m 2^{-\sigma} = q + f$ 拆分为三项数据：

$$q, \qquad g = [\,f \ge \tfrac12\,] = \text{bit } \sigma - 1 \text{ of } m, \qquad t = [\,f \notin \{0, \tfrac12\}\,] = [\,\nu_2(m) < \sigma - 1\,].$$

这就是经典的舍入位和粘滞位，由 test_bit 和 ctz 直接从系数中读出，无需构造任何移位后的副本。它们决定了每种舍入方向：记 $q_0$ 为 $q$ 的最低位，

方向	当以下条件成立时递增 $q$
RNE	$g \wedge (t \vee q_0)$
RNA	$g$
RZ	从不
RU	$(g \vee t) \wedge s = 0$
RD	$(g \vee t) \wedge s = 1$
RA	$g \vee t$

推导。 $f > \frac12 \iff g \wedge t$、$f = \frac12 \iff g \wedge \neg t$ 以及 $f > 0 \iff g \vee t$。RNE 在 $f > \frac12$ 时，或在 $f = \frac12$ 且 $q$ 为奇数时将绝对值向上舍入；即 $(g \wedge t) \vee (g \wedge \neg t \wedge q_0) = g \wedge (t \vee q_0)$。定向舍入各行恰在 $f > 0$ 且舍入方向对符号 $s$ 而言是远离零时将绝对值向上舍入。不精确性即 $g \vee t$。

原因。 由于 $\sigma$ 已包含次正规移位，极小的结果直接由 $m$ 一次舍入到网格 $\eta \mathbb{Z}$。若先舍入到 $p$ 位再舍入到次正规网格，就成了双重舍入：略高于粗网格中点的值可能被第一次舍入推到该中点上，再被第二次舍入朝错误方向舍入。从 $q$ 产生的进位（当 $q + 1 = 2^p$ 时）只会使 $\operatorname{top}$ 加一；结果随后重新规范化，并对舍入后的值施加下文的上溢检查。

由商和余数实现除法

问题。 满足 $a = c_a 2^{e_a}$ 的 $a/b$，$b = c_b 2^{e_b}$ 是有理数 $N/D \cdot 2^{e}$（$N = c_a$、$D = c_b$、$e = e_a - e_b$），其二进制展开通常是无限的。

选择。 首先确定精确的最高位指数：设 $k = \operatorname{bits}(N) - \operatorname{bits}(D)$，若 $N \ge D 2^{k}$ 则 $\lfloor \log_2 (N/D) \rfloor$ 为 $k$，否则为 $k - 1$，只需一次整数比较。这确定了目标指数 $\tau = \max(\operatorname{top} - p + 1,\ e_{\min} - p + 1)$，随后一次整数除法

$$N 2^{e - \tau} = q D + r, \qquad 0 \le r < D$$

直接给出 $q$，余数则给出舍入数据：$f = r/D$，因此

$$g = [\,2r \ge D\,], \qquad t = [\,r \ne 0 \wedge 2r \ne D\,].$$

不涉及对 $N/D$ 的任何近似；舍入由 $2r - D$ 的符号决定。当 $e - \tau < 0$ 时，移位被移到分母上；小于一个单位的商通过比较 $2N$ 与 $D 2^{\tau - e}$ 来判定，而无需构造它。

由整数平方根和中点测试实现平方根

问题。 除非 $c 2^{e}$ 是平方数，否则 $\sqrt{c 2^{e}}$ 是无理数。

选择。 最高位指数为 $\lfloor \operatorname{top}(x)/2 \rfloor$（向下取整除法），由此如上确定 $\tau$。将被开方数写作 $X = c\, 2^{e - 2\tau}$，从而 $\sqrt{x} = \sqrt{X}\, 2^{\tau}$。精确的整数平方根给出 $s = \lfloor \sqrt{X} \rfloor$ 及余数 $X - s^2$；当且仅当余数为零时平方根是精确的。否则，舍入数据来自中点 $s + \frac12$：

$$\sqrt{X} \gtrless s + \tfrac12 \iff X \gtrless \bigl(s + \tfrac12\bigr)^2 \iff 4X \gtrless (2s+1)^2,$$

这是一次精确的整数比较（当 $X$ 为分数时，将 2 的幂移到另一侧）。因此 $g = [4X \ge (2s+1)^2]$ 且 $t = \neg\text{exact} \wedge [4X \ne (2s+1)^2]$。当 $X$ 为整数时，右边为奇数而左边为偶数，因此平方根从不恰好是中点；这就是经典事实：$p$ 位数的 $\sqrt{x}$ 从不是 $(p+1)$ 位中点。⁴4 Muller 等，Handbook of Floating-Point Arithmetic，§5.3 与 §7.6。 相等分支仍然保留，因为位数多于上下文精度的操作数可能使 $\sqrt{X}$ 恰好成为中点（例如 $\sqrt{9/4} = 3/2$ 舍入到一位）。

指数范围、上溢与下溢

上溢依据舍入后的值判定：若舍入后的结果满足 $\operatorname{top} > e_{\max}$，则运算上溢，这就是 IEEE 754 的”舍入后”规则（第 7.4 条）。由舍入表可知，这发生在以下阈值处

$$\begin{aligned} \text{RNE, RNA:}\quad & |r| \ge 2^{e_{\max}}\bigl(2 - 2^{-p}\bigr) = \Omega + \tfrac12 \operatorname{ulp}(\Omega), \\ \text{RA, and RU for } r>0, \text{ RD for } r<0:\quad & |r| > \Omega, \\ \text{RZ, and RU for } r<0, \text{ RD for } r>0:\quad & \text{never to } \infty, \end{aligned}$$

因为 $(\Omega, \Omega + \frac12\operatorname{ulp})$ 中的值就近舍入时向下舍入到 $\Omega$，而平局值 $\Omega + \frac12\operatorname{ulp}$ 则舍入到偶数邻值 $2^{e_{\max}+1}$，它位于 $F$ 之外。上溢结果对前两组为 $\pm\infty$，对第三组为 $\pm\Omega$，并总是伴随 overflow 和 inexact。在 binary16 中，$\Omega = 65504$ 且 $\frac12\operatorname{ulp}(\Omega) = 16$：

///|
test "binary16 overflow threshold under nearest rounding" {
  let ctx = @bin_float.BinaryContext::binary16()
  let (below, below_flags) = @bin_float.BinFloat::from_int(65519).round_ctx(ctx)
  let (at, at_flags) = @bin_float.BinFloat::from_int(65520).round_ctx(ctx)
  inspect("\{below} \{below_flags.overflow()}", content="2047p5 false")
  inspect("\{at} \{at_flags.overflow()}", content="inf true")
}

微小性。 非零结果严格位于 $\pm 2^{e_{\min}}$ 之间时称为微小的。IEEE 754-2019（第 7.5 条）允许两种解读，由 TininessDetection 选择其一：

$$\text{before rounding: } |r| < 2^{e_{\min}}; \qquad \text{after rounding: } |\circ_{p,\infty}(r)| < 2^{e_{\min}},$$

其中 $\circ_{p,\infty}$ 在无界指数范围下舍入到 $p$ 位。终结器精确计算 $\operatorname{top}(r)$，并且对于舍入后规则，再仅按精度移位对同一绝对值做第二次拆分。两种解读仅对略低于 $2^{e_{\min}}$、且在 $p$ 位下向上舍入到它的 $r$ 有所不同：对于 binary16，$r = 2^{-14} - 2^{-27}$ 在舍入前是微小的，但在 11 位下舍入为 $2^{-14}$，后者并不微小。

下溢标志。 在默认异常处理下，仅当微小结果同时不精确时才引发下溢标志。终结器对精确结果完全不返回任何标志，因此精确的次正规数（例如如上所示的任意次正规差）不引发任何标志。binary16 乘积 $2^{-14} \times (1 - 2^{-11})$ 展示了完整规则：精确值 $2^{-14} - 2^{-25}$ 在两种解读下都是微小的，它恰好位于间距为 $2^{-24}$ 的两个次正规数的正中间，舍入到偶数邻值 $2^{-14}$（即最小正规数），并引发 underflow 和 inexact，尽管编码后的结果 0x0400 是正规数。

///|
test "underflow is raised for a tiny inexact result that rounds to normal" {
  let format = @bin_float.BinaryInterchangeFormat::Binary16
  let smallest_normal = @bin_float.BinaryInterchange::from_hex("0400", format)
    .unwrap()
    .to_bin_float()
  let below_one = @bin_float.BinaryInterchange::from_hex("3BFF", format)
    .unwrap()
    .to_bin_float()
  let (product, flags) = smallest_normal.mul_ctx(below_one, format.context())
  inspect(product.to_interchange(format).0.to_hex(), content="0400")
  inspect("\{flags.underflow()} \{flags.inexact()}", content="true true")
}

远低于范围。 绝对值必然小于 $\eta/2$ 的结果（由指数界判定，无需构造它，例如 $2^{-10^9} \cdot 2^{-10^9}$）根据舍入方向舍入为 $\pm 0$ 或 $\pm\eta$，并伴随 underflow 和 inexact。必然高于范围的结果则取上溢结果。

带符号零与 NaN

$a, b$ 符号相反时的精确零和 $a + b = 0$ 在除 RD 以外的所有方向下均为 $+0$，在 RD 下为 $-0$；$(-0) + (-0) = -0$（第 6.3 条）。积与商取两符号的异或。NaN 操作数产生第一个 NaN 操作数经静默化后的结果，保留其符号和载荷（第 6.2.3 条允许任一输入 NaN），并且当且仅当某个操作数是信号 NaN 或该运算本身无效（$\infty - \infty$、$0 \cdot \infty$、$0/0$、$\infty/\infty$、$\sqrt{x<0}$、$\operatorname{remainder}(\infty, y)$、$\operatorname{remainder}(x, 0)$）时引发 invalid_operation。标志是值：combine 是按位或，因此一次计算的标志构成一个交换幂等幺半群，可以按任意顺序累积，而全局粘滞寄存器无法为并发代码提供这一点。

加法中相距很远的操作数

问题。 $2^{10^9} + 2^{-10^9}$ 作为二进有理数是精确的，但构造它需要一个二十亿位的系数。

选择。 当两个操作数的首位指数相差超过 $p + 3$ 时，较小的操作数在位置 $\operatorname{exp}(\text{high}) - p - 3$ 处截断，其下的所有部分用一个粘滞位代替：截断后低位操作数的整数部分 $L$ 精确参与运算，若有任何部分被丢弃，则幅值变为以半个单位计的 $2(H \pm L) + 1$（减法时为 $2(H - L - 1) + 1$）。

为何对舍入而言是精确的。 结果满足 $\operatorname{top} \ge \operatorname{top}(\text{high}) - 1$，因此舍入位置至少为 $\operatorname{top}(\text{high}) - p$，舍入位至少比它低一位，而每个被丢弃的位都位于 $\operatorname{top}(\text{high}) - p - 3$ 或更低。因此被丢弃的部分既不改变 $q$ 也不改变 $g$，只影响 $t$ 是否被置位，而替代位恰在丢弃了非零内容时置位 $t$。对于减法，$H - (L + \varepsilon) = (H - L - 1) + (1 - \varepsilon)$ 且 $0 < 1 - \varepsilon < 1$，因此同样的替换也适用于借位形式。包括次正规移位在内的完整论证见下方附件。

融合乘加

fma_ctx 将乘积 $c_x c_y 2^{e_x + e_y}$ 精确地构造为一个精度等于其自身位长的值，并将其交给加法终结器，因此 $x y + z$ 只舍入一次（第 5.4.1 条）。与两次舍入的区别正是该运算的意义所在：对于 binary64 中的 $a = \operatorname{RN}(0.1)$，先 mul_ctx 再 sub_ctx 会丢失 $\operatorname{RN}(a \cdot a) - a \cdot a$（第二次运算看到的是两个相等的数），而 fma_ctx(a, a, -RN(a·a)) 精确地返回它，$-8.33\ldots \cdot 10^{-19}$。结果是精确的，这就是 Dekker 定理：在不发生下溢时，舍入乘积的误差本身属于 $F$。⁵5 T. J. Dekker，“A floating-point technique for extending the available precision”，Numerische Mathematik 18，1971；Muller 等，§4.4。 若乘积指数超出 Int 范围，则乘积要么必然上溢，要么小到相对非零加数而言只相当于一个粘滞位；代码在加数最低位以下 $p + 8$ 个位置处放置一个单独的位，根据上文关于远距操作数的论证，其舍入结果相同。

IEEE 余数是精确的

断言。 若 $x, y \in F$（精度 $p$ 相同，范围相同）且 $y \ne 0$，则满足 $n = \operatorname{RNE}_{\mathbb{Z}}(x/y)$ 的 $r = x - n y$ 属于 $F$。

证明。 记 $x = M_x 2^{q_x}$、$y = M_y 2^{q_y}$，其中 $|M_x|, |M_y| < 2^p$ 且 $q_x, q_y \ge e_{\min} - p + 1$。由 $n$ 的选取，$|r| \le |y|/2$。若 $n = 0$ 则 $r = x$。否则 $|x/y| \ge \frac12$，因此 $|r| \le |y|/2 \le |x|$。此时 $r$ 是 $2^{\min(q_x, q_y)}$ 的整数倍。

• 若 $q_x \ge q_y$：$r = k 2^{q_y}$，其中 $|k| 2^{q_y} \le |M_y| 2^{q_y}/2$，因此 $|k| < 2^{p-1}$。
• 若 $q_x < q_y$：$r = k 2^{q_x}$，其中 $|k| 2^{q_x} \le |x| = |M_x| 2^{q_x}$，因此 $|k| < 2^{p}$。

两种情形下均有 $|k| < 2^p$，指数至少为 $e_{\min} - p + 1$，且 $|r| \le \max(|x|, |y|) \le \Omega$，因此 $r \in F$。$\square$

实现从不构造 $n$，它可能有 $2^{30}$ 位。设 $X = |x| 2^{-m}$、$Y = |y| 2^{-m}$、$m = \min(q_x, q_y)$（均为整数），当 $q_x > q_y$ 时它通过 $2^{q_x - q_y}$ 的模幂计算 $X \bmod 2Y$。记 $X = Q (2Y) + R$、$0 \le R < 2Y$，得到 $\lfloor X/Y \rfloor = 2Q + [R \ge Y]$，因此仅凭 $R$ 即可同时得到 $X \bmod Y$ 和 $\lfloor X/Y \rfloor$ 的奇偶性，而这正是 $n$ 的偶数优先选择所需的全部信息。精确的 $r$ 随后经由通常的终结器处理；由上述断言，对上下文格式的操作数不会发生舍入。

邻值、缩放与整数值

next_up_ctx(x) 将正步长 $2^{\pi - 2}$ 加到 $x$ 上并朝 $+\infty$ 舍入，其中 $\pi = \min(e_{\min} - p,\ e(x),\ \operatorname{top}(x) - p)$。$x$ 旁 $F$ 中相邻点之间的每个间隙都至少为 $2^{\min(e_{\min} - p + 1,\ \operatorname{top}(x) - p)}$，且 $x$ 本身是 $2^{e(x)}$ 的倍数，因此对 $x \in F$ 有 $x < x + 2^{\pi - 2} < x^{+}$，再由 RU 的定义，结果是 $F$ 中高于 $x$ 的最小点。同样的论证对多于 $p$ 位的 $x$ 也成立，这就是此类操作数被接受的原因。这次内部加法的标志被丢弃，因为 nextUp 是静默的（第 5.3.1 条），即使它从 $\Omega$ 步进到 $+\infty$ 亦然。

scaleb_ctx(x, n) 是对 $(c, e + n)$ 应用终结器：在正规范围内精确，在范围之外以下溢和上溢正确舍入。logb_ctx 返回 $\operatorname{top}(x)$，它是精确的，并且对次正规的 $x$ 也正确，因为 $\operatorname{top}$ 是在整数系数上计算的。取整舍入使用相同的舍入位和粘滞位，取 $\sigma = -e$（即二进制小数点以下的位）；to_int_ctx 及其同类先舍入，再将整数与目标范围比较，报告 invalid_operation 而非返回由实现定义的哨兵值。

十进制转换

解析。 from_string_ctx 精确读取 $D \cdot 10^{k}$（$D$ 为不带末尾零的整数）。对于至多 $\max\bigl(400, \lfloor (3n + p)/2 \rfloor + 64\bigr)$ 的 $|k|$（$n$ 为数字位数），它被精确舍入：$k \ge 0$ 时由二进有理终结器处理 $D \cdot 5^{k} \cdot 2^{k}$，$k < 0$ 时由除法终结器处理 $D / 5^{|k|} \cdot 2^{k}$。超出该界后，它在工作精度 $w$ 下使用定向包络 $[\operatorname{RD}_w(D)\operatorname{RD}_w(10^k), \operatorname{RU}_w(D)\operatorname{RU}_w(10^k)]$，并将精度加倍，直到两端舍入到相同的值且标志相同。该界保证循环终止。某一方向的舍入断点是 $F$ 中的点（定向模式）或它们之间的中点（就近模式）；二者都是至多有 $p + 1$ 位有效位的二进有理数。对于 $k \ge 0$，$D 10^{k}$ 的奇数部分是 $5^{k}$ 的倍数，且超出该界后 $5^{k} > 2^{2.32 k} > 2^{p+2}$，因此该值不是断点。对于 $k < 0$，超出该界后 $5^{|k|} > 10^{n} > D$，因此 $5^{|k|} \nmid D$，而 $D/10^{|k|}$ 甚至不是二进有理数。不是断点的值到每个断点都有正的距离，而包络宽度随 $w$ 增大趋于零，因此某个 $w$ 必能完成认证。以 $\log_2 10$ 和安全余量估计、其二进制对数必然超出范围的值，无需任何算术即直接上溢或下溢，因此 1e100000000 不产生任何开销。

定点位数。 to_decimal_string_ctx(x, d) 需要满足 $E = \lfloor \log_{10}|x| \rfloor$ 的 $\operatorname{round}(x / 10^{E - d + 1})$。$E$ 从 $\lfloor \operatorname{top}(x) \log_{10} 2 \rfloor$ 开始，它与答案相差不超过一，再通过比较 $|x|$ 与 $10^{E+1}$ 加以校正，先使用幂的定向界，当它们跨越边界时再精确比较。当操作数至多有 $2^{20}$ 位时，商以整数除法精确构成，因此能识别平局和精确结果；否则就加宽定向包络，直到两端舍入到同一个整数，由于这样的商既不是整数也不是半整数，该过程必然终止。进位产生新的最高位数字（$9.99 \to 10.0$）时，$E$ 加一并重复。

最短输出。 设 $I(x)$ 为在该上下文中以 RNE 舍入到 $x$ 的实数集合；它是一个包含 $x$ 的区间。对每个位数 $n$，$x$ 旁的两个 $n$ 位十进制数（截断的和远离零舍入的）是候选值，若解析某候选值返回 $x$，即它位于 $I(x)$ 中，则接受该候选值。接受性关于 $n$ 是单调的：若 $n$ 位截断 $d_n$ 位于 $I(x)$ 中，则 $(n+1)$ 位截断满足 $d_n \le d_{n+1} \le x$（按绝对值），因此它也位于该区间内，上侧邻值同理。所以最小的可接受 $n$ 可以通过在 $[1, \lceil p \log_{10} 2 \rceil + 2]$（一个候选值总能被接受的上界）上二分查找得到。⁶6 位数 $\lceil p \log_{10} 2 \rceil + 1$ 足以保证往返转换（Matula 1968；Goldberg 1991，定理 15）；多出的一位是为二分查找上端留的余量。 若两个候选值都被接受，则选择较近的那个，再选偶数的那个，这就是该长度下最接近的十进制数。对于 binary64，在所有测试过的值上它都与宿主格式化器的输出一致。

经认证的初等函数

问题。 对于 $f = \exp, \ln, \sin, \ldots$，值 $f(x)$ 是超越数，必须在不知道其值的情况下正确舍入。

可选方案。 具有已证明误差界的固定多项式逼近（速度快，但绑定于单一精度）；Ziv 策略，即以先验误差界求值，在舍入有歧义时以更高精度重试；⁷7 A. Ziv，“Fast evaluation of elementary mathematical functions with correctly rounded last bit”，ACM TOMS 17(3)，1991；区间求值参见 W. Tucker，Validated Numerics，Princeton 2011，以及 F. Johansson，“Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”，IEEE Trans. Computers 66(8)，2017。 或区间求值。

选择。 采用误差界经计算而非估计的 Ziv 循环：每个初等函数在工作精度 $w$ 下求出一个包络 $[L, U] \ni f(x)$，其中每个内部运算对 $L$ 向下舍入，对 $U$ 向上舍入。若

$$\circ(L) = \circ(U) \quad\text{and}\quad \operatorname{flags}(L) = \operatorname{flags}(U),$$

则返回该公共值，否则增大 $w$。由 (R2) 可知这是可靠的： $L \le f(x) \le U$ 蕴含 $\circ(L) \le \circ(f(x)) \le \circ(U)$，因此两端相等必然得到 $\circ(f(x)) = \circ(L)$。标志也是一致的，因为上溢、微小性和不精确性在零的同一侧以同样方式单调；但该测试会显式比较它们，而不依赖这一点。

包络来自具有严格尾项界的级数以及单调的约简。对于 $[0, 1/8]$ 上的 $\exp$，各项 $t_k = x^k/k!$ 满足 $t_{k+1}/t_k = x/(k+1) \le 1/8$，因此在最后一个求和项 $t_n$ 之后

$$\sum_{j > n} t_j \le t_n \sum_{i \ge 1} 8^{-i} = \frac{t_n}{7} \le 2\, t_n,$$

一旦 $t_n < 2^{-(w+8)}$，代码即停止，并将 $2 t_n$（向上舍入）加到上和中；正项的下和本身已是下界。较大的参数先减半 $r$ 次，结果再平方 $r$ 次，平方在正数上是单调的，因此保持包络；$e^{-x} = 1/e^{x}$ 处理负参数。对于满足 $v \in [1, 2]$ 的 $\ln v$，级数为

$$\ln v = 2 \operatorname{artanh} z = 2 \sum_{k \ge 0} \frac{z^{2k+1}}{2k+1}, \qquad z = \frac{v - 1}{v + 1} \in \bigl[0, \tfrac13\bigr],$$

其尾项为 $\sum_{k \ge m} z^{2k+1}/(2k+1) \le \frac{z^{2m+1}}{2m+1} \cdot \frac{1}{1 - z^2}$，即代码所加的界。三角函数在 $w \ge p + \max(0, \operatorname{top}(x) + 1) + 96$ 位下用 $\pi/2$ 的包络（来自 $\pi = 4 \arctan 1$，而它本身由参数减半后的反正切级数包络）约简 $x$，使象限 $k = \operatorname{round}(x / (\pi/2))$ 在包络两端是同一个整数；否则以更高精度重试。这是以暴力提高精度而非存储 $2/\pi$ 表格的方式实现的 Payne–Hanek 思想；其代价随 $\log_2|x|$ 增长，因此需要超过 $10^6$ 位的输入会以 ResourceLimit 被拒绝，而不是运行数分钟。

预算。 循环从 $w_0 = p + 64$ 开始，步进 $w_{i+1} = w_i + \max(32, \lfloor w_i/2 \rfloor)$，至多尝试 12 次。对于 binary64，序列为 $117, 175, 262, \ldots, 10053$ 位。Ziv 关于终止性的论证是 $f(x)$ 不是 $\circ$ 的断点：由 Lindemann–Weierstrass 定理，除平凡例外外，$e^{x}$、$\ln x$、$\sin x$、$\cos x$、$\tan x$ 及其反函数在每个非零代数（特别是二进有理）参数处都取超越值，而断点是二进有理数。这些例外在循环之前被过滤：$e^0 = 1$、$\ln 1 = 0$、$\log_2 2^k = k$、整数 $n$ 的 $2^n$、$\sin(\pm 0)$、sinpi 和 cospi 的整数及半整数参数、$\operatorname{tanpi}(\pm 1/4) = \pm 1$、整数 $n$ 的 $10^n$、$\log_{10} 10^n = n$，等等。对于以 $\pi$ 缩放的函数，Niven 定理表明这些是仅有的二进有理结果。⁸8 I. Niven，Irrational Numbers，1956，推论 3.12：若 $r$ 为有理数且 $\sin(\pi r)$ 为有理数，则 $\sin(\pi r) \in \{0, \pm\frac12, \pm 1\}$；$\cos$ 同理，以及 $\tan(\pi r) \in \{0, \pm 1\}$。值 $\pm\frac12$ 需要分母为 6 或 3 的 $r$，而这不是二进有理数。 非断点到每个断点都有正的距离，因此足够大的 $w$ 必能完成认证。$w$ 必须多大，这就是制表者困境（table maker’s dilemma）：对于任意 $p$，尚无有用的先验界，因此预算是资源限制，而非正确性条件。预算耗尽时，try_* 形式报告一个带有阶段、原因和最后一个 $w$ 的 CertificationFailure，全函数形式则返回静默 NaN 并伴随 invalid_operation；两者都不会返回未经认证的值。固定版本的 MPFR 语料从未耗尽预算。当前分支上有一族例外没有被过滤：指数为 $1/2^k$ 以外的非整数、但结果仍是二进有理数的 pow，例如 $16^{3/4} = 8$。在就近舍入下，包络仍能认证正确的值，但会引发 inexact；在定向舍入下，循环无法完成认证，返回 CertificationFailure。

整数幂。 pow_int_ctx 在 $w = p + \operatorname{bits}(n) + \operatorname{bits}(p) + 4$ 位下以就近舍入通过加法链计算 $x^n$。每个链步 $a_k = a_i + a_j$ 将两个近似值相乘；若 $x^{a} \prod (1 + \delta)^{E(a)}$ 描述了误差结构，则 $E(a_k) = E(a_i) + E(a_j) + 1$ 且 $E(1) = 0$，由归纳得 $E(a) \le a - 1$。因此计算值为 $x^n (1 + \theta)$，用 Higham 的记号即 $|\theta| \le \gamma_{n-1} = (n-1)u_w/(1 - (n-1)u_w)$，它小于 $w$ 位结果的 $2^{\operatorname{bits}(n)}$ 个 ulp。代码使用半径 $2^{\operatorname{bits}(n) + 2}$ 个单位（对负幂再多一个因子 2，其倒数还会再增加 $n$ 个因子），构建区间，并在两端舍入相同时接受，依据同样的 (R2) 论证。若 12 次加倍仍无法认证，则构造精确幂 $c^n 2^{ne}$ 并舍入，因此结果在任何情况下都是正确舍入的。必然超出范围的幂会先依据经认证的 $\log_2$ 界判定，而精确值能放入 $p$ 位的幂会被精确计算，这也保证了 Ziv 路径只会遇到不精确的结果。

系数内核

问题。 精度的代价来自 $p$ 位数的整数乘法和除法，而 $p$ 的范围从 1 到 $2^{28}$。

选择。 BinCoeff 将至多 128 位的值内联存储，更大的值存储为小端序的 32 位 limb（在 JavaScript 上为宿主 bigint），并按较短操作数的 limb 长度 $n$ 进行分派：

乘积	Native	LLVM	Wasm、Wasm-GC
低于此值用竖式乘法	96	96	96
从此值起用 Karatsuba	96	96	96
从此值起用 Toom-3	2048	2048	4096
从此值起用双素数 NTT 乘法	2048	2048	4096
从此值起用 NTT 平方	768	768	3072
从此值起用递归平方	512	768	768

稀疏操作数（非零 limb 很少）使用稀疏乘积，具有 $m > 2n$ 个 limb 的操作数被切分为 $n$ limb 的块。除法使用单 limb 循环，除数少于 48 个 limb 时使用 Knuth 算法 D，从 48 起使用 Burnikel–Ziegler 递归，从 1024 起使用 Newton 倒数；GCD 在超过四个 limb 时从二进制（Stein）算法切换为 Lehmer 批处理。这些阈值由基准测试套件针对每个目标实测得出；它们是策略，而非语义。

为何能保持精确性。 竖式乘法、Karatsuba 和 Toom-3 计算的是整数多项式恒等式，例如

$$(a_1 B + a_0)(b_1 B + b_0) = a_1 b_1 B^2 + \bigl[(a_1 + a_0)(b_1 + b_0) - a_1 b_1 - a_0 b_0\bigr] B + a_0 b_0,$$

而 Toom-3（在 $0, 1, -1, 2, \infty$ 处求值）在插值时对已知为倍数的带符号中间量做精确的除以 2 和除以 3，因此它们都是精确的整数计算。NTT 是唯一的模运算步骤。它将每个操作数拆分为 16 位数字，因此数字卷积的每个系数至多为

$$\min(n_a, n_b)\,(2^{16} - 1)^2 < 2^{23} \cdot 2^{32} = 2^{55}$$

这对至多 $2^{23}$ 的变换长度成立。它对素数 $p_1 = 998244353 = 119 \cdot 2^{23} + 1$ 和 $p_2 = 754974721 = 45 \cdot 2^{24} + 1$ 取模计算卷积，两者都有 $2^{23}$ 次单位根，再用中国剩余定理重组，该重组在 $[0, p_1 p_2)$ 中唯一，其中 $p_1 p_2 \approx 2^{59.4} > 2^{55}$。因此重组得到的系数就是精确的整数。每次变换之前都会先做长度检查；更长的乘积使用长度可接受的重叠块，或回退到 Toom-3。除法路径按构造返回满足 $n = qd + r$ 和 $0 \le r < d$ 的 $(q, r)$；Newton 路径用余数校正其近似商，若需要超过两次校正则中止，因为这意味着存在缺陷，而非数值现象。由于所有路径计算的是同样的整数，算法的选择不会改变任何舍入结果、标志或编码。

compare 中 NaN 的排序

问题。 MoonBit 的 Compare trait 要求一个可供排序和有序映射依赖的三路比较。IEEE 比较是偏序：NaN 与任何值（包括它自身）都是无序的。

可选方案。 (1) 遇到 NaN 时中止，早期版本即是如此；这样一来，对任何可能含有 NaN 的数据排序都会崩溃。(2) 使用 IEEE totalOrder，它是全序，但区分 $-0 < +0$，并把负 NaN 排在 $-\infty$ 之下，因此 compare 在零上会与数值相等不一致。(3) 在数上保持数值顺序，并将所有 NaN 归为位于其上方的同一类。

选择。 方案 (3)。对数定义键 $\kappa(x) = (0, x)$，对 NaN 定义 $\kappa(\mathrm{NaN}) = (1, 0)$，按字典序排列；compare(x, y) 是 $\kappa(x)$ 与 $\kappa(y)$ 的比较，$-0$ 和 $+0$ 映射为同一个数。在全序集中比较键是自反、传递且完全的，因此 compare 是一个全预序；它不是反对称的（$-0$ 与 $+0$，或两个载荷不同的 NaN，比较结果相等但却是不同的值），而 Compare 并不要求反对称性。代价是在 < 下 nan > 1 为真，因此需要 IEEE 语义的代码必须使用 compare_checked（遇 NaN 报错）、compare_quiet / compare_signaling（四值，带标志）或 total_order。结构性的 == 仍是派生的 Eq，因为它是对每个方法（包括精度和载荷）都构成同余关系的唯一相等。

正确性 / 不变式

• 规范形式。 API 产生的每个有限值都满足 $c$ 为奇数或 $c = 0, e = 0$，且 $\operatorname{bits}(c) \le$ 为其精度；存储的指数从不饱和（饱和的指数会先被归类为上溢或下溢）。
• 正确舍入。 对每个算术运算、转换和初等函数以及每个上下文，返回的有限值都等于精确实数结果 $r$ 的 $\circ(r)$，并遵循上述范围规则。由 (R1)，每当 $r \in F$ 时有 $\circ(r) = r$ 且不引发任何标志；round_ctx 是幂等的。
• 标志。 inexact 当且仅当 $\circ(r) \ne r$；overflow 蕴含 inexact；underflow 当且仅当结果微小（按上下文规则）且不精确；division_by_zero 仅用于有限操作数产生精确无穷结果的情形；invalid_operation 当且仅当由非 NaN 操作数产生了静默 NaN，或消耗了信号 NaN。combine 满足结合律、交换律和幂等律。
• 误差模型。 因此，在正规范围内，就近舍入有 $|\circ(r) - r| \le u|r|$，定向舍入有 $< 2u|r|$，在 $2^{e_{\min}}$ 以下另有绝对项 $\eta/2$（相应地 $\eta$），而次正规的和与差是精确的。
• 单调性。 只要实函数在某个参数上单调，相应运算在该参数上也单调，因为它是 $\circ \circ f$ 且 $\circ$ 单调；特别地，RD 和 RU 的结果夹住精确值，ball_float 和 sqrt_bounds_for_precision 正依赖于这一点。
• 精确性定理。 remainder、正规范围内的 scaleb、copy_sign、neg、abs、logb、to_integral_* 以及解码都是精确的；fma(a, b, -RN(ab)) 在无下溢时是精确的。
• 复杂度。 加法与操作数长度呈线性关系（且由远距操作数规则，与指数差无关）；乘法遵循内核表格，从 $O(n^2)$ 到 $O(n \log n)$；在较大的 $n$ 下，除法和平方根的代价为常数次同等规模的乘法；remainder 的代价为 $O(\log(q_x - q_y))$ 次模乘；初等函数在工作精度 $w$ 下对其级数求值，且由于 $w$ 按几何级数增长，所有尝试的总代价不超过最后一次尝试的常数倍。

较长的证明（远距操作数加法规则、nextUp、余数约简、Ziv 接受测试以及 NTT 界）收录于附件中。

bin_float 的舍入与精确性证明

被否决的替代方案

• 对 from_double 以外的格式使用宿主 Double。 让 binary16、binary32 或 binary128 经由 Double 会导致双重舍入，在某些目标上丢失信号 NaN，而且根本无法表示 binary128。因此交换编码改为在 BinCoeff 位模式上完成。
• 固定数量的保护位。 三个保护位足以完成两个 $p$ 位操作数的加法，但不足以应对除法、平方根、从十进制的转换或宽于上下文的操作数。依据精确的整数数据（舍入位、粘滞位、余数符号、中点比较）进行判定，则可以用一个终结器处理所有这些情形。
• 先舍入到 $p$ 位，再舍入到次正规网格。 这种双重舍入会产生错误的次正规结果；终结器只做一次移位，移到两个位置中较粗的那个。
• 使用估计误差界的 Ziv 方法。 它需要对每个函数和每种约简单独做误差分析，而其中的错误会悄无声息地返回错误的最后一位。向外舍入的包络使误差界成为计算得出的量，代价是每个运算都要求值两次。
• 全局标志与舍入状态。 IEEE 754 将标志描述为粘滞的全局状态。返回的 BinaryFlags 值可以与纯代码、并发代码以及 @lf_arith 的 Result 风格组合，而 combine 可在需要时恢复粘滞行为。
• compare 遇到 NaN 时中止，或以 compare 作为 totalOrder。 参见compare 中 NaN 的排序。

边界

bin_float 有意不做以下事情：

• 提供区间或球算术：包络只在认证循环内部使用；ball_float 在 BinFloat 之上构建中点–半径算术；
• 实现 IEEE 754 替代异常处理（陷阱、替换）或粘滞的全局标志：标志是返回值；
• 实现十进制浮点（decimal、decimal_gda）或其上的非二进制 IEEE 运算；
• 为新生成的 NaN 承诺任何特定的 NaN 载荷（它使用载荷 0），或传播多于一个输入的载荷；
• 保证初等函数对每个输入都能在一定时间内完成：认证有预算，预算耗尽会被报告，而不会被隐藏；
• 暴露 limb 布局、阈值或变换参数：只要每个结果、标志和编码保持不变，它们可能随时更改而不另行通知；
• 声称超出符合性中所记录的有限语料之外的符合性。

Footnotes

1. IEEE Std 754-2019，IEEE Standard for Floating-Point Arithmetic：第 3 条（格式）、4.3（舍入方向属性）、5（运算）、6（无穷、NaN、带符号零）、7（默认异常处理）。 ↩

2. N. J. Higham，Accuracy and Stability of Numerical Algorithms，第 2 版，SIAM 2002，§2.2；D. Goldberg，“What every computer scientist should know about floating-point arithmetic”，ACM Computing Surveys 23(1)，1991。 ↩

3. J.-M. Muller 等，Handbook of Floating-Point Arithmetic，第 2 版，Birkhäuser 2018，§2.1 与 §4.3。 ↩

4. Muller 等，Handbook of Floating-Point Arithmetic，§5.3 与 §7.6。 ↩

5. T. J. Dekker，“A floating-point technique for extending the available precision”，Numerische Mathematik 18，1971；Muller 等，§4.4。 ↩

6. 位数 $\lceil p \log_{10} 2 \rceil + 1$ 足以保证往返转换（Matula 1968；Goldberg 1991，定理 15）；多出的一位是为二分查找上端留的余量。 ↩

7. A. Ziv，“Fast evaluation of elementary mathematical functions with correctly rounded last bit”，ACM TOMS 17(3)，1991；区间求值参见 W. Tucker，Validated Numerics，Princeton 2011，以及 F. Johansson，“Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”，IEEE Trans. Computers 66(8)，2017。 ↩

8. I. Niven，Irrational Numbers，1956，推论 3.12：若 $r$ 为有理数且 $\sin(\pi r)$ 为有理数，则 $\sin(\pi r) \in \{0, \pm\frac12, \pm 1\}$；$\cos$ 同理，以及 $\tan(\pi r) \in \{0, \pm 1\}$。值 $\pm\frac12$ 需要分母为 6 或 3 的 $r$，而这不是二进有理数。 ↩