arithmetic 設計

設計目標

線形代数のアルゴリズムには、環の演算子以外にもいくつかのスカラー演算が必要です。ピボット選択のための絶対値、ノルムのための平方根、順序づけのための比較、そして浮動小数点の結果を比較する手段です。arithmetic パッケージは、スカラー型が実際には満たさない法則を満たすと主張することなく、これらの演算に線形代数コード向けの名前を与えます。上流の語彙があればそれを再利用し、ない場合にだけ trait を追加します。

数学的背景

演算と構造

Field のような構造 trait は等式、つまり結合律・分配律・逆元を約束します。Sqrt のような演算 trait が約束するのは、関数 $\sqrt{\cdot} : T \to T$ が存在することだけです。この区別が重要なのは、組み込みの浮動小数点型は演算を実装していても、等式は近似的にしか満たさないからです。Double では

\mathrm{fl}\big(\mathrm{fl}(2^{53} + 1) - 2^{53}\big) = 0 \ne 1 = \mathrm{fl}\big(2^{53} - 2^{53}\big) + 1 ,

となるため、一般に $(a + b) - c \ne (a - c) + b$ です。このパッケージの trait はすべて演算 trait です。

部分演算とその全域化

除算と平方根は実数上の部分関数です。 $x / y$ は $\{(x, y) : y \ne 0\}$ 上で、 $\sqrt{x}$ は $\{x \ge 0\}$ 上で定義されます。IEEE 754 は特殊値（ $\pm\infty$ 、NaN）を返すことでこれらを全域化しますが、その値は黙って伝播します。検査付きの演算は代わりに、部分関数 $f : D \rightharpoonup C$ を次のように全域化します。

\hat f : X \to C + E, \qquad \hat f(x) = \begin{cases} \mathrm{Ok}\,(f(x)) & x \in D, \\ \mathrm{Err}\,(e(x)) & x \notin D, \end{cases}

ここで $E$ はエラー値の集合で、 $e(x)$ は $x$ が定義域の外にある理由を表します。CheckedDiv は浮動小数点オペランドに対して次の定義域分割を使います。

領域	IEEE の結果	検査付きの結果
$y \ne 0$ 、両方が無限大ではない	$\mathrm{fl}(x / y)$	同じ値の `Ok`
$x = y = 0$	NaN	`DomainError`
$x, y$ がともに無限大	NaN	`DomainError`
$x \ne 0$ , $y = 0$	$\pm\infty$	`DivisionByZero`

NaN を含む順序づけ

Double 上の Compare は全順序ではありません。NaN はどの値よりも小さくも、等しくも、大きくもないため、NaN を含むデータのソートやピボット選択は順序に依存した結果になります。CheckedCompare は比較を順序づけ可能な部分集合に制限し、その外では UnorderedComparison を報告するため、結果はその定義域上で全順序になります。

近似的な等しさ

ApproxEq は絶対的な規則 $a \approx b \iff |a - b| \le \varepsilon$ を使います。この関係は反射的かつ対称的ですが、推移的ではありません。 $|a - b| \le \varepsilon$ と $|b - c| \le \varepsilon$ から三角不等式で得られるのは

|a - c| \le |a - b| + |b - c| \le 2\varepsilon ,

だけであり、この限界は達成されます。 $a = 0$ , $b = \varepsilon$ , $c = 2\varepsilon$ とすると、 $a \approx b$ , $b \approx c$ ですが $a \not\approx c$ です。したがって approx_eq は同値関係ではなく、このパッケージもそのようには扱いません。

絶対的な許容誤差はスケール不変でもありません。浮動小数点数の間隔は大きさとともに広がり、 $x$ 付近で隣り合う Double の値はおよそ $2^{-52}|x|$ 離れています。 $10^{20}$ 付近ではその間隔はおよそ $1.6 \times 10^{4}$ なので、異なる 2 つの値が $10^{-12}$ 以内に収まることはありません。 $10^{-20}$ 付近ではどの 2 つの値も収まります。相対的な規則 $|a - b| \le \varepsilon \max(|a|, |b|)$ はスケールの問題を解決しますがゼロ付近で破綻するため、慎重なコードは両方を組み合わせます。このパッケージはその選択を呼び出し側に任せます。

設計上の判断

上流の名前を再利用する

スカラー trait Zero、One、Inverse、Conjugate と、Luna-Flow/arithmetic の解析的な trait は、再定義せずに pub using で再エクスポートしています。ローカルにコピーすると、互換性のない 2 つ目の Sqrt ができてしまい、上流の trait を実装したスカラー型がローカルの trait を満たさなくなります。pub using を使えば、@la_arithmetic.Sqrt は @lf_arith.Sqrt そのもの です。

小さなローカル trait

Abs、ApproxEq、CheckedDiv、CheckedSqrt、CheckedCompare が存在するのは、上流パッケージがこの形で提供していなかった時点で、線形代数コードがこれらの名前を必要としたからです。それぞれメソッドは 1 つなので、スカラー型はサポートする演算だけを選んで実装できます。検査付き trait は Float と Double について上流の検査付き trait に委譲するため、2 つの層はあらゆる入力で一致します。

2 進浮動小数点ではコンテキストを受け取るが無視する

checked_div と checked_sqrt は ArithmeticContext を受け取ります。これにより、1 つのシグネチャで固定精度と任意精度の両方のスカラー型に対応できます。Float と Double では精度はハードウェア形式で固定され、丸めは最近接偶数丸めなので、コンテキストは効果を持ちません。

固定の絶対許容誤差

ApproxEq は Double に $10^{-12}$ を、Float に $10^{-6}$ を使います。これは各形式の丸め単位 $u$ （ $2^{-53}$ と $2^{-24}$ ）のおよそ $9000u$ と $17u$ にあたるので、正規化されたベクトルの要素のような 1 程度の大きさの値に適しています。他のスケールでは、自分のコードで明示的な許容誤差を使って比較してください。

正しさと不変条件

Float と Double について、checked_div(x, y, ctx) が Ok(v) を返すのは、IEEE の除算が無効演算例外やゼロ除算例外によらない値を返すときに限ります。そのとき v は IEEE の商とビット単位で一致します。
checked_sqrt(x, ctx) は $x \ge 0$ と NaN に対して Ok(√x) を返し、 $x < 0$ に対して Err を返します。 $-0.0 \ge 0$ が成り立つので、checked_sqrt(-0.0) は Ok(-0.0) になることに注意してください。
checked_compare はその定義域上で反対称です。checked_compare(a, b) == Ok(k) ならば checked_compare(b, a) == Ok(-k) です。
approx_eq は NaN でない値について反射的かつ対称的で、NaN に対しては偽です。

却下した代替案

すべての演算をまとめたローカルの Real または Number trait。 アルゴリズムがどの演算を使うかが見えなくなり、特殊なスカラーにすべてを実装させることになります。
approx_eq を Eq の一部にする。 Eq は同値関係でなければなりませんが、近似的な等しさは推移的ではありません。
ApproxEq での相対許容誤差。 ゼロ付近でアプリケーションに依存する方針が必要になります。代わりに trait を単純に保ち、ドキュメント化しています。

境界

arithmetic はベクトル・行列・バックエンドの型を定義せず、このリポジトリの他のパッケージにも依存しません。代数的構造の trait は定義せず（それらは luna-generic から来ます）、行列アルゴリズムの許容誤差も選ばず（それは mutable の Tolerance trait の役割です）、任意精度演算も実装しません。