core チュートリアル

このチュートリアルでは、arithmetic のインストールから、必要なものを正確に述べる数値コードを書くところまで進む。解析的能力に対する汎用関数、不正な入力を処理できる値に変える検査付き演算、結果がどう丸められたかを報告するコンテキスト付き演算、そして同じトレイトに組み込む独自の型を扱う。どの例も _test.mbt ファイルに貼り付けて moon test で実行できるテストであり、期待される出力は inspect の呼び出しに書かれている。

クイックスタート

パッケージをモジュールに追加する：

moon add Luna-Flow/arithmetic@0.5.0

それを使うパッケージの moon.pkg でインポートする。Luna Flow のコードは別名 @lf_arith を使う：

import {
  "Luna-Flow/arithmetic" @lf_arith,
}

最小の実用的なプログラムは、一つの能力を求め、それを二つの数値型で使う：

fn[T : Add + Mul + @lf_arith.Sqrt] hypot(x : T, y : T) -> T {
  @lf_arith.Sqrt::sqrt(x * x + y * y)
}

test "quick start" {
  inspect(hypot(3.0, 4.0), content="5")
  inspect(hypot((6.0 : Float), (8.0 : Float)), content="10")
}

hypot は +、*、平方根を持つ任意の型で動き、それ以外では動かない。これがパッケージ全体の考え方である。アルゴリズムは自分が使う最小の能力の集合を名指しする。

日常的な作業

初等関数を汎用的に呼ぶ

非検査トレイト（Sqrt、Exponential、Logarithmic、Trigonometric など）は Self を返し、特殊な場合は型に任せる。コードを汎用に保つため、@lf_arith.Trigonometric::sin(x) のようにトレイト経由で呼ぶ：

fn[T : Add + Mul + @lf_arith.Trigonometric] sin_plus_cos_squared(x : T) -> T {
  let s = @lf_arith.Trigonometric::sin(x)
  let c = @lf_arith.Trigonometric::cos(x)
  s * s + c * c
}

test "pythagorean identity" {
  let v = sin_plus_cos_squared(0.7)
  inspect((v - 1.0).abs() < 1.0e-15, content="true")
  let pi : Double = @lf_arith.Constants::pi()
  inspect(@lf_arith.Logarithmic::ln(@lf_arith.Exponential::exp(pi)) == pi, content="true")
}

定数は引数をとらないので、上の pi のように欲しい型を注釈する。

不正な入力をエラー値に変える

検査付きトレイトは Result[_, ArithmeticError] を返す。負の数の平方根やゼロ除数が、後で NaN として見つかるものではなく呼び出し側が処理すべき場合であるときに使う。この二次方程式の解法は桁落ちのない形 $q = -\tfrac12\bigl(b + \operatorname{sign}(b)\sqrt{b^2 - 4ac}\bigr)$、$x_1 = q/a$、$x_2 = c/q$ を使う：

fn real_roots(
  a : Double,
  b : Double,
  c : Double,
) -> Result[(Double, Double), @lf_arith.ArithmeticError] {
  let ctx = @lf_arith.ArithmeticContext::new(53)
  let s = match @lf_arith.SqrtChecked::sqrt_checked(b * b - 4.0 * a * c, ctx) {
    Ok(s) => s
    Err(e) => return Err(e)
  }
  let q = if b >= 0.0 { -0.5 * (b + s) } else { -0.5 * (b - s) }
  let x1 = match @lf_arith.DivChecked::div_checked(q, a, ctx) {
    Ok(x) => x
    Err(e) => return Err(e)
  }
  @lf_arith.DivChecked::div_checked(c, q, ctx).map(x2 => (x1, x2))
}

test "real roots" {
  guard real_roots(1.0, -3.0, 2.0) is Ok((x1, x2)) else { fail("no roots") }
  inspect(x1, content="2")
  inspect(x2, content="1")
  guard real_roots(1.0, 0.0, 1.0) is Err(e) else { fail("expected an error") }
  inspect(e.is_domain_error(), content="true")
  guard real_roots(0.0, 2.0, -4.0) is Err(e) else { fail("expected an error") }
  inspect(e.is_division_by_zero(), content="true")
}

エラーの種別が何が悪かったかを呼び出し側に伝える。実数解がない場合は DomainError、$a = 0$ の退化した方程式は DivisionByZero である。

NaN かもしれない値を比較する

Double の < は NaN が関わると常に false を返すので、ソートや最大値の計算が黙って誤る。CompareChecked は順序付けられない場合をエラーにする：

fn checked_max(xs : Array[Double]) -> Result[Double, @lf_arith.ArithmeticError] {
  let mut best = xs[0]
  for x in xs {
    match @lf_arith.CompareChecked::compare_checked(x, best) {
      Ok(1) => best = x
      Ok(_) => ()
      Err(e) => return Err(e)
    }
  }
  Ok(best)
}

test "checked maximum" {
  inspect(checked_max([3.0, 7.5, -1.0]).unwrap(), content="7.5")
  let r = checked_max([3.0, 0.0 / 0.0, 7.5])
  inspect(r is Err(e) && e.is_unordered_comparison(), content="true")
}

複数のステップにわたって診断を集める

コンテキスト付き演算は ArithmeticOutcome、つまり値とその計算中に立った診断を返す。ステップをつなぐには、値を次へ渡し、診断を combine する。小さなヘルパーが任意のコンテキスト付き演算についてこれを行う：

fn[A, B] and_then(
  r : Result[@lf_arith.ArithmeticOutcome[A], @lf_arith.ArithmeticError],
  next : (A) -> Result[@lf_arith.ArithmeticOutcome[B], @lf_arith.ArithmeticError],
) -> Result[@lf_arith.ArithmeticOutcome[B], @lf_arith.ArithmeticError] {
  match r {
    Err(e) => Err(e)
    Ok(first) =>
      next(first.value).map(second => @lf_arith.ArithmeticOutcome::with_diagnostics(
        second.value,
        first.diagnostics.combine(second.diagnostics),
      ))
  }
}

test "diagnostics survive a chain" {
  let ctx = @lf_arith.ArithmeticContext::new(24)
  let embedded : Result[@lf_arith.ArithmeticOutcome[Float], _] = @lf_arith.IntegralContextual::from_int_contextual(
    16_777_217, ctx,
  )
  let root = and_then(embedded, x => @lf_arith.SqrtContextual::sqrt_contextual(x, ctx)).unwrap()
  inspect(root.value, content="4096")
  inspect(root.diagnostics.inexact, content="true")
}

$2^{24} + 1$ は Float に収まらないので、埋め込みはそれを $2^{24}$ に丸めて inexact を立てる。$2^{24}$ の平方根は正確だが、合成した診断は以前の丸めを覚えている。フラグは蓄積するだけである。

表現可能な数をたどる

AdjacentContextual は隣の表現可能な値を返す。ある数とその次の数との差は一最終桁単位（ulp）であり、その大きさで形式がどれだけ細かいかを示す：

fn ulp(x : Double) -> Double {
  let ctx = @lf_arith.ArithmeticContext::new(53)
  @lf_arith.AdjacentContextual::next_plus_contextual(x, ctx).unwrap().value - x
}

test "ulp grows with magnitude" {
  inspect(ulp(1.0), content="2.220446049250313e-16")
  inspect(ulp(1024.0), content="2.2737367544323206e-13")
  inspect(ulp(9007199254740992.0), content="2")
}

$2^{53}$ では連続する double は $2$ ずつ離れているので、それより大きい整数はもはやすべては表現できない。

証明の失敗を報告する

証明付きバックエンドは、結果を証明できなかったことを CertificationFailure エラーで報告する。定義域エラーとは分けて処理すること。入力は有効であり、予算を増やせば成功するかもしれない。

fn describe(err : @lf_arith.ArithmeticError) -> String {
  match err.certification_failure_detail() {
    Some(d) =>
      "\{d.operation()}: gave up after \{d.refinements()} refinements at \{d.work_precision()} bits"
    None => err.message
  }
}

test "describe errors" {
  let detail = @lf_arith.CertificationFailureDetail::new(
    "sinh",
    @lf_arith.CertificationStage::TargetRounding,
    @lf_arith.CertificationFailureReason::RefinementBudgetExhausted,
    53,
    1024,
    5,
  )
  inspect(
    describe(@lf_arith.ArithmeticError::certification_failure(detail)),
    content="sinh: gave up after 5 refinements at 1024 bits",
  )
  inspect(
    describe(@lf_arith.ArithmeticError::domain_error("negative input")),
    content="negative input",
  )
}

さらに進んで

独自の型にコンテキスト付きトレイトを実装する

能力は守れる場合にだけ実装すること。この固定小数点型は百分の一単位で値を保持する。二つの値を掛けると一万分の一単位になり、丸め戻す必要がある。実装はコンテキストの丸めモードに従い、丸めを診断で報告し、サポートしないモードは無視せずに拒否する：

struct Cents(Int) derive(Eq, Debug)

impl @lf_arith.MulContextual for Cents with mul_contextual(x, y, ctx) {
  let raw = x.0 * y.0
  let q = raw / 100
  let r = raw % 100
  if r == 0 {
    return Ok(@lf_arith.ArithmeticOutcome::exact(Cents(q)))
  }
  let away = if raw < 0 { q - 1 } else { q + 1 }
  let rounded = match ctx.rounding {
    TowardZero => q
    AwayFromZero => away
    ToNearestEven => {
      let twice = r.abs() * 2
      if twice > 100 || (twice == 100 && q % 2 != 0) { away } else { q }
    }
    _ =>
      return Err(
        @lf_arith.ArithmeticError::unsupported("Cents rounds only toward or away from zero, or to nearest"),
      )
  }
  let flags = @lf_arith.ArithmeticDiagnostics::new(inexact=true, rounded=true)
  Ok(@lf_arith.ArithmeticOutcome::with_diagnostics(Cents(rounded), flags))
}

test "fixed-point multiplication" {
  let nearest = @lf_arith.ArithmeticContext::new(2)
  let exact = @lf_arith.MulContextual::mul_contextual(Cents(150), Cents(150), nearest).unwrap()
  debug_inspect(exact.value, content="Cents(225)")
  inspect(exact.diagnostics.inexact, content="false")
  let tie_down = @lf_arith.MulContextual::mul_contextual(Cents(5), Cents(10), nearest).unwrap()
  debug_inspect(tie_down.value, content="Cents(0)")
  let tie_up = @lf_arith.MulContextual::mul_contextual(Cents(15), Cents(10), nearest).unwrap()
  debug_inspect(tie_up.value, content="Cents(2)")
  inspect(tie_up.diagnostics.rounded, content="true")
  let floor = @lf_arith.ArithmeticContext::new(2, rounding=@lf_arith.RoundingMode::TowardNegative)
  inspect(@lf_arith.MulContextual::mul_contextual(Cents(5), Cents(10), floor) is Err(_), content="true")
}

$0.05 \times 0.10 = 0.005$ と $0.15 \times 0.10 = 0.015$ はどちらもちょうど中間である。偶数への最近接丸めはそれらを、最後の桁が偶数の隣接値 $0.00$ と $0.02$ に送る。

包含区間を三通りの結果で比較する

区間やボールの型では、「$x < y$ か？」には三つの答えがある。許容されるすべての値で yes、すべての値で no、あるいは不明である。二つの関係トレイトでそれを計算できる：

struct Interval {
  lo : Double
  hi : Double
}

impl @lf_arith.DefinitelyLt for Interval with definitely_lt(x, y) { x.hi < y.lo }

impl @lf_arith.DefinitelyLe for Interval with definitely_le(x, y) { x.hi <= y.lo }

enum Truth {
  Yes
  No
  Unknown
} derive(Debug)

fn[X : @lf_arith.DefinitelyLt + @lf_arith.DefinitelyLe] less(x : X, y : X) -> Truth {
  if @lf_arith.DefinitelyLt::definitely_lt(x, y) {
    Yes
  } else if @lf_arith.DefinitelyLe::definitely_le(y, x) {
    No
  } else {
    Unknown
  }
}

test "three-valued comparison" {
  let x = Interval::{ lo: 1.0, hi: 2.0 }
  debug_inspect(less(x, Interval::{ lo: 3.0, hi: 4.0 }), content="Yes")
  debug_inspect(less(x, Interval::{ lo: 0.0, hi: 1.0 }), content="No")
  debug_inspect(less(x, Interval::{ lo: 1.5, hi: 2.5 }), content="Unknown")
}

Unknown という答えは失敗ではない。包含区間を狭めて（より高い精度で計算して）もう一度問えばよい。Yes や No は包含区間が縮んでも変わらない。その理由は設計ページで証明する。区間やボールのバックエンドは同じトレイトを実装するので、less はそれらでもそのまま動く。

代数構造と組み合わせる

arithmetic は環や体を定義しない。それは luna-generic の役割である。アルゴリズムが両方を必要とするときは、境界の中で二つを組み合わせる。例えば体の演算だけを使う $\sqrt{a}$ のニュートン反復を、Sqrt と 1 ulp 以内で照合する：

fn[T : @lf_alg.Field] newton_sqrt(a : T, x0 : T, steps : Int) -> T {
  let one : T = @lf_alg.One::one()
  let two = one + one
  let mut x = x0
  for _ in 0..<steps {
    x = (x + a / x) / two
  }
  x
}

test "newton agrees with sqrt" {
  let a = 2.0
  let diff = newton_sqrt(a, 1.0, 6) - @lf_arith.Sqrt::sqrt(a)
  inspect(diff.abs() <= 2.220446049250313e-16, content="true")
}

そのために @lf_arith の隣に "Luna-Flow/luna-generic" @lf_alg をインポートする。

性能のための層の選び方

非検査トレイトはアロケーションのない直接呼び出しにコンパイルされる。検査付きトレイトは結果を Result に包み、コンテキスト付きトレイトはさらに ArithmeticOutcome を作る。ネイティブスカラーの内側のループでは、入力を検査付きの呼び出しで一度検証し、ループ内では非検査トレイトを使う。汎用コードは型ごとに特殊化されるので、トレイト境界の実行時コストはない。

よくある落とし穴

• Float と Double の空の診断は正確を意味しない。 それらのコンテキスト付き算術はコンテキストを無視し、丸めを検出しない。add_contextual(0.1, 0.2, ctx) は空の診断とともに 0.30000000000000004 を返す。損失を検出するのは Float の整数埋め込みだけである。フラグが重要ならコンテキストに忠実なバックエンドを使うこと。
• コンテキストは要求である。 ArithmeticContext::new(16) は Double の算術を十進にはしない。可能ならどうすべきかをバックエンドに伝えるだけである。ArithmeticContext::new(0) は黙って精度 1 になり、e_min が e_max より大きいと中断する。
• 非検査の Power は負の整数指数で中断する。 @lf_arith.Power::pow(2, -1) は Int、Int16、Int64、BigInt で中断し、固定幅整数のべき乗はオーバーフロー時にラップアラウンドする。逆数には浮動小数点型で PowIntChecked を使うこと。
• 非常に小さい底と負の指数。 pow_int_checked(1.0e-200, -2, ctx) は DivisionByZero を返す。逆数を取る前に $x^2$ がゼロへアンダーフローするからである。
• epsilon_contextual は単位丸め誤差ではなく $\varepsilon$ である。 最近接丸めの誤差の上界は $u = \varepsilon/2$ である。
• definitely_lt が偽であることは「以上」を意味しない。 重なり合う包含区間ではどちらの向きも偽である。上の less のように、逆向きの問いも立てること。
• maybe_eq は等価性ではない。 包含区間が点を共有することを述べるだけである。異なる二つの値が重なり合う包含区間を持つことはありうる。
• Double の x.sqrt() はトレイトではなく core のメソッドである。 汎用コードではトレイト @lf_arith.Sqrt::sqrt(x) を呼ぶこと。Double では結果は一致するが、すべての T : Sqrt で動くのはトレイトの呼び出しだけである。

次のステップ

• core API は、すべてのトレイト、型、インスタンスを正確な意味とともに列挙する。
• core 設計 は三つの層、明示的なコンテキスト、包含区間の論理を、丸め誤差と正しさの導出とともに説明する。
• luna-generic は、これらの能力と組み合わせる代数トレイトを提供する。