bin_float の設計

bin_float は任意精度の IEEE 754 2 進浮動小数点算術を実装します。このページではその背後にある数学を説明します。値の集合、丸め関数とそれが満たす誤差モデル、各演算が厳密な整数データから正しく丸められた結果をどのように決定するか、指数範囲・極小性・ステータスフラグがどのように扱われるか、IEEE の剰余がなぜ厳密なのか、10 進変換と初等関数がどのように認証されるか、そして高速な整数カーネルがなぜ結果を変ええないのか、です。呼び出し可能な機能は API リファレンスに列挙されており、チュートリアルではその使い方を示しています。

設計目標

bin_float のすべての演算は、1 から $2^{28}$ ビットまでのすべての精度 $p$ と、$\pm(2^{30}-1)$ までのすべての指数範囲について、IEEE 754-2019 が正しく丸められた演算に要求する値、すなわち厳密な実数 $f(x)$ に対する $\circ(f(x))$ を、IEEE のステータスフラグとともに厳密に返します。同じコードが 2 種類の利用者に役立ちます。Double をはるかに超える 2 進有理数の値を必要とする任意精度数値計算と、非正規化数、5 つの丸め方向、両方の極小性規則を含む binary16、binary32、binary64、binary128 のビット単位で厳密なエミュレーションです。隠れた状態はありません。精度、丸め、範囲は不変な BinaryContext として渡され、フラグは値として返されます。

数学的背景

2 進有理数の値と保存される三つ組

有限の BinFloat は次の 2 進有理数を表します。

$$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 進有理数はちょうど 1 つのこのような形を持つ（$2^{\nu_2(c)}$ をくくり出す）ので、2 つの有限値が数値的に等しいのは、符号（非ゼロの値の場合）、係数、指数がすべて一致するときに限ります。先頭ビットの指数は次のとおりです。

$$\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 の 2 進形式

$k$ ビットの 2 進交換形式は、符号ビット、$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 の 6 つの丸め方向は、次の写像 $\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 が 10 進のコアと共有している GDA の「round-up」モードです。上記のすべての $\circ$ が持つ 2 つの性質が、このページの証明の大部分を支えます。

$$\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$ の 2 つの隣接点のうちの 1 つであり、その選択規則は $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, 2nd ed., 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 et al., Handbook of Floating-Point Arithmetic, 2nd ed., Birkhäuser 2018, §2.1 および §4.3。

正しい丸めはこのモデルより強い性質です。結果は $u$ 以内の何らかの点ではなく、ただ 1 つの点 $\circ(f(x))$ です。bin_float 全体がその点を返すように作られているので、上記のモデルは初等関数を含むすべての演算で成り立ちます。

設計上の判断

厳密なデータから一度だけ丸める

問題。 結果は厳密な実数 $r$ に対する $\circ(r)$ に等しくなければなりませんが、$r$ は $p$ よりはるかに多いビット（$p$ ビットの数 2 つの積は $2p$ ビット）や、無限のビット（商、平方根、$e^x$）を必要とすることがあります。

選択肢。 FMA を持たないハードウェアのように、より広い形式で計算してから丸め直す方法、固定数のガードビットを保持する方法、あるいは厳密な情報から丸めを決定する方法があります。

選択。 すべての演算は $r$ の厳密な記述を計算し、1 つの最終処理関数を呼び出します。2 進有理数の結果（和、差、積、fma、scaleb、remainder、変換）では、記述は $r = \pm m 2^{e}$ を満たす厳密な整数の絶対値 $m$ と指数 $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$ と分解し、3 つのデータを得ます。

$$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$ ビットに丸めてから非正規化数の格子へ丸めると二重丸めになります。粗い格子の中点のすぐ上にある値が、最初の丸めでその中点に押し出され、2 回目の丸めで誤った向きに丸められることがあるからです。$q$ からの桁上がり（$q + 1 = 2^p$ のとき）は $\operatorname{top}$ を 1 だけ上げるだけです。結果は再正規化され、丸められた値に対して後述のオーバーフロー判定が適用されます。

商と剰余による除算

問題。 $a = c_a 2^{e_a}$、$b = c_b 2^{e_b}$ として、$a/b$ は有理数 $N/D \cdot 2^{e}$（$N = c_a$、$D = c_b$、$e = e_a - e_b$）であり、その 2 進展開は通常無限です。

選択。 まず厳密な先頭指数を求めます。$k = \operatorname{bits}(N) - \operatorname{bits}(D)$ とすると、$\lfloor \log_2 (N/D) \rfloor$ は $N \ge D 2^{k}$ ならば $k$、そうでなければ $k - 1$ であり、整数比較 1 回で決まります。これにより目標の指数 $\tau = \max(\operatorname{top} - p + 1,\ e_{\min} - p + 1)$ が定まり、続いて 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$ のときはシフトを分母に移し、1 単位未満の商は $D 2^{\tau - e}$ を作ることなく $2N$ と比較して判定します。

整数平方根と中点判定による平方根

問題。 $\sqrt{c 2^{e}}$ は、$c 2^{e}$ が平方数でない限り無理数です。

選択。 先頭指数は $\lfloor \operatorname{top}(x)/2 \rfloor$（床除算）であり、これにより上と同様に $\tau$ が定まります。被開平数を $X = c\, 2^{e - 2\tau}$ と書くと、$\sqrt{x} = \sqrt{X}\, 2^{\tau}$ です。厳密な整数平方根は剰余 $X - s^2$ とともに $s = \lfloor \sqrt{X} \rfloor$ を与え、根が厳密であるのは剰余がゼロのときに限ります。そうでなければ、丸めのデータは中点 $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 et al., Handbook of Floating-Point Arithmetic, §5.3 および §7.6。 それでも等号の分岐は残してあります。コンテキストの精度より多くのビットを持つオペランドでは、$\sqrt{X}$ がちょうど中点になりうるからです（たとえば 1 ビットに丸める $\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$ の外にあるからです。オーバーフローした結果は、最初の 2 つのグループでは $\pm\infty$、3 つ目のグループでは $\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")
}

極小性（tininess）。 非ゼロの結果が $\pm 2^{e_{\min}}$ の間に真に含まれるとき、その結果は極小です。IEEE 754-2019（第 7.5 節）は 2 通りの解釈を認めており、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 回目の分解を行います。2 つの解釈が異なるのは、$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 つの非正規化数のちょうど中間にあり、偶数側の隣接点である最小の正規化数 $2^{-14}$ に丸められ、エンコードされた結果 0x0400 は正規化数であるにもかかわらず、underflow と inexact を発生させます。

///|
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$、それ以外のすべての方向では $+0$ です。また $(-0) + (-0) = -0$ です（第 6.3 節）。積と商は符号の排他的論理和をとります。NaN のオペランドは、最初の NaN オペランドを符号とペイロードを保ったまま quiet にしたものを返し（第 6.2.3 節は任意の入力 NaN を許容しています）、invalid_operation は、オペランドが signaling NaN であるか、演算そのものが不正である場合（$\infty - \infty$、$0 \cdot \infty$、$0/0$、$\infty/\infty$、$\sqrt{x<0}$、$\operatorname{remainder}(\infty, y)$、$\operatorname{remainder}(x, 0)$）に限り発生します。フラグは値です。combine はビット単位の OR なので、計算のフラグは可換で冪等なモノイドをなし、任意の順序で蓄積できます。これは大域的なスティッキーレジスタでは並行コードに提供できない性質です。

加算における大きく離れたオペランド

問題。 $2^{10^9} + 2^{-10^9}$ は 2 進有理数として厳密ですが、それを作るには 20 億ビットの係数が必要です。

選択。 先頭の指数の差が $p + 3$ を超える場合、小さい方のオペランドは位置 $\operatorname{exp}(\text{high}) - p - 3$ で切り捨てられ、それより下はすべて 1 つのスティッキービットに置き換えられます。切り捨てられた下位オペランドの整数部 $L$ は厳密に入り、何かが捨てられた場合、大きさは半単位で $2(H \pm L) + 1$（減算では $2(H - L - 1) + 1$）になります。

丸めにとって厳密である理由。 結果は $\operatorname{top} \ge \operatorname{top}(\text{high}) - 1$ を満たすので、丸め位置は少なくとも $\operatorname{top}(\text{high}) - p$ であり、ラウンドビットはそれより少なくとも 1 つ下にあります。一方、捨てられるビットはすべて $\operatorname{top}(\text{high}) - p - 3$ 以下にあります。したがって捨てられた部分は $q$ も $g$ も変えず、$t$ が立つかどうかだけに影響し、代わりのビットは非ゼロの何かが捨てられたときに限り $t$ を立てます。減算では、$0 < 1 - \varepsilon < 1$ として $H - (L + \varepsilon) = (H - L - 1) + (1 - \varepsilon)$ なので、借りを伴う形にも同じ置き換えが適用できます。非正規化数のシフトを含む完全な議論は、以下の添付資料にあります。

融合積和演算

fma_ctx は積 $c_x c_y 2^{e_x + e_y}$ を、精度がそれ自身のビット長に等しい値として厳密に作り、加算の最終処理関数に渡します。したがって $x y + z$ は一度だけ丸められます（第 5.4.1 節）。2 回の丸めとの違いこそがこの演算の要点です。binary64 で $a = \operatorname{RN}(0.1)$ とすると、mul_ctx に続けて sub_ctx を行うと $\operatorname{RN}(a \cdot a) - a \cdot a$ は失われます（2 番目の演算には等しい 2 つの数が見えるだけです）が、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 et al., §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$ に属する。

証明. $|M_x|, |M_y| < 2^p$、$q_x, q_y \ge e_{\min} - p + 1$ として $x = M_x 2^{q_x}$、$y = M_y 2^{q_y}$ と書く。$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$ の場合：$|k| 2^{q_y} \le |M_y| 2^{q_y}/2$ として $r = k 2^{q_y}$ なので、$|k| < 2^{p-1}$。
• $q_x < q_y$ の場合：$|k| 2^{q_x} \le |x| = |M_x| 2^{q_x}$ として $r = k 2^{q_x}$ なので、$|k| < 2^{p}$。

どちらの場合も $|k| < 2^p$ であり、指数は少なくとも $e_{\min} - p + 1$ であり、$|r| \le \max(|x|, |y|) \le \Omega$ なので、$r \in F$ である。$\square$

実装は、$2^{30}$ ビットにもなりうる $n$ を決して作りません。$m = \min(q_x, q_y)$ として $X = |x| 2^{-m}$、$Y = |y| 2^{-m}$（整数）とし、$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) は、$\pi = \min(e_{\min} - p,\ e(x),\ \operatorname{top}(x) - p)$ として正のステップ $2^{\pi - 2}$ を $x$ に加え、$+\infty$ 方向に丸めます。$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 の定義により結果は $x$ より大きい $F$ の最小の点になります。同じ議論は $p$ ビットを超える $x$ に対しても成り立ち、そのようなオペランドが受け付けられるのはこのためです。この内部加算のフラグは捨てられます。nextUp は $\Omega$ から $+\infty$ へ進む場合でも quiet（第 5.3.1 節）だからです。

scaleb_ctx(x, n) は $(c, e + n)$ に最終処理関数を適用したものです。正規化範囲では厳密であり、その外ではアンダーフローとオーバーフローを伴って正しく丸められます。logb_ctx は $\operatorname{top}(x)$ を返します。$\operatorname{top}$ は整数の係数に対して計算されるので、これは非正規化数の $x$ に対しても厳密かつ正しい値です。整数への丸めは $\sigma = -e$（2 進小数点より下のビット）として同じラウンドビットとスティッキービットを用います。to_int_ctx とその仲間は、まず丸めてから整数を対象の範囲と比較し、実装定義の番兵値を返す代わりに invalid_operation を報告します。

10 進変換

パース。 from_string_ctx は $D \cdot 10^{k}$（$D$ は末尾にゼロを持たない整数）を厳密に読み取ります。$|k|$ が $\max\bigl(400, \lfloor (3n + p)/2 \rfloor + 64\bigr)$（$n$ は桁数）以下であれば厳密に丸められます。$k \ge 0$ では 2 進有理数の最終処理関数によって $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)]$ を用い、両端が同じフラグとともに同じ値に丸められるまで作業精度を 2 倍にしていきます。上限があることでこのループは停止します。ある方向の丸めの切り替わり点は、$F$ の点（方向付きモード）またはそれらの中点（最近接モード）であり、いずれも有効ビット数が高々 $p + 1$ の 2 進有理数です。$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|}$ は 2 進有理数ですらありません。切り替わり点でない値はすべての切り替わり点から正の距離を持ち、包含の幅は $w$ が大きくなるにつれて 0 に近づくので、ある $w$ がそれを認証します。2 進対数が確実に範囲外である値は、$\log_2 10$ と安全余裕を用いて見積もられ、一切の算術なしにオーバーフローまたはアンダーフローするので、1e100000000 にはコストがかかりません。

固定桁数。 to_decimal_string_ctx(x, d) には、$E = \lfloor \log_{10}|x| \rfloor$ として $\operatorname{round}(x / 10^{E - d + 1})$ が必要です。$E$ は答えとの差が 1 以内である $\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$ の隣にある 2 つの $n$ 桁の 10 進数（切り捨てたものとゼロから遠ざかる向きに丸めたもの）を候補とし、候補をパースして $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, Theorem 15）。追加の 1 桁は二分探索の上端のための余裕です。 両方の候補が受理された場合は、近い方を、次いで偶数の方を選びます。これはその桁数での最も近い 10 進数です。binary64 では、テストしたすべての値でホストの書式化関数と同じ結果を再現します。

認証付き初等関数

問題。 $f = \exp, \ln, \sin, \ldots$ について、値 $f(x)$ は超越数であり、それを知らないまま正しく丸めなければなりません。

選択肢。 誤差限界が証明された固定の多項式近似（高速だが 1 つの精度に縛られる）、事前の誤差限界を用いて評価し、丸めが曖昧なときはより高い精度で再試行する 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 ループです。すべての初等関数は包含区間 $[L, U] \ni f(x)$ を評価し、その内部のすべての演算は、作業精度 $w$ で $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$ 回 2 乗されます。2 乗は正の数の上で単調なので包含が保たれます。負の引数は $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$ およびそれらの逆関数は、自明な例外を除くすべての非ゼロの代数的（特に 2 進有理数の）引数で超越数になりますが、切り替わり点は 2 進有理数です。例外はループの前に除外されます。$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 の定理によりこれらが唯一の 2 進有理数の結果であることが示されます。⁸8 I. Niven, Irrational Numbers, 1956, Corollary 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$ が必要であり、これは 2 進有理数ではない。 切り替わり点でない値はすべての切り替わり点から正の距離を持つので、十分大きな $w$ がそれを認証します。$w$ がどれだけ大きくなければならないかはテーブルメーカーのジレンマであり、任意の $p$ に対して有用な事前の上界は知られていません。したがって予算はリソースの上限であって、正しさの条件ではありません。予算が尽きると、try_* 形式は段階、理由、最後の $w$ とともに CertificationFailure を報告し、total 形式は invalid_operation を伴う quiet NaN を返します。どちらも認証されていない値を返すことはありません。固定された MPFR コーパスで予算が尽きることはありません。現在のブランチでは、例外の 1 つの族が除外されていません。$1/2^k$ 以外の非整数の指数を持ちながら結果が 2 進有理数になる 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$ は 2 つの近似値を掛け合わせます。$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$ です。したがって計算値は、Higham の記法で $|\theta| \le \gamma_{n-1} = (n-1)u_w/(1 - (n-1)u_w)$ を満たす $x^n (1 + \theta)$ であり、これは $w$ ビットの結果の最終桁の単位で $2^{\operatorname{bits}(n)}$ 未満です。コードは半径 $2^{\operatorname{bits}(n) + 2}$ 単位（負の冪では、逆数がさらに $n$ 個の因子を加えるので、もう 1 つ 2 の因子を加えます）を用いて区間を構築し、同じ (R2) の議論により、両端が同じように丸められるときに受理します。12 回の倍増で認証できない場合は、厳密な冪 $c^n 2^{ne}$ を作って丸めるので、結果はどの場合でも正しく丸められます。確実に範囲外となる冪は、まず認証付きの $\log_2$ の上下界から判定されます。また厳密な値が $p$ ビットに収まる冪は厳密に計算されるので、Ziv の経路が不正確な結果しか扱わないことも保証されます。

係数カーネル

問題。 精度のコストは $p$ ビットの数の整数乗算と整数除算であり、$p$ は 1 から $2^{28}$ までの範囲をとります。

選択。 BinCoeff は 128 ビットまではインラインで保存し、それより大きな値はリトルエンディアンの 32 ビットのリム（JavaScript ではホストの bigint）として保存し、短い方のオペランドのリム数 $n$ に基づいて処理を振り分けます。

積	Native	LLVM	Wasm、Wasm-GC
筆算法（未満）	96	96	96
Karatsuba（以上）	96	96	96
Toom-3（以上）	2048	2048	4096
2 素数 NTT 乗算（以上）	2048	2048	4096
NTT 2 乗（以上）	768	768	3072
再帰的 2 乗（以上）	512	768	768

疎なオペランド（非ゼロのリムが少ないもの）には疎な積を用い、$m > 2n$ リムのオペランドは $n$ リムのブロックに分割します。除算は、1 リムのループ、除数が 48 リム未満では Knuth のアルゴリズム D、48 以上では Burnikel–Ziegler の再帰、1024 以上では Newton 法による逆数を用います。GCD は 4 リムを超えると 2 進（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 だけです。NTT は各オペランドを 16 ビットの桁に分割するので、桁の畳み込みの各係数は高々

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

です（変換長が $2^{23}$ までの場合）。NTT は素数 $p_1 = 998244353 = 119 \cdot 2^{23} + 1$ と $p_2 = 754974721 = 45 \cdot 2^{24} + 1$ を法として畳み込みを計算します。どちらも $2^{23}$ 乗根の 1 を持ちます。そして中国剰余定理によって再結合しますが、これは $p_1 p_2 \approx 2^{59.4} > 2^{55}$ なので $[0, p_1 p_2)$ において一意です。したがって再結合された係数は厳密な整数です。長さの検査はすべての変換に先立って行われ、より長い積は許容される長さの重なり合うブロックを用いるか、Toom-3 にフォールバックします。除算の経路は、構成上 $n = qd + r$ かつ $0 \le r < d$ を満たす $(q, r)$ を返します。Newton の経路は近似的な商を剰余で補正し、2 回を超える補正が必要になる場合は中断します。そのような事態は数値的な現象ではなくバグを示すものだからです。すべての経路が同じ整数を計算するので、アルゴリズムの選択によって丸められた結果、フラグ、エンコーディングが変わることはありません。

compare における NaN の順序付け

問題。 MoonBit の Compare トレイトは、ソートや順序付きマップが依拠できる 3 通りの比較を要求します。IEEE の比較は半順序であり、NaN は自分自身を含むあらゆるものと順序付けられません。

選択肢。 (1) 以前のバージョンのように NaN で中断する。この場合、NaN を含みうるデータのソートはすべてクラッシュになります。(2) IEEE の totalOrder を用いる。これは全順序ですが $-0 < +0$ を区別し、負の NaN を $-\infty$ より下に置くので、compare がゼロについて数値的な等価性と食い違ってしまいます。(3) 数については数値的な順序を保ち、すべての NaN をその上の 1 つのクラスに置く。

選択。 選択肢 (3) です。数に対してキー $\kappa(x) = (0, x)$、$\kappa(\mathrm{NaN}) = (1, 0)$ を定義し、辞書式に順序付けます。compare(x, y) は $\kappa(x)$ と $\kappa(y)$ の比較であり、$-0$ と $+0$ は同じ数に対応付けられます。全順序集合におけるキーの比較は反射的、推移的、全域的なので、compare は全前順序です。反対称ではありません（$-0$ と $+0$、あるいはペイロードの異なる 2 つの NaN は等しいと比較されますが、異なる値です）が、Compare はそれを要求しません。代償として < のもとで nan > 1 が真になるので、IEEE の意味論を必要とするコードは compare_checked（NaN でエラー）、compare_quiet / compare_signaling（4 値、フラグ付き）、または 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 でないオペランドから quiet NaN が生成されたか、signaling NaN が消費されたときに限り立ちます。combine は結合的、可換、冪等です。
• 誤差モデル。 したがって正規化範囲では、最近接丸めで $|\circ(r) - r| \le u|r|$、方向付き丸めで $< 2u|r|$ が成り立ち、$2^{e_{\min}}$ 未満では絶対誤差項 $\eta/2$（方向付き丸めでは $\eta$）となり、非正規化数の和と差は厳密です。
• 単調性。 各演算は $\circ$ が単調な $\circ \circ f$ なので、実関数が単調である引数について単調です。特に 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 経由で扱うと二重丸めが起き、一部のターゲットでは signaling NaN が失われ、binary128 はまったく表現できません。代わりに交換形式のエンコーディングは BinCoeff のビットパターン上で行います。
• 固定数のガードビット。 2 つの $p$ ビットのオペランドの加算には 3 つのガードビットで十分ですが、除算、平方根、10 進からの変換、コンテキストより広いオペランドには不十分です。厳密な整数データ（ラウンドビット、スティッキービット、剰余の符号、中点との比較）から判定する方法なら、1 つの最終処理関数でそのすべてに対応できます。
• $p$ ビットに丸めてから非正規化数の格子に丸めること。 この二重丸めは誤った非正規化数の結果を生みます。最終処理関数は、2 つの位置のうち粗い方へ一度だけシフトします。
• 見積もった誤差限界による Ziv。 関数ごと、還元ごとに個別の誤差解析が必要になり、その誤りは黙って誤った最終ビットを返すことにつながります。外向きに丸めた包含区間を用いれば、すべての演算を 2 回評価するコストと引き換えに、限界は計算される量になります。
• 大域的なフラグと丸め状態。 IEEE 754 はフラグをスティッキーな大域状態として記述しています。返される BinaryFlags の値は、純粋なコード、並行コード、@lf_arith の Result スタイルと組み合わせられ、必要な場合は combine によってスティッキーな振る舞いを再現できます。
• NaN で中断する compare、あるいは totalOrder としての compare。 compare における NaN の順序付けを参照してください。

境界

bin_float は意図的に以下を行いません。

• 区間演算やボール演算を提供すること。包含区間は認証ループの内部でのみ用いられます。ball_float は BinFloat の上に中点・半径演算を構築しています。
• IEEE 754 の代替例外処理（トラップ、置換）やスティッキーな大域フラグを実装すること。フラグは返り値です。
• 10 進浮動小数点（decimal、decimal_gda）や、それに対する 2 進以外の IEEE 演算を実装すること。
• 新たに生成された NaN に特定のペイロードを約束すること（ペイロード 0 を用います）、あるいは複数の入力のペイロードを伝播すること。
• すべての入力について、初等関数が一定時間内に完了することを保証すること。認証には予算があり、予算が尽きた場合は隠されずに報告されます。
• リムの配置、しきい値、変換のパラメータを公開すること。すべての結果、フラグ、エンコーディングが変わらない限り、これらは予告なく変更される可能性があります。
• 適合性に記録された有限のコーパスを超える適合性を主張すること。

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, 2nd ed., 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 et al., Handbook of Floating-Point Arithmetic, 2nd ed., Birkhäuser 2018, §2.1 および §4.3。 ↩

4. Muller et al., 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 et al., §4.4。 ↩

6. 往復変換には桁数 $\lceil p \log_{10} 2 \rceil + 1$ で十分です（Matula 1968; Goldberg 1991, Theorem 15）。追加の 1 桁は二分探索の上端のための余裕です。 ↩

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, Corollary 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$ が必要であり、これは 2 進有理数ではない。 ↩