decimal の設計

このページでは Luna-Flow/floating/decimal の算術モデルを説明します。10 進浮動小数点数とは何か、コホートと推奨指数がどのように情報を運ぶか、交換形式のエンコーディングが数字をどのようにビットに詰め込むか、すべての結果がどのようにちょうど 1 回だけ丸められるか、そこからどのような誤差限界が従うか、指数範囲がどのように強制されるか、そして初等関数がどのように認証されるか、です。各関数の仕様は decimal API に、使用例は decimal チュートリアルにあります。

設計目標

decimal は IEEE 754-2019¹1 IEEE Std 754-2019, Standard for Floating-Point Arithmetic、3.3–3.5 節（10 進形式とエンコーディング）、4 節（属性と丸め）、5 節（演算）、7 節（例外）、9 節（推奨演算）。M. F. Cowlishaw による General Decimal Arithmetic 仕様（バージョン 1.70）は同じモデルを任意精度の形で与えています。本ページではその用語 coefficient（係数）、adjusted exponent（調整指数）、Etiny、clamp を使います。の 10 進算術を任意の精度で実装しており、次の 3 つの性質を持ちます。

すべてのコンテキスト演算は正しく丸められます。 結果は、厳密な数学的結果を、選択された方向に、コンテキストの精度と指数範囲へ 1 回だけ丸めたものです。これは基本演算でも初等関数でも同じです。
何も黙って失われません。 結果の指数（その量子）、ゼロの符号、NaN のペイロード、およびすべての例外条件は、返される値または返される DecimalFlags の一部になります。
隠れた状態はありません。 精度、丸め、指数範囲、フラグは、Luna-Flow のどこでもそうであるように、明示的に渡され返される通常の不変な値です。

数学的背景

10 進浮動小数点数

精度 $p$ 、調整指数の範囲 $[e_{\min}, e_{\max}]$ の 10 進浮動小数点形式は、次の数の集合です。

x = (-1)^s \cdot c \cdot 10^{q}, \qquad s \in \{0,1\},\quad c \in \mathbb{Z},\ 0 \le c < 10^{p},\quad E_{\text{tiny}} \le q \le E_{\text{top}},

これに $\pm\infty$ と NaN を加えます。ここで $c$ は係数、 $q$ は指数または量子であり、

E_{\text{tiny}} = e_{\min} - p + 1, \qquad E_{\text{top}} = e_{\max} - p + 1 .

ゼロでない $x$ の調整指数は $\operatorname{adj}(x) = q + \operatorname{digits}(c) - 1 = \lfloor \log_{10} |x| \rfloor$ です。これは $x$ を科学的記数法 $d_0.d_1d_2\ldots \times 10^{\operatorname{adj}(x)}$ で書いたときの指数です。ゼロでない $x$ は $\operatorname{adj}(x) \ge e_{\min}$ のとき正規化数、そうでないとき非正規化数です。最小の正の非正規化数は $10^{E_{\text{tiny}}}$ 、最大の有限値は次のとおりです。

N_{\max} = (10^{p} - 1)\cdot 10^{E_{\text{top}}} = 10^{e_{\max}+1} - 10^{e_{\max}-p+1}.

DecimalContext が保持するのは、ちょうど $p$ 、 $e_{\min}$ 、 $e_{\max}$ 、丸めモード、clamp（ $E_{\text{top}}$ を超える指数を許すかどうか。クランプを参照）、および極小性（tininess）の規則です。交換形式は次のとおりです。

形式	$p$	$e_{\max}$	$e_{\min}$	$E_{\text{tiny}}$	$E_{\text{top}}$	バイアス $=-E_{\text{tiny}}$	指数の個数 $E_{\text{top}}-E_{\text{tiny}}+1$
decimal32	7	96	$-95$	$-101$	90	101	$192 = 3\cdot 2^{6}$
decimal64	16	384	$-383$	$-398$	369	398	$768 = 3\cdot 2^{8}$
decimal128	34	6144	$-6143$	$-6176$	6111	6176	$12288 = 3\cdot 2^{12}$

IEEE 754 は $e_{\min} = 1 - e_{\max}$ と定めているので、指数の個数は $E_{\text{top}} - E_{\text{tiny}} + 1 = e_{\max} - e_{\min} + 1 = 2e_{\max}$ です。標準は $e_{\max} = 3 \cdot 2^{w-1}$ を選んでこの個数を $3 \cdot 2^{w}$ にしています。これはちょうど、値 $0,1,2$ をとる先頭 2 ビットの指数ビットと、さらに $w$ ビットで符号化できる個数です（エンコーディングを参照）。 $q$ を非負の格納指数に変換するバイアスは $-E_{\text{tiny}} = p - e_{\min} - 1$ で、decimal64 では $16 + 383 - 1 = 398$ です。

なぜ 10 進か

既約分数で表した有理数 $n/m$ が基数 $b$ で有限展開を持つのは、 $m$ のすべての素因数が $b$ を割り切るとき、かつそのときに限ります。 $b = 2$ で許される分母は 2 のべきだけで、 $b = 10$ では $2^{i}5^{j}$ です。したがって、すべての 2 進浮動小数点数は有限の 10 進展開を持ちますが、 $0.1 = 1/(2\cdot 5)$ は 2 進では有限展開を持ちません。これに最も近い Double は次のとおりです。

0.1000000000000000055511151231257827021181583404541015625 = \frac{3602879701896397}{2^{55}} .

したがって、10 進で定義される量（価格、利率、計測値、プロトコルのフィールド）は 10 進浮動小数点で厳密に表現され、10 進の丸めは人や規則が指定する小数位で行われます。その代償は、より大きな wobble（後述）と、より高価な数字の算術です。

コホートと量子

写像 $(s, c, q) \mapsto (-1)^s c\,10^{q}$ は単射ではありません。同じゼロでない値のすべての表現は、その値のコホートをなします。 $c$ が $d$ 桁で末尾に $k$ 個のゼロを持つとき、指数範囲を無視すれば、メンバーは $-k \le j \le p - d$ に対する $(c\cdot 10^{j},\, q - j)$ であり、コホートは $p - d + k + 1$ 個のメンバーを持ちます。たとえば decimal32 における $1000$ （ $c = 1$ 、 $q = 3$ 、 $d = 1$ 、 $k = 0$ ）は、 $1\text{E+}3, 10\text{E+}2, \ldots, 1000000\text{E-}3$ の 7 つのメンバーを持ちます。ゼロはすべての指数に対してメンバーを持ちます。

コホートのメンバーは、数値的な等価性では表せない情報を運びます。12.30 は小数 2 桁を、1.2E+3 は有効数字 2 桁を表します。そのため IEEE 754 は各演算について推奨指数を定めており、厳密な結果は、指数がそれに最も近いメンバーで返されます。推奨指数は、厳密な結果が自然に属する位置から決まります。

\begin{aligned} c_a 10^{q_a} \pm c_b 10^{q_b} &= \bigl(c_a 10^{q_a - m} \pm c_b 10^{q_b - m}\bigr)\,10^{m}, & m &= \min(q_a, q_b),\\ c_a 10^{q_a} \cdot c_b 10^{q_b} &= (c_a c_b)\,10^{q_a + q_b},\\ c_a 10^{q_a} / c_b 10^{q_b} &= (c_a / c_b)\,10^{q_a - q_b},\\ \sqrt{c\,10^{q}} &= \sqrt{c\,10^{q - 2\lfloor q/2\rfloor}}\;10^{\lfloor q/2 \rfloor},\\ x\cdot y + z &: \ \min(q_x + q_y,\ q_z). \end{aligned}

和と積では括弧内の係数が整数なので、厳密な係数が $p$ 桁に収まる限り推奨指数が達成されます：1.20 + 3.40 = 4.60、1.25 × 2.50 = 3.1250。商では、 $c_a / c_b$ が有限の 10 進展開を持つ場合にのみ厳密な結果が存在し、そのとき $p$ 桁が許す範囲で $q_a - q_b$ に向けて移動されるので、2.400 / 1.2 = 2.00 となります。厳密でない結果は常に $p$ 桁すべてを使い、これは指数が最小のメンバーです。quantize は指数を明示的な引数とし、reduce_ctx/normalized は指数が最大のメンバーを選びます。

///|
test "design: preferred exponents" {
  let ctx = @decimal.DecimalContext::decimal64()
  let d = fn(s : String) { @decimal.Decimal::from_string(s).unwrap() }
  inspect(d("1.20").add_ctx(d("3.40"), ctx).0, content="4.60")
  inspect(d("1.25").mul_ctx(d("2.50"), ctx).0, content="3.1250")
  inspect(d("2.400").div_ctx(d("1.2"), ctx).0, content="2.00")
  inspect(d("0.0400").sqrt_ctx(ctx).0, content="0.20")
  inspect(d("1.5").fma_ctx(d("2.0"), d("0.25"), ctx).0, content="3.25")
}

丸めの方向

$x > 0$ を厳密な値とし、目標指数を $t$ （ $p$ 桁を残す指数、または極小な結果では $E_{\text{tiny}}$ ）とします。次のように書きます。

x \cdot 10^{-t} = c + f, \qquad c \in \mathbb{Z}_{\ge 0},\ 0 \le f < 1 .

各丸め方向は $c$ または $c + 1$ （に $10^{t}$ を掛けたもの）を返します。どちらを選ぶかは $f$ 、 $c$ の最後の桁、および符号に依存します。

モード	IEEE での名前	$c+1$ を返す条件（ $f > 0$ のとき）
`Down`	roundTowardZero	しない
`Up`	—	常に
`Ceiling`	roundTowardPositive	$x > 0$
`Floor`	roundTowardNegative	$x < 0$
`HalfUp`	roundTiesToAway	$f \ge \tfrac12$
`HalfDown`	—	$f > \tfrac12$
`HalfEven`	roundTiesToEven	$f > \tfrac12$ 、または $f = \tfrac12$ かつ $c$ が奇数
`ZeroFiveUp`	—	$c \bmod 5 = 0$

負の $x$ に対しては、同じ規則を $|x|$ に適用してから符号を戻すので、Ceiling と Floor が入れ替わります。すべてのモード $\circ$ は単調です： $x \le y \Rightarrow \circ(x) \le \circ(y)$ 。初等関数の認証が依拠しているのはこの単調性です。

誤差モデル

$x$ を正規化数の範囲にあるとし、 $10^{e} \le |x| < 10^{e+1}$ とします。この 10 進の桁区間（decade）における表現可能な数は、最終桁の 1 単位 $\operatorname{ulp}(x) = 10^{e - p + 1}$ の間隔で並んでいます。最近接への丸めの誤差はその半分以下なので、

\begin{aligned} \frac{|\operatorname{fl}(x) - x|}{|x|} \le \frac{\tfrac12\, 10^{e-p+1}}{10^{e}} = \tfrac12\, 10^{1-p} =: u , \end{aligned} \qquad\text{equivalently}\qquad \operatorname{fl}(x) = x(1 + \delta),\ |\delta| \le u .

この限界は decade の下端付近で達成されます。上端付近 $|x| \approx 10^{e+1}$ では、同じ絶対誤差は相対誤差にしてわずか $\tfrac12 10^{-p}$ です。したがって、1 つの decade 内での最悪と最良の相対誤差の比は次のとおりです。

\frac{\tfrac12\,10^{1-p}}{\tfrac12\,10^{-p}} = 10 = \beta ,

これが基数 $\beta = 10$ の wobble です。²2 Goldberg, “What every computer scientist should know about floating-point arithmetic”, ACM Computing Surveys 23(1), 1991, §1.2; Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.1–2.2. 2 進では wobble は 2 なので、同じ記憶容量では 10 進の方が最悪ケースの相対誤差がわずかに大きくなります。decimal64 は $u = \tfrac12 10^{-15} = 5\cdot 10^{-16}$ 、binary64 は $u = 2^{-53} \approx 1.1 \cdot 10^{-16}$ です。方向付き丸めのモードでは $|\delta| < 10^{1-p} = 2u$ です。マシンイプシロン、すなわち 1 から次に大きい数までの距離は $10^{1-p}$ です。epsilon_contextual はこれを返し、decimal64 での next_plus(1) は 1.000000000000001 です。1 未満では間隔が 10 分の 1 になります：next_minus(1) は 0.9999999999999999 です。

非正規化数の結果では間隔が固定の $10^{E_{\text{tiny}}}$ なので、誤差限界は絶対誤差 $|\operatorname{fl}(x) - x| \le \tfrac12 10^{E_{\text{tiny}}}$ になります。すべてのコンテキスト演算は正しく丸められるので、標準モデル $\operatorname{fl}(a \circ b) = (a \circ b)(1+\delta)$ 、 $|\delta| \le u$ は、 $+,-,\times,/,\sqrt{\ }$ 、fma、および結果が正規化数であるすべての初等関数について成り立ちます。したがって Higham の古典的な前進誤差解析と後退誤差解析が $u = \tfrac12 10^{1-p}$ でそのまま適用できます。

設計上の判断

1 つの表現、多数のコンテキスト

問題。 アプリケーションは、型の間で変換することなく、decimal32/64/128 と任意精度、IEEE の意味論と GDA 互換性を必要とします。

選択肢。 形式ごとに 1 つの型（ハードウェアのように）。形式をコンテキストで運ぶ 1 つの任意精度型。形式でパラメータ化された型。

選択。 符号、任意長の係数、指数、クラス、NaN の種類を示すビット、作業精度フィールドを保持する 1 つの Decimal 型と、別個の不変な DecimalContext。decimal64 の計算は DecimalContext::decimal64() のもとでの計算であり、交換形式のエンコーダはエンコードの前に形式のコンテキストを適用します。これにより算術が 1 か所にまとまり、呼び出し側は 50 桁で作業して最後に decimal64 へ丸めることができ、General Decimal Arithmetic のモデルとも一致します。値の中の精度フィールドは、コンテキストを持たない演算子と変換のためだけに使われます。

1 か所で 1 回だけ丸める

問題。 二重丸め（すでに丸められた値を再び丸めること）は、正しく丸められた結果を変えてしまうことがあります。

選択。 すべての有限なコンテキスト結果は 1 つの最終化ルーチンを通ります。このルーチンは厳密な結果 $(s, C, Q)$ （ $C$ は $p$ よりはるかに長いこともあります）を受け取り、 $p$ 桁への丸め、オーバーフローの検査、非正規化数のグリッドへの丸め、非正規化数／アンダーフローのフラグ設定、フォールドダウンをこの順に行います。演算のカーネルは厳密な整数を計算するだけで、丸めは一切行いません。例外は厳密な結果が無限桁になる演算（除算、平方根、初等関数）で、これらについては後述します。いずれも最後の桁は厳密な情報から決定します。

丸め桁とスティッキー情報の求め方

$C$ （ $D$ 桁の非負整数）を $D - s$ 桁に丸めるには、10 のべきで割ります。

C = Q \cdot 10^{s} + R, \qquad 0 \le R < 10^{s},

そして $2R$ と $10^{s}$ を比較して $Q$ と $Q+1$ のどちらにするかを決めます。この 1 回の比較が、古典的な丸め桁とスティッキービットの情報をちょうど担っています。丸め桁 $r \in \{0,\ldots,9\}$ と残り $0 \le R' < 10^{s-1}$ を用いて $R = r\,10^{s-1} + R'$ と書くと、

2R - 10^{s} = 2(r - 5)\,10^{s-1} + 2R' ,

そして $0 \le 2R' < 2\cdot 10^{s-1}$ なので、

\begin{aligned} r \ge 6 &\implies 2R - 10^{s} \ge 2\cdot 10^{s-1} > 0,\\ r = 5 &\implies \operatorname{sign}(2R - 10^{s}) = \operatorname{sign}(R'),\\ r \le 4 &\implies 2R - 10^{s} \le -2\cdot 10^{s-1} + 2R' < 0 . \end{aligned}

したがって $2R > 10^{s}$ 、 $2R = 10^{s}$ 、 $2R < 10^{s}$ はそれぞれ「中点より上、中点ちょうど、中点より下」を意味し、half 系のモードに必要なのはこれだけです。 $R \ne 0$ は方向付き丸めのモードが必要とするスティッキー情報であり、ZeroFiveUp と HalfEven に必要なのは $Q$ の最後の桁だけです。 $s \ge D$ のときは係数全体が捨てられて $Q = 0$ となり、中点に達しうるのは $s = D$ の場合だけです。このときコードは $C$ と $5\cdot 10^{D-1}$ を比較します。 $10^{p}-1$ を $10^{p}$ に変える繰り上がりは、10 による 1 回の厳密な除算と指数のインクリメントで取り除かれます。

rounded は $s > 0$ のとき常に、inexact は $R \ne 0$ のとき常に立てられます。丸めの前に末尾のゼロを落とすので、ゼロを捨てるだけなら rounded だけが立ちます。

除算

2 つの有限でゼロでない値の商 $a/b$ は、次の 3 つの厳密な経路のいずれかを、この順に試します。

除数が 10 のべきの場合。 商は指数をずらした $a$ です。
有限小数になる商。 $c_a / c_b$ を $g = \gcd(c_a, c_b)$ で約分して $c_a'/c_b'$ とします。商が有限の 10 進展開を持つのは $c_b' = 2^{i}5^{j}$ のとき、かつそのときに限ります。 $k = \max(i, j)$ とすると、 $\frac{c_a'}{2^{i}5^{j}} = \frac{c_a'\,2^{k-i}\,5^{k-j}}{10^{k}},$ したがって厳密な結果は指数 $q_a - q_b - k$ の整数 $c_a' 2^{k-i}5^{k-j}$ であり、他の厳密な結果と同様に最終化ルーチンに渡されます。
有限小数にならない商。 結果は必ず厳密ではありません。 $a_c = \operatorname{adj}(c_a/c_b)$ とします。これは $c_a$ と、適切な 10 のべきでスケールした $c_b$ とを比較することで厳密に求まります。分子（または分母）を $10^{p-1-a_c}$ でスケールすると、整数商がちょうど $p$ 桁になります： $N = Q D + R$ 。インクリメントの判定は $2R$ と $D$ を比較します。これは上と同じ中点テストで、今度は厳密な剰余がスティッキー情報になります。したがって商は 1 回だけ丸められます。

3 番目の経路は、拡張コンテキストで結果が正規化数の場合に使われます。結果が非正規化数の場合やサブセットコンテキストでは、コードは $p + \operatorname{digits}(c_b) + 2$ 桁を計算してコンテキストのモードで丸め、さらに精度と $E_{\text{tiny}}$ へもう一度丸めます。コンテキストを持たない演算子 / も、同じガード付きの方式を HalfEven で使います。最近接への丸めを 2 回行うことは常に正しいとは限りません。1 回目の丸めが 2 回目の丸めの中点にちょうど一致すると、厳密な値ではすでに決まっていたケースを 2 回目の丸めのタイ規則が決めてしまうからです。したがってこの方式は、そのような中点付近の場合を除いてすべての場合に正確です。演算子 / では、たとえば 5 桁での $15/83294$ でこれが現れます（0.00018009 ではなく 0.00018008）。結果が正規化数の場合の div_ctx にはこの弱点はありません。

ZeroFiveUp はまさにこのような 2 段階の方式を安全にするために存在します。 $x$ をまず ZeroFiveUp で $p + k$ 桁（ $k \ge 1$ ）に丸め、次に任意のモード $\circ$ で $p$ 桁に丸めると、結果は $\circ(x)$ に等しくなります。厳密でない ZeroFiveUp の結果は 0 と 5 以外の桁で終わるので、 $p$ 桁の数にも $p$ 桁の中点にもならず、すべての $p$ 桁の数と中点に対して $x$ と同じ側にあるからです。証明は添付資料にあります。³3 decimal の丸めに関する証明には、二重丸めの補題、オーバーフローの表、フォールドダウンの限界、認証の補題、NTT の限界の完全な証明が含まれています。

平方根

目標指数を $t = \max(E_{\text{tiny}},\ \lfloor \operatorname{adj}(x)/2 \rfloor - p + 1)$ とし、オペランドを整数 $M = c\,10^{q - 2t}$ にスケールします（負の分岐で $q - 2t$ が奇数なら先に 10 を掛けます）。整数平方根 $r = \lfloor \sqrt{M} \rfloor$ は整数上の Newton 反復で計算します。

a_{k+1} = \left\lfloor \frac{a_k + \lfloor M / a_k \rfloor}{2} \right\rfloor , \qquad a_0 = 10^{\lceil (\operatorname{digits}(M)+1)/2 \rceil} > \sqrt{M},

この反復は $a_k > \lfloor\sqrt M\rfloor$ である間は狭義に減少し（相加相乗平均の不等式 $\tfrac12(a + M/a) \ge \sqrt{M}$ により、また $a_{k+1} < a_k$ となるのは $a_k^2 > M$ のとき、かつそのときに限るため）、 $\lfloor \sqrt M \rfloor$ で停止します。剰余 $M - r^2$ がスティッキー情報です。中点テストは $\sqrt{M} \gtrless r + \tfrac12 \iff 4M \gtrless (2r+1)^{2}$ であり、 $4M$ は偶数で $(2r+1)^2$ は奇数なので等号は起こりえません。平方根がちょうど中点になることはないので、HalfEven、HalfUp、HalfDown は平方根では一致します。非正規化数の目標指数で直接丸めることで、極小な平方根の二重丸めを避けています。厳密な平方根は最初に検出され（指数を偶数にした後で約分した係数が完全平方数である場合）、推奨指数 $\lfloor q/2 \rfloor$ で返されます。

指数範囲

最終化ルーチンは $p$ 桁の丸められた係数 $C'$ と指数 $Q'$ 、すなわち指数範囲を無制限として丸めた値に対して処理を行います。

オーバーフロー

$\operatorname{adj}(C' 10^{Q'}) > e_{\max}$ であれば結果はオーバーフローします。IEEE 754 §7.4 は、返される結果を、同じ $p$ で指数が無制限の形式で厳密な値を丸め、それを飽和させたものと定義しています。絶対値を決して増やさない方向は有限の範囲を出られず、増やしうる方向は無限大になります。したがって $N_{\max} = (10^{p}-1)10^{E_{\text{top}}}$ とすると、

モード	$x > 0$	$x < 0$
`HalfEven`, `HalfUp`, `HalfDown`, `Up`	$+\infty$	$-\infty$
`Down`	$+N_{\max}$	$-N_{\max}$
`Ceiling`	$+\infty$	$-N_{\max}$
`Floor`	$+N_{\max}$	$-\infty$
`ZeroFiveUp`	$+N_{\max}$	$-N_{\max}$

half 系のモードが無限大になるのは、オーバーフローする厳密な値は少なくとも $N_{\max} + \tfrac12 10^{E_{\text{top}}}$ であり（それより小さいものは有限値に丸められてオーバーフローしません）、これは $N_{\max}$ と次の 10 のべきとの中点以上だからです。ZeroFiveUp が飽和するのは、 $N_{\max}$ の最後の桁が 9 だからです。すべてのオーバーフローは overflow、inexact、rounded を立てます。

///|
test "design: overflow depends on the rounding direction" {
  let ctx = @decimal.DecimalContext::decimal64()
  let big = @decimal.Decimal::from_string("9E+384").unwrap()
  let ten = @decimal.Decimal::from_int(10)
  inspect(big.mul_ctx(ten, ctx).0, content="inf")
  let down = ctx.with_rounding(@def.RoundingMode::TowardZero)
  let (sat, flags) = big.mul_ctx(ten, down)
  inspect(sat, content="9.999999999999999E+384")
  inspect(flags.overflow && flags.inexact, content="true")
}

非正規化数、極小性、アンダーフロー

厳密な結果が $E_{\text{tiny}}$ 未満の指数を必要とする場合、それは非正規化数のグリッドに丸められます。シフト量は $s = E_{\text{tiny}} - Q$ となり、保持する桁数が $p$ 未満の状態で上記の丸め規則が適用されます。結果は、その調整指数が $e_{\min}$ 未満のとき極小です。調整指数は厳密な値で測る（BeforeRounding）か、指数無制限で $p$ 桁に丸めた値で測ります（AfterRounding、デフォルト）。2 つの規則が異なるのは、 $10^{e_{\min}}$ のすぐ下にあってそれへ切り上げられる値の場合だけです。極小な結果は subnormal を立て、極小かつ厳密でない結果は、IEEE 754 §7.5 がデフォルトの例外処理として要求するとおり underflow も立てます。ゼロに丸められた結果は指数 $E_{\text{tiny}}$ と clamped を得ます。

クランプ

clamp が有効な場合（すべての交換形式）、エンコーディングには $E_{\text{top}} - E_{\text{tiny}} + 1$ 個の指数を入れる余地しかないため、 $E_{\text{top}}$ を超える指数は表現できません。 $Q' > E_{\text{top}}$ でオーバーフローしなかった結果はフォールドダウンされます。係数に $10^{Q' - E_{\text{top}}}$ を掛け、指数を $E_{\text{top}}$ に設定し、clamped を立てます。これに $p$ 桁を超える桁が必要になることはありません。

\operatorname{digits}(C') + (Q' - E_{\text{top}}) = \bigl(\operatorname{adj} - Q' + 1\bigr) + Q' - (e_{\max} - p + 1) = \operatorname{adj} - e_{\max} + p \le p ,

ここで $\operatorname{adj} \le e_{\max}$ を用いました。値は変わらず、コホートのメンバーだけが変わります。decimal32 では、1E+96 は 1000000E+90 として格納されます。

///|
test "design: fold-down in decimal32" {
  let ctx = @decimal.DecimalContext::decimal32()
  let (x, flags) = @decimal.Decimal::from_string_ctx("1E+96", ctx)
  inspect(x.coefficient(), content="1000000")
  inspect(x.exponent10(), content="90")
  inspect(flags.clamped, content="true")
}

ゼロには埋める桁がないので、その指数は単に $[E_{\text{tiny}}, E_{\text{top}}]$ （clamp なしでは $[E_{\text{tiny}}, e_{\max}]$ ）にクランプされ、変化したときに clamped が立てられます。

戻り値としてのフラグ

問題。 IEEE 754 のステータスフラグは、ハードウェアではスティッキーなプロセス状態です。GDA はさらにトラップを加えます。隠れた状態は、意味論を明示的にするという Luna-Flow の規則と衝突し、並行なコードや合成的なコードを壊れやすくします。

選択。 各コンテキスト演算はそれ自身の DecimalFlags を返します。combine はフィールドごとの OR なので、フラグ集合は単位元 DecimalFlags::new() を持つ可換で冪等なモノイドをなします。パイプライン全体でどのようにグループ化して蓄積しても同じ集合が得られ、これはまさに状態を持たないスティッキーフラグの意味論です。decimal_checked はこの蓄積をパッケージ化し、decimal_gda はスティッキーなステータスとトラップを別のモデルとして実装します。IEEE の 5 つの例外は invalid_operation、division_by_zero、overflow、underflow、inexact に対応し、GDA の条件（rounded、subnormal、clamped、lost_digits、conversion_syntax、division_impossible、division_undefined、invalid_context）はそれらを細分化したものです。

quantize と same-quantum

x.quantize(y) は、指数がちょうど $t = q_y$ である $x$ の値を返します。

\operatorname{quantize}(c\,10^{q},\, t) = \begin{cases} c\,10^{q-t} \cdot 10^{t} & q \ge t \text{ (exact padding)},\\ \circ\!\left(c\,10^{q - t}\right)\cdot 10^{t} & q < t \text{ (rounding)} . \end{cases}

結果はその指数で表現可能でなければなりません。すなわち新しい係数は高々 $p$ 桁で、 $t$ は $[E_{\text{tiny}}, e_{\max}]$ に含まれ、結果の調整指数は $e_{\max}$ を超えてはなりません。そうでなければ演算は無効です。指数こそが契約であるため（セント単位に量子化した金額は小数 2 桁でなければならない）、別の指数で代用することは決してありません。係数を切り上げると桁が増えることがあるため（精度 2 桁で $9.99 \to 10.0$ ）、桁数の検査は丸めの後に行います。same_quantum は述語 $q_x = q_y$ です（2 つの無限大または 2 つの NaN に対しては真）。指数が一致していなければならない値を組み合わせる前に使うべき判定です。

交換エンコーディング

3 つの形式はすべて同じレイアウトを共有します。符号ビット、5 ビットの組み合わせフィールド $G$ 、 $w$ ビットの指数継続ビット、 $10J$ ビットの後続仮数です。ここで $p = 3J + 1$ です。

形式	$w$	$J$	$1 + 5 + w + 10J$
decimal32	6	2	32
decimal64	8	5	64
decimal128	12	11	128

バイアス付き指数 $E = q + \text{bias}$ は $w + 2$ ビットで、その上位 2 ビットは 00、01、10 の値しかとりません。これが上で導いた $3\cdot 2^{w}$ という個数です。

DPD。 組み合わせフィールドは指数の上位 2 ビットと先頭の桁 $d_0$ を保持します。 $G = ab\,cde$ で $ab \ne 11$ なら、指数ビットは $ab$ 、 $d_0 = cde \in [0,7]$ です。 $G = 11\,cd\,e$ で $cd \ne 11$ なら、指数ビットは $cd$ 、 $d_0 = 8 + e$ です。 $G = 11110$ は無限大、 $G = 11111$ は NaN で、次のビットが signaling と quiet を区別します。残りの $3J$ 桁は、densely packed decimal により $J$ 個のデクレット（3 桁を 10 ビットで表す）に格納されます。⁴4 M. F. Cowlishaw, “Densely packed decimal encoding”, IEE Proceedings — Computers and Digital Techniques 149(3), 2002. IEEE 754-2019 §3.5.2 がエンコーディングの表を与えています。コードはそれらをブール式として実装しており、適合性コーパスは 1024 個すべてのデクレットを検査します。 3 桁は 1000 通りの値、10 ビットは 1024 通りの符号を持ち、効率は $\log_2 1000 / 10 = 99.66\%$ です。0–7 の桁（3 ビット）を小さい桁、8 または 9 の桁（1 ビット）を大きい桁と呼びます。デクレット $pqr\,stu\,v\,wxy$ は、3 桁とも小さい場合に $v = 0$ を使い、そのとき桁は $pqr$ 、 $stu$ 、 $wxy$ にそのまま格納されます。 $v = 1$ は少なくとも 1 つの大きい桁があることを示し、 $wx$ 、次いで $st$ がそれがどれかを示します。大きい桁の個数で数えると、

\underbrace{8^3}_{v=0} = 512,\quad \underbrace{3\cdot 2\cdot 8^2}_{v=1,\ wx\ne 11} = 384,\quad \underbrace{3\cdot 2^2\cdot 8}_{wx=11,\ st\ne 11} = 96,\quad \underbrace{2^3}_{wx=11,\ st=11} = 8,

そして $512 + 384 + 96 + 8 = 1000$ です。最初の 3 つの場合はちょうど 512、384、96 個の符号を使います。最後の場合は 8 個の値に対して 32 個の符号があります。 $r$ 、 $u$ 、 $y$ が 3 桁の下位ビットを運び、 $p$ 、 $q$ は無視されるので、24 個の符号は冗長です。3 桁とも大きい各値には 4 つのエンコーディングがあり、そのうち $pq = 00$ のものが正準です。デコードは 4 つすべてを受け付け、canonical() はそれらを書き換えます。たとえば 125 はデクレット 0010100101 = 0x0A5（3 桁とも小さい）であり、999 は 0011111111 = 0x0FF で、0x1FF、0x2FF、0x3FF とも書けます。

BID。 係数は 2 進整数として格納されます。符号の後の 2 ビットが 11 でなければ、続く $w+2$ ビットがバイアス付き指数、残りの $10J + 3$ ビットが係数です。そうでなければ指数は 11 の後に続き、係数は $2^{10J+3}$ に残りの $10J+1$ ビットを加えたものです（“100” が暗黙に補われます）。 $10^{7} - 1 < 2^{24}$ 、 $10^{16}-1 < 2^{54}$ 、 $10^{34}-1 < 2^{114}$ なので、すべての係数が収まります。係数が $\ge 10^{p}$ のエンコーディングは非正準であり、ゼロにデコードされます。

///|
test "design: redundant DPD declets decode and canonicalize" {
  let fmt = @decimal.DecimalInterchangeFormat::Decimal64
  let canonical = @decimal.DecimalInterchange::from_hex("#22300000000004FF", fmt).unwrap()
  let redundant = @decimal.DecimalInterchange::from_hex("#22300000000007FF", fmt).unwrap()
  inspect(canonical.to_decimal(), content="19.99")
  inspect(redundant.to_decimal(), content="19.99")
  inspect(redundant.is_canonical(), content="false")
  inspect(redundant.canonical().to_hex(), content="#22300000000004FF")
}

DecimalInterchange は生のビットを保持するので、非正準な入力は呼び出し側が正準化を決めるまでそのまま残ります。算術は常に Decimal 上で行われます。

認証付き初等関数

問題。 超越関数 $f$ に対して、ほとんどすべての 10 進数 $x$ で $f(x)$ は無理数になるため、近似することしかできません。正しく丸められた結果を得るには、丸めを決定できるほど良い近似が必要です（table maker’s dilemma）。

選択肢。 事前の誤差限界を伴う固定精度の評価（高速だが、正しさはその限界がすべての関数と引数で正しいことに依存する）。厳密な誤差限界を伴う Ziv の適応的戦略。⁵5 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, ch. 12; 包含区間のモデルについては J. van der Hoeven, “Ball arithmetic”, 2009.

選択。 厳密な包含区間に対する Ziv ループ。演算 $f$ と入力 $x$ について：

$x$ を 2 進の区間に厳密に変換します： $\underline{x} = \nabla_w(x)$ 、 $\overline{x} = \Delta_w(x)$ 。これは $w$ ビットでの TowardNegative と TowardPositive による to_bin_float であり、 $x \in [\underline{x}, \overline{x}]$ となります。
ball_float を使ってその区間上で $f$ を評価します。これは入力区間のすべての点について $[L, U] \ni f(x)$ となる区間を返します。
$L$ と $U$ （2 進有理数なので厳密な 10 進数）を厳密に Decimal に変換し、両方を目標コンテキストで丸めます。両者が同じ表現（compare_total で等しい）と同じフラグを与えれば、それを返します。
そうでなければ $w \leftarrow w + \max(32, \lfloor w/2 \rfloor)$ と増やして繰り返します。これを最大 12 回行い、その後は認証の失敗を報告します。

丸めは単調なので、この受理判定は健全です。 $L \le f(x) \le U$ ならば $\circ(L) \le \circ(f(x)) \le \circ(U)$ であり、外側の 2 つが同じ表現なら、中央のものも同じです。フラグも同様に引き継がれます。ただし $f(x)$ 自身が表現可能でない場合に限ります。そのとき一方の端点は共通の結果と異なるので両方とも厳密ではなく、ゼロを含まない区間上では overflow、subnormal、underflow は $|f(x)|$ について単調だからです。したがって表現可能な結果はループの前に処理されます（後述）。

初期作業精度は $w_0 = \max(128,\ 4\max(D, p) + 64)$ ビットです。ここで $D$ は入力の桁数です。10 進 1 桁には $\log_2 10 \approx 3.32 < 4$ ビットが必要なので、 $4\max(D,p)$ ビットで入力と目標を余裕をもって表現でき、さらに 64 ビットで包含区間による損失を賄います。スケジュールは 1 ステップごとにおよそ $3/2$ 倍で増え、12 ステップ後の予算は約 $w_0 \cdot (3/2)^{12} \approx 130\,w_0$ ビットです。

2 種類の入力は一致判定を決して通過しないので、ループの前に決定されます。

厳密な結果。 $f(x)$ が表現可能な 10 進数であれば、方向付き丸めのモードでは $L < f(x) < U$ がいつまでも異なる隣接値に丸められます。コードは厳密なケース（exp(0)、ln(1)、 $\log_{10} 10^{k}$ 、sinpi/cospi/tanpi の整数および半整数の引数、exp2/exp10 の整数の引数、 $x^{1/2}$ 、整数べき、…）を検出します。ball_float は $\log_2 8 = 3$ 、 $\sqrt[3]{8} = 2$ 、 $\operatorname{hypot}(3, 4) = 5$ のような厳密な 2 進の結果に対して点区間を返します。どちらでも検出されない厳密な結果（たとえば power_ctx による $4^{1.5} = 8$ ）は認識されず、精緻化の予算をすべて使い切ります。
範囲外の結果。 $[0, \text{tiny}]$ のような包含区間では、端点が異なるフラグで丸められます。 $\ge 10^{e_{\max}+2}$ の値はすべて同じようにオーバーフローし、 $\le 10^{E_{\text{tiny}}-2} < \tfrac12 10^{E_{\text{tiny}}}$ の非ゼロ値はすべて同じように（モードと符号のみに応じてゼロまたは最小の非正規化数に）丸められるため、遠い端点はそうした代表値に置き換えられます。2 進のテストでは $3.322 > \log_2 10$ を用いるので、この置き換えは保守的です。exp では、 $x \log_{10} e \ge 3 \cdot 0.434\,(e_{\max}+1) > e_{\max}+1$ であるため $x \ge 3(e_{\max}+1)$ はオーバーフローし、同じ評価により $x \le -3\,|E_{\text{tiny}} - 1|$ は $10^{E_{\text{tiny}}-1}$ 未満へアンダーフローします。power は同じ目的で、方向付き丸めの 128 ビット演算により $y \log_2 x$ を評価します。

初等関数は、 $p$ または $|e|$ が 999,999 を超えるコンテキストを拒否します（invalid_context）。これにより端点の厳密な変換の大きさが抑えられます。

係数カーネル

問題。 10 進の位置での丸めには、10 のべきによる高速な除算と高速な桁数計算が必要です。大きな精度には、2 乗未満の計算量の乗算と除算が必要です。

選択。 パッケージプライベートの DecCoeff は、 $10^{9}$ 未満のインラインの UInt か、基数 $10^{9}$ のリムからなるリトルエンディアンの配列（先頭にゼロのリムを持たず、桁数をキャッシュする）のどちらかです。基数 $10^9$ は $2^{30}$ 未満で最大の 10 のべきなので、リム同士の積は $10^{18} < 2^{60}$ 未満になります。9 桁の倍数だけシフトするのはリムの移動で済み、digits10 はリムの個数と最上位リムの桁数から求まります。BigInt が現れるのは公開境界だけです。

乗算のディスパッチ：

形状	アルゴリズム	コスト
それぞれ 1 リム	インライン	$O(1)$
ゼロのリムが多い（ゼロでないリムが $n_a' n_b' \cdot 4 < n_a n_b$ ）	疎な積	$O(n_a' n_b')$
$n_{\text{long}} > 2 n_{\text{short}}$	短い方の長さのブロックに均等分割	$\frac{n_{\text{long}}}{n_{\text{short}}} M(n_{\text{short}})$
小さい	Comba の列方式 / 筆算	$O(n^2)$
Karatsuba の閾値以上	Karatsuba	$O(n^{\log_2 3}) = O(n^{1.585})$
Toom-3 の閾値以上	Toom-3	$O(n^{\log_3 5}) = O(n^{1.465})$
NTT の閾値以上	2 素数 NTT	$O(n \log n)$

Comba カーネルは、 $\min(n_a,n_b) \le 18$ のときに限り、 $n$ 個のリム積の列を UInt64 に累積します。列の和に入ってくる繰り上がりを加えても高々 $18(10^{9}-1)^{2} + 2\cdot 10^{10} < 1.8\cdot 10^{19} < 2^{64}$ ですが、19 個の積では $2^{64} \approx 1.845\cdot 10^{19}$ を超えうるからです。

NTT は係数を基数 $10^4$ の桁に分割し、素数 $p_1 = 998\,244\,353 = 119\cdot 2^{23} + 1$ と $p_2 = 754\,974\,721 = 45 \cdot 2^{24} + 1$ を法として畳み込みます。どちらも $2^{23}$ 乗根の 1 の原始根を持つので、長さ $2^{23}$ までの変換が存在します。畳み込みの係数は $10^{4}$ 未満の桁同士の積を高々 $m = \min(n_a, n_b)$ 個足したものなので $m \cdot 9999^{2}$ 未満であり、中国剰余定理によって剰余から厳密に復元されます。

z = r_1 + p_1 \bigl((r_2 - r_1)\, p_1^{-1} \bmod p_2\bigr),

これは $m \cdot 9999^{2} < p_1 p_2 \approx 7.54 \cdot 10^{17}$ 、すなわち $m < 7.5 \cdot 10^{9}$ である限り成り立ち、変換の上限未満では常に真です。基数 $10^{9}$ の桁では $m\cdot 10^{18} < p_1p_2$ が必要になり、これはどの $m \ge 1$ でも不可能です。NTT がより小さい桁を使うのはこのためです。限界が成り立たない場合、カーネルは Toom-3 にフォールバックします。

除算は、1 リムによる除算、Knuth の Algorithm D、⁶6 D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., §4.3.1 (Algorithm D) and §4.3.3; R. Brent and P. Zimmermann, Modern Computer Arithmetic, Cambridge 2010, §1.3–1.4 and §2.4 (NTT); C. Burnikel and J. Ziegler, “Fast recursive division”, MPI-I-98-1-022, 1998. Burnikel–Ziegler の再帰的除算、または Newton の逆数除算を使います。Newton 法は $S = B^{n+1}$ に対して $r_{k+1} = \lfloor r_k (2S - d r_k)/S \rfloor$ により $r \approx S/d$ を計算します。 $r_k = S/d - \varepsilon_k$ とすると $S/d - r_{k+1} \approx d\,\varepsilon_k^{2}/S$ となるので、正しいリムの数は 1 ステップごとに倍になります。反復は下から単調に増加しなければなりません。そうならない場合、あるいは最後の商の補正に $2n+8$ ステップより多くかかる場合、ルーチンは Burnikel–Ziegler にフォールバックし、それ自体もサポートしない形状では Algorithm D にフォールバックします。

切り替え点はターゲットごとに計測され、ターゲット固有のファイルに保存されています。

ターゲット	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

（リムは 9 桁）。native では NTT の閾値は変換長にも依存します。乗算では 1,728、2,816、4,608、7,680、その後 8,192 リム、2 乗では 640、1,040、1,824、3,648、7,296、その後 8,192 リムです。また Burnikel–Ziegler に入る点は、ブロック長が大きくなると 5,120 リム、10,240 リムに移ります。これらはディスパッチの境界であり、コストを変えるだけで結果を変えることはありません。native の Newton 経路は実装・テスト済みですが、native での計測では切り替え点が見られないため無効になっています。

正しさ／不変条件

表現。 有限の Decimal は $c \ge 0$ を持ちます。DecCoeff のリムは正準です（先頭にゼロのリムがなく、桁数が正確）。ゼロの符号は符号ビットに保持され、coefficient() が符号を持つことはありません。
1 回だけの丸め。 +、-、×、fma、sqrt、quantize、変換、および（結果が正規化数の場合の）/ のすべてのコンテキスト結果は、厳密な結果を 1 回だけ丸めたものです。初等関数の結果は返される場合は常に正しく丸められており、失敗は近似されることなく報告されます。
誤差限界。 これらの演算の結果が正規化数であれば、 $\operatorname{fl}(x) = x(1+\delta)$ であり、half 系のモードでは $|\delta| \le \tfrac12 10^{1-p}$ 、方向付き丸めのモードでは $|\delta| < 10^{1-p}$ です。
厳密さが見える。 inexact は返される値が厳密な結果と異なるとき、かつそのときに限り立てられます。rounded は桁が落とされたときに常に立てられます。
コホートの保存。 収まる厳密な結果は推奨指数で返されます。quantize は指数 $q_y$ を返すか、失敗するかのどちらかです。
フラグ。 combine は結合的、可換、冪等で、単位元は new() です。
順序。 compare は全前順序です（NaN 同士は等しく、すべての数より大きい。 $-0 = +0$ ）。compare_total は表現上の全順序であり、NaN でない値については compare を細分化します。
エンコーディング。 正準なビットをデコードしてからエンコードすると恒等写像になります。形式に収まる値をエンコードしてからデコードすると、コホート、ゼロの符号、NaN のペイロード（DPD）を含めて恒等写像になります。
計算量。 比較、加算、シフト、1 リムによる除算はリム数について $O(n)$ です。乗算と除算はディスパッチ表に従います。

中点テスト、ZeroFiveUp の二重丸めの補題、オーバーフローの表、フォールドダウンの限界、認証の補題、カーネルの限界の証明は、添付資料にまとめられています。

decimal の丸めに関する証明

適合性のページには、有限の証拠（固定された IEEE コーパス、デクレットの網羅的検査、MPFR で認証された初等関数の行、4 つのターゲット）が記録されています。

却下した代替案

2 進の BigInt 係数。 10 進の位置での丸めには、毎回の演算で $10^{s}$ による除算と桁数の計算が必要です。2 進の係数ではどちらも高価ですが、基数 $10^{9}$ のリムならリムの移動と 1 回の小さな除算で済みます。
暗黙のコンテキストとスティッキーフラグ。 明示性と合成可能性のために不採用としました。フラグのモノイドが同じ情報を与えます。
NaN で中断する compare。 NaN を含むデータに対して、ソートや汎用の Compare コードが中断してしまいます。現在の compare は全前順序であり、IEEE の意味論は compare_checked、compare_ctx、compare_signal_ctx、compare_total を通じて利用できます。
解析的な誤差限界を持つ固定精度の超越関数カーネル。 小さな精度では高速ですが、関数ごと・精度ごとに正しさの証明が必要になります。包含区間のループは構成上正しく、その失敗モードは明示的です。
IEEE と GDA を 1 つのパッケージにすること。 スティッキーなステータス、トラップ、トラップの優先順位はすべての演算の型を変えてしまいます。これらは decimal_gda にあります。
すべての結果を正規化すること。 コホートを失うと 12.30 と 12.3 が区別できなくなり、量子に依存するプロトコルが壊れます。

境界

decimal は意図的に次のことを行いません。

スティッキーなステータスやトラップを保持すること（decimal_checked または decimal_gda を使ってください）。
コンテキストを持たない演算子をコンテキストへ丸めること。* は厳密で、+ と / はオペランドの精度にのみ丸め、いずれも指数範囲を適用しません。
初等関数の認証が成功することを保証すること。精緻化の予算（最後のステップでは非常に長い数を扱い、長い時間がかかることがあります）を使い切ると、CertificationFailure（try_*_ctx）または invalid_operation を伴う NaN（*_ctx）を報告します。
精度または指数が 999,999 を超えるコンテキストで初等関数を評価すること。
2 進との変換を通して NaN のペイロードを保持すること、あるいは値の精度が形式の精度と異なる BID の NaN ペイロードを保持すること。
係数の表現、カーネルの選択、閾値を公開すること。
適合性のページにある有限の証拠を超えて適合性を主張すること。

IEEE Std 754-2019, Standard for Floating-Point Arithmetic、3.3–3.5 節（10 進形式とエンコーディング）、4 節（属性と丸め）、5 節（演算）、7 節（例外）、9 節（推奨演算）。M. F. Cowlishaw による General Decimal Arithmetic 仕様（バージョン 1.70）は同じモデルを任意精度の形で与えています。本ページではその用語 coefficient（係数）、adjusted exponent（調整指数）、Etiny、clamp を使います。 ↩
Goldberg, “What every computer scientist should know about floating-point arithmetic”, ACM Computing Surveys 23(1), 1991, §1.2; Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §2.1–2.2. ↩
decimal の丸めに関する証明には、二重丸めの補題、オーバーフローの表、フォールドダウンの限界、認証の補題、NTT の限界の完全な証明が含まれています。 ↩
M. F. Cowlishaw, “Densely packed decimal encoding”, IEE Proceedings — Computers and Digital Techniques 149(3), 2002. IEEE 754-2019 §3.5.2 がエンコーディングの表を与えています。コードはそれらをブール式として実装しており、適合性コーパスは 1024 個すべてのデクレットを検査します。 ↩
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, ch. 12; 包含区間のモデルについては J. van der Hoeven, “Ball arithmetic”, 2009. ↩
D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., §4.3.1 (Algorithm D) and §4.3.3; R. Brent and P. Zimmermann, Modern Computer Arithmetic, Cambridge 2010, §1.3–1.4 and §2.4 (NTT); C. Burnikel and J. Ziegler, “Fast recursive division”, MPI-I-98-1-022, 1998. ↩