ball_float チュートリアル

このチュートリアルでは、保証付きの包含区間を使った計算方法を学びます。不確かな実数を含むことがわかっている区間を構築し、それを算術と初等関数に通し、起こりうるすべての厳密な結果を依然として含む上下界を読み取ります。各節は出力付きの完全な例です。保証の背後にある数学は設計ページに、すべての項目は API リファレンスにあります。

クイックスタート

モジュールを追加し、パッケージをインポートします(端点の型を提供する bin_float も併せて)。

moon add Luna-Flow/floating@0.8.0
import {
  "Luna-Flow/floating/bin_float",
  "Luna-Flow/floating/ball_float",
}

BallFloat は、二つの端点が BinFloat の値である閉区間 [x‾,x‾][\underline{x}, \overline{x}] です。このページの例では次のヘルパーで区間を出力します。下端は下向きに、上端は上向きに丸めて書き出すので、出力された十進の区間は格納されている区間を依然として含みます。

///|
fn show(x : @ball_float.BallFloat, digits : Int) -> String {
  if x.is_empty() {
    return "[empty]"
  }
  let down = @bin_float.BinaryContext::unbounded(
    x.precision(),
    rounding=@bin_float.BinaryRoundingMode::RoundTowardNegative,
  )
  let up = @bin_float.BinaryContext::unbounded(
    x.precision(),
    rounding=@bin_float.BinaryRoundingMode::RoundTowardPositive,
  )
  let lo = x.lower_bound().to_decimal_string_ctx(digits, down).0
  let hi = x.upper_bound().to_decimal_string_ctx(digits, up).0
  "[" + lo + ", " + hi + "]"
}

最小限の有用なプログラムは、1/31/3 を 53 ビットで包含します。

///|
test "quick start: enclose one third" {
  let one = @ball_float.BallFloat::from_int(1, precision=53)
  let three = @ball_float.BallFloat::from_int(3, precision=53)
  let third = one / three
  inspect(show(third, 20), content="[3.3333333333333331482e-1, 3.3333333333333337035e-1]")
  // The endpoints are neighbouring 53-bit numbers: the width is one unit
  // in the last place, 2^-54.
  inspect(third.width().to_string(), content="1p-54")
}

53 ビットの BinFloat で 1/31/3 に等しいものはないので、結果はそれを囲む最小の 53 ビットの区間になります。下端は 1/31/3 を下向きに丸めたもの、上端は 1/31/3 を上向きに丸めたものです。

日常的なタスク

区間を構築する

値について何がわかっているかに応じてコンストラクタを選びます。

///|
test "constructors" {
  let two = @bin_float.BinFloat::from_int(2, precision=53)
  let five = @bin_float.BinFloat::from_int(5, precision=53)
  // A known point: the singleton {2}.
  let point = @ball_float.BallFloat::exact(two)
  inspect(point.is_singleton(), content="true")
  // Known bounds: [2, 5].
  let range = @ball_float.BallFloat::from_bounds(two, five)
  inspect(show(range, 3), content="[2.00e+0, 5.00e+0]")
  // A measurement: 5 +/- 2, stored as the endpoints [3, 7].
  let measured = @ball_float.BallFloat::new(five, two)
  inspect(show(measured, 3), content="[3.00e+0, 7.00e+0]")
  inspect(measured.center().to_string(), content="5p0")
  inspect(measured.radius().to_string(), content="1p1")
  // Untrusted bounds: reversed endpoints are an error value, not an abort.
  let reversed = @ball_float.BallFloat::try_from_bounds(five, two)
  inspect(reversed is Err(_), content="true")
}

from_bounds と exact は不正な入力(NaN の端点、逆転した上下界、有限でない点)で異常終了します。それらの try_ 形式は同じ条件を ArithmeticError として返します。new(center, radius) はコンストラクタとしての見え方にすぎず、格納される値は常に端点の組であり、center() / radius() はそこから中心–半径形式を再計算します。実数直線全体と空集合は BallFloat::whole() と BallFloat::empty() です。

十進定数を包含する

BallFloat::from_double(0.1) は 0.10.1 に最も近い binary64 の数を包含しますが、それは 0.10.1 ではありません。十進の値そのものを包含するには、十進文字列を一度は下向きに、一度は上向きに丸め、その二つの結果を上下界として使います。

///|
test "enclose decimal 0.1" {
  let down = @bin_float.BinaryContext::unbounded(
    53,
    rounding=@bin_float.BinaryRoundingMode::RoundTowardNegative,
  )
  let up = @bin_float.BinaryContext::unbounded(
    53,
    rounding=@bin_float.BinaryRoundingMode::RoundTowardPositive,
  )
  let lo = @bin_float.BinFloat::from_string_ctx("0.1", down).unwrap().0
  let hi = @bin_float.BinFloat::from_string_ctx("0.1", up).unwrap().0
  let tenth = @ball_float.BallFloat::from_bounds(lo, hi)
  inspect(show(tenth, 20), content="[9.9999999999999991673e-2, 1.0000000000000000556e-1]")
  // The binary64 double 0.1 is one of the two endpoints, so the singleton
  // built from it misses the other half of the uncertainty.
  let double_tenth = @ball_float.BallFloat::from_double(0.1)
  inspect(double_tenth.is_singleton(), content="true")
  inspect(double_tenth.subset(tenth), content="true")
}

十進テキストや丸められた計算からデータが入ってくるすべての境界でこれを行ってください。そうすれば区間の領域の内部では保証が自動的に維持されます。

計算して結果を読み取る

算術演算子は、オペランドの点に演算を適用したすべての結果を含む区間を返します。

///|
test "arithmetic on intervals" {
  let x = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(1, precision=53),
    @bin_float.BinFloat::from_int(2, precision=53),
  )
  let y = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(-3, precision=53),
    @bin_float.BinFloat::from_int(5, precision=53),
  )
  inspect(show(x + y, 3), content="[-2.00e+0, 7.00e+0]")
  inspect(show(x - y, 3), content="[-4.00e+0, 5.00e+0]")
  inspect(show(x * y, 3), content="[-6.00e+0, 1.00e+1]")
  // y contains 0 in its interior, so x / y is the whole real line.
  inspect(show(x / y, 3), content="[-inf, inf]")
  // y / x is bounded because x stays away from 0.
  inspect(show(y / x, 3), content="[-3.00e+0, 5.00e+0]")
}

結果は、端点には lower_bound() / upper_bound()、上向きに丸めた x‾−x‾\overline{x} - \underline{x} には width()、代表点には midpoint()、形状には is_bounded() / is_entire() / is_empty() を使って読み取ります。

区間を比較する

区間はその内部にある未知の一点を表すので、「x<yx < y か?」には三通りの答えがあります。確実にそうである、そうかもしれない、確実にそうでない、です。関係演算はそのどれかを示します。

///|
test "relations" {
  let a = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(9, precision=53),
    @bin_float.BinFloat::from_int(11, precision=53),
  )
  let b = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(10, precision=53),
    @bin_float.BinFloat::from_int(16, precision=53),
  )
  let c = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(12, precision=53),
    @bin_float.BinFloat::from_int(13, precision=53),
  )
  inspect(a.definitely_lt(c), content="true") // every point of a < every point of c
  inspect(a.definitely_lt(b), content="false") // not certain ...
  inspect(a.maybe_eq(b), content="true") // ... because they share points
  inspect(c.subset(b), content="true")
  inspect(a.contains(@bin_float.BinFloat::from_int(10, precision=53)), content="true")
  debug_inspect(a.overlap_state(b), content="OverlapsState")
}

contains は点を判定し、subset、interior、set_equal は集合を比較します。definitely_lt、definitely_le、definitely_gt は、すべての点の組について順序が成り立つ場合にのみ成り立ちます。maybe_eq(overlaps と同じ)は、ある組が等しい可能性がある場合に成り立ちます。overlap_state は二つの区間の IEEE 1788 による完全な分類を与えます。

初等関数を評価する

関数は、引数のうち ff が定義されているすべての ξ\xi について f(ξ)f(\xi) を結果が含むように評価されます。区間内部の極値や極も考慮されます。

///|
test "elementary functions" {
  let zero_to_four = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(0, precision=53),
    @bin_float.BinFloat::from_int(4, precision=53),
  )
  // sin reaches its maximum 1 at pi/2, which lies inside [0, 4].
  inspect(show(zero_to_four.sin_interval(), 6), content="[-7.56803e-1, 1.00000e+0]")
  inspect(show(zero_to_four.exp_interval(), 6), content="[1.00000e+0, 5.45982e+1]")
  // ln is undefined at 0; the result encloses ln over (0, 4].
  inspect(show(zero_to_four.ln_interval(), 6), content="[-inf, 1.38630e+0]")
  // [1, 2] contains the pole pi/2 of tan.
  let one_to_two = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(1, precision=53),
    @bin_float.BinFloat::from_int(2, precision=53),
  )
  inspect(one_to_two.tan_interval().is_entire(), content="true")
  // Outside the domain the result is the empty set.
  let negative = zero_to_four.neg() - one_to_two
  inspect(negative.sqrt_interval().is_empty(), content="true")
}

すべての *_interval 関数には、常に有効な包含区間を返す全域的な形式があり、多くには try_*_interval 形式もあります。両者の違いが現れるのは、認証付きの評価が精度予算を使い切った場合だけです。そのとき全域的な形式は安全な範囲(例えば sin_interval では [−1,1][-1, 1])に広げ、try_ 形式は失敗を説明する ArithmeticError を返します。

目標の形式に丸める

BallContext は二進形式(精度と指数範囲)を記述します。apply_ctx は区間をその形式へ外向きに丸め、何が起きたかを BallFlags で報告します。*_ctx 演算子は計算した後にコンテキストを適用します。

///|
test "round into binary32" {
  let ctx = @ball_float.BallContext::binary32()
  let x = @ball_float.BallFloat::from_int(1, precision=64)
  let y = @ball_float.BallFloat::from_int(3, precision=64)
  let (q, flags) = x.div_ctx(y, ctx)
  inspect(q.precision(), content="24")
  inspect(show(q, 9), content="[3.33333313e-1, 3.33333344e-1]")
  inspect(flags.inexact(), content="true")
  inspect(flags.overflow(), content="false")
  // A bound beyond the binary32 range becomes infinite (or the largest
  // finite value, on the side where that is still an enclosure).
  let big = @ball_float.BallFloat::from_double(1.0e300)
  let (clamped, big_flags) = big.apply_ctx(ctx)
  inspect(show(clamped, 9), content="[3.40282346e+38, inf]")
  inspect(big_flags.overflow(), content="true")
}

さらに進んで

装飾

BallFloatDecorated は区間と IEEE 1788 の装飾を組にしたもので、装飾はその区間を生成した関数評価について何がわかっているかを記録します。com(定義済み・連続・有界)、dac(定義済み・連続)、def(定義済み)、trv(何もわからない)、ill(特殊値 NaI、「区間ではない」)です。

///|
test "decorated evaluation" {
  let x = @ball_float.BallFloatDecorated::new(
    @ball_float.BallFloat::from_bounds(
      @bin_float.BinFloat::from_int(-1, precision=53),
      @bin_float.BinFloat::from_int(4, precision=53),
    ),
  )
  inspect(x.decoration(), content="com")
  // sqrt is only defined on part of [-1, 4]: the bare result [0, 2] is
  // correct, but the decoration drops to trv.
  let root = x.sqrt_interval()
  inspect(show(root.interval(), 3), content="[0.00e+0, 2.00e+0]")
  inspect(root.decoration(), content="trv")
  // Decorations only go down: later operations cannot restore com.
  inspect((root + x).decoration(), content="trv")
  // NaI absorbs everything and is different from the empty set.
  let nai = @ball_float.BallFloatDecorated::nai()
  inspect((nai + x).is_nai(), content="true")
  inspect((nai + x).is_empty(), content="false")
}

結果が入力全体で定義済みかつ連続な関数に裏付けられているかどうかを呼び出し側が知る必要がある場合に装飾を使ってください。これは、区間の結果がブラウワーの不動点定理のような存在定理を証明できる条件です。包含区間だけが重要な場合は、素の BallFloat の方が簡単です。

包含トレイト上のジェネリックなコード

BallFloat は Luna-Flow/arithmetic の包含関係(Contains、Overlaps、DefinitelyLt、DefinitelyLe、MaybeEq)を実装しているので、それらのトレイトに対して書かれたコードは任意の包含型で動作します。

///|
fn[T : @lf_arith.DefinitelyLt] certainly_increasing(xs : Array[T]) -> Bool {
  for i in 1..<xs.length() {
    if !@lf_arith.DefinitelyLt::definitely_lt(xs[i - 1], xs[i]) {
      return false
    }
  }
  true
}

///|
test "generic enclosure code" {
  let xs = [1, 3, 5].map(fn(n) {
    @ball_float.BallFloat::from_int(n, precision=53) /
    @ball_float.BallFloat::from_int(3, precision=53)
  })
  inspect(certainly_increasing(xs), content="true")
  // The trait form of contains tests set inclusion, unlike the
  // point-taking method BallFloat::contains.
  inspect(@lf_arith.Contains::contains(xs[2], xs[2]), content="true")
}

checked の能力トレイト DivChecked、PowNatChecked、PowIntChecked は @lf_arith.ArithmeticContext を受け取り、その precision が結果の精度を決めます。BallFloat ではこれらが失敗することはありません。ゼロを含む区間による除算は非有界な包含区間を返します。

checked パイプライン

計算が失敗しうる複数のステップ(不正な構築、認証されない初等関数)からなる場合、関連パッケージ ball_float_checked が最初のエラーを演算の連鎖に通して伝えるので、失敗の判定は最後に一度だけ行えば済みます。

精度とコスト

各区間はビット単位の作業 precision を持ちます。二項演算の結果には二つのオペランドのうち大きい方の精度が使われます。算術のコストは、その精度での端点演算数回分です。初等関数は約 p+64p + 64 から p+192p + 192 ビットで認証付きの級数を評価し、難しい入力では最大 12 回まで精緻化することがあるため、乗算 1 回の数倍のコストがかかります。精度を上げると幅のうち丸めに由来する部分は狭まりますが、入力の幅に由来する部分は狭まりません。

よくある落とし穴

x - x はゼロではない。 区間算術は二つのオペランドを独立した未知数として扱うので、[1,2]−[1,2]=[−1,1][1, 2] - [1, 2] = [-1, 1] となります。x * x も同様で、x がゼロを含むときは x.square() より広くなります。

///|
test "dependency" {
  let x = @ball_float.BallFloat::from_bounds(
    @bin_float.BinFloat::from_int(-1, precision=53),
    @bin_float.BinFloat::from_int(2, precision=53),
  )
  inspect(show(x - x, 3), content="[-3.00e+0, 3.00e+0]")
  inspect(show(x * x, 3), content="[-2.00e+0, 4.00e+0]")
  inspect(show(x.square(), 3), content="[0.00e+0, 4.00e+0]")
  inspect(show(x.pown(2), 3), content="[0.00e+0, 4.00e+0]")
}

各不確かな量が一度だけ現れるように式を書き換え、手書きの積よりも square、pown、初等関数を優先してください。

== は集合の等価性ではない。 == は精度タグを含む格納された表現を比較します。集合の比較には set_equal を使ってください。

contains は点を受け取る。 メソッド BallFloat::contains は BinFloat の点を判定します。集合の包含は subset(またはトレイトメソッド @lf_arith.Contains::contains。こちらは「部分集合として含む」を意味します)です。

区間は順序付けられていない。 区間には全順序がありません。less(IEEE 1788 の集合関係)を数の比較であるかのように使ってソートしてはいけません。

非有界な区間には中心がない。 center()、radius()、midpoint() は空集合と片側だけ有界な区間で異常終了します(midpoint() は実数直線全体に対しては 0 を返します)。先に is_bounded() を確認するか、代わりに +∞+\infty を返す radius_extended() を使ってください。

精度より幅の広い整数。 現在 from_int(n, precision=p) と from_coefficient は、区間を構築する前に n を max⁡(p,8)\max(p, 8) ビットの最も近い値に丸めるため、それより幅の広い整数では結果が n を含まないことがあります。少なくとも整数のビット長以上の精度を使ってください(デフォルトの 16 ビットでカバーされるのは ∣n∣<216|n| < 2^{16} だけです)。

Double からの入力。 from_double(x) は x の二進値を厳密に包含しますが、x がどの十進数を近似したものかを知ることはできません。十進データは十進定数を包含するのようにして包含してください。

次のステップ

  • 設計ページでは、これらの結果がなぜ包含区間になるのかを証明し、端点の公式を導出し、装飾と認証付き初等関数を説明しています。
  • API リファレンスは、すべての項目をその特殊ケースとともに文書化しています。
  • ball_float_checked は失敗しうる区間演算を合成します。bin_float は端点の算術です。
  • 適合性のページでは、IEEE 1788 のテストコーパスが何をカバーしているかを述べています。