ball_float の設計

ball_float は単一の近似値ではなく、実数の集合を用いて計算します。このページでは数学的な契約を述べ、コードが用いる公式を導き、その背後にある選択を説明し、パッケージが行わないことを列挙します。各関数は API リファレンスに記載されており、チュートリアルではその使い方を示しています。

設計目標

浮動小数点計算は真の結果に近い数を返しますが、その誤差は呼び出し側が別途見積もらなければなりません。ball_float は任意の作業精度で、真の結果を含むことが証明可能な区間を返すため、誤差限界が値の一部になります。ここから 3 つの要件が導かれます。

1. 緊密さより包含を優先。 すべての演算は、オペランドの厳密な像の上位集合を返します。より広い結果は許容されますが、取りうる値を取りこぼす結果はバグです。
2. 明示的な意味論。 精度、対象形式、フラグ、装飾は受け渡しされる値であり、プロセス全体の大域状態ではありません（ハードウェアの丸めモード切り替えは行いません）。
3. 標準の語彙。 集合、関係、装飾は IEEE 1788-2015¹1 IEEE Std 1788-2015, IEEE Standard for Interval Arithmetic。このパッケージはその集合ベースのフレーバーに従います。 に従うため、結果をそのテストコーパスと照合できます（適合性を参照）。

数学的背景

区間と包含性

区間とは、$\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 ,$$

また、その包（hull）$\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, Thm. 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 数なので、厳密な端点の値は丸める必要があります。$p$ ビットの 2 進数の集合を $F_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 \}.$$

これらの定義から直ちに 3 つの性質が導かれます。

$$\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）。

丸めによって加わる幅は、端点ごとに高々 1 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 は 2 つの端点を保存する（中点と半径ではなく端点を参照）ので、各演算は端点に関する公式になります。

和と差。 $\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$、残りの 2 つは $\ge 0$ です。最小値は前者の中に、最大値は後者の中にあります。非有界な区間では、4 つの積すべてを $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})$ です。コードは 2 回丸める代わりに、選ばれた端点の商を方向付き除算で直接評価します。$0 \in \boldsymbol{y}$ のときは次のとおりです。

• $\boldsymbol{y} = \{0\}$：許容される $\eta$ がないので、結果は $\emptyset$ です。
• $\underline{y} < 0 < \overline{y}$、$\boldsymbol{x} \ne \{0\}$：$\eta$ は両側から 0 に近づくので、$\xi/\eta$ は両方向に非有界であり、結果は Entire です。
• $\overline{y} > 0$ かつ $\underline{x} > 0$ で $\boldsymbol{y} = [0, \overline{y}]$ の場合：$\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$ を返し、これらは 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 のもう 1 つの定理）。同じ効果により、$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 および初等関数 — を提供しており、その結果は独立な因子の積ではなく、（丸めを除いて）真の値域の包になります。また区間値は加法について群をなしません。一般に $(\boldsymbol{x} + \boldsymbol{y}) - \boldsymbol{y} \ne \boldsymbol{x}$ なので、加算を打ち消す演算は cancel_minus です。

初等関数：単調性、臨界点、極

連続関数では区間上の値域は区間であり、その端点は $\boldsymbol{x}$ の端点またはその内部の臨界点で達成されます。パッケージは 3 つのパターンを用います。

• 単調関数（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$ では入れ替え）なので、端点だけを評価し、下端は切り下げ、上端は切り上げます。区間内部にある定義域の境界は、そこでの極限に置き換えられます（$\xi \downarrow 0$ のとき $\ln \xi \to -\infty$）。

• 極値が既知の関数。 cosh は 0 で最小値 1 をとります。偶数指数の pown は 0 で最小値 0 をとります。sin は $k \equiv 1 \pmod 4$ である $\xi = k\pi/2$ で $+1$ を、$k \equiv 3$ で $-1$ をとり、cos は $k \equiv 0$ で $+1$ を、$k \equiv 2$ で $-1$ をとります。連続する 2 つの臨界点の間ではどちらの関数も単調なので、$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}$$

  コードは $K$ を $[\lceil q^-(\underline{x}) \rceil, \lfloor q^+(\overline{x}) \rfloor]$ で包含します。ここで $q^- \le 2\xi/\pi \le q^+$ は $\pi$ の認証済み包含区間 $[\pi^-, \pi^+]$ から計算されます。この集合は大きすぎることしかありえません。余分な臨界点は結果を $\pm 1$ まで広げるだけですが、臨界点を見落とせば包含性が崩れます。要素が 4 つ以上あれば、すべての剰余類が現れます。sinpi、cospi、tanpi は厳密な臨界点 $k/2$ を用い、これは $\pi$ を一切近似せずに 2 進有理数の端点から計算されます。

• 極。 $\tan$ は連続する極 $(2k+1)\pi/2$ の間で連続かつ増加です。$K$ が奇数の添字を含めば区間は極を含む可能性があり、結果は Entire です。そうでなければ $[\tan\underline{x}, \tan\overline{x}]$ です。極の位置が厳密な tanpi では、端点にある極は代わりに片側が非有界な結果を与えます。負の冪と除算は、前節の規則によって 0 にある極を扱います。

定義域 $\xi > 0$（および $\xi = 0$、$\eta > 0$）上の pow_interval(x, y) は、$\ln \xi$ と $\eta$ の符号を固定すると $\xi^\eta$ が各引数について単調であることを利用します。したがってボックス上の極値は 4 つの角の中にあり、それに加えて、ボックスが $\xi = 1$ または $\eta = 0$（単調性の向きが変わる場所）をまたぐときは値 1、$\xi \to 0$ のときは極限 0 と $+\infty$ が候補になります。atan2 も同様に扱われ、軸との交差が追加の候補となり、ボックスが負の $\xi$ 軸上の分枝切断をまたぐときは結果は $[-\pi, \pi]$ になります。

超越関数の端点の認証付き評価

$e^{\underline{x}}$ は 2 進有理数ではないので、端点には厳密な値の有理数による包含 $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$ とし、テイラー級数を方向付き丸めで総和し、その結果を区間演算で $k$ 回 2 乗します（$e^{a} = (e^{a/2^k})^{2^k}$）。2 乗するたびに相対幅が 2 倍になるので、作業精度は $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$ であり、$z = (m-1)/(m+1) \in [0, 1/3]$ として $\ln m = 2\operatorname{artanh} z = 2\sum_k z^{2k+1}/(2k+1)$ です。$\ln 2$ は $m = 2$ における同じ級数です。連続する項は少なくとも $z^2 \le 1/9$ の比で縮小するので、2 倍した後の省略された剰余は、最初に省略された項 $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]$（有理数区間として）はテイラー級数に渡され、象限によって $(\sin r, \cos r)$ が $(\sin a, \cos a)$ に対応付けられます。初期作業精度は $p + 96$ に $|a|$ の整数部のビット数を加えたもので、巨大な引数を還元しても十分なビットが残るようにしています。³3 Payne–Hanek 還元を使えば余分なビットは不要になりますが、このパッケージでは代わりに作業精度を引き上げ、total 形式と try_ 形式の項で述べるリソース上限によってそれを制限しています。

三角関数と逆正接のカーネルは、Ziv 方式の受理判定を追加します。⁴4 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, §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 の total 形式（まずこちらを試みます）は、端点の計算を bin_float の認証付き try_*_ctx 関数に委ね、$-\infty$ 方向と $+\infty$ 方向に丸める無制限のコンテキストで評価します。total の双曲線関数と asin/acos は、その定義式を 64 から 192 ビットの追加精度で区間演算により評価します。これは基本定理により正当です。

IEEE 1788 の装飾モデル

単なる区間の結果は、値がどこにあるかを示すだけで、関数が定義されていたかどうかは示しません。IEEE 1788 は、ボックス $\boldsymbol{x}$ 上で $f$ を評価した各結果に装飾（decoration）を付加します。

装飾	性質 $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}$ と全順序付けられます。装飾付きの入力 $(\boldsymbol{y}_j, d_j)$ に適用された装飾付き演算 $g$ は次を返します。

$$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$ が $\boldsymbol{x}$ 上の $f_j$ によって性質 $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 の有界性は最終結果の性質であり、最終結果に対して検査されます。帰納法により、式全体の装飾は入力ボックス上の式についての真の命題になります。これは区間による存在証明が必要とするものです。たとえば装飾が dac 以上で $F(\boldsymbol{x}) \subseteq \boldsymbol{x}$ ならば、関数は $\boldsymbol{x}$ 上で連続かつそれをそれ自身の中へ写すので、Brouwer の定理により $\boldsymbol{x}$ 内に不動点が存在します。装飾がなければ、$\sqrt{[-1, 4]} = [0, 2]$ は $\sqrt{\cdot}$ が $[-1, 4]$ 上で定義されていると誤って示唆してしまいます。

パッケージは、API リファレンスに列挙された定義域判定によって、オペランドから $d_g$ を計算します。0 を含む区間による除算、$\xi \le 0$ に達する対数、0 未満での sqrt などは trv を与えます。分枝切断をまたぐ atan2 は def（定義されているが不連続）を与え、上側から切断に接する場合は dac を与えます。集合演算（intersection、convex_hull、cancel_*）は点関数ではないので、常に trv を与えます。結果は正準化されます。空の結果は常に trv であり（空集合の逆像上の $f$ については何も主張できません）、非有界な結果の com は dac になります（オーバーフローはこのようにして報告されます）。不正な装飾付き構築の結果である NaI はすべての演算を吸収し、有効な集合である $\emptyset$ とは区別されます。

関係：確実に（certainly）と可能性として（possibly）

区間は未知の点を表すので、2 つの区間の比較は量化された問いになります。2 つの量化子から有用な関係が得られます。

$$\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}$$

1 行目について、「$\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$ を破り、端点の判定も偽になります。2 行目は 2 つの区間の共通部分が空でないことです。「可能性として小さい」は「確実に小さくない」の否定なので、2 つの族は双対です。definitely_lt(x, y) が偽になるのは、ある組が $\xi \ge \eta$ を満たすときに限ります。空のオペランドでは全称命題は空虚に真となり、呼び出し側が空の包含区間から何でも証明できてしまいます。そのため definitely_* の関係は空のオペランドに対して偽を返します。一方、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$ で、両方の端点の組を比較することに帰着します）も、同様に端点の比較によって評価されます。これらの関係はいずれも全順序ではありません。またトレイトメソッド @lf_arith.Contains::contains(x, y) は集合の包含 $\boldsymbol{y} \subseteq \boldsymbol{x}$ であり、arithmetic のトレイトにおける包含区間としての解釈と一致します。

設計上の判断

中点と半径ではなく端点

問題。 ボールは端点 $[\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 回あたり最大 $1.5$ 倍まで幅が増大しうることを示しました。⁶6 S. M. Rump, “Fast and parallel interval arithmetic”, BIT 39(3), 1999。

選択肢。 (a) 低精度の半径を伴う中点・半径による保存、(b) 完全な精度の端点による保存、(c) その両方。

採用：(b)。 IEEE 1788 は端点に基づいて定義されており、片側非有界な集合や空集合には中点・半径形式が存在せず、端点の公式における包は丸める前は厳密なので、inf–sup の結果は精度が許す限り緊密です。代償として両方の端点が $p$ ビットを持つため、高精度の広い区間では無駄なビットを多く保存することになります。$c \pm r$ で考える利用者のために中点・半径のビュー（new、center、radius、with_precision）を残しており、その厳密な変換は上記で導いています。

厳密な候補と、一度だけの方向付き丸め

問題。 端点の候補は、各ステップで方向付き丸めを行って計算することも、厳密に計算して一度だけ丸めることもできます。

採用。 BinFloat の端点の和と積は厳密に作られ（係数はそれを保持できるまで大きくなります）、quantize_interval が最後に $\operatorname{RD}_p$/$\operatorname{RU}_p$ を一度だけ適用します。性質 (iii) により、結果は厳密な包の最も緊密な $p$ ビット包含区間になります。厳密な結果が 2 進有理数にならない除算と平方根は、精度 $p$ の方向付き丸めで直接計算されます。この規則には、次に述べる 2 つの安全策が必要です。大きく離れた加数の桁合わせの制限と、指数範囲の端での外向きクランプです。

離れた加数：端点の和を精度で抑える

問題。 $A = 2^{10^9}$ と $s = 2^{-10^9}$ を厳密に加算すると、2 つの係数の桁を合わせて 20 億ビットの数を作ることになります（区間の加算 1 回で約 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$ なので、$A$ と $A \pm 2^c$ の間に真に挟まる $p$ ビット数は存在しません。$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$ のとき、方向付き丸めでは 2 回丸めても害はありません。

$$\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 にはなりません。フラグは返されるのであり、大域的に保存されることはありません。

total 形式と try_ 形式

問題。 認証付き評価は失敗することがあります。精緻化の予算が尽きることもあれば、引数が大きすぎて還元できないこともあります。中断すれば長い計算が台無しになり、黙って広い区間を返せば、緊密さの保証が失われたことが隠されてしまいます。

採用：両方。 total 形式は有効なフォールバックの包含区間を返し（sin/cos には $[-1, 1]$、tan には Entire、atan2 には $[-\pi, \pi]$、pow、hypot、rootn には合成した式、expm1 には値域の範囲 $[-1, +\infty)$）、包含性を保ちます。try_ 形式は認証の詳細（演算、段階、理由、精度）を伴う ArithmeticError を返します。三角関数の還元には上限があります。大きい方の端点の絶対値が $2^{\max(65536, 4p)+1}$ 以上になると、還元にはそれ以上のビット数が必要になるため、total 形式は直ちにフォールバックを返し、try_ 形式はリソース上限を報告します。

装飾は別の型に

問題。 装飾には演算ごとにフィールド 1 つと最小値の計算 1 回のコストがかかり、ほとんどの利用者はそれを必要としません。

採用。 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 は厳密な包の外向き丸めを返すので、各端点は最適値から 1 ulp 以内です。三角関数と逆正接の端点は、認証された場合は正しい方向付き丸めです。その他の初等関数は数 ulp 以内であり、フォールバックは緊密ではありません。

装飾。 上記の帰納法の議論により、結果の装飾は評価された関数についての真の命題です。

計算量。 算術は高々 4 つの端点の積または 2 つの和を計算し、乗算コストを $M(p)$ とすると $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 進で約 $2^{16}$ 桁以上離れていて、新しい精度が代理値を厳密に保存できるほど大きい（約 65536 ビットを超える）場合、結果は小さい方の端点を失うことがあります。
• 負の次数を持つ装飾付き rootn は、引数に 0 が含まれるときに装飾を trv に下げません。
• midpoint_ctx はコンテキストの $e_{\max}$ を適用せず、overflow を発生させることもありません。

却下した代替案

• 中点・半径による保存（Arb 方式）：半径は安価ですが、積を過大評価し、片側非有界な集合を表現できません。上記を参照してください。
• 丸めモードを切り替えたハードウェアの Double 端点： 精度が固定され、プロセス全体の大域状態となり、しかも丸めモードの制御はすべての MoonBit ターゲットで利用できるわけではありません。
• 最近接に丸めて 1 ulp 膨らませる（「イプシロン膨張」）：より単純ですが、方向付き丸めより端点ごとに約 1 ulp 広くなり、最近接丸めの結果が 1 ulp 以内であることがわかっている演算でしか有効でないため、初等関数のカーネルの大半が除外されます。
• 認証できない初等関数で中断すること： 難しい引数が 1 つあるだけで計算全体が止まってしまいます。失敗を知る必要がある呼び出し側には try_ 形式があります。
• ゼロを含む区間による除算をエラーとして扱うこと： IEEE 1788 はこれらの結果を集合（Entire、片側非有界、または空）として定義しており、拡張された除算こそが区間ニュートン法を機能させるものです。

境界

このパッケージは意図的に以下を行いません。

• IEEE 1788 の逆演算（sqrRev、mulRevToPair、…）や 2 出力の除算を提供すること。
• 緊密な結果を約束すること。フォールバックと依存性問題によって包含区間は任意に広がりえます。
• 区間上の全順序を定義すること、あるいは Eq を集合の等価性として扱うこと。
• 10 進データを独自に外向き変換すること。from_double と exact は 2 進の値を包含するものであり、10 進リテラルを包含するのは呼び出し側の役割です（チュートリアルを参照）。
• 10 進の端点による区間、複素ボール、区間ベクトルや区間行列、テイラーモデルを実装すること。
• ステータスを大域フラグで報告すること。フラグは *_ctx 呼び出しから返され、装飾は値とともに受け渡されます。

Footnotes

1. IEEE Std 1788-2015, IEEE Standard for Interval Arithmetic。このパッケージはその集合ベースのフレーバーに従います。 ↩

2. R. E. Moore, Interval Analysis, Prentice-Hall, 1966; および Moore, Kearfott, Cloud, Introduction to Interval Analysis, SIAM, 2009, Thm. 5.1。 ↩

3. Payne–Hanek 還元を使えば余分なビットは不要になりますが、このパッケージでは代わりに作業精度を引き上げ、total 形式と try_ 形式の項で述べるリソース上限によってそれを制限しています。 ↩

4. 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, §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。 ↩