ball_float 设计

ball_float 以实数集合而非单个近似值进行计算。本页陈述其数学契约，推导代码所用的公式，解释其背后的取舍，并列出本包不做的事情。API 参考 说明了每个函数；教程 展示了它们的用法。

设计目标

浮点计算返回一个接近真实结果的数，其误差需要调用者另行估计。ball_float 则在任意工作精度下返回一个可证明包含真实结果的区间，从而使误差界成为值本身的一部分。由此得出三项要求：

1. 包含性先于紧致性。 每个运算都返回其操作数精确像的超集。结果偏宽是可以接受的；结果遗漏某个可能的值则是缺陷。
2. 显式语义。 精度、目标格式、标志和装饰都是作为值传入和返回的，从不是进程全局状态（不切换硬件舍入模式）。
3. 标准词汇。 集合、关系和装饰遵循 IEEE 1788-2015¹1 IEEE Std 1788-2015，IEEE Standard for Interval Arithmetic。本包遵循其基于集合的风格（set-based flavor）。 ，因此结果可以用其测试语料进行检验（参见符合性）。

数学背景

区间与包含性质

区间是满足 $\underline{x} \in \mathbb{R} \cup \{-\infty\}$、$\overline{x} \in \mathbb{R} \cup \{+\infty\}$ 和 $\underline{x} \le \overline{x}$ 的集合 $\boldsymbol{x} = [\underline{x}, \overline{x}] = \{\xi \in \mathbb{R} : \underline{x} \le \xi \le \overline{x}\}$，或空集 $\emptyset$。按照 IEEE 1788 基于集合的模型，无穷端点仅表示集合无界；元素总是实数。整条实数线 $(-\infty, +\infty)$ 称为 Entire。

对于定义域为 $D_f \subseteq \mathbb{R}^n$ 的实函数 $f$，其在一个盒子上的值域为

$$f(\boldsymbol{x}) = \{ f(\xi) : \xi \in \boldsymbol{x} \cap D_f \}, \qquad \boldsymbol{x} = \boldsymbol{x}_1 \times \cdots \times \boldsymbol{x}_n ,$$

而它的凸包 $\operatorname{hull} f(\boldsymbol{x})$ 是包含该值域的最小区间。$f$ 的区间扩展是对每个盒子都满足 $f(\boldsymbol{x}) \subseteq F(\boldsymbol{x})$ 的映射 $F$。$D_f$ 之外的点会被忽略而不报告为错误：$\sqrt{[-1, 4]} = [0, 2]$ 且 $\sqrt{[-2, -1]} = \emptyset$。是否发生了定义域违例由装饰另行报告。

BallFloat 的每个公开运算都是区间扩展。复合运算于是给出支撑本包的结论：

区间算术基本定理。²2 R. E. Moore，Interval Analysis，Prentice-Hall，1966；以及 Moore、Kearfott、Cloud，Introduction to Interval Analysis，SIAM，2009，定理 5.1。 若对表达式中的每个运算都使用区间扩展求值，则结果包含该表达式在输入盒子中所有有定义之点处的值。

证明。 对表达式归纳。变量 $\xi_i$ 求值为 $\boldsymbol{x}_i \ni \xi_i$。对于 $e = g(e_1, \dots, e_k)$，设 $\xi$ 是 $e$ 有定义的点；则每个 $e_j$ 都在 $\xi$ 处有定义，且由归纳假设 $y_j = e_j(\xi) \in E_j$，即 $e_j$ 的区间值。由于 $g$ 在 $(y_1, \dots, y_k)$ 处有定义且 $G$ 扩展了 $g$，

$$\begin{aligned} e(\xi) = g(y_1, \dots, y_k) &\in g(E_1 \times \cdots \times E_k) && y_j \in E_j \\ &\subseteq G(E_1, \dots, E_k) = E(\boldsymbol{x}) && G \text{ is an interval extension.} \end{aligned}$$

$\square$ 完整论证（含下文各引理）见附件：

ball_float 的包含性证明

向外舍入

端点是具有 $p$ 位有效位的 BinFloat 数，因此精确的端点值必须经过舍入。记 $F_p$ 为 $p$ 位二进制数的集合，并且

$$\operatorname{RD}_p(a) = \max\{ f \in F_p : f \le a \}, \qquad \operatorname{RU}_p(a) = \min\{ f \in F_p : f \ge a \}.$$

由这些定义可直接得到三条性质：

$$\begin{aligned} &\text{(i)}\quad \operatorname{RD}_p(a) \le a \le \operatorname{RU}_p(a), \\ &\text{(ii)}\quad a \le b \implies \operatorname{RD}_p(a) \le \operatorname{RD}_p(b) \text{ and } \operatorname{RU}_p(a) \le \operatorname{RU}_p(b), \\ &\text{(iii)}\quad \min_{s \in S} \operatorname{RD}_p(s) = \operatorname{RD}_p(\min S), \quad \max_{s \in S} \operatorname{RU}_p(s) = \operatorname{RU}_p(\max S). \end{aligned}$$

(ii) 成立是因为对 $\operatorname{RD}_p(b)$ 取最大值的集合包含对 $\operatorname{RD}_p(a)$ 取最大值的集合；(iii) 由 (ii) 推出，因为在 $S$ 上的最小值在 $\min S$ 处取得。因此，若对集合 $A$ 有 $L \le \inf A$ 且 $U \ge \sup A$，则

$$A \subseteq [L, U] \subseteq [\operatorname{RD}_p(L), \operatorname{RU}_p(U)].$$

这是舍入进入本包的唯一途径：每个运算计算精确值域的下界 $L$ 和上界 $U$，并存储 $\operatorname{RD}_p(L)$ 和 $\operatorname{RU}_p(U)$。性质 (iii) 意味着候选值在取最小值之前还是之后舍入并无区别；代码只对极值候选舍入一次（quantize_interval）。

舍入增加的宽度在每个端点上至多为一个 ulp。对于精确值域为 $[S_\ell, S_u]$ 的和：

$$\begin{aligned} \operatorname{RU}_p(S_u) - \operatorname{RD}_p(S_\ell) &\le (S_u + \operatorname{ulp}_p(S_u)) - (S_\ell - \operatorname{ulp}_p(S_\ell)) \\ &\le w(\boldsymbol{x}) + w(\boldsymbol{y}) + 2^{1-p}\bigl(|S_u| + |S_\ell|\bigr), \end{aligned}$$

这里用到了 $S_u - S_\ell = w(\boldsymbol{x}) + w(\boldsymbol{y})$ 和 $\operatorname{ulp}_p(a) \le 2^{1-p}|a|$。因此宽度在计算过程中以加性方式增长，每次舍入再加上相对量 $2^{1-p}$。

端点公式

BallFloat 存储两个端点（参见存储端点，而非中点与半径），因此每个运算都是关于端点的公式。

和与差。 $\xi + \eta$ 对两个参数均递增，$\xi - \eta$ 对 $\xi$ 递增、对 $\eta$ 递减，因此盒子上的极值出现在使每个参数朝正确方向取极值的角点处：

$$\begin{aligned} \boldsymbol{x} + \boldsymbol{y} &= [\underline{x} + \underline{y},\ \overline{x} + \overline{y}], & \boldsymbol{x} - \boldsymbol{y} &= [\underline{x} - \overline{y},\ \overline{x} - \underline{y}]. \end{aligned}$$

积。 对固定的 $\eta$，$\xi \mapsto \xi\eta$ 是仿射的，因此其在 $\boldsymbol{x}$ 上的极值位于 $\underline{x}$ 或 $\overline{x}$；对 $\eta$ 运用同样的论证可得

$$\boldsymbol{x}\boldsymbol{y} = [\min S, \max S], \qquad S = \{\underline{x}\,\underline{y},\ \underline{x}\,\overline{y},\ \overline{x}\,\underline{y},\ \overline{x}\,\overline{y}\}.$$

符号模式决定了 $S$ 中哪些元素可能取到极值。若 $\underline{x}, \underline{y} \ge 0$，乘积对两个参数均递增，因此值域为 $[\underline{x}\,\underline{y}, \overline{x}\,\overline{y}]$；其余单一符号的情形由 $\xi\eta = (-\xi)(-\eta) = -((-\xi)\eta)$ 推出，由此得到 multiplication_bounds 所用的表格：

$\boldsymbol{x}$	$\boldsymbol{y}$	下界	上界
$\ge 0$	$\ge 0$	$\underline{x}\,\underline{y}$	$\overline{x}\,\overline{y}$
$\le 0$	$\le 0$	$\overline{x}\,\overline{y}$	$\underline{x}\,\underline{y}$
$\ge 0$	$\le 0$	$\overline{x}\,\underline{y}$	$\underline{x}\,\overline{y}$
$\le 0$	$\ge 0$	$\underline{x}\,\overline{y}$	$\overline{x}\,\underline{y}$
$\ni 0$	$\ni 0$	$\min(\underline{x}\,\overline{y}, \overline{x}\,\underline{y})$	$\max(\underline{x}\,\underline{y}, \overline{x}\,\overline{y})$

在最后一行中两个区间都跨越 0，因此 $\underline{x}\,\overline{y}$ 和 $\overline{x}\,\underline{y}$ 为 $\le 0$，另外两个为 $\ge 0$；最小值在前者之中，最大值在后者之中。对于无界区间，四个乘积都按 $0 \cdot \infty := 0$ 构成：无论另一个因子为何，靠近零端点的实数点给出的乘积都接近 0。

除法。 IEEE 1788 定义了 $\boldsymbol{x}/\boldsymbol{y} = \operatorname{hull}\{\xi/\eta : \xi \in \boldsymbol{x}, \eta \in \boldsymbol{y}, \eta \ne 0\}$。当 $0 \notin \boldsymbol{y}$ 时，$1/\eta$ 在 $\boldsymbol{y}$ 上连续且递减，因此 $1/\boldsymbol{y} = [1/\overline{y}, 1/\underline{y}]$ 且 $\boldsymbol{x}/\boldsymbol{y} = \boldsymbol{x} \cdot (1/\boldsymbol{y})$；代码直接用定向除法计算所选的端点商，而不是舍入两次。当 $0 \in \boldsymbol{y}$ 时：

• $\boldsymbol{y} = \{0\}$：不存在可取的 $\eta$，结果为 $\emptyset$。
• $\underline{y} < 0 < \overline{y}$，$\boldsymbol{x} \ne \{0\}$：$\eta$ 从两侧趋于 0，因此 $\xi/\eta$ 在两个方向上都无界：Entire。
• $\boldsymbol{y} = [0, \overline{y}]$ 且 $\overline{y} > 0$、$\underline{x} > 0$：当 $\eta \downarrow 0$ 时 $\xi/\eta \to +\infty$，而最小的商为 $\underline{x}/\overline{y}$，因此结果为 $[\underline{x}/\overline{y}, +\infty)$；其他符号情形与之镜像对称。若 $\boldsymbol{x}$ 跨越 0，则两侧都无界：Entire。
• $\boldsymbol{x} = \{0\}$ 且 $\boldsymbol{y} \ne \{0\}$：每个可取的商都为 0。

从中点–半径到端点

用户常以 $c \pm r$ 的形式了解一个值。BallFloat::new(c, r, precision=p) 将中心就近舍入为 $\tilde c = \operatorname{RN}_p(c)$，并存储

$$[\tilde c - R,\ \tilde c + R], \qquad R = \operatorname{RU}_p\bigl(\operatorname{RU}_p(r) + \operatorname{RU}_p(|c - \tilde c|)\bigr).$$

断言： $[c - r, c + r] \subseteq [\tilde c - R, \tilde c + R]$。对于 $|\xi - c| \le r$，

$$\begin{aligned} |\xi - \tilde c| &\le |\xi - c| + |c - \tilde c| && \text{triangle inequality} \\ &\le \operatorname{RU}_p(r) + \operatorname{RU}_p(|c - \tilde c|) && \text{property (i)} \\ &\le R && \text{property (i) again.} \end{aligned}$$

随后精确地构成端点 $\tilde c \pm R$（它们可能多于 $p$ 位）。with_precision 以当前的中心和半径使用相同的构造，中心的舍入模式由调用者选择。反向视图是精确的：center 返回 $(\underline{x} + \overline{x})/2$，radius 返回 $(\overline{x} - \underline{x})/2$，它们都是二进有理数，无需舍入（半径仅在低于指数范围下溢时才向上舍入），因此 $[\text{center} - \text{radius}, \text{center} + \text{radius}]$ 恰好就是所存储的区间。Show 打印的正是这一对值。

依赖问题

基本定理把变量的每次出现都视为独立的点。对于 $\boldsymbol{x} = [1, 2]$，

$$\boldsymbol{x} - \boldsymbol{x} = [1 - 2,\ 2 - 1] = [-1, 1] \supsetneq \{0\} = \{\xi - \xi : \xi \in \boldsymbol{x}\}.$$

结果是正确的（它包含 0），但并不紧致：区间减法是 $(\xi, \eta) \mapsto \xi - \eta$ 的扩展，而盒子 $\boldsymbol{x} \times \boldsymbol{x}$ 同时包含 $(1, 2)$ 和 $(2, 1)$。一般地 $f(\boldsymbol{x}) \subseteq F(\boldsymbol{x})$，且仅当每个变量在表达式中只出现一次时等号成立（Moore 的另一条定理）。同样的效应解释了对 $0 \in \operatorname{int}\boldsymbol{x}$ 有 $\boldsymbol{x}\boldsymbol{x} \supsetneq \boldsymbol{x}^2$（$[-1, 2]\cdot[-1, 2] = [-2, 4]$ 但 $[-1, 2]^2 = [0, 4]$），以及次分配律 $\boldsymbol{x}(\boldsymbol{y} + \boldsymbol{z}) \subseteq \boldsymbol{x}\boldsymbol{y} + \boldsymbol{x}\boldsymbol{z}$。因此本包提供单次出现的运算——square、pown、fma、hypot、cancel_minus 以及各初等函数——其结果是真实值域的凸包（至多相差舍入），而不是独立因子的乘积。区间值在加法下也不构成群：cancel_minus 才是撤销加法的运算，因为一般而言 $(\boldsymbol{x} + \boldsymbol{y}) - \boldsymbol{y} \ne \boldsymbol{x}$。

初等函数：单调性、临界点与极点

对于连续函数，其在区间上的值域是一个区间，值域的端点在 $\boldsymbol{x}$ 的端点或其内部的临界点处取得。本包使用三种模式。

• 单调函数（exp、exp2、exp10、expm1、ln、log2、log10、log1p、sqrt、sinh、tanh、asinh、acosh、atanh、asin、atan，以及递减的 acos）：在 $\boldsymbol{x} \cap D_f$ 上的值域为 $[f(\underline{x}), f(\overline{x})]$（对于递减的 $f$ 则交换），因此只需对端点求值，下端点向下舍入，上端点向上舍入。区间内部的定义域边界以该处的极限代替（$\ln \xi \to -\infty$ 当 $\xi \downarrow 0$）。

• 已知极值的函数。 cosh 在 0 处取最小值 1。偶数指数的 pown 在 0 处取最小值 0。sin 在 $\xi = k\pi/2$ 处取 $+1$（其中 $k \equiv 1 \pmod 4$），在 $k \equiv 3$ 处取 $-1$；cos 在 $k \equiv 0$ 处取 $+1$，在 $k \equiv 2$ 处取 $-1$。在两个相邻临界点之间这两个函数都是单调的，因此对 $K = \{k \in \mathbb{Z} : k\pi/2 \in \boldsymbol{x}\}$ 有

  $$\sin(\boldsymbol{x}) \subseteq [m, M], \qquad m = \begin{cases} -1 & \exists k \in K,\ k \equiv 3 \\ \min(\sin\underline{x}, \sin\overline{x}) & \text{otherwise,} \end{cases} \qquad M = \begin{cases} 1 & \exists k \in K,\ k \equiv 1 \\ \max(\sin\underline{x}, \sin\overline{x}) & \text{otherwise.} \end{cases}$$

  代码用 $[\lceil q^-(\underline{x}) \rceil, \lfloor q^+(\overline{x}) \rfloor]$ 包络 $K$，其中 $q^- \le 2\xi/\pi \le q^+$ 由 $\pi$ 的认证包络 $[\pi^-, \pi^+]$ 计算得到。这个集合只可能偏大：多出一个伪临界点会使结果扩大到 $\pm 1$，而遗漏一个则会破坏包含性。当它有四个或更多元素时，所有余数都会出现。sinpi、cospi 和 tanpi 使用精确的临界点 $k/2$，它们由二进有理端点计算得到，不涉及对 $\pi$ 的任何近似。

• 极点。 $\tan$ 在相邻极点 $(2k+1)\pi/2$ 之间连续且递增。若 $K$ 含有奇数下标，则区间可能包含极点，结果为 Entire；否则结果为 $[\tan\underline{x}, \tan\overline{x}]$。对于极点位置精确已知的 tanpi，位于端点上的极点改为给出半无界结果。负幂和除法按前几节的规则处理 0 处的极点。

pow_interval(x, y) 在定义域 $\xi > 0$ 上（以及 $\xi = 0$、$\eta > 0$）利用了如下事实：在 $\ln \xi$ 和 $\eta$ 的符号固定时，$\xi^\eta$ 对每个参数都是单调的，因此它在盒子上的极值位于四个角点之中，另外当盒子跨越 $\xi = 1$ 或 $\eta = 0$（单调方向在此改变）时加上值 1，当 $\xi \to 0$ 时加上极限 0 和 $+\infty$。atan2 也以同样方式处理，额外将与坐标轴的交点作为候选，当盒子跨越负 $\xi$ 轴上的分支切割时结果为 $[-\pi, \pi]$。

超越函数端点的认证求值

$e^{\underline{x}}$ 不是二进有理数，因此端点需要精确值的一个有理包络 $L \le f(a) \le U$，之后由性质 (i)，$\operatorname{RD}_p(L)$ 就是一个有效的下端点。本包的内核在工作精度 $w > p$ 下，用精确有理算术和定向 BinFloat 运算构建 $L$ 和 $U$：

• exp。 对于 $|a| < 2^{e+1}$，将参数减半 $k = e + 4$ 次使得 $0 \le |a|/2^k \le 1/8$，以定向舍入对 Taylor 级数求和，再以区间算术将结果平方 $k$ 次，$e^{a} = (e^{a/2^k})^{2^k}$。每次平方使相对宽度加倍，因此工作精度为 $w = p + 64 + 2k$。项 $t_n$ 之后的级数尾项利用 $t_{j+1}/t_j = a'/(j+1) \le 1/16$ 加以界定： $$\sum_{j > n} t_j \le t_{n+1}\sum_{i \ge 0} 16^{-i} = \tfrac{16}{15} t_{n+1} \le 2 t_{n+1}.$$ 负参数使用 $e^{-a} = 1/e^{a}$。对于 $|a| \ge 2^{30}$，在整个 BinFloat 指数范围内 $|a| > (e_{\max}+1)\ln 2$，结果直接为 $[\text{largest finite}, +\infty)$ 或 $[0, \text{smallest positive}]$。
• ln。 对于满足 $m \in [1, 2)$ 的 $a = m \cdot 2^{e}$，$\ln a = \ln m + e \ln 2$，且 $\ln m = 2\operatorname{artanh} z = 2\sum_k z^{2k+1}/(2k+1)$，其中 $z = (m-1)/(m+1) \in [0, 1/3]$；$\ln 2$ 是在 $m = 2$ 处的同一级数。相邻各项至少按 $z^2 \le 1/9$ 缩小，因此对第一个省略项 $t'$，加倍后省略的尾项至多为 $\tfrac{9}{4} t' \le 3t'$。
• π。 使用 Machin 公式 $\pi = 16\arctan\tfrac15 - 4\arctan\tfrac1{239}$ 及交错级数，其截断误差以第一个省略项为界。
• sin、cos。 参数按象限 $q = \lfloor 2a/\pi + 1/2 \rfloor$ 约简，该象限同时用 $\pi^-$ 和 $\pi^+$ 计算；若两者不一致，则提高工作精度。约简后的参数 $r = a - q\pi/2 \in [-\pi/4, \pi/4]$（作为有理区间）送入 Taylor 级数，象限则将 $(\sin r, \cos r)$ 映射到 $(\sin a, \cos a)$。初始工作精度为 $p + 96$ 加上 $|a|$ 的整数位数，使得约简巨大参数时仍保留足够的位数。³3 Payne–Hanek 约简可以省去这些额外位数；本包则改为提高工作精度，并以在全函数形式和 try_ 形式下所述的资源上限加以封顶。

三角函数和反正切内核增加了 Ziv 风格的接受测试。⁴4 A. Ziv，“Fast evaluation of elementary mathematical functions with correctly rounded last bit”，ACM TOMS 17(3)，1991；J.-M. Muller 等，Handbook of Floating-Point Arithmetic，第 2 版，Birkhäuser，2018，§10。 由于 $\operatorname{RD}_p$ 是单调的，

$$\operatorname{RD}_p(L) = \operatorname{RD}_p(U) \implies \operatorname{RD}_p(L) \le \operatorname{RD}_p(f(a)) \le \operatorname{RD}_p(U) = \operatorname{RD}_p(L),$$

因此当 $L$ 和 $U$ 的 $\operatorname{RD}$ 与 $\operatorname{RU}$ 都一致时，端点就是 $f(a)$ 的正确定向舍入。否则工作精度 $w$ 增长到 $w + \max(32, w/2)$，至多 12 次（CertifiedRefinementBudget）。exp 和 log 内核跳过该测试并保留 $\operatorname{RD}_p(L)$、$\operatorname{RU}_p(U)$，它们总是有效的，且在 64 个保护位下几乎总是紧致的。try_ 形式，以及 expm1、log1p、sinpi、cospi、tanpi、pow_interval、hypot 和 atan2 的全函数形式（它们会先尝试前者），将端点委托给 bin_float 中经认证的 try_*_ctx 函数，在分别朝 $-\infty$ 和 $+\infty$ 舍入的无界上下文中求值。全函数形式的双曲函数以及 asin/acos 在多出 64 到 192 位的精度下以区间算术计算其定义公式，由基本定理可知这是有效的。

IEEE 1788 装饰模型

单纯的区间结果只说明值位于何处，而不说明函数是否有定义。IEEE 1788 为在盒子 $\boldsymbol{x}$ 上对 $f$ 求值的每个结果附加一个装饰：

装饰	性质 $p_d(f, \boldsymbol{x})$
com	$\boldsymbol{x} \subseteq D_f$，$f$ 在 $\boldsymbol{x}$ 上连续，且结果有界
dac	$\boldsymbol{x} \subseteq D_f$ 且 $f \vert_{\boldsymbol{x}}$ 连续
def	$\boldsymbol{x} \subseteq D_f$
trv	恒为真
ill	该值为 NaI，不是区间

装饰按强度全序排列，$\text{com} > \text{dac} > \text{def} > \text{trv} > \text{ill}$，因为每条性质都蕴含下一条。将装饰运算 $g$ 作用于装饰输入 $(\boldsymbol{y}_j, d_j)$ 时返回

$$d = \min(d_1, \dots, d_k, d_g), \qquad d_g = \text{strongest } d \text{ with } p_d(g, \boldsymbol{y}_1 \times \cdots \times \boldsymbol{y}_k).$$

为何取最小值是可靠的。 假设 $\boldsymbol{y}_j$ 由 $f_j$ 在 $\boldsymbol{x}$ 上以性质 $p_{d_j}$ 产生，且 $g$ 在 $\boldsymbol{y}_j$ 的盒子上具有 $p_{d_g}$。若 $d \ge \text{def}$，则每个 $f_j$ 都在 $\boldsymbol{x}$ 上有定义且取值于 $\boldsymbol{y}_j$（包含性质），而 $g$ 在这些值上有定义，因此 $g \circ (f_1, \dots, f_k)$ 在 $\boldsymbol{x}$ 上有定义。若 $d \ge \text{dac}$，连续性同理成立，因为连续函数的复合仍是连续的。com 所要求的有界性是最终结果的性质，在最终结果上检查。由归纳，整个表达式的装饰是关于该表达式在输入盒子上的一个真命题。这正是区间存在性证明所需要的：例如，若 $F(\boldsymbol{x}) \subseteq \boldsymbol{x}$ 且装饰至少为 dac，则函数在 $\boldsymbol{x}$ 上连续并将其映入自身，于是由 Brouwer 定理可知在 $\boldsymbol{x}$ 中存在不动点。若没有装饰，$\sqrt{[-1, 4]} = [0, 2]$ 会错误地暗示 $\sqrt{\cdot}$ 在 $[-1, 4]$ 上有定义。

本包通过 API 参考 中列出的定义域测试，由操作数计算 $d_g$：除以包含 0 的区间、对数到达 $\xi \le 0$、sqrt 低于 0 等情形给出 trv；atan2 跨越其分支切割时给出 def（有定义但不连续），从上方触及切割则给出 dac。集合运算（intersection、convex_hull、cancel_*）不是逐点函数，总是给出 trv。结果会被规范化：空结果总是 trv（对空集的原像，无法就 $f$ 断言任何事），无界结果上的 com 变为 dac（上溢正是这样报告的）。NaI 是无效装饰构造的结果，它吸收所有运算，并且不同于 $\emptyset$——后者是一个有效的集合。

关系：必然与可能

区间代表一个未知的点，因此两个区间的比较是一个带量词的问题。两种量词给出有用的关系：

$$\begin{aligned} \forall \xi \in \boldsymbol{x}, \forall \eta \in \boldsymbol{y} : \xi < \eta &\iff \sup\boldsymbol{x} < \inf\boldsymbol{y} \iff \overline{x} < \underline{y} && \texttt{definitely\_lt} \\ \exists \xi \in \boldsymbol{x}, \exists \eta \in \boldsymbol{y} : \xi = \eta &\iff \boldsymbol{x} \cap \boldsymbol{y} \ne \emptyset \iff \underline{x} \le \overline{y} \wedge \underline{y} \le \overline{x} && \texttt{maybe\_eq} \end{aligned}$$

对于第一行，"$\Leftarrow$" 即 $\xi \le \overline{x} < \underline{y} \le \eta$。对于 "$\Rightarrow$"：若 $\overline{x}$ 和 $\underline{y}$ 有限，则它们是元素，因此 $\overline{x} < \underline{y}$；若 $\overline{x} = +\infty$ 或 $\underline{y} = -\infty$，则足够大的 $\xi$ 或足够小的 $\eta$ 会违反 $\xi < \eta$，而端点测试同样为假。第二行是两个区间的非空交集。“可能小于” 是 “必然不小于” 的否定，因此这两族关系是对偶的：definitely_lt(x, y) 为假当且仅当存在某对值满足 $\xi \ge \eta$。对于空操作数，全称命题会空真成立，这将使调用者能从空包络证明任何结论；因此 definitely_* 关系对空操作数返回 false，而 IEEE 1788 的关系 precedes（$\forall\xi\,\forall\eta: \xi \le \eta$）则被定义为空真。

集合关系（subset、interior、disjoint、set_equal）和 IEEE 1788 序（less：$\forall\xi\,\exists\eta\, \xi \le \eta$ 与 $\forall\eta\,\exists\xi\, \xi \le \eta$，归结为比较两对端点）同样通过端点比较求值。这些关系都不是全序，trait 方法 @lf_arith.Contains::contains(x, y) 是集合包含 $\boldsymbol{y} \subseteq \boldsymbol{x}$，与 arithmetic trait 的包络解读一致。

设计决策

存储端点，而非中点与半径

问题。 球可以存储为端点 $[\underline{x}, \overline{x}]$（inf–sup），也可以像 Arb 那样存储为中点与半径 $\langle m, r\rangle$。⁵5 J. van der Hoeven，“Ball arithmetic”，2009；F. Johansson，“Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”，IEEE Trans. Computers 66(8)，2017。

中点–半径算术。 设 $x = m_x + \delta_x$，则 $|\delta_x| \le r_x$，$y$ 同理：

$$\begin{aligned} x + y &= (m_x + m_y) + (\delta_x + \delta_y), & |\delta_x + \delta_y| &\le r_x + r_y, \\ xy - m_x m_y &= m_x\delta_y + m_y\delta_x + \delta_x\delta_y, & |xy - m_x m_y| &\le |m_x| r_y + |m_y| r_x + r_x r_y . \end{aligned}$$

对于舍入后的中点 $m = \operatorname{RN}_p(m_x \circ m_y)$，需加上舍入误差 $|m - m_x \circ m_y| \le 2^{-p}|m|$，因此

$$r_{x+y} = \operatorname{RU}\bigl(r_x + r_y + 2^{-p}|m|\bigr), \qquad r_{xy} = \operatorname{RU}\bigl(|m_x| r_y + |m_y| r_x + r_x r_y + 2^{-p}|m|\bigr).$$

这些运算代价低廉（半径只需少量位数），但乘积的半径会高估：对 $\boldsymbol{x} = \boldsymbol{y} = \langle 1, 1\rangle = [0, 2]$ 它给出 $\langle 1, 3\rangle = [-2, 4]$，而精确乘积为 $[0, 4]$。Rump 证明了与 inf–sup 相比，每次乘法宽度可能增大至多 $1.5$ 倍。⁶6 S. M. Rump，“Fast and parallel interval arithmetic”，BIT 39(3)，1999。

可选方案。 (a) 中点–半径存储，半径采用低精度；(b) 端点存储，端点采用全精度；(c) 两者兼有。

选择：(b)。 IEEE 1788 是基于端点定义的，半无界集合和空集没有中点–半径形式，并且端点公式中的凸包在舍入前是精确的，因此 inf–sup 结果在精度允许的范围内尽可能紧致。代价是两个端点都携带 $p$ 位，因此高精度下的宽区间会存储大量无用的位。中点–半径视图（new、center、radius、with_precision）予以保留，供以 $c \pm r$ 方式推理的用户使用，其精确转换已在上文推导。

精确候选值，一次定向舍入

问题。 端点候选值既可以在每一步都用定向舍入计算，也可以精确计算后只舍入一次。

选择。 BinFloat 端点的和与积精确构成（系数会增长以容纳它们），quantize_interval 在最后施加一次 $\operatorname{RD}_p$/$\operatorname{RU}_p$；由性质 (iii)，结果是精确凸包的最紧 $p$ 位包络。除法和平方根的精确结果不是二进有理数，因此直接在精度 $p$ 下以定向舍入计算。该规则需要两项保护措施，如下所述：对相距很远的加数做有界对齐，以及在指数范围边缘向外钳制。

相距很远的加数：按精度界定端点和

问题。 精确地将 $A = 2^{10^9}$ 与 $s = 2^{-10^9}$ 相加会对齐两个系数，构造出一个二十亿位的数（一次区间加法约需 750 MB，见 issue #24）。

选择。 设 $t(v)$ 为 $v$ 最高位的指数，$e(A)$ 为较大加数 $A$ 最低位的指数，$M = \max(65536, \operatorname{prec}(A), \operatorname{prec}(s))$ 且 $c = \min(e(A), t(A) - M)$。若 $t(s) < c - 2$，当 $s$ 将和推向所计算端点的舍入方向时，小加数被替换为粘滞替身 $s' = \operatorname{sign}(s)\,2^{c-2}$，否则替换为 $s' = 0$。于是：

• 对任意精度均可靠。 $|s| < 2^{t(s)+1} \le 2^{c-2}$。当 $s > 0$ 时，上端点得到 $A + 2^{c-2} > A + s$；当 $s < 0$ 时得到 $A > A + s$。下端点对称。因此定向和始终位于精确和所要求的一侧，而这正是包含性质所需的全部。
• 当 $p \le M$ 时不损失紧致性。 $A$ 附近的 $p$ 位数是 $2^{t(A) - p}$ 的倍数，且 $t(A) - p \ge t(A) - M \ge c$，因此没有 $p$ 位数严格位于 $A$ 与 $A \pm 2^c$ 之间；$A + s$ 和 $A + s'$ 落在同一间隙中，舍入到同一端点（附件中的引理 7）。

现在对齐的代价至多比 $A$ 的宽度多出约 $M$ 位，与指数差无关。该替身也用于 center 内部的就近舍入求和，在那里它并不精确；这正是正确性中所述局限的来源。

在指数范围边缘向外钳制

BinFloat 具有有限（但非常宽）的指数范围。超出该范围的精确候选值无法存储，因此精确辅助函数会采用所构建端点的方向：下端点朝 $-\infty$ 钳制，上端点朝 $+\infty$ 钳制，半径向上舍入，指数之和以 64 位计算并饱和。例如，这能使 $\exp([10^9, 10^9])$ 保持为包络 $[\text{largest finite}, +\infty)$，而不会坍缩。

上下文与双重舍入

BallContext 携带目标精度 $q$ 和指数范围；apply_ctx 将端点向外舍入到其中，*_ctx 运算先在操作数的精度 $p$ 下计算，再应用上下文。当 $q \le p$ 时，对定向舍入而言舍入两次是无害的：

$$\operatorname{RD}_q(\operatorname{RD}_p(a)) = \operatorname{RD}_q(a) \quad (q \le p),$$

因为 $F_q \subseteq F_p$（$q$ 位有效数就是用零填充的 $p$ 位有效数）：$a$ 以下的每个 $f \in F_q$ 都属于 $F_p$，因此也在 $\operatorname{RD}_p(a)$ 以下，反之亦然。所以只要操作数至少与上下文一样精确，x.add_ctx(y, ctx) 就等于将精确和一次性向外舍入到上下文中。上溢按包络方向处理：超出范围的上端点变为 $+\infty$，正的下端点变为最大有限数（仍低于真实值）。下溢在次正规数网格上向外舍入，因此极小的正上界变为最小次正规数，而绝不会变为 0。标志作为返回值给出，从不全局存储。

全函数形式与 try_ 形式

问题。 认证求值可能失败：细化预算可能耗尽，或参数可能大到无法约简。中止会打断长时间的计算；而悄无声息地返回宽区间则会掩盖紧致性保证的丧失。

选择。 两者兼备。全函数形式返回一个有效的后备包络（sin/cos 返回 $[-1, 1]$，tan 返回 Entire，atan2 返回 $[-\pi, \pi]$，pow、hypot 和 rootn 返回复合公式，expm1 返回值域界 $[-1, +\infty)$），从而保持包含性质；try_ 形式返回带有认证细节（运算、阶段、原因、精度）的 ArithmeticError。三角函数约简设有上限：当较大端点的绝对值至少为 $2^{\max(65536, 4p)+1}$ 时，约简需要超过这么多位，因此全函数形式立即返回后备值，try_ 形式则报告资源限制。

装饰放在独立的类型中

问题。 装饰在每次运算中都要多耗费一个字段和一次取最小值，而大多数用户并不需要它们。

选择。 BallFloat 是裸集合；BallFloatDecorated 用装饰和 NaI 状态将其包装。裸运算保持廉价、简单，且装饰类型不会被意外地与裸值混用。

精度作为标签

每个区间都携带一个精度标签；二元运算使用较大的标签。这使一连串运算无需上下文参数即可保持在最精确输入的精度上，而集合值类型要与运算符（x + y）组合正需要这一点。需要强制指定特定格式时，使用 BallContext。

正确性 / 不变式

表示不变式。 非空的 BallFloat 具有非 NaN 端点、$\underline{x} \le \overline{x}$、$\underline{x} \ne +\infty$、$\overline{x} \ne -\infty$ 以及精度 $\ge 1$；每条构造路径都以 store_interval 结束，它检查上述条件，不满足则中止。空区间是一个标志，其端点为 $(+\infty, -\infty)$。

包含性。 BallFloat 的每个公开运算 $F$ 都满足 $f(\boldsymbol{x}) \subseteq F(\boldsymbol{x})$：对算术运算依据端点公式和推论 2，对初等函数依据临界点分析和经认证的端点包络，对未认证的情形依据后备值（例外列于已知局限中）。由基本定理，每个复合运算也满足这一点。

紧致性。 基本算术、square、pown、fma、abs、minimum、maximum、集合运算和 sqrt_interval 返回精确凸包的向外舍入，因此每个端点与最优值相差不超过一个 ulp。三角函数和反正切的端点在通过认证时是正确的定向舍入；其他初等函数相差在几个 ulp 以内；后备值不紧致。

装饰。 由上文的归纳论证，结果的装饰是关于所求值函数的一个真命题。

复杂度。 算术运算至多执行四次端点乘积或两次求和，对于 $p$ 位端点、乘法代价为 $M(p)$ 时为 $O(M(p))$，求和另需对齐至多约 $65536 + p$ 位。初等函数在工作精度 $w = p + O(1)$ 下对 $O(w)$ 项级数求和（三角函数约简还需加上参数的整数位数），并最多进行 12 次细化，使 $w$ 按几何级数增长。

证据。 包测试检查定向端点、相距很远的加数与指数范围的情形、装饰和关系；固定版本的 ITF1788 语料在 strict 模式下运行 4,656 个用例（参见符合性）。

已知局限

以下输入目前会破坏包含性质或装饰规则；它们已被报告待修复，并记录于此，以便调用者避开。

• from_int(n, precision=p) 和 from_coefficient 在构建单点区间之前先将整数在 $\max(p, 8)$ 位下就近舍入，因此 from_int(257, precision=8) 为 $\{256\}$。
• with_precision（因而 normalized 也）从 center() 重建有界区间，而后者使用带就近舍入的远加数替身。当端点相距超过约 $2^{16}$ 个二进制数量级，且新精度大到足以精确存储该替身（超过约 65536 位）时，结果可能丢失较小的端点。
• 负次数的装饰 rootn 在参数包含 0 时不会将装饰降为 trv。
• midpoint_ctx 不应用上下文的 $e_{\max}$，且从不引发 overflow。

被否决的替代方案

• 中点–半径存储（Arb 风格）：半径更廉价，但会高估乘积，且无法表示半无界集合；参见上文。
• 使用硬件 Double 端点并切换舍入模式： 精度固定、依赖进程全局状态，且并非每个 MoonBit 目标都支持舍入模式控制。
• 就近舍入后膨胀一个 ulp（“epsilon 膨胀”）：更简单，但每个端点比定向舍入宽约一个 ulp，且只对就近舍入结果已知误差在一个 ulp 以内的运算有效，这排除了大多数初等函数内核。
• 对未认证的初等函数中止： 单个困难参数就会使整个计算停止；需要获知失败的调用者可使用 try_ 形式。
• 将除以含零区间视为错误： IEEE 1788 将这些结果定义为集合（Entire、半无界或空集），而扩展除法正是区间 Newton 方法得以奏效的关键。

边界

本包有意不做以下事情：

• 提供 IEEE 1788 逆运算（sqrRev、mulRevToPair、…）或双输出除法；
• 承诺紧致的结果：后备值和依赖问题可能使包络任意变宽；
• 在区间上定义全序，或将 Eq 视为集合相等；
• 自行将十进制数据向外转换：from_double 和 exact 包络的是二进制值，包络十进制字面量是调用者的职责（参见教程）；
• 实现以十进制数为端点的区间、复球、区间向量或矩阵，或 Taylor 模型；
• 通过全局标志报告状态：标志由 *_ctx 调用返回，装饰随值一起传递。

Footnotes

1. IEEE Std 1788-2015，IEEE Standard for Interval Arithmetic。本包遵循其基于集合的风格（set-based flavor）。 ↩

2. R. E. Moore，Interval Analysis，Prentice-Hall，1966；以及 Moore、Kearfott、Cloud，Introduction to Interval Analysis，SIAM，2009，定理 5.1。 ↩

3. Payne–Hanek 约简可以省去这些额外位数；本包则改为提高工作精度，并以在全函数形式和 try_ 形式下所述的资源上限加以封顶。 ↩

4. A. Ziv，“Fast evaluation of elementary mathematical functions with correctly rounded last bit”，ACM TOMS 17(3)，1991；J.-M. Muller 等，Handbook of Floating-Point Arithmetic，第 2 版，Birkhäuser，2018，§10。 ↩

5. J. van der Hoeven，“Ball arithmetic”，2009；F. Johansson，“Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”，IEEE Trans. Computers 66(8)，2017。 ↩

6. S. M. Rump，“Fast and parallel interval arithmetic”，BIT 39(3)，1999。 ↩