RFC 0001：Luna-Flow 数値微積分の基盤アーキテクチャ

ステータス：実験的検証
対象バージョン：次期メジャーバージョン
範囲：アーキテクチャ、アルゴリズムの選定、API の実現可能性。現行の本番実装は置き換えない

1. 決定の要約

次世代の calculus-numerical は、GSL の Double API と実装構造を手本とするのをやめ、Luna-Flow の代数トレイト、checked arithmetic、明示的な数値コンテキストの上に構築する。

正式な移行時に、現行の basic、deriv、diff、integration API は legacy/* に移し、メジャーバージョン 1 つ分の互換期間を設ける。
新 API の第 1 段階では、同一型のスカラー関数 (T) -> T のみをサポートする。座標の型、戻り値の型、誤差スカラーは当面分離しない。
Double、Float はコンテキストなしの高速実装を使い、Decimal、BinFloat は明示的な context adapter を使う。
完全なエントリポイントはコンテキスト、許容誤差、リソース上限を必ず受け取る。簡便なエントリポイントは文書化されたデフォルトコンテキストを使ってよい。Decimal のデフォルトは Decimal128 とする。
数値アルゴリズムは、意味の不透明な多値タプルではなく、名前付きの結果と診断情報を返す。
estimated_error は経験的な誤差推定を表す。certified_bound を返せるのは、区間演算または BallFloat に基づく実装だけである。

2. 標準と独立実装の原則

2.1 数値セマンティクス

2 進および 10 進浮動小数点の分類、丸め、特殊値は IEEE 754-2019 に従う。
Decimal のコンテキスト、フラグ、cohort、例外のセマンティクスは General Decimal Arithmetic と decNumber のテスト仕様に従う。
任意精度の参照値は MPFR または Arb で生成する。特定のマシン Double 実装を正しさの基準としない。

2.2 アルゴリズムの資料

実装は論文と公開されたアルゴリズム記述に基づく。成熟したライブラリは、インターフェースと挙動の相互検証にのみ参照する。

分野	第 1 段階の候補	主な根拠	選定理由
数値微分	Fornberg の重み	Fornberg, 1988/1998	任意の節点・任意の微分階数に対応し、手書きのステンシル定数を避けられる
適応的微分	Richardson/Ridders	Richardson 補外と Ridders 法	説明可能なステップ幅の列と誤差推定
固定求積	Gauss–Kronrod 15/21	Gauss–Kronrod 則と QUADPACK の文献	1 回のサンプリングで入れ子の誤差推定も得られる
適応積分	誤差最大の区間を優先して細分	QUADPACK のアルゴリズム解説	成熟したワークスペースと停止戦略
高精度積分	Tanh–Sinh	高橋・森の二重指数関数型公式	端点特異性と高精度 Decimal に適する
滑らかな関数の積分	Clenshaw–Curtis	Chebyshev 展開の文献	サンプルを再利用でき、高精度へ自然に拡張できる
有限和	pairwise、Kahan、Neumaier	浮動小数点誤差解析の文献	性能と安定性の段階を提供する
無限級数	Wynn epsilon、Euler 変換	収束加速の文献	交代級数や収束の遅い一般的な場合をカバーする
べき級数	打ち切り係数の代数	形式的べき級数の標準アルゴリズム	評価、微分、積分、合成の基礎をサポートする

GSL は GPL ソフトウェアである。本プロジェクトはその公開された挙動と出典論文を比較してよいが、ソースコードをコピー、翻訳、または構造的に書き換えてはならない。

3. 依存関係の境界

依存関係は一方向に保たなければならない。

luna-generic -> arithmetic -> floating
                         \-> calculus-numerical
floating ----------------> calculus-numerical
luna-poly ---------------> calculus-numerical/power_series adapters
autodiff ----------------> optional calculus bridge

この段階では arithmetic と floating を変更しない。通常の代数演算は既存のトレイトを直接組み合わせる。calculus の内部では、欠けていてコンテキストに依存する能力に限って極小の provider を使う。検証が安定した後、別の RFC を提出し、汎用的な next_up、数値フォーマット、コンテキストの能力を arithmetic に移し、floating が Decimal/BinFloat のインスタンスを実装する。

floating が calculus に依存してはならない。既存の Decimal 型のコンテキスト能力は calculus が所有する provider 値で提供する。四則演算、比較、ゼロ、1 などコンテキストに依存しない能力は引き続きトレイト合成を使う。

4. 型システムの決定

4.1 第 1 段階ではスカラーと値域を統一する

安定したエントリポイントは (T) -> T を採用する。現在の MoonBit/Luna-Flow エコシステムには、次の関係を表現する成熟した公開能力がない。

V + V -> V
S * V -> V
norm(V) -> E

linear-algebra の現在のベクトルのスケーリングも同じ要素型を使っている。S/V/E の 3 つの型を持つ API を公開すると、呼び出し側に大きな関数テーブルを渡させることになり、積分のホットパスに間接呼び出しが入るおそれがある。そのため、ベクトル、行列、異種の値域に対する積分は別の RFC に先送りする。

4.2 トレイト合成と最小限の provider

すべての初等関数を含む RealNumber トレイトは新設せず、既存のトレイト演算を関数テーブルに包み直すこともしない。アルゴリズムは必要な最小限の能力をトレイト合成で受け取る。

コアのスカラー：Zero、One、IntegralHomomorphism、Add、Sub、Mul、Div、Compare などの既存トレイトを直接組み合わせる。
数値フォーマット：epsilon、隣接値、最小正規数、最大有限値、値の分類。
オプションの能力：コンテキスト対応の sqrt、整数べき、パース。

既存のトレイトでは正しく表現できないコンテキスト依存の操作に限って極小の provider を使う。例えば (Context) -> T のコンテキスト付きの 1 や、(T, Context) -> T の next-up である。安定した後、これらの能力は arithmetic に移す。

4.3 コンテキストと epsilon

Float と Double のフォーマット定数は型によって決まる。Decimal の epsilon と範囲は精度、指数範囲、丸めコンテキストによって共に決まるため、引数なしの Decimal epsilon() を提供してはならない。

Decimal の作業用 epsilon は次のように定義する。

next_plus(1, context) - 1

アルゴリズムのエントリポイントは派生環境を 1 回だけ構築し、epsilon、sqrt(epsilon) とよく使う係数をキャッシュする。ホットループの中でコンテキストの構築、規則テーブルの変換、フォーマット定数の計算を繰り返してはならない。

5. パッケージと公開 API の草案

core             context, tolerance, status, diagnostics, experimental adapter
differentiation  Fornberg, Richardson/Ridders
integration      quadrature rule, adaptive controller, workspace
series           compensated sum, convergence acceleration
power_series     truncated series and calculus operations
legacy           current public API compatibility layer

提案する結果の形：

pub struct IntegrationResult[T] {
  value : T
  estimated_error : T
  evaluations : Int
  intervals : Int
  status : NumericalStatus
}

pub struct IntegrationOptions[T] {
  absolute_tolerance : T
  relative_tolerance : T
  max_evaluations : Int
  max_intervals : Int
}

NumericalStatus は少なくとも、完了、リソース枯渇、無効な許容誤差、丸めによる停滞、特異点の疑い、非有限の入力・出力、context arithmetic failure を区別する。下層の Decimal フラグは、単一の整数エラーコードに圧縮するのではなく、診断情報に集約するべきである。

6. アルゴリズムの境界

6.1 微分

規則層は節点と微分階数から Fornberg の重みを生成し、規則オブジェクトはキャッシュできる。
実行層はサンプリングと重み付けだけを担当する。
制御層は Richardson/Ridders でステップ幅を調整し、打ち切り誤差と丸め誤差を報告する。
autodiff は独立したままとする。calculus はオプションのブリッジとアルゴリズム選択の説明だけを提供し、二重数を重複して実装しない。

6.2 積分

QuadratureRule[T] は出典のある節点、Gauss の重み、拡張の重みを保持する。
固定規則と適応コントローラは分離する。コントローラは規則係数の出典を知らない。
適応ワークスペースは局所的な変更を許すが、可変状態を公開しない。
Double/Float の規則は具体的な静的テーブルを使う。Decimal/BinFloat の規則はエントリポイントでコンテキストに合わせて 1 回だけ変換し、その呼び出しの間キャッシュする。
Tanh–Sinh は固定精度のテーブルに依存せず、高精度および端点特異性のある場合の第一の補完手段とする。

6.3 級数

有限列には naive、pairwise、Kahan/Neumaier を提供し、1 つの sum が精度のコストを暗黙に選ぶことはしない。
無限級数では許容誤差と最大項数を明示しなければならず、使用した項数と終了理由を返す。
収束加速は合成可能な戦略とし、数列の生成器には混ぜない。
PowerSeries[T] は打ち切り次数を明示的に保持する。luna-poly との変換は明示的な adapter とし、打ち切り級数を通常の多項式と同一視しない。

7. Legacy の移行

正式な移行は 2 つのバージョンに分けて行う。

新パッケージのリリース時に、旧エントリポイントは legacy/* に転送し、ドキュメントで非推奨と明記して項目ごとの移行表を提供する。
次のメジャーバージョンで legacy を削除する。重大な正しさの問題は引き続き修正するが、legacy に新しいアルゴリズムや数値型は追加しない。

この RFC の段階では既存のファイルを移動しない。新 API が検証される前に互換性の変化を生まないためである。

8. 検証の基準

各分野ではまず Double/Decimal の縦断的なスライスを 1 つ提供してから、アルゴリズムの数を増やす。
Double adapter のプロトタイプは、直接実装に対するオーバーヘッドを 15% 以内とすることを目標とし、ホットループで項目ごとに割り当ててはならない。
Decimal は全体を通して同じ明示的コンテキストを使う。テストは inexact、rounded、underflow、overflow、invalid operation を網羅しなければならない。
解析的な関数、難しい関数、高精度オラクルを階層に分けてテストし、実際の誤差と報告された誤差を別々に記録する。
性質テストは線形性、区間の反転、定数関数、スケーリング関係、許容誤差に対する単調性、級数の打ち切りの一貫性を網羅する。

9. 今後の RFC

安定した数値フォーマットの能力を arithmetic に移す。
スカラーと値域の分離、およびベクトル・行列・複素数の積分。
BallFloat/区間演算による certified integration。
Fourier 級数と FFT エコシステムとの境界。
ODE/PDE ソルバー。この基盤能力の段階には含まれない。

参考文献

B. Fornberg, “Generation of Finite Difference Formulas on Arbitrarily Spaced Grids”, 1988.
B. Fornberg, “Calculation of Weights in Finite Difference Formulas”, 1998.
R. Piessens et al., QUADPACK: A Subroutine Package for Automatic Integration, 1983.
H. Takahasi and M. Mori, “Double Exponential Formulas for Numerical Integration”, 1974.
L. N. Trefethen, “Is Gauss Quadrature Better than Clenshaw–Curtis?”, 2008.
P. Wynn, “On a Device for Computing the e_m(S_n) Transformation”, 1956.
IEEE 754-2019, Standard for Floating-Point Arithmetic.
General Decimal Arithmetic Specification and decNumber test suite.