Luna Flow
日本語

bin_float チュートリアル

このチュートリアルでは、floating の任意精度二進浮動小数点型である BinFloat を使った計算方法を学びます。厳密な値を構築し、binary32 のような IEEE 754 形式でステータスフラグとともに計算を実行し、交換形式のビット列と十進テキストを読み書きし、正しく丸められた初等関数を呼び出します。すべての例はコンパイル可能で、出力を示しています。丸め規則の背後にある数学は設計ページに、すべての関数は API リファレンスにあります。

クイックスタート

モジュールを追加し、moon.pkg でパッケージをインポートします。

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

最小限の有用なプログラムは、1 を 3 で割る計算を 53 ビットで一度、binary32 形式で一度行い、両方の結果を、読み戻すと厳密に同じ値になる最短の十進数として出力します。

///|
test "one third, two precisions" {
  let one = @bin_float.BinFloat::one()
  let three = @bin_float.BinFloat::from_int(3)
  inspect((one / three).to_shortest_string(), content="0.3333333333333333")
  let binary32 = @bin_float.BinaryContext::binary32()
  let (third, flags) = one.div_ctx(three, binary32)
  inspect(third.to_shortest_string_ctx(binary32), content="0.33333334")
  inspect(flags.inexact(), content="true")
}

演算子 / はオペランドの精度（デフォルトでは 53 ビット）で最近接偶数丸めを行います。div_ctx は binary32 のコンテキストの下で丸め、IEEE フラグも返します。ここでは inexact が、三分の一は二進数ではないことを示しています。

日常的なタスク

厳密な値を構築する

有限の BinFloat は、整数係数 cc と指数 ee による c⋅2ec \cdot 2^{e} です。make はその両方と精度を受け取ります。to_string は格納された形 c p e を出力し、これは厳密です。

///|
test "exact dyadic values" {
  let x = @bin_float.BinFloat::make(
    @bin_float.BinCoeff::from_uint64(3UL),
    -1,
    32,
  )
  inspect(x, content="3p-1")
  inspect(x.to_shortest_string(), content="1.5")
  // Powers of two move into the exponent: 12 = 3 * 2^2.
  let twelve = @bin_float.BinFloat::from_int(12)
  inspect("\{twelve.coefficient()} \{twelve.exponent2()}", content="3 2")
}

十進の入力は from_string で読み込みます。これは十進値そのものを丸めます。ホストの Double を経由して変換すると、53 ビット以外の精度を求める場合に丸めが 2 回起こり、from_double はホストの二進近似をそのまま忠実に保持します。

///|
test "decimal input" {
  let tenth = @bin_float.BinFloat::from_string("0.1").unwrap()
  inspect(tenth, content="3602879701896397p-55")
  let tenth_200 = @bin_float.BinFloat::from_string("0.1", precision=200).unwrap()
  inspect(tenth_200.precision(), content="200")
  let host = @bin_float.BinFloat::from_double(0.1)
  // Same 53-bit value, but a 200-bit parse is much closer to one tenth.
  inspect(host.compare(tenth), content="0")
  inspect(host.compare(tenth_200), content="1")
}

IEEE 形式で計算してフラグを保持する

BinaryContext は演算の精度、指数範囲、丸め方向、極小性（tininess）の規則を固定します。*_ctx メソッドはその一回の演算の値とフラグを返します。フラグを蓄積するには自分で組み合わせてください。

///|
test "a binary16 computation with accumulated flags" {
  let half = @bin_float.BinaryContext::binary16()
  let x = @bin_float.BinFloat::from_int(60000)
  let (sum, add_flags) = x.add_ctx(x, half)
  let (product, mul_flags) = @bin_float.BinFloat::from_string("0.001")
    .unwrap()
    .mul_ctx(@bin_float.BinFloat::from_string("0.0001").unwrap(), half)
  let all = add_flags.combine(mul_flags)
  inspect(sum, content="inf")
  inspect(product.to_shortest_string_ctx(half), content="1e-7")
  inspect(
    "\{all.overflow()} \{all.underflow()} \{all.inexact()}",
    content="true true true",
  )
}

120000120000 は binary16 の最大値 6550465504 を超えるので、和は無限大にオーバーフローします。積 10−710^{-7} は binary16 の最小の正規化数 2−14≈6.1⋅10−52^{-14} \approx 6.1 \cdot 10^{-5} より小さいので、非正規化数の格子上で丸められ、アンダーフローを発生させます。丸め方向によってオーバーフローの結果は変わります。ゼロ方向への丸めでは、結果は最大の有限値にとどまります。

///|
test "overflow depends on the rounding direction" {
  let toward_zero = @bin_float.BinaryContext::binary16(
    rounding=@bin_float.BinaryRoundingMode::RoundTowardZero,
  )
  let x = @bin_float.BinFloat::from_int(60000)
  let (sum, flags) = x.add_ctx(x, toward_zero)
  // 2047 * 2^5 = 65504, the largest finite binary16 value.
  inspect(sum, content="2047p5")
  inspect(flags.overflow(), content="true")
}

交換形式のビット列を読み書きする

ファイルやプロトコルが binary16、binary32、binary64、binary128 のビットパターンを運ぶ場合は BinaryInterchange を使います。ホストの Double を決して経由しないので、signaling NaN、ペイロード、binary128 がそのまま保たれます。

///|
test "binary32 bits in and out" {
  let format = @bin_float.BinaryInterchangeFormat::Binary32
  let value = @bin_float.BinaryInterchange::from_hex("3FC00000", format)
    .unwrap()
    .to_bin_float()
  inspect(value, content="3p-1")
  let (doubled, flags) = value.add_ctx(value, format.context())
  let (bits, encode_flags) = doubled.to_interchange(format)
  inspect(bits.to_hex(), content="40400000")
  inspect(flags.combine(encode_flags).to_testfloat_bits(), content="0")
}

エンコードは形式へ丸めるので、フラグも返します。デコードは常に厳密です。

十進テキストを出力する

to_shortest_string_ctx は、ある形式で同じ値として読み戻せる最短の十進数を与えます。binary64 ではホストのフォーマッタと一致します。to_decimal_string_ctx は、コンテキストの方向に丸めた固定桁数の有効数字を与え、桁が捨てられたかどうかを報告します。

///|
test "decimal output" {
  let ctx = @bin_float.BinaryContext::binary64()
  let tenth = @bin_float.BinFloat::from_string("0.1").unwrap()
  inspect(tenth.to_shortest_string_ctx(ctx), content="0.1")
  let (exact, exact_flags) = tenth.to_decimal_string_ctx(55, ctx)
  inspect(exact, content="1.000000000000000055511151231257827021181583404541015625e-1")
  inspect(exact_flags.inexact(), content="false")
  let up = @bin_float.BinaryContext::binary64(
    rounding=@bin_float.BinaryRoundingMode::RoundTowardPositive,
  )
  let (short, _) = tenth.to_decimal_string_ctx(3, up)
  inspect(short, content="1.01e-1")
}

すべての二進数は有限の十進展開を持つので、55 桁あれば 0.1 の binary64 の値を厳密に出力できます。

初等関数を呼び出す

各初等関数には、素の形式、*_ctx 形式、try_*_ctx 形式があります。三つとも正しく丸められた結果を返します。違いは、定義域エラーや精緻化予算の枯渇をどう報告するかだけです。素の形式と *_ctx 形式は invalid_operation とともに quiet NaN を返し、try_*_ctx は何が問題だったかを示す @lf_arith.ArithmeticError を返します。

///|
test "certified elementary functions" {
  let ctx = @bin_float.BinaryContext::binary64()
  let two = @bin_float.BinFloat::from_int(2)
  inspect(two.ln().to_shortest_string(), content="0.6931471805599453")
  let (ln2, _) = two.ln_ctx(@bin_float.BinaryContext::binary128())
  inspect(ln2.to_shortest_string(), content="0.6931471805599453094172321214581766")
  match @bin_float.BinFloat::from_int(-2).try_ln_ctx(ctx) {
    Ok(_) => fail("ln(-2) is not real")
    Err(error) => inspect(error.message, content="ln requires a positive value")
  }
  let (sine, flags) = @bin_float.BinFloat::from_string("0.5")
    .unwrap()
    .sinpi_ctx(ctx)
  inspect("\{sine} \{flags.inexact()}", content="1p0 false")
}

sinpi(0.5) は厳密に 1 であり、関数はそれを知っています。厳密な結果は inexact フラグを発生させません。

さらに進んで

方向付き丸めで実数を包含する

負の無限大方向と正の無限大方向への丸めによって、厳密な結果を挟み込めます。これが区間ライブラリなしで保証付きの包含を構築する方法です。

///|
test "an enclosure of the square root of two" {
  let two = @bin_float.BinFloat::from_int(2)
  let (low, high) = @bin_float.sqrt_bounds_for_precision(two, 20).unwrap()
  inspect("\{low} \{high}", content="741455p-19 46341p-15")
  let down = @bin_float.BinaryContext::unbounded(
    20,
    rounding=@bin_float.BinaryRoundingMode::RoundTowardNegative,
  )
  let (low_exp, _) = two.exp_ctx(down)
  let up = @bin_float.BinaryContext::unbounded(
    20,
    rounding=@bin_float.BinaryRoundingMode::RoundTowardPositive,
  )
  let (high_exp, _) = two.exp_ctx(up)
  inspect(low_exp.compare(high_exp), content="-1")
}

ball_float はこの考え方を中心–半径算術としてパッケージ化しています。ball_float チュートリアルを参照してください。

融合演算と厳密な IEEE 演算を使う

fma は xy+zx y + z を一度だけ丸めるので、積の丸め誤差を復元できます。remainder は同じ精度のオペランドに対して常に厳密です。

///|
test "fma recovers the rounding error of a product" {
  let ctx = @bin_float.BinaryContext::binary64()
  let a = @bin_float.BinFloat::from_string("0.1").unwrap()
  let (p, _) = a.mul_ctx(a, ctx)
  let (error, flags) = a.fma_ctx(a, p.neg(), ctx)
  inspect(error.to_shortest_string(), content="-8.326672684688674e-19")
  // The error of a rounded product is exactly representable (Dekker).
  inspect(flags.inexact(), content="false")
  let r = @bin_float.BinFloat::one().remainder(a)
  inspect(r.to_hex(), content="-0x1p-54")
}

二つの積の誤差項は補償付きアルゴリズム（Kahan の総和法、double-double 算術）の基礎です。それが厳密である理由は設計ページで示しています。

アンダーフローと極小性に注意する

IEEE 754 では、極小性を丸めの前に検出するか後に検出するかを実装が選べ、その選択によって最小の正規化数のすぐ下の値に対するアンダーフローフラグが変わります。BinaryContext で選択できます。

///|
test "tininess before and after rounding" {
  // 2^-14 - 2^-27 lies below the smallest binary16 normal 2^-14 but
  // rounds up to it.
  let x = @bin_float.BinFloat::make(
    @bin_float.BinCoeff::from_uint64(8191UL),
    -27,
    13,
  )
  let after = @bin_float.BinaryContext::binary16()
  let before = @bin_float.BinaryContext::binary16(
    tininess=@bin_float.TininessDetection::BeforeRounding,
  )
  let (a, a_flags) = x.round_ctx(after)
  let (b, b_flags) = x.round_ctx(before)
  inspect("\{a} \{a_flags.underflow()}", content="1p-14 false")
  inspect("\{b} \{b_flags.underflow()}", content="1p-14 true")
}

値は同じで、違うのはフラグだけです。ハードウェアも異なり、x86 は丸めの後に、ARM は丸めの前に極小性を検出します。TestFloat は両方の規則を検査します。

ジェネリックなコードを書く

BinFloat は、floating のすべてのスカラーコアが共有する @def.Floating トレイトと、@lf_arith の checked トレイトおよび contextual トレイトを実装しています。それらのトレイトに対して書かれたコードは BinFloat、Decimal、その他のコアで動作します。

///|
fn[F : @lf_arith.DivContextual] reciprocal(
  x : F,
  ctx : @lf_arith.ArithmeticContext,
  one : F,
) -> Result[F, @lf_arith.ArithmeticError] {
  one.div_contextual(x, ctx).map(fn(outcome) { outcome.value })
}

///|
test "generic reciprocal on BinFloat" {
  let ctx = @lf_arith.ArithmeticContext::new(24, e_min=-126, e_max=127)
  let r = reciprocal(
    @bin_float.BinFloat::from_int(3),
    ctx,
    @bin_float.BinFloat::one(),
  ).unwrap()
  inspect(r.to_shortest_string(), content="0.33333334")
  inspect(
    reciprocal(@bin_float.BinFloat::zero(), ctx, @bin_float.BinFloat::one()) is Err(_),
    content="true",
  )
}

contextual トレイトは IEEE の division_by_zero フラグと invalid_operation フラグをエラーに変え、その他のフラグを診断情報として報告します。最初のエラーで停止するパイプライン全体には bin_float_checked を使ってください。

性能

精度にかかる時間は、おおよそそのビット数の整数乗算と同程度です。係数カーネルは筆算法から Karatsuba、Toom-3、数論変換へと自動的に切り替わるので、数千ビットでも実用的です。初等関数は結果が確実になるまで作業精度を精緻化し、ほぼすべての入力で最初の試行が成功します。精度を繰り返し広げたり狭めたりするのではなく、アルゴリズムごとに一つの精度を選んでください。測定結果は性能ページに記録されています。

よくある落とし穴

  • == は表現を比較する。 BinFloat は Eq を導出しているので、53 ビットの 1 と 24 ビットの 1 は == ではなく、−0≠+0-0 \ne +0 であり、NaN は同一の NaN と等しくなります。数値としての等価性には compare(x, y) == 0 を、IEEE の述語には equal_quiet を使ってください。
  • compare は NaN を最後に置く。 異常終了することはなく、すべての NaN はすべての NaN と等しく、すべての数より大きいと判定されるので、nan > x は真です。NaN を順序付けられないものとして扱う必要がある場合は compare_checked、compare_quiet、less_quiet を使ってください。
  • 素の演算子はあなたの形式を無視する。 x + y は大きい方のオペランドの精度と約 2±2302^{\pm 2^{30}} の指数範囲で動作し、binary32 の限界でオーバーフローすることはありません。形式のコンテキストとともに add_ctx を使ってください。
  • from_double(0.1) は十分の一ではない。 それはホストがすでに丸めた binary64 の値です。十進テキストは from_string でパースしてください。
  • to_string は十進ではない。 3602879701896397p-55 と出力します。人が読むためには to_shortest_string または to_decimal_string_ctx を使ってください。
  • from_hex には十六進の小数点がない。 0x3p-1 は 3⋅2−13 \cdot 2^{-1} を意味し、0x1.8p0 は拒否されます。
  • with_precision はフラグを隠す。 精度の縮小が inexact だったか、オーバーフローしたかを知る必要がある場合は round_ctx を使ってください。
  • 既知の不具合。 現在のブランチでは、acos は NaN、無限大、または ∣x∣>1|x| > 1 に対して終わりなく再帰し、atan2 は IEEE の二つの特殊ケースで符号付きゼロを誤って扱い、pow は絶対値が 2312^{31} 以上の整数指数を持つ負の底や、ほとんどの非整数指数を持つ底 −0-0 を拒否し、163/4=816^{3/4} = 8 のような厳密な結果に inexact フラグを立てます。詳細は API リファレンスにあります。

次のステップ