core 設計

設計目標

arithmetic は代数構造と具体的な数との間の層である。luna-generic は型が何であるか（環、体）を述べ、arithmetic は型がどの解析的演算をできるか、そしてそれをどれだけ誠実に行えるかを述べる。平方根が黙って NaN を返してよいのか、拒否した引数を報告するのか、さらに指定された精度のもとで結果がどう丸められたかまで報告するのか、である。汎用アルゴリズムは依存する能力を正確に述べ、数値バックエンドは意味論を守れる能力だけを実装する。

パッケージは語彙（トレイト、コンテキスト、診断、エラー）と、ネイティブな Float、Double、整数型のための基本的なインスタンスを提供する。正しく丸められ、コンテキストに忠実で、証明付きの算術は数値バックエンドにあり、それらがこれらのトレイトを実装する。

数学的背景

浮動小数点形式

基数 $\beta$、精度 $p$、指数範囲 $[e_{\min}, e_{\max}]$ の浮動小数点形式は、次の有限集合である

$$\mathbb{F} = \{0\} \cup \{\, \pm m \cdot \beta^{\,e-p+1} \;:\; m \in \mathbb{Z},\ \beta^{p-1} \le m < \beta^{p},\ e_{\min} \le e \le e_{\max} \,\} \cup \{\, \pm m \cdot \beta^{\,e_{\min}-p+1} \;:\; 0 < m < \beta^{p-1} \,\}.$$

二つ目の集合は正規化数で、先頭の桁は $\beta^{e}$ の位にある（$e$ は調整済み指数）。三つ目は $\beta^{e_{\min}}$ 未満の非正規化数である。IEEE 754 は $\mathbb{F}$ を符号付きゼロとともに $\overline{\mathbb{F}} = \mathbb{F} \cup \{-\infty, +\infty, \mathrm{NaN}\}$ へ拡張する。FpClass は次の写像である

$$\operatorname{class} : \overline{\mathbb{F}} \to \{\texttt{Finite}, \texttt{Infinity}, \texttt{NaN}\}, \qquad \operatorname{class}(x) = \begin{cases} \texttt{Finite} & x \in \mathbb{F},\\ \texttt{Infinity} & x = \pm\infty,\\ \texttt{NaN} & x = \mathrm{NaN}. \end{cases}$$

形式	$\beta$	$p$	$e_{\min}$	$e_{\max}$	出所
binary32	2	24	$-126$	$127$	Float
binary64	2	53	$-1022$	$1023$	Double
decimal32	10	7	$-95$	$96$	ArithmeticContext::decimal32
decimal64	10	16	$-383$	$384$	ArithmeticContext::decimal64
decimal128	10	34	$-6143$	$6144$	ArithmeticContext::decimal128

ArithmeticContext は $p$ を precision として、調整済み指数の上下限を e_min と e_max として保存する。基数はコンテキストではなくバックエンドの型の性質である。三つの十進プリセットはいずれも $e_{\min} = 1 - e_{\max}$ を満たす。これは、最小の正規化数の逆数 $\beta^{-e_{\min}} = \beta^{e_{\max}-1}$ が有限範囲に収まるようにする IEEE 754 の関係である。

マシンイプシロン。 $\mathbb{F}$ における $1$ の次の値は $1 + \beta^{1-p}$ である。$1 = \beta^{p-1} \cdot \beta^{\,0-p+1}$ と書くと $m = \beta^{p-1}$、$e = 0$ であり、次の仮数 $m + 1$ から次が得られる

$$(\beta^{p-1} + 1)\,\beta^{1-p} = 1 + \beta^{1-p} =: 1 + \varepsilon .$$

より一般に、二進桁区間（binade）$[\beta^{e}, \beta^{e+1})$ では連続する数は $\beta^{e-p+1} = \varepsilon\,\beta^{e}$ だけ離れている。epsilon_contextual は $\varepsilon$ を返す。Float では $2^{-23}$、Double では $2^{-52}$ である。

丸め

丸め関数 $\circ : \mathbb{R} \to \overline{\mathbb{F}}$ は実数を表現可能な数へ写す。$a < x < b$ が $x \notin \mathbb{F}$ の隣接値であるとき：

$$\begin{aligned} \operatorname{RN}(x) &= \text{the nearer of } a, b \text{ (ties to the even significand)} && \texttt{ToNearestEven}\\ \operatorname{RZ}(x) &= \text{the one of smaller magnitude} && \texttt{TowardZero}\\ \operatorname{RU}(x) &= b && \texttt{TowardPositive}\\ \operatorname{RD}(x) &= a && \texttt{TowardNegative}\\ \operatorname{RA}(x) &= \text{the one of larger magnitude} && \texttt{AwayFromZero} \end{aligned}$$

そしてどのモードも $x \in \mathbb{F}$ なら $x$ 自身を返す。どのモードも単調、すなわち $x \le y \Rightarrow \circ(x) \le \circ(y)$ であり、$\mathbb{F}$ 上で冪等である。後の証明の議論はこの二つの性質だけを使う。

相対誤差。 $x$ が正規範囲にあり、$e_{\min} \le e \le e_{\max}$ で $\beta^{e} \le |x| < \beta^{e+1}$ とする。$x$ の隣接値は $\varepsilon\beta^{e}$ 離れているので、有向モードは $x$ を間隔一つ未満だけ動かし、最近接丸めは高々その半分だけ動かす：

$$\begin{aligned} |\operatorname{RN}(x) - x| &\le \tfrac12\,\varepsilon\,\beta^{e} \le \tfrac12\,\varepsilon\,|x|,\\ |\circ_{\text{dir}}(x) - x| &< \varepsilon\,\beta^{e} \le \varepsilon\,|x|. \end{aligned}$$

したがって $|\delta| \le u$ として $\circ(x) = x(1 + \delta)$ であり、ここで単位丸め誤差は

$$u = \begin{cases} \tfrac12\,\beta^{1-p} & \texttt{ToNearestEven},\\ \beta^{1-p} & \text{directed modes}. \end{cases}$$

形式	ToNearestEven での $u$
binary32	$2^{-24} \approx 5.96 \times 10^{-8}$
binary64	$2^{-53} \approx 1.11 \times 10^{-16}$
decimal32	$5 \times 10^{-7}$
decimal64	$5 \times 10^{-16}$
decimal128	$5 \times 10^{-34}$

$\beta^{e_{\min}}$ 未満では間隔はもう縮まらず、相対誤差の上界は絶対誤差の上界 $|\operatorname{RN}(x) - x| \le \tfrac12\,\beta^{e_{\min}-p+1}$ になる。最大の有限値を超えると、$\circ(x)$ はモードに応じて $\pm\infty$ か最大の有限値になる。これらが ArithmeticDiagnostics の subnormal、underflow、overflow の条件である。

丸め誤差の標準モデル

IEEE 754 は $+, -, \times, /$ と $\sqrt{\ }$ が正しく丸められること、すなわち計算結果が正確な結果を丸めたものであることを要求する。上の上界から、オペランドが $\mathbb{F}$ にあり結果が正規範囲にあるとき、

$$\operatorname{fl}(x \circ y) = (x \circ y)(1 + \delta), \qquad |\delta| \le u, \qquad \circ \in \{+, -, \times, /\}.$$

これはコンテキスト付き算術トレイトの契約である。コンテキストに忠実な AddContextual はコンテキストの $p$ と丸めモードについての $\operatorname{fl}(x + y)$ を返し、$\delta \ne 0$ のときにちょうど inexact を立てる。これがコンテキストが重要な理由でもある。decimal64 と ToNearestEven のもとでは、アルゴリズムはバックエンドの型によらず演算ごとに $|\delta| \le 5 \times 10^{-16}$ を前提にできる。

組み込みの Float と Double のインスタンスは、ハードウェアの演算が正しく丸められるので、最近接丸めのもとで自身の固定形式についてこのモデルを満たす。しかしコンテキストの精度とモードを無視し、$\delta \ne 0$ を検出しない。初等関数は Kaida-Amethyst/math から来ており正しい丸めは保証されないので、それらについてはモデルは $u$ の未規定の小さな倍数でしか成り立たない。

隣接値と IEEE エンコーディング

AdjacentContextual は次を計算する

$$\operatorname{succ}(x) = \min\{\, y \in \overline{\mathbb{F}} : y > x \,\}, \qquad \operatorname{pred}(x) = \max\{\, y \in \overline{\mathbb{F}} : y < x \,\}.$$

Float と Double のインスタンスはこれらをビットパターン上で計算する。符号 $s$、バイアス付き指数フィールド $E$、仮数フィールド $F$ を持つ IEEE 二進値は、$k$ ビット形式で符号なし整数 $\operatorname{bits}(x) = s \cdot 2^{k-1} + E \cdot 2^{p-1} + F$ として格納される。非負の値の上でこの写像は順序を保つ：

• $E$ を固定すると、値 $2^{E - \text{bias}}(1 + F\,2^{1-p})$（$E = 0$ のときは $2^{1-\text{bias}}\,F\,2^{1-p}$）は $F$ について狭義単調増加である。
• フィールド $E$ を持つ最大の値は $2^{E-\text{bias}}(2 - 2^{1-p}) < 2^{E+1-\text{bias}}$ で、右辺はフィールド $E + 1$ を持つ最小の値である。非正規化数（$E = 0$）はすべて最小の正規化数 $2^{1-\text{bias}}$ より小さい。
• $+\infty$ は $E = 2^{k-p} - 1$、$F = 0$ であり、すべての有限値のパターンより上にある。

エンコーディングは $(E, F)$ について辞書式であり、$(E, F)$ の辞書式順序は $E \cdot 2^{p-1} + F$ の整数順序と一致するので、$0 \le x < y$ なら $\operatorname{bits}(x) < \operatorname{bits}(y)$ であり、連続する値の間に入るパターンは存在しない。したがって

$$\operatorname{succ}(x) = \begin{cases} \operatorname{bits}^{-1}(\operatorname{bits}(x) + 1) & x > 0,\\ \operatorname{bits}^{-1}(\operatorname{bits}(x) - 1) & x < 0 \quad (\text{since } \operatorname{succ}(x) = -\operatorname{pred}(|x|)),\\ \text{smallest positive subnormal} & x = \pm 0, \end{cases}$$

$\operatorname{pred}$ も対称的である。ゼロの場合を分けるのは、$+0$ と $-0$ が等しい値で異なるパターンを持つからである。この式は API に挙げた境界の場合を再現する。最大の有限値のパターンに 1 を足すと $+\infty$ のパターンになり、$\operatorname{succ}(1) - 1 = \varepsilon$ である。$\operatorname{succ}(x)$ は定義により $\overline{\mathbb{F}}$ に属するので丸めは起きず、これらのインスタンスの空の診断は単に未検出なのではなく正確である。

包含区間と三値比較

包含区間 $X \subseteq \mathbb{R}$ は、$x \in X$ を満たすと分かっている未知の実数 $x$ を表す。区間 $[a, b]$ やボール $B(m, r) = [m - r, m + r]$ がその例である。包含区間の関係トレイトは、許容されるすべての組について量化することで、包含区間だけから未知の値についての問いに答える：

$$\begin{aligned} \texttt{definitely\_lt}(X, Y) &\iff \forall x \in X,\ \forall y \in Y:\ x < y,\\ \texttt{definitely\_le}(X, Y) &\iff \forall x \in X,\ \forall y \in Y:\ x \le y,\\ \texttt{maybe\_eq}(X, Y) &\iff \exists x \in X,\ \exists y \in Y:\ x = y \iff X \cap Y \ne \emptyset,\\ \texttt{overlaps}(X, Y) &\iff X \cap Y \ne \emptyset,\\ \texttt{contains}(X, Y) &\iff Y \subseteq X. \end{aligned}$$

区間の式。 $X = [a, b]$ と $Y = [c, d]$ を空でない区間とする。

$$\texttt{definitely\_lt}(X, Y) \iff b < c .$$

($\Rightarrow$) $x = b \in X$、$y = c \in Y$ を取る。($\Leftarrow$) 任意の $x \in X$、$y \in Y$ について $x \le b < c \le y$。$\le$ で同じ議論をすると $\texttt{definitely\_le}(X, Y) \iff b \le c$ を得る。可能的な関係は存在量化によるものである：

$$\exists x \in X,\ \exists y \in Y:\ x < y \iff a < d .$$

($\Rightarrow$) $a \le x < y \le d$。($\Leftarrow$) $x = a$、$y = d$ を取る。これは引数を入れ替えた確定的関係の否定なので、専用のトレイトは要らない：

$$\neg\,\texttt{definitely\_le}(Y, X) \iff \neg\,\forall y, x:\ y \le x \iff \exists x, y:\ x < y \iff a < d .$$

最後に $X \cap Y \ne \emptyset \iff a \le d \wedge c \le b$ である。両方が成り立てば $\max(a, c) \le \min(b, d)$ が共通点であり、逆に共通点 $z$ があれば $a \le z \le d$ かつ $c \le z \le b$ である。ボールについては端点を代入すると $\texttt{definitely\_lt}(B(m_1, r_1), B(m_2, r_2)) \iff m_1 + r_1 < m_2 - r_2$ を得る。

三値の真理値。 包含区間だけからは、「$x < y$」は三つの真理値のいずれかをとる：

$$[\![\, x < y \,]\!] = \begin{cases} \mathsf{T} & \texttt{definitely\_lt}(X, Y),\\ \mathsf{F} & \texttt{definitely\_le}(Y, X),\\ \mathsf{U} & \text{otherwise}. \end{cases}$$

この値は well-defined である。$\mathsf{T}$ と $\mathsf{F}$ が同時に成り立つには $b < c$ かつ $d \le a$ が必要で、$a \le b < c \le d \le a$ となり矛盾する。同様に $[\![\, x = y \,]\!]$ は $\neg\,\texttt{maybe\_eq}(X, Y)$ のとき $\mathsf{F}$、それ以外は $\mathsf{U}$ である（両方の包含区間が同じ一点である場合を除く）。複合条件は Kleene の強三値論理で組み合わせる：

$p$	$q$	$\neg p$	$p \wedge q$	$p \vee q$
$\mathsf{T}$	$\mathsf{T}$	$\mathsf{F}$	$\mathsf{T}$	$\mathsf{T}$
$\mathsf{T}$	$\mathsf{U}$	$\mathsf{F}$	$\mathsf{U}$	$\mathsf{T}$
$\mathsf{T}$	$\mathsf{F}$	$\mathsf{F}$	$\mathsf{F}$	$\mathsf{T}$
$\mathsf{U}$	$\mathsf{T}$	$\mathsf{U}$	$\mathsf{U}$	$\mathsf{T}$
$\mathsf{U}$	$\mathsf{U}$	$\mathsf{U}$	$\mathsf{U}$	$\mathsf{U}$
$\mathsf{U}$	$\mathsf{F}$	$\mathsf{U}$	$\mathsf{F}$	$\mathsf{U}$
$\mathsf{F}$	$\mathsf{T}$	$\mathsf{T}$	$\mathsf{F}$	$\mathsf{T}$
$\mathsf{F}$	$\mathsf{U}$	$\mathsf{T}$	$\mathsf{F}$	$\mathsf{U}$
$\mathsf{F}$	$\mathsf{F}$	$\mathsf{T}$	$\mathsf{F}$	$\mathsf{F}$

$\mathsf{U}$ を「許容される値のうちあるものでは真、あるものでは偽」と読むと、各欄はそれらすべてについて成り立つ最も強い主張である。例えば $\mathsf{F} \wedge \mathsf{U} = \mathsf{F}$ となるのは、偽の連言肢を含む連言はもう一方が何であれ偽だからである。¹1 S. C. Kleene, Introduction to Metamathematics, 1952, §64。三通りの結果を持つ区間比較は R. E. Moore, Interval Analysis, 1966 にさかのぼる。

証明の段階

証明付きバックエンドは、超越関数 $f$ を目標精度 $p_t$ で正しく丸めた $\circ_{p_t}(f(x))$ を、CertificationStage の値を段階とするパイプラインで計算する。例として $f = \exp$ の場合：

1. RangeReduction：$k \in \mathbb{Z}$、$|r| \le \tfrac12 \ln 2$ として $x = k \ln 2 + r$ と書き、$\exp(x) = 2^{k}\exp(r)$ とする。還元後の引数 $r$ 自体も包含区間で囲む必要があり、そのために $\ln 2$ がおよそ $\log_2 |x|$ ビット余分に必要になる。

2. SeriesEvaluation：$\sum_{j<N} r^{j}/j!$ を計算し、残りの項を上から抑える。$|r| < N + 1$ のとき、

   $$\Bigl|\sum_{j \ge N} \frac{r^{j}}{j!}\Bigr| \le \frac{|r|^{N}}{N!}\sum_{i \ge 0}\Bigl(\frac{|r|}{N+1}\Bigr)^{i} = \frac{|r|^{N}}{N!}\cdot\frac{1}{1 - |r|/(N+1)},$$

   ここで $\frac{N!}{(N+i)!} \le (N+1)^{-i}$ を用いた。

3. EnclosurePropagation：打ち切り誤差の上界と、作業精度 $p_w > p_t$ でのすべての丸め誤差を残りの演算を通して伝播し、包含区間 $[\ell, h] \ni f(x)$ を得る。

4. TargetRounding：$\circ_{p_t}(\ell) = \circ_{p_t}(h)$ ならその値が答えである。単調性により

   $$\ell \le f(x) \le h \;\Longrightarrow\; \circ_{p_t}(\ell) \le \circ_{p_t}(f(x)) \le \circ_{p_t}(h) = \circ_{p_t}(\ell).$$

   そうでなければ $[\ell, h]$ は丸めの境界をまたいでいる。バックエンドは $p_w$ を上げて繰り返す。これが Ziv の戦略である。²2 A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991。$f(x)$ が境界にどれだけ近づきうるかは表作成者のジレンマ（table maker’s dilemma）と呼ばれる。ほとんどの関数では $p_w$ についての有用な事前の上界が知られておらず、そのためループには予算が必要である。

ループが諦めると、バックエンドは段階、理由、$p_t$（target_precision）、最後の $p_w$（work_precision）、精度を上げた回数（refinements）を持つ ArithmeticError::certification_failure を返す。TargetRounding で想定される理由は RefinementBudgetExhausted であり、その他の理由はより前の段階に属する。このパッケージは語彙だけを定義し、何も評価しない。

設計判断

一つのシグネチャではなく三つの層

問題。 内側のループで Double の平方根を使うなら fn sqrt(Double) -> Double が望ましく、負の引数に NaN を返すことも許容できる。十進バックエンドには精度と丸めモードが必要で、結果が丸められたかどうかを報告しなければならない。一つのシグネチャで両方に応えようとすると、ネイティブの呼び出しすべてに Result とコンテキストを強いるか、十進の呼び出し側が必要とする情報を捨てるかのどちらかになる。

選択肢。 (a) 非検査トレイトのみ。(b) コンテキスト付きトレイトのみ。(c) 独立した三つの層。

選択。 (c)。層は情報を段階的に多く運び、各結果型は次の結果型へ埋め込まれる：

$$\underbrace{T}_{\text{unchecked}} \;\xrightarrow{\;v \,\mapsto\, \mathrm{Ok}(v)\;}\; \underbrace{\mathrm{Result}[T, E]}_{\text{checked}} \;\xrightarrow{\;\mathrm{Ok}(v) \,\mapsto\, \mathrm{Ok}(v, \mathbf{0})\;}\; \underbrace{\mathrm{Result}[(T, D), E]}_{\text{contextual}},$$

ここで $D$ は診断の集合、$\mathbf{0}$ はその空の値である。組み込みのコンテキスト付きアダプタは、Float の整数埋め込みを除き、まさに検査付きまたは非検査の結果にこれらの埋め込みを適用したものである。境界（bound）はアルゴリズムが何を扱うかを述べる。T : Sqrt はバックエンド自身の振る舞いを受け入れ、T : SqrtChecked は拒否を処理し、T : SqrtContextual はコンテキストと診断を必要とする。層は上位トレイトでつながっていないので、型は守れる層だけを実装できる。区間型は、コンテキストなしの Sqrt を持つふりをせずに DivChecked と包含区間の関係を実装できる。

コンテキストと診断は明示的な値である

問題。 IEEE 754 は丸め方向とステータスフラグを実行環境の属性として記述する。C はそれらを <fenv.h> で公開し、Python の decimal はスレッドローカルな現在のコンテキストを持つ。どちらも隠れた状態である。a + b の結果が引数でないものに依存し、ある計算で立ったフラグが次の計算でも立ったままになる。

選択肢。 (a) グローバルな可変コンテキスト。(b) スレッドまたはタスクにローカルなコンテキスト。(c) コンテキストを引数にし、フラグを戻り値に入れる。

選択。 (c)。すべてのコンテキスト付き演算は次の関数である

$$\operatorname{op} : T^{n} \times \texttt{ArithmeticContext} \to \mathrm{Result}[(T, D), E],$$

したがって等しい入力は等しい出力を与え、計算の診断はそれが組み合わせた演算の診断そのものである。これはどの MoonBit ターゲットでも同じように動き（どれもスレッドローカル記憶域を必要としない）、並行利用を安全にし、テストが演算の入力全体を記述できるようにする。代償は冗長さであり、バックエンドは独自のヘルパーでそれを減らせる。

エラー、診断、証明の失敗は別物である

問題。 IEEE 754 には五つの例外がある。無効演算、ゼロ除算、オーバーフロー、アンダーフロー、不正確である。存在しない結果を表すものもあれば、存在するが丸められた結果を表すものもある。

選択。 パッケージは値が返されるかどうかでそれらを分ける：

状況	経路	IEEE 754 での対応
意味のある値がない（定義域外、不定形）	Err, DomainError	無効演算
極：ゼロでない有限値をゼロで割る	Err, DivisionByZero	ゼロ除算
正確な値と異なる $\overline{\mathbb{F}}$ の値	inexact、rounded 付きの Ok	不正確
有限範囲を超える、または正規範囲を下回る値	overflow、underflow、subnormal 付きの Ok	オーバーフロー、アンダーフロー
結果を証明できなかった有効な入力	Err, CertificationFailure	なし

フラグがエラーを隠してはならず、フラグを運ぶためにエラーをでっち上げてもならない。$\operatorname{fl}(10^{300} \times 10^{300}) = +\infty$ は IEEE として正しい答えなので、失敗ではなく overflow 付きの値である。証明の失敗はエラーだが定義域エラーではない。入力は有効で、予算を増やせば成功するかもしれないので、再試行するかどうかを決めるのに必要なデータを運ぶ。パッケージは再試行の方針を課さない。

包含区間の関係は順序ではない

問題。 区間は順序付けられているように見え、Compare を実装すれば汎用のソートコードに渡せるようになる。

選択。 五つの別々の関係トレイト。Compare は全順序を約束するが、空でない区間上の definitely_lt は狭義半順序にすぎない。非反射的（$[a, b]$ について $b < a$ は成り立たない）かつ推移的である：

$$b_1 < c_2 \;\wedge\; b_2 < c_3 \;\Longrightarrow\; b_1 < c_2 \le b_2 < c_3 \;\Longrightarrow\; b_1 < c_3 ,$$

しかし全順序ではない。重なり合う $X, Y$ では $\texttt{definitely\_lt}(X, Y)$ も $\texttt{definitely\_lt}(Y, X)$ も成り立たず、両者が等しいわけでもない。Compare のインスタンスはそこで $<, =, >$ のどれかを答えなければならず、未知の値について偽の主張をすることになる。

ネイティブスカラーは守れるものだけを実装する

問題。 Float と Double は、コンテキストを無視すればすべてのトレイトを実装できてしまう。

選択。 結果が意味を持つ場合にだけ能力を実装する：

• すべての非検査トレイトと Power。これらはバックエンド以上のことを約束しない。
• 検査付きトレイト。追加の約束は不正な引数を拒否することだけである。
• コンテキスト付きの四則演算、絶対値、平方根、指数関数、整数埋め込み、隣接値、形式の問い合わせ。汎用のコンテキスト付きコードがネイティブスカラーでも動くようにするためのアダプタである。
• ConstantsContextual と HyperbolicContextual は実装しない。これらの目的は任意の精度に従い意味のある診断を伴う結果であり、固定精度のライブラリ関数にはそれを提供できない。

アダプタは、固定形式によって正確になる場所（隣接値、Double の整数埋め込み）では誠実であり、損失の検出が安価な場所（後述の Float の整数埋め込み）では損失を検出する。算術のアダプタは丸めを検出しない。その空の診断は「検出していない」を意味する。API ページ は使用箇所でこのことを述べている。

トレイトごとに一つの能力、Real はない

Real : Field + Sqrt + Exponential + Trigonometric + Compare のようなトレイトは便利だが、重要な違いを隠してしまう。区間型には全順序がなく、十進型には安価な sin がなく、整数型には Power はあっても Sqrt はない。このパッケージのトレイトはどれも一つの能力であり、アルゴリズムは斜辺の計算なら T : Add + Mul + Sqrt のように、使う境界を組み合わせる。Radical は唯一の連言であり、平方根と立方根はよく一緒に必要になるからである。

三つの浮動小数点クラス

IEEE 754 の class は十のクラス（シグナリング NaN と静かな NaN、負と正の無限大、正規化数、非正規化数、ゼロ）を区別する。FpClass が三つにとどめるのは、それが汎用コードの分岐する場合だからである。有限値はさらに演算に入れられ、無限大は有効な極限であり、NaN は無効である。符号、ゼロ、非正規化かどうかは形式自身の API で調べる。

IntegralContextual は Int だけを埋め込む

どのバックエンドも MoonBit の Int を受け取れ、ループカウンタや添字は Int である。BigInt の埋め込みにはすべてのバックエンドで任意精度の丸めが必要になる。共通の場合の実装を安価に保つため、それは別の能力に任せる。

Power は一つのシグネチャを保つ

Power::pow(Self, Self) は浮動小数点型でも整数型でも同じなので、汎用のべき乗は型の族に依存しない。代償は指数の型が底の型になることで、符号付き整数と BigInt のインスタンスは負の指数で中断しなければならない。一般に $x^{-n} \notin \mathbb{Z}$ だからである。定義された失敗が必要なコードは PowNatChecked（指数は UInt）か PowIntChecked（指数は Int）を使う。

正しさと不変条件

コンテキストの不変条件

ArithmeticContext::new は $p \ge 1$ をクランプで、$e_{\min} \le e_{\max}$（両方あるとき）を中断で保証し、フィールドはパッケージ外から読み取り専用である。したがってすべてのコンテキスト値は両方を満たし、バックエンドは再検査する必要がない。

combine の法則

ArithmeticDiagnostics はブール束 $D = \{0, 1\}^{6}$ であり、combine は成分ごとの $\vee$ である。各成分がブール代数の法則を満たすので、すべての $d_1, d_2, d_3 \in D$ について：

$$\begin{aligned} (d_1 \vee d_2) \vee d_3 &= d_1 \vee (d_2 \vee d_3) && \text{associativity}\\ d_1 \vee d_2 &= d_2 \vee d_1 && \text{commutativity}\\ d \vee d &= d && \text{idempotence}\\ d \vee \mathbf{0} &= d && \mathbf{0} = \texttt{empty()} \end{aligned}$$

したがって $(D, \vee, \mathbf{0})$ は可換かつ冪等なモノイド、すなわち最小元を持つ結び半束である。計算の診断は各ステップの診断の結びであり、評価の順序やまとめ方に依存せず、一度立ったフラグは合成によって消えることがない。

コンテキスト付き演算の連結

二つのコンテキスト付き演算 $f : A \to \mathrm{Result}[(B, D), E]$ と $g : B \to \mathrm{Result}[(C, D), E]$ を合成すると次を得る

$$(g \mathbin{\bar\circ} f)(a) = \begin{cases} \mathrm{Err}(e) & f(a) = \mathrm{Err}(e),\\ \mathrm{Err}(e) & f(a) = \mathrm{Ok}(b, d_1),\ g(b) = \mathrm{Err}(e),\\ \mathrm{Ok}(c, d_1 \vee d_2) & f(a) = \mathrm{Ok}(b, d_1),\ g(b) = \mathrm{Ok}(c, d_2). \end{cases}$$

これはエラーモナドの上に積んだ $(D, \vee, \mathbf{0})$ 上の writer モナドであり、上のモノイド法則こそが $\bar\circ$ を結合的にし、ArithmeticOutcome::exact をその単位元にする。³3 $\bar\circ$ の結合性は、診断上の $\vee$ の結合性と値上の関数合成の結合性に帰着し、単位法則は $d \vee \mathbf{0} = d$ に帰着する。E. Moggi, “Notions of computation and monads”, 1991 を参照。 パッケージは合成子ではなく部品（exact、with_diagnostics、combine）を提供する。チュートリアルに短いヘルパーを示す。

Float への整数埋め込み

binary32 は $p = 24$ である。$|n| \le 2^{24}$ の整数 $n$ は有効ビットが高々 24 ビット（あるいは 2 のべき $2^{24}$ そのもの）なので表現可能である。$2^{24} + 1$ は有効ビットが 25 ビット必要なので表現できず、偶数への最近接丸めで $2^{24}$ になる。Float のインスタンスは幅の広い整数の比較なしに損失を検出する。binary64 は $p = 53 > 31$ なので Double::from_int はすべての Int で正確であり、binary32 の値はすべて binary64 の値なので拡張も正確である。したがって

$$\operatorname{double}(\operatorname{RN}_{32}(n)) = \operatorname{double}(n) \iff \operatorname{RN}_{32}(n) = n ,$$

であり、インスタンスは変換で情報が失われたときにちょうど inexact と rounded を設定する。$|n| \le 2^{31} < 2^{128}$ なのでオーバーフローは起こらない。

二進べき乗法の誤差の上界

Float と Double の PowNatChecked は次のループで $x^{n}$ を計算する

$$\textit{acc} \leftarrow 1,\ \textit{f} \leftarrow x,\ k \leftarrow n;\quad \text{while } k > 0:\ \text{if } k \text{ odd}: \textit{acc} \leftarrow \textit{acc}\cdot \textit{f};\ k \leftarrow \lfloor k/2 \rfloor;\ \text{if } k > 0: \textit{f} \leftarrow \textit{f}^{2}.$$

正しさ。 正確な算術では、各反復の開始時に $\textit{acc}\cdot \textit{f}^{\,k} = x^{n}$ が成り立つ。初期状態で成り立ち、$k = 2j + 1$ なら $\textit{acc}\,\textit{f}\cdot(\textit{f}^{2})^{j} = \textit{acc}\,\textit{f}^{\,k}$、$k = 2j$ なら $\textit{acc}\,(\textit{f}^{2})^{j} = \textit{acc}\,\textit{f}^{\,k}$ である。$k = 0$ で不変条件から $\textit{acc} = x^{n}$ を得る。ループは $\lfloor \log_2 n \rfloor + 1$ 回まわり、$\lfloor\log_2 n\rfloor$ 回の二乗と $\operatorname{popcount}(n)$ 回の積を行うが、その最初の積（$1 \cdot \textit{f}$）は正確である。整数の Power インスタンスも $\mathbb{Z}/2^{k}$ で同じ不変条件を使い、そこではすべてのステップが正確である。

丸め誤差。 $x^{m}$ を近似する各計算値 $q$ に、$|\delta_i| \le u$ かつ $\sum_i k_i \le c(q)$ として $q = x^{m}\prod_i (1 + \delta_i)^{k_i}$ となるような誤差の個数 $c(q)$ を割り当てる。すると $c(x) = 0$ であり、$q_1 \approx x^{m_1}$ と $q_2 \approx x^{m_2}$ の丸めを伴う積一回で $c \le c(q_1) + c(q_2) + 1$ となる。帰納法により $c(q) \le m - 1$：

$$c(q_1 q_2) \le (m_1 - 1) + (m_2 - 1) + 1 = (m_1 + m_2) - 1,$$

であり、二乗は $q_1 = q_2$ の場合で、共有された誤差が二度数えられる。したがってオーバーフローとアンダーフローがなければ、

$$\operatorname{fl}(x^{n}) = x^{n}(1 + \theta_{n-1}), \qquad |\theta_{n-1}| \le (1 + u)^{n-1} - 1 \le \gamma_{n-1} := \frac{(n-1)u}{1 - (n-1)u}.$$

最後の段階は、$ku < 1$ のときの標準的な補題 $|\prod_{i=1}^{k}(1+\delta_i)^{\pm 1} - 1| \le \gamma_k$ である。⁴4 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1 および §3.1。 負の指数の PowIntChecked はもう一度除算を行い、$(1 + \delta)/(1 + \theta_{n-1})$ から $|\theta| \le \gamma_{n}$ を得る。二進べき乗法は $n - 1$ 回の逐次乗算の最悪時の上界を改善しない。乗算の回数を $n - 1$ から $O(\log n)$ に減らすだけである。

検査付き除算

Float と Double の DivChecked インスタンスはすべてのゼロ除数を拒否し、その理由を種別で示す。$0/0$ と $\infty/\infty$ は不定形である。$t \to 0$ や $t \to \infty$ のとき極限 $\lim (\lambda t)/t = \lambda$ はあらゆる値をとるので意味のある商はなく、種別は DomainError である。$x \ne 0$ なら $t \to 0$ で $x/t$ は発散する（極）ので、種別は DivisionByZero である。これは IEEE 754 の無効演算とゼロ除算の例外に対応するが、より厳格である。IEEE は $\infty/0$ には $\pm\infty$ を、NaN$/0$ には NaN を黙って返すが、検査付きインスタンスは DivisionByZero を返す。NaN の被除数とゼロでない除数は Ok(NaN) として伝播し、SqrtChecked も同様に NaN を通す。検査付き演算は不正な引数を拒否するのであって、以前の不正な結果を報告し直すわけではない。

包含区間の関係の健全性と単調性

健全性。 $x \in X$、$y \in Y$ かつ $\texttt{definitely\_lt}(X, Y)$ なら $x < y$ である。この関係は $(x, y)$ を含む $X \times Y$ 全体についての全称命題だからである。双対的に、$\neg\,\texttt{maybe\_eq}(X, Y)$ なら $x \ne y$ である。

精緻化に対する単調性。 $X' \subseteq X$ かつ $Y' \subseteq Y$ なら

$$\texttt{definitely\_lt}(X, Y) \Rightarrow \texttt{definitely\_lt}(X', Y'), \qquad \texttt{maybe\_eq}(X', Y') \Rightarrow \texttt{maybe\_eq}(X, Y),$$

である。全称命題は定義域を縮めても成り立ち、存在命題は定義域を広げても成り立つからである。三値で言えば、包含区間を精緻化すると $\mathsf{U}$ が $\mathsf{T}$ や $\mathsf{F}$ に変わることはあるが、$\mathsf{T}$ が $\mathsf{F}$ に変わることはない。これが、上の目標丸め段階のような「決まるまで精緻化する」ループを正しくする。一度下した判断は有効なままである。そのようなループは Contains、すなわち $\texttt{contains}(X, X')$ で、精緻化した包含区間 $X'$ が古いものの内側にあることを確かめる。

トレイトは空の包含区間の規約を定めず、各バックエンドが自身の規約を文書化する。量化の読み方では確定的関係は空の引数について空虚に成り立つので、$\mathsf{T}$ が情報の欠如から生じないようにしたいバックエンドは、代わりにそれらに false を返す。

退けた代替案

• スティッキーなフラグを持つグローバルまたはスレッドローカルなコンテキスト（C の <fenv.h> や Python の decimal のようなもの）。明示的な値の節で述べた理由で退けた。結果が隠れた状態に依存し、計算の間でフラグが漏れる。
• Real や Number のような上位トレイト。 正確な型、近似的な型、包含区間値の型の違いを隠すので退けた。
• 包含区間に対する Compare。 確定的な順序は全順序ではないので退けた。
• 比較のための三値の結果型（True | False | Unknown）。上で示したとおり、二つのブール値の射影 definitely_* と maybe_eq から再構成でき、通常の if と組み合わせられる。別の型にすると、一方向だけを問う場合でもすべての呼び出し側に $\mathsf{U}$ の処理を強いることになる。
• definitely_eq 関係。 退化していない包含区間では常に偽なので、等しい一点どうしの判定にしかならない。
• エラーとして報告される診断。 不正確な、あるいはオーバーフローした IEEE の結果は正しい答えであり、それを Err にすると丸めを伴うすべての演算が失敗することになるので退けた。
• 検査付き結果のための Option や raise。 Option は理由を失う。MoonBit の raise は失敗を戻り値の型の外に置く。Luna-Flow はリポジトリ全体で構造化エラーを伴う Result を使う。
• コンテキストを無視してネイティブスカラーに ConstantsContextual と HyperbolicContextual を実装すること。 これらのトレイトはコンテキストに忠実な結果を約束するためにあるので退けた。

範囲外のこと

• このパッケージは任意精度、十進、区間、ボールの算術を実装せず、証明付きの関数評価も行わない。それらのバックエンドが実装するトレイトを定義する。
• 組み込みの Float と Double のインスタンスは、コンテキストの精度、丸めモード、指数範囲に従わず、その算術アダプタは丸め、オーバーフロー、アンダーフローを検出しない。
• Float と Double について正しく丸められた初等関数は約束しない。それらは Kaida-Amethyst/math から来ている。
• 代数構造（Ring、Field など）は定義しない。それは luna-generic の役割である。ベクトル、行列、複素数、多項式も定義しない。
• 提供される各インスタンスが受け継ぐもの以上に、非検査トレイトの分岐切断や特殊値の規約を選ぶことはしない。
• 証明の失敗に対する再試行や精度引き上げの方針は定めない。

Footnotes

1. S. C. Kleene, Introduction to Metamathematics, 1952, §64。三通りの結果を持つ区間比較は R. E. Moore, Interval Analysis, 1966 にさかのぼる。 ↩

2. A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991。$f(x)$ が境界にどれだけ近づきうるかは表作成者のジレンマ（table maker’s dilemma）と呼ばれる。ほとんどの関数では $p_w$ についての有用な事前の上界が知られておらず、そのためループには予算が必要である。 ↩

3. $\bar\circ$ の結合性は、診断上の $\vee$ の結合性と値上の関数合成の結合性に帰着し、単位法則は $d \vee \mathbf{0} = d$ に帰着する。E. Moggi, “Notions of computation and monads”, 1991 を参照。 ↩

4. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1 および §3.1。 ↩