internal の設計

設計目標

2 進・10 進・区間の各コアは、いずれも難しいステップを正確な整数演算に帰着させます。指数の整列、末尾のゼロの除去、丸めモード付きの除算、2 つの 2 進有理数（dyadic）の間への比の包含です。internal はこれらのステップを、不変条件を明記した小さな関数として一度だけ実装します。これにより、すべてのコアが同じ規則で丸め、consistency テストで各ヘルパーを BigInt オラクルに対して検査できます。internal であるため、コアは公開互換性の約束なしにヘルパーをまとめて変更できます。

数学的背景

剰余付き整数除算

$n \ge 0$ と $d > 0$ に対し、 $n = qd + r$ かつ $0 \le r < d$ を満たす整数 $q, r$ が一意に存在します（ユークリッド除算）。このとき $q = \lfloor n/d \rfloor$ であり、 $n/d$ が整数であるのはちょうど $r = 0$ のときです。小数部分は $r/d \in [0, 1)$ であり、

\frac{r}{d} \lesseqgtr \frac12 \iff 2r \lesseqgtr d .

整数への丸め

実数 $x$ と丸め方向に対し、 $\circ(x)$ は $\lfloor x \rfloor$ と $\lceil x \rceil$ の中から選ばれる整数です。ゼロ方向、 $+\infty$ 方向、 $-\infty$ 方向、ゼロから離れる方向、または最近接（タイは偶数の候補へ）のいずれかです。¹1 IEEE 754-2019 の第 4.3 節は丸め方向属性を定義しています。整数の場合は、表現可能な数の集合を $\mathbb{Z}$ とした同じ定義です。

2 進有理数による包含区間

2 進有理数（dyadic number）とは、 $n \in \mathbb{Z}$ 、 $s \ge 0$ に対する $n \cdot 2^{-s}$ です。実数 $x$ とスケール $s$ に対し、 $x$ を囲むスケール $s$ の 2 進有理数は $\lfloor x 2^{s} \rfloor 2^{-s} \le x \le \lceil x 2^{s} \rceil 2^{-s}$ であり、幅が高々 $2^{-s}$ の区間をなします。

設計上の判断

絶対値と符号を別々に渡す

round_positive_div(n, d, negative, mode) は $n \ge 0$ と、真の商の符号をフラグとして受け取ります。符号と絶対値による表現はコアが係数を格納する方法であり、絶対値の方向付き丸めは符号に依存します。 $-2.5$ を $-\infty$ 方向に丸めると絶対値は増加します。符号を別の引数として保持することで、規則を 1 つの表にまとめられ、言語によって規約が異なる負の剰余を避けられます。

共有キャッシュによる 10 と 5 のべき乗

10 進スケーリングでは $10^{k}$ を、2 進–10 進変換では $5^{k}$ を使います（ $10^{k} = 5^{k} 2^{k}$ であり、因子 $2^{k}$ はシフトだからです）。キャッシュは 64 ビットに収まる 19 個のべき乗から始まり、必要に応じて $k = 4096$ まで拡張されます。それより大きなべき乗は呼び出しのたびに直接計算されるため、極端な指数によってキャッシュが際限なく大きくなることはありません。

文字列を使わない桁数の計算

digits10 は、 $b$ を $|x|$ のビット長として、推定値 $d_0 = \lfloor b \log_{10} 2 \rfloor + 1$ から始め、10 のべき乗との比較によってこれを補正します。 $2^{b-1} \le |x| < 2^{b}$ であるため、真の桁数 $d$ は次を満たします。

\lfloor (b-1) \log_{10} 2 \rfloor + 1 \;\le\; d \;\le\; \lfloor b \log_{10} 2 \rfloor + 1 = d_0 ,

そして $\log_{10} 2 < 1$ なので、2 つの境界の差は高々 1 です。したがって下方向の補正は高々 1 回しか必要なく（上方向のループは推定値の浮動小数点誤差に対する保護です）、コストは定数回の BigInt 比較と 1 つの 10 のべき乗です。

飽和する指数の解析

split_decimal_string は、桁を読み取る間に、記述された指数を $\pm 1\,500\,000\,000$ で打ち切ります。それを超える指数はサポートされるどの指数範囲からも大きく外れているため、値はすでにオーバーフローまたはアンダーフローしています。打ち切ることで、観測可能な結果を何も変えずに演算を Int の範囲に保てます。

正準な有理数

ExactRat::new は $\gcd(n, d)$ で割り、分母を正にします。正準形があれば、導出された構造的等価性が有理数としての等価性そのものになります。これは、semantic が 2 進表現と 10 進表現をまたいで値を比較するために必要なことです。

認証付き初等関数は、作業精度 $p_k$ で包含区間を評価し、両端が同じ目標の数に丸められる場合に結果を受け入れます（Ziv の戦略²2 A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991. ）。そうでなければ、より高い精度で再試行します。予算は次のスケジュールを定めます。

p_{k+1} = p_k + \max\bigl(32, \lfloor p_k / 2 \rfloor\bigr), \qquad k < L ,

既定では精密化の回数は $L = 12$ です。 $p_k \ge 64$ ではこれは $p_{k+1} = \lfloor 3 p_k / 2 \rfloor$ であり、比 $3/2$ の幾何級数的な増加なので、 $p_L \approx p_0 (3/2)^{L}$ となります（ $L = 12$ では約 $130\,p_0$ ）。精度 $p$ での 1 回の試行のコストが $C(p) \ge c\,p$ （少なくとも線形）であれば、全試行の総コストは最後の試行によって支配されます。

\sum_{k=0}^{m} C(p_k) \le C(p_m) \sum_{j \ge 0} (2/3)^{j} = 3\, C(p_m) \quad \text{when } C(p_k) \le (2/3)^{m-k} C(p_m),

これは、幾何級数的な段階において $\alpha \ge 1$ の $C(p) = c\,p^{\alpha}$ に対し、スケジュール中の床関数の分を除いて成り立ちます。予算の枯渇は、認証されていない値を返す代わりに、目標精度と最終的な作業精度を記録する certified_failure によって報告されます。

正しさ／不変条件

丸め表。 $n = qd + r$ 、 $0 \le r < d$ 、 $x = (-1)^{\sigma} n/d$ とします。このとき $|\circ(x)| \in \{q, q+1\}$ であり、 $q + 1$ となるのは $r > 0$ の場合に限られます。round_positive_div は $|\circ(x)|$ を返します。

\begin{aligned} \text{toward zero:} &\quad |\circ(x)| = q, \\ \text{toward } +\infty: &\quad |\circ(x)| = q + [r > 0 \wedge \sigma = 0], \\ \text{toward } -\infty: &\quad |\circ(x)| = q + [r > 0 \wedge \sigma = 1], \\ \text{away from zero:} &\quad |\circ(x)| = q + [r > 0], \\ \text{nearest even:} &\quad |\circ(x)| = q + [2r > d \vee (2r = d \wedge q \text{ odd})] . \end{aligned}

証明。 $|x| = q + r/d$ は $[q, q+1)$ にあります。ゼロ方向は小さい方の絶対値をとります。 $+\infty$ 方向は、正の非整数では大きい方の絶対値を、負の非整数では小さい方をとります。 $-\infty$ 方向はその鏡像です。ゼロから離れる方向は、任意の非整数で大きい方の絶対値をとります。最近接では、 $q$ までの距離は $r/d$ 、 $q+1$ までの距離は $1 - r/d$ なので、 $q + 1$ の方が近いのは $2r > d$ のときかつそのときに限り、 $2r = d$ がタイです。タイは $q, q+1$ のうち偶数の方に解決され、それが $q + 1$ になるのは $q$ が奇数のときかつそのときに限ります。 $\square$

round_shift(m, s, …) は $d = 2^{s}$ 、 $q = m \gg s$ 、 $r = m - (q \ll s)$ の場合であり、同じ表が適用されます。consistency テストは、タイと方向付きの場合について、両方の関数をこれらの式に対して検査します。

因子の除去。 remove_factor2(sig, e) は $t = \operatorname{ctz}(|sig|)$ として $(sig / 2^{t}, e + t)$ を返すので、 $(sig / 2^{t}) \cdot 2^{e+t} = sig \cdot 2^{e}$ であり、新しい仮数は奇数になります。remove_factor10 と trim_trailing_decimal_zeros も同様に、1 ステップごとに 10 の因子を 1 つずつ除いて $c \cdot 10^{e}$ を保存します。後者は max_drop ステップで停止します。

包含。 certified_dyadic_fraction(n, d, s) は $\ell = \lfloor n 2^{s} / d \rfloor 2^{-s}$ と $u = \lceil n 2^{s} / d \rceil 2^{-s}$ を返します。 $y = n 2^{s}/d$ に対する $\lfloor y \rfloor \le y \le \lceil y \rceil$ から $\ell \le n/d \le u$ と $u - \ell \le 2^{-s}$ が得られ、 $\ell = u$ となるのは $y \in \mathbb{Z}$ のときかつそのときに限ります。負の比の床関数は $-\lceil |n| / d \rceil$ として計算されるため、包含はどちらの符号でも正しくなります。certified_dyadic_div は符号を分子に移して $a / b = (n_a 2^{s_b}) / (n_b 2^{s_a})$ と書き換えるので、同じ限界を継承します。round_down と round_up はより小さいスケールでの同じ床関数と天井関数なので、 $\text{round\_down}(x, s) \le x \le \text{round\_up}(x, s)$ が成り立ちます。

予算。 精度の列は狭義単調増加であり（各ステップで少なくとも 32 ビット増える）、精密化は高々 limit 回しか行われないため、available() を検査するすべての精密化ループは停止します。

中断の契約。 負の指数での pow2、pow5、pow10、 $n < 0$ または $d \le 0$ での round_positive_div、 $d = 0$ での ExactRat::new、および負のスケールでの 2 進有理数のコンストラクタは中断（abort）します。これらはコア内部のプログラマのエラーであり、ユーザー入力から到達することはありません。

却下した代替案

浮動小数点による丸め。 比を丸めるために Double に変換すると、53 ビットを超える正確性が失われます。すべてのヘルパーは BigInt の範囲にとどまります。
符号付きの剰余付き除算。 切り捨て除算と床除算は負のオペランドで異なります。符号と絶対値の表現はこの曖昧さを避けます。
正規化された CertifiedDyadic。 演算のたびに正規化すると末尾ゼロの走査のコストがかかります。包含区間は compare で比較され、これは正準形を必要としません。
非有界なべき乗キャッシュ。 病的な指数によって、巨大な BigInt 値がプロセスの存続期間中保持されてしまいます。

境界

浮動小数点形式、コンテキスト、フラグはありません。コアがこれらのヘルパーの上にそれらを構築します。
split_decimal_string には 10 進文字列の整形も特殊値（inf、nan）もありません。
公開の安定性はありません。このパッケージは Luna-Flow/floating の内部からのみインポートできます。
べき乗キャッシュはプロセス全体の可変状態であり、パッケージ内の唯一の状態です。これが結果を変えることはありません。

IEEE 754-2019 の第 4.3 節は丸め方向属性を定義しています。整数の場合は、表現可能な数の集合を $\mathbb{Z}$ とした同じ定義です。 ↩
A. Ziv, “Fast evaluation of elementary mathematical functions with correctly rounded last bit”, ACM TOMS 17(3), 1991. ↩