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 是一个闭区间 [x‾,x‾][\underline{x}, \overline{x}],其两个端点都是 BinFloat 值。本页的示例用下面这个辅助函数打印区间:它把下端点向下舍入、上端点向上舍入后输出,因此打印出的十进制区间仍然包含所存储的区间。

///|
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 + "]"
}

最小的实用程序以 53 位精度包络 1/31/3:

///|
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() 读取端点,用 width() 读取向上舍入的 x‾−x‾\overline{x} - \underline{x},用 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")
}

当调用方必须知道结果是否由一个在整个输入上都有定义且连续的函数所支撑时,就使用装饰——这正是区间结果能够证明诸如 Brouwer 不动点定理之类存在性定理的条件。如果只关心包络,直接使用 BallFloat 更简单。

基于包络 trait 的泛型代码

BallFloat 实现了 Luna-Flow/arithmetic 的包络关系(Contains、Overlaps、DefinitelyLt、DefinitelyLe、MaybeEq),因此针对这些 trait 编写的代码适用于任何包络类型。

///|
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 次,因此其代价是一次乘法的数倍。提高精度会缩小宽度中由舍入造成的部分,但不会缩小来自输入宽度的部分。

常见陷阱

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(或 trait 方法 @lf_arith.Contains::contains,其含义是“作为子集包含”)。

区间之间没有序。 区间上不存在全序;切勿用 less(一种 IEEE 1788 集合关系)对区间排序,把它当作数的比较。

无界区间没有中心。 center()、radius() 和 midpoint() 在空集和半无界区间上会中止(midpoint() 对整条实轴返回 0)。请先检查 is_bounded(),或使用 radius_extended(),它在这些情形下返回 +∞+\infty。

比精度更宽的整数。 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 测试语料覆盖的范围。