decimal_gda の設計

このページでは、decimal_gda が実装する数学と、パッケージがこのように作られている理由を説明します。使い方はチュートリアルに、すべての名前の一覧は API ページにあります。

設計目標

decimal_gda は、M. F. Cowlishaw による General Decimal Arithmetic Specification（バージョン 1.70）¹1 M. F. Cowlishaw, General Decimal Arithmetic Specification, version 1.70 (2009), https://speleotrove.com/decimal/decarith.html。テストスイートは同じ著者による dectest コレクションのバージョン 2.62 です。を純粋な MoonBit の値として実装します。GDA の演算は、数値以上のものを生成するよう規定されています。すなわち結果、条件の集合、コンテキストのスティッキーなステータスの更新、そして場合によっては、定義済みの結果を利用可能なまま制御を移すトラップです。目標は、そのすべてを正確かつ観測可能な形でモデル化することであり、その結果として

各演算は、（オペランド、コンテキスト）から（結果、立てられた条件、次のコンテキスト、トラップの判定）への全域関数であり、隠れた状態を持ちません。
結果は仕様のものとビット単位で一致し、固定された公式テストスイートに対して検査されています（適合性を参照）。
パッケージは IEEE 754 のパッケージ decimal から独立しているので、一方の契約の変更がもう一方に漏れ出すことはありません。

数学的背景

数、コホート、調整指数

有限の GDA 数は、符号 $s \in \{0, 1\}$ 、係数 $c \in \mathbb{N}$ 、指数 $e \in \mathbb{Z}$ からなる三つ組 $(s, c, e)$ であり、次の値を表します。

v(s, c, e) = (-1)^s \cdot c \cdot 10^{e}.

写像 $v$ は単射ではありません： $(0, 250, -2)$ と $(0, 25, -1)$ はどちらも $2.5$ を表します。同じ値を持つ三つ組はコホートをなします。GDA は値だけでなく三つ組を保持します。指数が意味（量子 $10^e$ ：「2.50」はセント単位まで計測された）を持つからです。ゼロは符号を持つので、 $(1, 0, e)$ は $-0$ です。 $c > 0$ に対し、 $c$ の 10 進桁数を $d(c)$ と書き、

\hat e = e + d(c) - 1

を調整指数、すなわち先頭の桁の指数とします： $10^{\hat e} \le |v| < 10^{\hat e + 1}$ 。有限の数のほかに、 $\pm\infty$ と、符号とペイロード $c$ を持つ quiet NaN および signaling NaN があります。

コンテキストと表現可能な集合

コンテキストは、精度 $p \ge 1$ 、丸めモード、指数範囲 $e_{\min} \le e_{\max}$ 、クランプビットを定めます。有限の数がそのコンテキストで表現可能であるとは、

d(c) \le p, \qquad \hat e \le e_{\max}, \qquad e \ge E_{\mathrm{tiny}} \;\text{ where }\; E_{\mathrm{tiny}} = e_{\min} - p + 1,

かつ、クランプが有効なら $e \le E_{\mathrm{top}} = e_{\max} - p + 1$ でもあることです。値 $E_{\mathrm{tiny}}$ は最小単位の指数です。 $\hat e = e_{\min}$ で $p$ 桁の数は $e = e_{\min} - p + 1$ を持ち、 $10^{e_{\min}}$ 未満の数（非正規化数の範囲）はその単位を保ったまま先頭の桁を失っていきます。 $E_{\mathrm{top}}$ は $\hat e = e_{\max}$ で $p$ 桁の最大の数の指数です。クランプによって指数の集合はちょうど IEEE の交換形式が符号化できるものになります。decimal32/64/128 のプリセットがクランプするのはこのためです。

算術の規則：厳密な結果を求め、1 回だけ丸める

GDA のすべての算術演算は、同じ 2 段階の規則で定義されています。

厳密な数学的結果 $x$ を計算し、 $x$ を表す三つ組のうち、指数が演算の理想指数 $e^\ast$ に最も近いものを選ぶ。
その三つ組が表現可能でなければ、それをコンテキストに丸め、何が起きたかを表す条件を立てる。

ステップ 1 は厳密なので、ステップ 2 は 1 回の丸めであり、以下のすべての限界は 1 回の丸めについての限界です。パッケージはこの規則を文字どおりに守っています。各演算はオペランドから厳密な係数と指数を構築し、その後、共通の最終化ルーチン（normalize_decimal_parts_ctx_repr に続いて finalize_finite_ctx）を 1 回呼び出します。これが、有限の結果を丸めたり、精度や範囲に関する条件を立てたりする唯一のコードです。

理想指数

理想指数とは、オペランドにある情報だけを使って厳密な結果を表す、最も粗い量子です。導出は簡単です。

\begin{aligned} x_1 \pm x_2 &= \bigl(c_1 10^{e_1 - m} \pm c_2 10^{e_2 - m}\bigr) 10^{m}, & m &= \min(e_1, e_2),\\ x_1 \times x_2 &= (c_1 c_2)\, 10^{e_1 + e_2}, & &\\ x_1 \div x_2 &= \frac{c_1}{c_2}\, 10^{e_1 - e_2}, & &\\ \sqrt{x} &= \sqrt{c\,10^{e - 2\lfloor e/2 \rfloor}}\;10^{\lfloor e/2 \rfloor}. & & \end{aligned}

加算では、 $m$ は両方のオペランドが整数となる最大の指数なので、和が量子の整数倍であることが保証される最大の指数です。このとき括弧内は厳密な整数係数になります。乗算では、係数 $c_1 c_2$ は指数 $e_1 + e_2$ で整数であり、一般にそれより大きな指数では整数になりません。除算では商 $c_1 / c_2$ が整数とは限らないので、理想指数 $e_1 - e_2$ は到達可能な場合にだけ使われます。厳密な商は、高々 $p$ 桁の整数係数を保つ範囲で $e_1 - e_2$ に最も近い指数で書かれ、厳密でない商は $p$ 桁すべてを使います。平方根では、 $\lfloor e/2 \rfloor$ は、その 2 倍が $e$ 以下で最大の偶数となる指数です。表はパッケージが使う規則をまとめたものです。

演算	理想指数 $e^\ast$
`add`、`subtract`、`fma`（和の部分）	$\min(e_1, e_2)$
`multiply`、`fma`（積の部分）	$e_1 + e_2$
`divide`	厳密なら $e_1 - e_2$ 、そうでなければ $p$ 桁
`sqrt`	厳密なら $\lfloor e/2 \rfloor$ 、そうでなければ $p$ 桁
`remainder`, `remainder_near`	$\min(e_1, e_2)$
`divide_integer`, `to_integral_*`	整数値には $\max(e, 0)$ 、商には 0
`quantize`, `rescale`	要求された指数
`scaleb`	$e + n$
整数 $n$ での `power`、厳密な場合	$n e$
`exp`、`ln`、`log10`、整数でない `power`	$p$ 桁（`exp(0)`、`ln(1)`、`log10(10^k)` を除き、常に厳密でない）
`plus`, `minus`, `abs`, `apply`	$e$

ゼロの結果には先頭の桁がないので、その指数は単に理想指数を $[E_{\mathrm{tiny}}, e_{\max}]$ （または $E_{\mathrm{top}}$ ）にクランプしたものです。厳密にゼロとなる和の符号は、両方のオペランドが負である場合、またはモードが Floor でどちらかが負である場合を除いて $+$ です。これは、 $-\infty$ 方向への丸めの場合を除いて $x - x = +0$ とする IEEE の規則の 10 進版です。

///|
test "ideal exponents" {
  let ctx = @decimal_gda.GdaContext::decimal64()
  let d = (s : String) => @decimal_gda.Decimal::from_string(s).unwrap()
  inspect(@decimal_gda.add(d("1.30"), d("1.2"), ctx).value(), content="2.50")
  inspect(@decimal_gda.multiply(d("1.30"), d("1.2"), ctx).value(), content="1.560")
  inspect(@decimal_gda.divide(d("2.40"), d("2"), ctx).value(), content="1.20")
  inspect(@decimal_gda.divide(d("1"), d("4"), ctx).value(), content="0.25")
  inspect(@decimal_gda.sqrt(d("1.00"), ctx).value(), content="1.0")
  inspect(@decimal_gda.subtract(d("1.0"), d("1.00"), ctx).value(), content="0.00")
}

精度への丸め

厳密な結果を $(-1)^s c\,10^{e}$ 、 $d(c) = p + k$ 、 $k \ge 1$ とします。下位 $k$ 桁を割り落として $c = q\,10^{k} + r$ 、 $0 \le r < 10^k$ とします。すべての丸めモードは、インクリメント $\delta \in \{0, 1\}$ を用いて $(-1)^s (q + \delta)\,10^{e+k}$ を返し、8 つのモードは $\delta$ だけが異なります（should_increment_decimal_repr）。

\begin{aligned} \delta_{\mathrm{Down}} &= 0, & \delta_{\mathrm{Up}} &= [r > 0],\\ \delta_{\mathrm{Ceiling}} &= [r > 0 \wedge s = 0], & \delta_{\mathrm{Floor}} &= [r > 0 \wedge s = 1],\\ \delta_{\mathrm{HalfUp}} &= [2r \ge 10^k], & \delta_{\mathrm{HalfDown}} &= [2r > 10^k],\\ \delta_{\mathrm{HalfEven}} &= [2r > 10^k \vee (2r = 10^k \wedge q \text{ odd})], & \delta_{\mathrm{ZeroFiveUp}} &= [r > 0 \wedge q \equiv 0 \pmod 5]. \end{aligned}

ZeroFiveUp は「保持される最後の桁が 0 または 5 でない限り、ゼロ方向へ丸める」と読みます。その目的は、後でどのモードでより少ない桁へ丸めても、最初の丸めの影響を受けないようにすることです。 $2r$ と $10^k$ の比較は 10 進のリム上で厳密に行われます（gda_coeff_div_pow10_round_info_repr は商、 $r > 0$ かどうか、 $2r - 10^k$ の符号を返します）。 $q + \delta = 10^p$ であれば、結果は $10^{p-1} \cdot 10^{e+k+1}$ に再正規化されます。丸めの前に、最終化ルーチンはまず $c$ の末尾のゼロを落とそうとします。厳密な係数が長すぎる原因がゼロだけであれば、結果は厳密であり、Rounded だけが立てられます。

$u = 10^{e+k}$ を結果 $\hat x$ の最終桁の単位とします。half 系のモードでは $|\delta\,10^k - r| \le 10^k/2$ 、それ以外では $< 10^k$ なので、

|\hat x - x| \le \tfrac12 u \;\;(\text{half modes}), \qquad |\hat x - x| < u \;\;(\text{others}),

そして $c \ge 10^{p+k-1}$ から $|x| \ge 10^{p-1} u$ が得られるので、

\frac{|\hat x - x|}{|x|} \le \tfrac12\,10^{1-p} \;\;(\text{half modes}), \qquad \frac{|\hat x - x|}{|x|} < 10^{1-p} \;\;(\text{others}).

これは標準モデル $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \varepsilon)$ 、 $|\varepsilon| \le \mathbf u$ における 10 進の単位丸め誤差 $\mathbf u = \frac12 10^{1-p}$ です。²2 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.2. 10 進の wobble は 2 進より大きく、相対間隔は 1 つの decade の中で 10 倍変化します（2 進では 2 倍）（Goldberg 1991, §1.2）。丸めは、 $k \ge 1$ 桁が取り除かれたときには常に Rounded を、取り除かれた桁にゼロでないものがあったとき（ $r > 0$ ）には常に Inexact を立てます。

この限界と、このページにある他の丸めに関する事実の完全な証明は、添付資料にあります。

decimal_gda の丸めに関する証明

非正規化数の結果と $E_{\mathrm{tiny}}$

厳密な結果が $\hat e < e_{\min}$ を持つとき、それは極小です。極小な結果が必要とする単位は高々 $10^{E_{\mathrm{tiny}}}$ なので、丸めの位置は「 $p$ 桁を残す」ではなく「 $10^{E_{\mathrm{tiny}}}$ 以上の桁を残す」になります。厳密な指数が $e < E_{\mathrm{tiny}}$ のとき、最終化ルーチンは $k = E_{\mathrm{tiny}} - e$ 桁を取り除きます（すべての桁を取り除くこともあります）。このとき絶対誤差は非正規化数の単位で抑えられます。

|\hat x - x| \le \tfrac12\,10^{E_{\mathrm{tiny}}} \;\;(\text{half modes}),

しかし相対誤差はもはや $\mathbf u$ では抑えられません。 $|x|$ が $10^{e_{\min}}$ から $10^{E_{\mathrm{tiny}}}$ へ小さくなるにつれて、精度は $p$ 桁から 1 桁へと徐々に下がります。条件がこれを記録します。Subnormal はすべての極小な結果で立てられ（GDA、およびこのパッケージの GDA 関数は、厳密な $\hat e$ から丸めの前に極小性を判定します）、Underflow は極小な結果が厳密でもない場合に、Clamped はゼロに丸められた場合に立てられます（そのときゼロは指数 $E_{\mathrm{tiny}}$ をとります）。ステータスを持たない層では、IEEE 流の用途のために丸め後の判定（DecimalTininessDetection::AfterRounding）も提供されています。これはどの結果を極小とみなすかを変えるだけで、値は決して変えません。

クランプ

clamp が有効な場合、 $e > E_{\mathrm{top}}$ だが $\hat e \le e_{\max}$ である結果は、値としては表現可能ですが三つ組としては表現できません。 $\hat e \le e_{\max}$ は $d(c) + e - 1 \le e_{\max}$ を意味するので、係数に $e - E_{\mathrm{top}}$ 個のゼロを付け足すと

d\bigl(c\,10^{e - E_{\mathrm{top}}}\bigr) = d(c) + e - E_{\mathrm{top}} \le e_{\max} + 1 - E_{\mathrm{top}} = p,

となり、付け足した三つ組は高々 $p$ 桁で、値も同じです。パッケージはまさにこれを行い、Clamped を立てます。値は変わらず、コホートのメンバーだけが異なります。ゼロも同様に、指数を $[E_{\mathrm{tiny}}, E_{\mathrm{top}}]$ の中へ移すことでクランプされます。

オーバーフロー

丸めた結果が $\hat e > e_{\max}$ を持つとき、演算はオーバーフローします。Overflow、Inexact、Rounded が立てられ、結果は $\pm\infty$ か、同じ符号の最大の有限数 $N_{\max} = (10^p - 1) \cdot 10^{E_{\mathrm{top}}}$ のどちらかになります。どちらになるかは、 $\infty$ を $N_{\max}$ の先にある表現可能な数とみなし、モード自身の方向を適用することで決まります。half 系のモードと Up はゼロから遠ざかる方向に動くので $\infty$ に達し、Down と ZeroFiveUp はゼロ方向に動くので（ZeroFiveUp がゼロから遠ざかる方向に丸めるのは最後の桁が 0 または 5 の場合だけで、 $N_{\max}$ の最後の桁は 9 です） $N_{\max}$ で止まります。Ceiling は正の結果に対して $+\infty$ 、負の結果に対して $-N_{\max}$ を与え、Floor はその鏡像になります（overflow_to_infinity）。

モード	正のオーバーフロー	負のオーバーフロー
`HalfEven`, `HalfUp`, `HalfDown`, `Up`	$+\infty$	$-\infty$
`Down`, `ZeroFiveUp`	$+N_{\max}$	$-N_{\max}$
`Ceiling`	$+\infty$	$-N_{\max}$
`Floor`	$+N_{\max}$	$-\infty$

条件、シグナル、トラップ

GDA は、条件（何が起きたか：ゼロ除算、結果が丸められた、変換の入力が不正だった、…）、シグナル（条件が引き起こす名前付きのイベント）、トラップ（ユーザーが計算を中断するよう求めたシグナル）を区別します。パッケージは条件ごとに 1 つのフラグを保持し、1 つの規則で条件を GDA のシグナルに対応付けます。4 つの詳細な無効条件 ConversionSyntax、DivisionImpossible、DivisionUndefined、InvalidContext はすべて InvalidOperation をシグナルします。形式的には、 $\mathcal S$ を 13 個のフラグ、 $\iota \subset \mathcal S$ を無効系統の 5 個のフラグとします。フラグ集合 $R$ とシグナル $\sigma$ に対して、

\sigma \in^\ast R \iff \begin{cases} R \cap \iota \neq \emptyset & \sigma = \mathrm{InvalidOperation},\\ \sigma \in R & \text{otherwise,} \end{cases}

これが GdaFlags::contains です。このとき演算は次の状態機械になります。

(\vec x, C) \;\longmapsto\; \begin{cases} \mathrm{Completed}(v, C, \emptyset) & R = \emptyset,\\ \mathrm{Completed}(v, C', R) & R \neq \emptyset,\ \tau = \bot,\\ \mathrm{Trapped}(\tau, v, C', R) & \tau \neq \bot, \end{cases}

ここで $(v, R) = F(\vec x, \pi(C))$ は、オペランドとポリシー $\pi(C)$ （精度、丸め、指数範囲、クランプ、extended）だけから計算される結果と立てられたフラグであり、

C'.\mathrm{status} = C.\mathrm{status} \cup R \cup \bigl\{\mathrm{InvalidOperation} \mid R \cap \iota \neq \emptyset\bigr\}, \qquad C'.\pi = C.\pi,\quad C'.\mathrm{traps} = C.\mathrm{traps},

$\tau$ は、優先順位リスト InvalidOperation, DivisionByZero, DivisionUndefined, DivisionImpossible, InvalidContext, ConversionSyntax, Overflow, Underflow, Subnormal, Inexact, Rounded, Clamped, LostDigits のうち、 $\tau \in^\ast R$ かつ $\tau \in C.\mathrm{traps}$ を満たす最初のシグナル、そのようなものがなければ $\bot$ です（complete_gda と trapped_signal）。

定義から直ちに 3 つの性質が従い、ユーザーはこれらに依拠しています。

値はステータスにもトラップにも依存しない。 $v$ と $R$ は $(\vec x, \pi(C))$ の関数であり、ステータスとトラップは $C'$ と $\tau$ にしか関与しません。したがってトラップを有効にしても結果が変わることはなく、Trapped は定義済みの結果を運ぶことができます。
ステータスは単調かつ冪等である。 $C.\mathrm{status} \subseteq C'.\mathrm{status}$ であり、 $\cup$ は結合的、可換、冪等なので、受け渡された一連の演算の後のステータスは、立てられたすべての集合の和集合であり、列のグループ化の仕方に依存しません。
トラップの選択は決定的である。 優先順位リストはシグナル上の全順序なので、条件がどの順序で検出されたかにかかわらず、1 つの演算が選ぶトラップは高々 1 つです。

///|
test "status is the union of the raised sets" {
  let ctx = @decimal_gda.context(precision=3)
  let d = (s : String) => @decimal_gda.Decimal::from_string(s).unwrap()
  let a = @decimal_gda.divide(d("1"), d("3"), ctx) // Inexact, Rounded
  let b = @decimal_gda.divide(d("1"), d("0"), a.next_context()) // DivisionByZero
  let expected = a.raised().combine(b.raised())
  inspect(b.next_context().status() == expected, content="true")
  // Trapping changes the variant, never the value.
  let trapping = ctx.trap(Inexact)
  let t = @decimal_gda.divide(d("1"), d("3"), trapping)
  inspect(t.value() == a.value(), content="true")
  inspect(t is @decimal_gda.GdaOutcome::Trapped(Inexact, _, _, _), content="true")
}

設計上の判断

ステータスはグローバル状態ではなく `GdaOutcome` で受け渡す

問題。 GDA 仕様はコンテキストを、演算がステータスフラグを設定する可変オブジェクトとして記述しています。ほとんどの実装（decNumber、Python の decimal）はスレッドごとにカレントコンテキストを保持します。

選択肢。 (a) 参照で渡される可変なコンテキスト。(b) スレッドローカルまたはグローバルなカレントコンテキスト。(c) すべての結果とともに返される不変なコンテキスト。

選択：(c)。 すべての GDA 関数は GdaOutcome を返し、呼び出し側は next_context() を次の演算に渡します。理由は上記の性質です。 $F$ はポリシーしか読まないので、状態機械は純粋な数値計算の部分と純粋な記録管理の部分に分解され、どちらも参照透過です。式は、結果やフラグを変えることなく、再評価したり、メモ化したり、別のスレッドで実行したりできます。テストランナーはコンテキストのスナップショットを取り、そのもとで同じ演算を何度も実行できます。frontend/gda_expr の .decTest ランナーはまさにそれを行っています。また MoonBit にはスレッドローカルストレージがないので、(b) はプロセス全体のグローバル状態を意味することになり、Luna-Flow はそれを排除しています。代償は明示的な受け渡しですが、線形なパイプラインについては decimal_gda_checked パッケージがそれを不要にします。何も立てられなかった場合は入力のコンテキストがそのまま返されるので、厳密な演算が新しいコンテキストを割り当てることはありません。

トラップの定義済みの結果を保持する

GDA ではトラップは制御を移しますが、仕様は演算が返したはずの結果をそれでも定義しています。トラップをエラー（Result::Err、raise）として表現すると、その値が捨てられてしまいます。Trapped は値、次のコンテキスト、立てられた集合を保持するので、呼び出し側はそれを検査したり、ログに記録したり、そこから再開したりできます。性質 1 により、それはトラップがない場合と同じ値であることが保証されます。

詳細な無効条件を保持する

仕様は ConversionSyntax、DivisionImpossible、DivisionUndefined、InvalidContext を InvalidOperation を通じて報告します。パッケージはこれらを別々のフラグとして保持し、かつそれぞれについて contains(InvalidOperation) を真にし、いずれかが立てられたときは常にステータスに invalid_operation を設定します。8 つの GDA シグナルしか知らないプログラムにはちょうど GDA の振る舞いが見え、一方でテストハーネスや診断では不正なリテラルと $0/0$ を区別できます。InvalidOperation を優先順位リストの先頭に置くことで、仕様が要求するとおり、InvalidOperation のトラップが 4 つすべてを捕捉します。

1 つの最終化ルーチン、特殊値を先に

各演算はまず特殊なケース（NaN の伝播、無効演算、無限大、厳密なゼロ）を決定し、次に厳密な有限の係数と指数を構築し、その後で初めて最終化ルーチンを呼び出します。係数のアルゴリズムがフラグを設定することはありません。これにより、条件は厳密な結果とポリシーの関数となり、どの乗算・除算カーネルが選ばれたかに依存しません。また、標準に関わる振る舞いに触れずにカーネルを調整できます。

`decimal` から独立した GDA パッケージ

IEEE 754-2008 の 10 進算術は GDA から発展したものなので、両者はほとんどの有限の結果で一致しますが、契約は異なります。

観点	`decimal_gda` (GDA 1.70)	`decimal` (IEEE 754-2019)
精度	コンテキストごとに任意の $p \ge 1$	形式の $p$ （または選択した値）
丸めモード	`HalfUp`、`HalfDown`、`ZeroFiveUp` を含む 8 種類	IEEE の属性
条件	コンテキスト内のスティッキーなステータス、定義済みの結果を伴うトラップ	値とともに返される演算ごとのフラグ
詳細な無効条件、`LostDigits`、サブセット算術	あり	いいえ
極小性	丸め前	選択可能
初等関数	`sqrt`, `exp`, `ln`, `log10`, `power`	IEEE の推奨演算の集合
交換形式	DPD	DPD と BID

共通のコアを持つと、両方の状態モデルの和を抱え込まなければならず、一方の標準のために行った変更が、もう一方の結果を黙って変えてしまうおそれがあります。そのためパッケージは独自の値の型、係数カーネル、コンテキスト、最終化ルーチン、DPD コーデックを持ち、本番用の依存関係スキャンによって decimal を決してインポートしないことを検査しています。係数の閾値は現在 decimal のものと等しいですが、これは計測による偶然の一致であり、共有された依存関係ではありません。

ステータスを持たない層（DecimalContext、DecimalFlags、*_ctx メソッド）は、パッケージ自身のエンジンを公開したものです。IEEE のように演算ごとにフラグを返すのでアダプタにとって便利です（Luna-Flow/arithmetic のトレイト実装はこれを使っています）が、実装しているのは GDA の算術です。IEEE 754 の契約は decimal にあります。

等価性が証明された高速経路

2 種類の近道が、観測可能な結果を何も変えずに汎用の機構を回避します。

小さな厳密な整数。 parse、add、subtract、multiply、fma はまず次の条件を確認します。すべてのオペランドが指数 0 で係数が $10^{18}$ 未満の整数であること、厳密な結果が $10^{18}$ 未満で高々 $p$ 桁のゼロでない整数であること、その調整指数が $[e_{\min}, e_{\max}]$ にあること、クランプによって指数 0 が許されること、そしてコンテキストが extended であることです。これらの述語のもとでは、汎用の経路は丸めを行わず、条件を立てず、指数 0（これらの演算すべての理想指数）を返すので、近道は同じ $v$ で Completed(v, C, none) を返します。ゼロの結果は、その符号が丸めモードに依存するため除外されます。

吸収される加数。 $x_1$ と、 $\hat e_2 < \hat e_1 - p + 1$ を満たすはるかに小さい $x_2$ を加えるとき、厳密な和を作る必要はありません。 $x_2$ のすべての桁は $x_1$ の保持される最後の桁より下にあるので、 $x_2$ が丸めに影響するのは、その符号と半単位との比較を通じてだけです。

\hat x = \mathrm{round}\bigl(x_1 + x_2\bigr) = \bigl(q + \delta(\operatorname{sign} x_2,\; |x_2| \lessgtr \tfrac12 u)\bigr) u,

これは、上の丸めの表で $r$ を同じ比較クラスのスティッキーな代表値に置き換えたものです。exact_base_small_addend_result はまさにそれを評価し、compare_magnitude_to_half_ulp が $x_2$ と $u/2$ を厳密に比較します。 $x_1$ が 10 のべきで $x_2$ が逆の符号を持つ特別な場合（ $x_1$ のすぐ下の単位が 10 分の 1 になる）は除外され、汎用の経路をとります。

除算

divide は、 $0$ 、 $\infty$ 、10 のべきによる厳密な除算、小さな厳密な除数による除算を別々に扱います。それ以外の場合は、被除数を $t = p + d(c_2) + 2$ として $10^{t}$ でスケールし、 $Q = c_1 10^{t} / c_2 \ge 10^{p+2}$ が少なくとも $p + 3$ 桁の整数部を持つようにします。そして $Q$ をコンテキストのモードで整数 $q$ に丸め、さらに $q$ を $p$ 桁に丸めます。商が有限小数になるかどうかは厳密に判定され（has_finite_decimal_expansion_repr： $c_1/c_2$ が有限小数になるのは $c_2 / \gcd(c_1, c_2)$ が 2 と 5 以外の素因数を持たないとき、かつそのときに限る）、それによって厳密な商には理想指数のコホートが選ばれ、それ以外には $p$ 桁の係数で Inexact が強制されます。

方向付き丸めのモードでは、切り捨ては合成できるので、2 回の連続した丸めは 1 回の丸めと等価です： $\lfloor \lfloor y/10^j \rfloor / 10^m \rfloor = \lfloor y/10^{j+m} \rfloor$ 。half 系のモードでは、1 回目の丸めがタイを作り出す場合、すなわち $Q \ne q$ でありながら $q$ の捨てられる桁がちょうど $50\cdots0$ になる場合を除いて等価です。

警告

このケースは現在のブランチで発生します。精度 1 で $1/2222 = 0.00045004\ldots$ の場合、1 回目の丸めで $q = 4500$ （単位 $10^{-7}$ ）となり、2 回目の丸めはタイとみなして最近接偶数丸めで 0.0004 に切り下げますが、GDA は 0.0005 を要求します。 $c_1 < 20$ 、 $c_2 < 3000$ 、 $p \le 3$ を走査したところ、このような商が 12 個見つかりました。固定されたテストスイートにはこのケースは含まれていません。最終化ルーチンが他の箇所で行っているように、2 回目の丸めのために $q$ とともにスティッキービット（ $Q \notin \mathbb Z$ ）を保持すれば、この欠陥は解消されます。これが唯一の失敗モードであることの証明は添付資料にあります。

平方根

sqrt はまず厳密な平方根を試みます。末尾のゼロを取り除き、指数を偶数にし、Newton 反復 $a \leftarrow \lfloor (a + \lfloor c/a \rfloor)/2 \rfloor$ によって $s^2 + \rho = c$ となる整数平方根 $(s, \rho)$ を求め、 $\rho = 0$ なら受理し、その後 $p$ 桁の範囲で理想指数 $\lfloor e/2 \rfloor$ に向けてゼロを付け足します。そうでない場合は、最終的な位置で直接丸めます。平方根の指数を $f = \max(E_{\mathrm{tiny}}, \lfloor \hat e/2 \rfloor - p + 1)$ として $s = \lfloor \sqrt{c\,10^{e - 2f}} \rfloor$ を計算し、被開平数と中点の 2 乗を比較してインクリメントを決めます。これは整数で厳密に行えます。

\sqrt{N} \gtrless s + \tfrac12 \iff 4N \gtrless (2s + 1)^2 .

（まず $p$ 桁で、次に $E_{\mathrm{tiny}}$ で丸めるのではなく） $f$ で 1 回だけ丸めることで、極小な平方根の二重丸めを避けています。厳密な平方根は先に処理されているので、この比較が判定するのは無理数の平方根だけであり、それが中点に等しくなることはありえません。GDA 関数は常に最近接偶数丸めを行います。

整数べき

整数の指数 $n$ に対しては、結果は GDA 仕様の規定どおりに計算されます。厳密なべきが $p$ 桁に収まれば、指数 $n e$ で厳密に返されます。そうでなければ、作業精度 $w = p + d(|n|) + 2$ （サブセットコンテキストでは 1 桁少ない）で 2 進べき乗法を実行し、各積の後で最近接偶数丸めを行います。負の $n$ は $1/x$ から始め、最後の積をコンテキストへ丸めます。 $|n|$ 個の因子はそれぞれ高々 $|n| - 1$ 回の丸めを経るので、作業結果は $x^n(1 + \theta)$ となり、次が成り立ちます。²2 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.2. 10 進の wobble は 2 進より大きく、相対間隔は 1 つの decade の中で 10 倍変化します（2 進では 2 倍）（Goldberg 1991, §1.2）。

|\theta| \le (1 + \mathbf u_w)^{|n|-1} - 1 \approx |n|\,\mathbf u_w < 10^{d(|n|)} \cdot \tfrac12\,10^{1 - p - d(|n|) - 2} = \tfrac12\,10^{-1-p},

これは結果の最終桁の単位で高々 $0.05$ です。最後の丸めを加えると、厳密でない整数べきは、half 系のモードで最終桁の単位の $0.55$ 以内になります。これは正しい丸めに近いものの正しい丸めではありません。仕様も整数べきについてそれ以上は要求していません。結果が必ずオーバーフローまたはアンダーフローするほど大きな指数は、結果の調整指数に対する限界 $n(\hat e + 1) - 1$ と $n \hat e$ から事前に検出されます。

正しく丸められた `exp`、`ln`、`log10` と整数でない `power`

これらの関数は、認証付きの区間演算と Ziv 流の精緻化ループによって評価されます。³3 A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991. ボール演算については J. van der Hoeven, “Ball arithmetic”, 2009、および F. Johansson, “Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”, IEEE Trans. Computers 66(8), 2017 を参照してください。 $f(x)$ について：

まず厳密なケース。 $\exp 0 = 1$ 、 $\ln 1 = 0$ 、 $\log_{10} 10^k = k$ 、特殊なオペランド、定義域エラー、そして power が厳密に決定できるケース（sqrt による $x^{1/2}$ 、10 のべき、 $1^y$ 、確実なオーバーフローまたはアンダーフロー）はループに到達しません。
入力を包含する。 $x$ は方向付き丸め（ $-\infty$ 方向と $+\infty$ 方向への to_bin_float）によって $w$ ビットの 2 進ボール $[x^-, x^+]$ に変換されるので、厳密に $x \in [x^-, x^+]$ となります。
出力を包含する。 ball_float はボール上で $f$ を評価し、 $[L, U] \ni f(x)$ を返します（log10 では $\ln$ を $\ln 10$ のキャッシュされた包含区間で割り、power では $\ln x$ の包含区間を通じて $\exp(y \ln x)$ を計算します）。
候補を認証する。 10 進の近似値 $\tilde y$ （ボールの中心、または 10 進級数による評価）を $r$ に丸めます。 $(a, b)$ を、 $r$ に丸められる実数からなる $r$ のまわりの開区間とします。half 系のモードでは隣接値との 2 つの中点、方向付き丸めのモードでは $(r, r^+)$ または $(r^-, r)$ です。その端点は厳密な 10 進数であり、それぞれを 2 進で包含し、次が成り立つとき結果を受理します。 $a \le a^+ < L \le f(x) \le U < b^- \le b ,$ これにより $\mathrm{round}(f(x)) = r$ が証明されます。2 つ目の判定では、厳密な 2 進有理数の端点 $L$ と $U$ の両方をコンテキストへ丸めます。すべての丸めモードは単調なので、 $\mathrm{round}(L) = \mathrm{round}(U) = r$ からも $\mathrm{round}(f(x)) = r$ が証明されます。
精緻化する。 そうでなければ作業精度を $w \leftarrow w + \max(32, \lfloor w/2 \rfloor)$ と増やします。これを最大 12 回行います。

開始精度は $D = \max(d(c), p)$ として $w_0 = \max(128, 4D + 64)$ ビットです。10 進 1 桁あたり 4 ビットは $\log_2 10 \approx 3.32$ を上回るので、入力は 10 進の情報を失わずに運ばれ、64 ビットのガードビットが残ります。 $p \le 64$ の広いクランプなしのコンテキストで $\hat e = 0$ の引数に対しては、約 $\frac{10}{3} D + 12$ ビットを使う安価な初回試行を行い、認証できなければ $w_0$ にフォールバックします。

$f(x)$ 自身が表現可能な数でも中点でもなければ、包含区間が 1 点に縮むのでループは停止します。exp、ln、log10 では、ステップ 1 の後これは常に成り立ちます。Lindemann–Weierstrass の定理により、有理数 $x \ne 0$ に対して $e^x$ は超越数であり、したがって有理数 $x \ne 1$ に対して $\ln x$ は無理数です。また $\log_{10} x$ が有理数になるのは 10 の整数べきの場合だけです。表現可能な数と中点は有理数です。それでも予算を使い切った場合（たとえば厳密な値が中点で、ステップ 1 で捕捉されない power）、try_* メソッドは認証の失敗を返し、GDA 関数は InvalidOperation を伴う NaN を返します。証明されていない結果を返すことは決してありません。

GDA 関数 exp、ln、log10、sqrt は、丸めを HalfEven にしたコンテキストのコピーを渡します。仕様がこれらの関数を、コンテキストのモードにかかわらず最近接偶数丸めで正しく丸められるものとして定義しているからです。power は、仕様とそのテストベクタが要求するとおりコンテキストのモードを使い、整数でないべきについては、たまたま厳密であっても $p$ 桁で Inexact として報告します（decimal_power_noninteger_gda_result）。さらに 2 つの GDA の規則を文字どおりに実装しています。これらの関数は $p$ 、 $e_{\max}$ 、 $-e_{\min}$ が高々 999,999 のコンテキストでのみ定義され（math_context_is_restricted、それ以外では InvalidContext を立てます）、サブセット算術では ln は古典的な参照アルゴリズムの結果を再現します。これは正しく丸められた結果より最終桁の単位で 1 大きくなることがあります（decimal_ln_subset_result）。従来のサブセットのテストベクタがその結果を記録しているためです。

係数の表現とカーネル

係数は、 $10^{18}$ 未満の値に対する Small(UInt64) か、桁数をキャッシュした基数 $10^9$ のリムの配列のどちらかで、どちらも永続的です。演算がオペランドのリムを書き換えることはありません。10 進のリムにより、桁数の計算、末尾のゼロの除去、半単位との比較、ZeroFiveUp の最終桁の判定が、2 進から 10 進への変換なしに定数時間または線形時間で行えます。乗算と除算は、リム数に応じて、筆算、Karatsuba、Toom-3、2 法 NTT、Knuth の Algorithm D、Burnikel–Ziegler、Newton の逆数除算から選択します。閾値はターゲットごとに次のとおりです。

ターゲット	Karatsuba mul/square	Toom-3	first NTT mul/square	Burnikel–Ziegler	Newton
native	96 / 48	1,152	1,728 / 640	2,816 以上	無効
LLVM	96 / 96	2,048	4,096 / 2,048	2,048	4,096
Wasm、Wasm-GC、JS	96 / 96	4,096	8,192 / 4,096	2,048	4,096

閾値は計測に基づく性能上の方針であり、意味論ではありません。どのカーネルも厳密な積または商を返すので、選択によって結果やフラグが変わることはありません。計測結果は性能を参照してください。

正しさ／不変条件

1 回だけの丸め。 整数の power と上記の除算のケースを除くすべての有限の結果は、厳密な結果を 1 回だけ丸めたものです。したがって half 系のモードでは $|\hat x - x| \le \frac12 u$ 、それ以外では $< u$ であり、非正規化数の範囲外では $|\hat x - x|/|x| \le \frac12 10^{1-p}$ です。
理想のコホート。 厳密な結果は、その理想指数に最も近い表現可能な指数で返されるので、厳密なデータに対する演算は量子を保ちます（2.50 × 3 = 7.50）。
ステータスからの独立性。 $F(\vec x, C)$ はポリシーを通じてのみ $C$ に依存します。結果がスティッキーなステータスやトラップに依存することはありません。
ステータスの代数。 ステータスは増えるだけであり、受け渡された一連の演算の後のステータスは、立てられた集合の和集合です。
トラップの決定性。 1 つの演算で発火するトラップは高々 1 つで、固定された全順序によって選ばれます。Trapped は Completed が運ぶはずの値と同じ値を運びます。
認証付きの初等関数。 exp、ln、log10、整数でない power の有限で厳密でない結果は、それが正しく丸められた値であることの証明（添付資料の補題 3 または 4）とともにしか返されません。証明できなければ NaN と InvalidOperation になります。
全順序。 compare_total は三つ組上の全順序（符号、次にクラス、次に値、次に指数、次にペイロード）であり、Decimal::compare は NaN がすべての数より上の 1 つのクラスをなす全前順序です。
証拠。 固定された official テストスイートでは、正当な実行可能スカラー行 64,986/64,986 が、従来の official0 スイートでは 16,124/16,124 が合格しています（適合性）。これらは有限の主張です。上記の除算の欠陥と、API ページに記載されている to-integral の相違は、固定されたどの行でもカバーされていません。

却下した代替案

可変またはグローバルなカレントコンテキスト。 結果とフラグが評価順序と隠れた状態に依存してしまうため不採用としました（最初の設計判断を参照）。
エラーとしてのトラップ。 GDA で定義された結果が失われるため不採用としました。
decimal とエンジンを共有すること。 一方の標準の変更がもう一方の結果を変えうるため不採用としました。計測された重複は、偶発的な意味論上の結合よりも安上がりです。
固定数のガード桁を持つ 2 進浮動小数点による初等関数。 固定数のガード桁では正しい丸めを保証できないため不採用としました（table-maker’s dilemma）。精緻化を伴う認証は正しい丸めを保証し、できない場合は目に見える形で失敗します。
すべての結果を正規化すること。 量子は GDA の数の一部なので不採用としました。正規化は reduce として明示的に利用できます。

境界

パッケージは GDA の演算だけを実装し、それ以外は何も実装しません。仕様外の三角関数、双曲線関数、その他の関数はなく、ステータスを持たない層のわずかな追加機能を除いて IEEE 754 の演算もありません。
.decTest ファイルの解析、テストスイートの実行、ファイルの読み込みは行いません。それらは frontend/gda_expr とリポジトリのツールの役割です。
BID の交換形式は提供しません。DPD のみです。
整数べきについて GDA の要求を超える正しい丸めは提供しません。また現在のブランチでは、上記の half 系モードでの除算のケースについても正しい丸めを提供しません。
可変またはグローバルなコンテキストは提供せず、コンテキストは同一性を持ちません。同じフィールドを持つ 2 つのコンテキストは交換可能です。

M. F. Cowlishaw, General Decimal Arithmetic Specification, version 1.70 (2009), https://speleotrove.com/decimal/decarith.html。テストスイートは同じ著者による dectest コレクションのバージョン 2.62 です。 ↩
N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.2. 10 進の wobble は 2 進より大きく、相対間隔は 1 つの decade の中で 10 倍変化します（2 進では 2 倍）（Goldberg 1991, §1.2）。 ↩ ↩²
A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991. ボール演算については J. van der Hoeven, “Ball arithmetic”, 2009、および F. Johansson, “Arb: efficient arbitrary-precision midpoint-radius interval arithmetic”, IEEE Trans. Computers 66(8), 2017 を参照してください。 ↩