Luna Flow
日本語

dzmingli_vs_floating の設計

設計目標

このパッケージは DzmingLi/decimal@0.2.2 と Luna-Flow/floating/decimal_gda@0.7.1 について二つの問いに答えます。同じ厳密な十進結果を計算するか、そして係数が数千桁に伸びたときどれだけ速く計算するか。速度は両方の答えが正しい場合にだけ報告するので、設計の中心には、ライブラリ同士が一致しているだけでは騙されない正しさの検査があります。

API ページが各項目を挙げ、チュートリアルがそれらを実行します。計測値は性能の章にあります。

数学的背景

十進値

有限の十進数とは組 (c,s)∈Z×Z(c, s) \in \mathbb{Z} \times \mathbb{Z} で、次を表します。

v(c,s)=c⋅10−s.v(c, s) = c \cdot 10^{-s}.

写像 vv は単射ではありません。(c,s)(c, s) と (10c,s+1)(10c, s + 1) は同じ数を表します。GDA 算術はこうした組を区別します（指数は結果の一部です）が、ベンチマークが比べるのは数です。そのため各同値類の正準な代表が必要で、このパッケージは s≥0s \ge 0 かつ cc の末尾の 0 が最も少ない組を使います。

c≠0c \ne 0 について ∣c∣|c| の十進桁数を d(c)d(c) と書くと、

10d(c)−1≤∣c∣<10d(c).10^{d(c) - 1} \le |c| < 10^{d(c)} .

working_precision は dd を ∣c∣|c| の十進文字列の長さとして計算するので、d(0)=1d(0) = 1 です。

精度と 0 方向への丸め

精度 pp の GDA コンテキストは、∣c∣<10p|c| < 10^{p} かつ ee が指数範囲内にある数 c⋅10ec \cdot 10^{e} をちょうど表現します。演算は厳密な結果 xx を計算し、それをこうした数に丸めます。丸めモード Down（0 方向）では結果は次のとおりです。

round⁡p(x)=sgn⁡(x) ⌊∣x∣⋅10 p−E(x)⌋⋅10 E(x)−p,E(x)=⌊log⁡10∣x∣⌋+1,\operatorname{round}_p(x) = \operatorname{sgn}(x)\, \bigl\lfloor |x| \cdot 10^{\,p - E(x)} \bigr\rfloor \cdot 10^{\,E(x) - p}, \qquad E(x) = \lfloor \log_{10} |x| \rfloor + 1 ,

したがって ∣round⁡p(x)∣≤∣x∣|\operatorname{round}_p(x)| \le |x| で、等号は xx の有効桁数が pp 以下のときに限り成り立ちます。小数部を固定の kk 桁で同様に切り捨てる操作を次のように書きます。

trunc⁡k(x)=sgn⁡(x) ⌊∣x∣⋅10k⌋⋅10−k.\operatorname{trunc}_k(x) = \operatorname{sgn}(x)\,\bigl\lfloor |x| \cdot 10^{k} \bigr\rfloor \cdot 10^{-k} .

参照オラクルを用いた差分テスト

I1,I2I_1, I_2 を検査対象の実装、OO をオラクル、κ\kappa を等価が判定可能な集合への正準化写像とします。入力 xx に対する実装 kk の判定は次のとおりです。

valid⁡k(x)  ⟺  κ(Ik(x))=κ(O(x)).\operatorname{valid}_k(x) \iff \kappa(I_k(x)) = \kappa(O(x)) .

ペアごとの差分テストは κ(I1(x))=κ(I2(x))\kappa(I_1(x)) = \kappa(I_2(x)) だけを検査します。この検査は両実装が同じ誤った結果を返す共通モードの誤りを見逃し、失敗してもどちらが間違っているかを示しません。オラクルがあれば各実装を個別に判定でき、二つの不一致も判定結果で説明できます。

設計上の決定

三つの実装、計測するのはそのうち二つ

課題。 ベンチマークは「速い」と「速いが間違っている」を区別しなければなりません。選択肢。 ライブラリ間のペアごとの一致、一方のライブラリを他方の参照とする、独立したオラクル。決定。 計時の外で全データセットに対して実行する、BigInt 上の独立した厳密なオラクル（ValidationCoverage::EveryDataset）。理由。 オラクルは一行ずつ確認できるほど小さく、どちらのライブラリともコードを共有せず、4,096 桁以降の DzmingLi の 108 件の失敗を見つけました。ペアごとの一致では「ライブラリが食い違う」としか報告できなかったはずです。あるサイズが速度比のプロットに入るのは、そのサイズで両実装のすべての検証が通った場合だけです。

中立な表現と許容誤差ゼロ

課題。 二つのライブラリは結果の表記（指数表記、末尾の 0、0 の符号）が異なり、近似的な比較では本物の誤りが隠れてしまいます。決定。 すべての結果を DecimalValue に解析し、canonical_string の等価、つまり許容誤差ちょうど 0 で比較します。理由。 後述の精度契約のもとでは、ベンチマーク対象のすべての演算の結果が厳密に表現できるので、正しい実装はそれを一桁ずつ再現しなければなりません。健全性の議論は正しさと不変条件にあります。代償として、比較は指数のコホートとステータスフラグを無視します。それらは公式の decTest 監査（tools/run_dzmingli_dectest_audit.sh）がカバーします。

厳密有限オラクル

オラクルは係数に対する整数演算だけを使います。

加算と減算は両オペランドを s=max⁡(sℓ,sr)s = \max(s_\ell, s_r) に揃えます：

cℓ10−sℓ±cr10−sr=(cℓ10s−sℓ±cr10s−sr) 10−s.c_\ell 10^{-s_\ell} \pm c_r 10^{-s_r} = \bigl(c_\ell 10^{s - s_\ell} \pm c_r 10^{s - s_r}\bigr)\, 10^{-s} .

乗算は係数を掛け、スケールを足します。除算は商を整数の分数として書き、約分します。

cℓ10−sℓcr10−sr=NM,N=cℓ10sr, M=cr10sℓ,N′=Ng, M′=∣M∣g, g=gcd⁡(N,M),\frac{c_\ell 10^{-s_\ell}}{c_r 10^{-s_r}} = \frac{N}{M},\quad N = c_\ell 10^{s_r},\ M = c_r 10^{s_\ell},\quad N' = \frac{N}{g},\ M' = \frac{|M|}{g},\ g = \gcd(N, M),

そして M′M' が 22 と 55 以外の素因数を持たない場合にだけ受け入れます。次に M′∣N′10kM' \mid N' 10^{k} を満たす最小の kk を求め、(N′10k/M′, k)(N' 10^{k} / M',\ k) を返します。整数除算、剰余、べき乗、完全平方数の平方根、FMA、単項演算はこれらの規則から導かれます。一覧はAPI ページの表にあります。

切り捨て系の規則（DivideInteger、Quantize、Rescale、ToIntegralExact、ToIntegralValue）が 0 方向に切り捨てるのは、BigInt の除算がそうだからです。GDA の quantize、rescale、to_integral_* はコンテキストのモードで丸めるので、両方のコンテキストで Down を使います。他のモードではこれらの演算でオラクルが誤ります。なお公開されたフィクスチャでは、これらの演算はいずれにせよ厳密です（次の決定を参照）。

精度契約

課題。 すべての演算は厳密な結果を保持できる精度 pp で実行しなければなりません。そうすれば丸めは決して起きず、許容誤差 0 も公平になります。大きな pp はただではありません。GDA の除算の仕事量は精度とともに増えうるので、過大な pp は結果に不要な仕事まで計時しかねません。決定。 pp はフィクスチャごとにオペランドから計算し、厳密な結果の桁数に対する最小の単純な上界にガード桁を足したもの、すなわち working_precision(op, ℓ, r) + 2 とし、Power と Fma は上書きします。理由。 以下の導出は、この上界がランナーの生成するすべてのオペランドクラスで成り立ち、さらに各オペランドも厳密に保持できるので、pp での解析が丸めを起こさないことを示します。

オペランドの桁数を dℓ,dr,dtd_\ell, d_r, d_t、m=max⁡(dℓ,dr)+∣sℓ−sr∣m = \max(d_\ell, d_r) + |s_\ell - s_r| と書きます。ランナーが使うオペランドは次のとおりです。Add、Subtract、Multiply、Fma、Quantize、Compare はスケールがプロファイル (0,2)(0, 2)、(8,18)(8, 18)、(28,0)(28, 0) から来る生成オペランドを取ります。Divide、DivideInteger、Remainder はスケール 00、88、2828 の生成オペランドを 22、88、2525 で割ります。Power は二乗（指数 22）、SquareRoot は生成した整数の二乗を取ります。Quantize は量子 10−sℓ10^{-s_\ell}、Rescale は指数 −sℓ-s_\ell、ScaleB はシフト 00 を使い、残りの単項演算は整数を取ります。

Add/Subtract:∣cℓ10s−sℓ±cr10s−sr∣<10m+10m≤10m+1⇒ ≤m+1 digits, p=m+4Multiply:∣cℓcr∣<10dℓ+dr⇒ ≤dℓ+dr, p=dℓ+dr+4Power (n=2):∣cℓ2∣<102dℓ⇒ ≤2dℓ, p=2dℓ+4Fma:∣cℓcr10σ−sℓ−sr+ct10σ−st∣<10m′+1, m′=max⁡(dℓ+dr,dt)+∣sℓ+sr−st∣⇒ p=dℓ+dr+dt+∣sℓ+sr−st∣+4≥m′+1Divide (b∈{2,8,25}):12=5⋅10−1, 18=125⋅10−3, 125=4⋅10−2⇒ ≤dℓ+3, p=dℓ+dr+4≥dℓ+5DivideInteger:∣trunc⁡(a/b)∣≤∣a∣<10dℓ−sℓ⇒ ≤dℓRemainder:∣r∣≤∣a∣, scale(r)=sℓ⇒ ≤dℓSquareRoot:cℓ=g2, d(g)≤⌈dℓ/2⌉⇒ p=dℓ+4Unary group:∣result coefficient∣≤∣cℓ∣⇒ ≤dℓ, p=dℓ+4Compare:result∈{−1,0,1}⇒ p=max⁡(dℓ,dr)+3\begin{aligned} \textbf{Add/Subtract:}\quad & |c_\ell 10^{s-s_\ell} \pm c_r 10^{s-s_r}| < 10^{m} + 10^{m} \le 10^{m+1} &&\Rightarrow\ \le m + 1 \text{ digits},\ p = m + 4 \\ \textbf{Multiply:}\quad & |c_\ell c_r| < 10^{d_\ell + d_r} &&\Rightarrow\ \le d_\ell + d_r,\ p = d_\ell + d_r + 4 \\ \textbf{Power}\ (n = 2):\quad & |c_\ell^2| < 10^{2 d_\ell} &&\Rightarrow\ \le 2 d_\ell,\ p = 2 d_\ell + 4 \\ \textbf{Fma:}\quad & |c_\ell c_r 10^{\sigma - s_\ell - s_r} + c_t 10^{\sigma - s_t}| < 10^{m'+1},\ m' = \max(d_\ell + d_r, d_t) + |s_\ell + s_r - s_t| &&\Rightarrow\ p = d_\ell + d_r + d_t + |s_\ell + s_r - s_t| + 4 \ge m' + 1 \\ \textbf{Divide}\ (b \in \{2, 8, 25\}):\quad & \tfrac{1}{2} = 5 \cdot 10^{-1},\ \tfrac{1}{8} = 125 \cdot 10^{-3},\ \tfrac{1}{25} = 4 \cdot 10^{-2} &&\Rightarrow\ \le d_\ell + 3,\ p = d_\ell + d_r + 4 \ge d_\ell + 5 \\ \textbf{DivideInteger:}\quad & |\operatorname{trunc}(a/b)| \le |a| < 10^{d_\ell - s_\ell} &&\Rightarrow\ \le d_\ell \\ \textbf{Remainder:}\quad & |r| \le |a|,\ \text{scale}(r) = s_\ell &&\Rightarrow\ \le d_\ell \\ \textbf{SquareRoot:}\quad & c_\ell = g^2,\ d(g) \le \lceil d_\ell / 2 \rceil &&\Rightarrow\ p = d_\ell + 4 \\ \textbf{Unary group:}\quad & |\text{result coefficient}| \le |c_\ell| &&\Rightarrow\ \le d_\ell,\ p = d_\ell + 4 \\ \textbf{Compare:}\quad & \text{result} \in \{-1, 0, 1\} &&\Rightarrow\ p = \max(d_\ell, d_r) + 3 \end{aligned}

表のどの pp も max⁡(dℓ,dr,dt)\max(d_\ell, d_r, d_t) 以上なので、オペランド自体は厳密に解析されます。Power と Fma を上書きするのは、working_precision だけでは dℓ+dr+2d_\ell + d_r + 2、指数 22 では dℓ+3d_\ell + 3 となり、二乗を保持するには小さすぎるからです。

この契約が証明されているのはこれらのオペランドクラスについてであり、任意の入力についてではありません。一般の有限小数を与える除数 q=2αq = 2^{\alpha} では商の係数は cℓ⋅5αc_\ell \cdot 5^{\alpha} で約 dℓ+0.699αd_\ell + 0.699\alpha 桁ですが、上界は dr≈0.301α+1d_r \approx 0.301\alpha + 1 しか増えません。α=13\alpha = 13 で上界は破れます。1/8192=0.00012207031251 / 8192 = 0.0001220703125 は 1010 桁必要ですが、p=1+4+4=9p = 1 + 4 + 4 = 9 です。このようなフィクスチャが黙って受け入れられることはなく、次節で示すとおりオラクルとの比較で拒否されます。

計時範囲

OperationOnly（arithmetic_only）は計時前に解析したオペランドに対する公開演算一回を計時します。FullPath（full_path）は加えて準備済みのコンテキストで正準オペランド文字列を解析します。これはテキストを受け取る呼び出し側が払うコストです。両範囲は同じフィンガープリントの同じデータセットを実行するので、数値は同じ入力を表します。コンテキストの構築、正準化、検証、レポートはどちらの範囲にも含まれません。

対応のある統計

Mare Mark はデータセット jj、繰り返し rr、ブロック bb ごとに、実装ごとの較正済みレイテンシを一つ記録します。パッケージは (j,r,b)(j, r, b) を共有する DzmingLi と GDA のサンプルを対にし、差をとります。

Δi=tiGDA−tiDZ.\Delta_i = t^{\mathrm{GDA}}_i - t^{\mathrm{DZ}}_i .

サンプルを t=μimpl+βb+εt = \mu_{\mathrm{impl}} + \beta_b + \varepsilon（βb\beta_b はブロックに共通のドリフト。周波数の変化やキャッシュ状態など）とモデル化すると、Δi=μGDA−μDZ+(ε−ε′)\Delta_i = \mu_{\mathrm{GDA}} - \mu_{\mathrm{DZ}} + (\varepsilon - \varepsilon') となり、ブロック効果は打ち消されます。BalancedBlocks は先に走る実装を交互に入れ替えるので、順序の効果も平均すれば打ち消されます。報告される量は次のとおりです。

δ=100⋅median⁡iΔimedian⁡itiDZ %,speedup=median⁡itiGDAmedian⁡itiDZ,\delta = 100 \cdot \frac{\operatorname{median}_i \Delta_i}{\operatorname{median}_i t^{\mathrm{DZ}}_i}\ \%, \qquad \text{speedup} = \frac{\operatorname{median}_i t^{\mathrm{GDA}}_i}{\operatorname{median}_i t^{\mathrm{DZ}}_i},

判定は δ≤−3\delta \le -3 なら gda_faster、δ≥3\delta \ge 3 なら dzmingli_faster、それ以外は equivalent です。中央値の破綻点は 50 % なので、Mare Mark が報告はするが除去しない外れ値によって大きく動くことはありません。3 %3\,\% のしきい値は実用上の有意性の目安であって仮説検定ではなく、レポートには信頼区間はありません。11 Mare Mark の compare_paired は Δi\Delta_i の四分位範囲を区間として保存します。これは対応差のばらつきを表すもので、中央値の不確かさではありません。

サイズごとに 3 データセット、それぞれ確認用の繰り返し 20 回なので、有効なサイズには 60 組があります。

演算ごとのサイズ上限

DzmingLi の digit_count は bit_length * 30103 を 32 ビットの Int で評価します。積は次の条件でオーバーフローします。

bit_length>231−130103≈71 337⟺digits≳71 337⋅log⁡102≈21 475.\text{bit\_length} > \frac{2^{31} - 1}{30103} \approx 71\,337 \quad\Longleftrightarrow\quad \text{digits} \gtrsim 71\,337 \cdot \log_{10} 2 \approx 21\,475 .

nn 桁のオペランド二つの積は最大 2n2n 桁なので、乗算、FMA、二乗は n≈10 738n \approx 10\,738 から限界を超えうります。そのためスケーリング用ランナーはすべての演算を 10,000 桁で止め、16,384 桁と 20,000 桁では add、subtract、divide、compare だけを実行します。これらはオーバーフロー条件から得た解析的な上界で、32,768 桁と 65,536 桁での実行では中断が再現しました。

正しさと不変条件

正準形。 c≠0c \ne 0 について、normalize は s′≥0s' \ge 0 かつ（s′=0s' = 0 または 10∤c′10 \nmid c'）を満たす (c′,s′)(c', s') を返し、各ステップは (10q,s)(10q, s) を (q,s−1)(q, s - 1) に置き換えるだけなので v(c′,s′)=v(c,s)v(c', s') = v(c, s) です。このような二つの組が同じ数を表すのは、両者が等しいときに限ります：

c110−s1=c210−s2, s1<s2 ⇒ c2=c110 s2−s1 ⇒ 10∣c2 and s2>0,\begin{aligned} c_1 10^{-s_1} = c_2 10^{-s_2},\ s_1 < s_2 &\ \Rightarrow\ c_2 = c_1 10^{\,s_2 - s_1} \\ &\ \Rightarrow\ 10 \mid c_2 \text{ and } s_2 > 0 , \end{aligned}

これは (c2,s2)(c_2, s_2) の正規形に矛盾し、s1=s2s_1 = s_2 なら c1=c2c_1 = c_2 となります。canonical_string は符号、∣c′∣|c'| の各桁、小数点の位置 s′s' を書きますが、これらはすべて数によって決まるので、文字列の等価は数の等価です。

オラクルの除算。 N′/M′N'/M' は既約分数です。M′∣N′10kM' \mid N' 10^{k} なら gcd⁡(M′,N′)=1\gcd(M', N') = 1 より M′∣10k=2k5kM' \mid 10^{k} = 2^k 5^k なので、M′M' は 22 と 55 以外の素因数を持ちません。逆に M′=2α5βM' = 2^{\alpha} 5^{\beta} なら、そのような最小の kk は max⁡(α,β)\max(\alpha, \beta) です。オラクルはループの前に因数分解を確認するのでループは停止し、返すスケールは最短の厳密なスケールです。

整数平方根。 オラクルは次を反復します。

xk+1=⌊xk+⌊n/xk⌋2⌋,x0=10⌈D/2⌉>n,x_{k+1} = \Bigl\lfloor \frac{x_k + \lfloor n / x_k \rfloor}{2} \Bigr\rfloor, \qquad x_0 = 10^{\lceil D/2 \rceil} > \sqrt{n},

ここで DD は nn の桁数です。相加相乗平均の不等式 (x+n/x)/2≥n(x + n/x)/2 \ge \sqrt{n} より、各反復値は ⌊n⌋\lfloor \sqrt{n} \rfloor 以上です。xk>⌊n⌋x_k > \lfloor \sqrt{n} \rfloor の間は n/xk<xkn / x_k < x_k なので xk+1<xkx_{k+1} < x_k です。列は ⌊n⌋\lfloor \sqrt{n} \rfloor に達するまで狭義単調に減少し、そこでループが止まります。その後オラクルは x2=nx^2 = n を要求するので、平方数でない入力は丸めた根を出すのではなく中断します。

許容誤差 0 の健全性。 xx を厳密な結果、pp をフィクスチャの精度とします。

  1. xx の有効桁数が pp 以下なら、規格に準拠した GDA 演算は xx に等しい数を返します。xx はどの丸めモードでも自分自身が正しく丸めた値だからです。よって κ(I(x))=κ(O(x))\kappa(I(x)) = \kappa(O(x)) で、正しい実装が拒否されることはありません。
  2. xx が pp 桁を超えるなら、0 方向への丸めは ∣round⁡p(x)∣<∣x∣|\operatorname{round}_p(x)| < |x| を与えるので正準文字列が異なり、検証は失敗します。精度不足が受け入れられることはありません。
  3. 正準文字列の等価は数の等価なので、誤った結果が受け入れられることはありません。

精度契約は生成されるすべてのオペランドクラスについて (1) の前提を確立します。それらのクラスの外でも (2) は成り立つので、検査はフェイルクローズドになります。

決定性。 オペランドは Mare Mark のシード付き derive_seed(seed, "<operation>:<digits>", profile) から導かれ、時計は使わないので、同じシードは同じコーパスと同じフィンガープリントを再現します。generate_decimal は数字 1..91..9 だけを使うので、生成された係数は要求どおりの桁数を持ち、末尾に 0 を持ちません。

オラクルのコスト。 nn 桁の係数では、加算は 10∣sℓ−sr∣10^{|s_\ell - s_r|} を掛ける桁揃えが、乗算は BigInt の積一回が、除算は GCD 一回と 1010 倍 max⁡(α,β)\max(\alpha, \beta) 回（M′=2α5βM' = 2^{\alpha} 5^{\beta}）が支配的です。いずれも計時領域の外で実行されます。

却下した代替案

  • ペアごとの一致だけ。 却下：不一致の原因を特定できず、共通の誤りを見逃します。
  • 一方のライブラリを他方のオラクルにする。 却下：GDA ライブラリ自体が検査対象の一つであり、ベンチマークが見つけた DzmingLi の失敗を GDA の失敗と区別できなくなります。
  • 許容誤差付きの比較（ulp や相対誤差）。却下：ベンチマーク対象の結果はすべて十進で厳密なので、0 でない許容誤差は誤りを隠すだけです。
  • 固定の大きな精度（10510^5 桁など）。却下：除算の仕事量を結果に必要な分以上に増やしうるうえ、計時を任意の定数に結び付けてしまいます。
  • 実行ごとにランダムなオペランド。 却下：結果は再現可能でなければならず、フィンガープリントも計時範囲をまたいで安定している必要があります。
  • exp、ln、log10 のスケーリングベンチマーク。 却下：恒等的なフィクスチャではアルゴリズムを動かせず、一般の結果には独立した高精度の超越関数オラクルが必要です。これらは decTest 監査だけでカバーします。

境界

  • このパッケージが検査するのは数値だけです。指数のコホート、結果の末尾の 0、ステータスフラグ、NaN、無限大、符号付きゼロは比較の対象外で、decTest 監査がカバーします。
  • 除算は有限小数になる除数 22、88、2525 だけで試験します。循環小数の商や大きな分母はベンチマークせず、オラクルも循環小数の商を拒否します。
  • Parse と Format は恒等パスです。GDA の書式化は試験せず、DzmingLi の 329 件の toSci 失敗は decTest 監査が別に記録しています。
  • OperandShape はラベルにすぎません。コーパスはサイズごとに三つの決定的なスケールプロファイルで、ランダムな分布ではありません。
  • 精度契約が証明されているのは上に挙げた生成オペランドクラスについてであり、任意の入力についてではありません。
  • 実行ファイルは native ターゲットでだけ動作し、他のターゲットではメッセージを表示して終了します。
  • DzmingLi/decimal@0.2.2 は非推奨となり moonbit-community/decimal に移行しています。ベンチマークは歴史的な比較のためにこれを固定しています。

Footnotes

  1. Mare Mark の compare_paired は Δi\Delta_i の四分位範囲を区間として保存します。これは対応差のばらつきを表すもので、中央値の不確かさではありません。 ↩