bin_float 教程

本教程教你使用 BinFloat——floating 的任意精度二进制浮点类型——进行计算:构造精确值,在 binary32 等 IEEE 754 格式下连同其状态标志执行计算,读写交换位模式和十进制文本,以及调用正确舍入的初等函数。每个示例都能编译并展示其输出。舍入规则背后的数学见设计页面,所有函数都列在 API 参考中。

快速入门

在你的 moon.pkg 中添加模块并导入该包:

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

最小的实用程序计算一除以三,一次用 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 是 c⋅2ec \cdot 2^{e},其中系数 cc 为整数,ee 为指数。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 读取十进制输入,它对十进制值本身进行舍入。当你想要 53 位以外的精度时,经由宿主 Double 转换会舍入两次,而 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,因此信号 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 的静默 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) 恰好等于一,函数也知道这一点:精确结果不会引发 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 允许实现在舍入之前或之后检测微小性(tininess),这一选择会改变略小于最小正规数的值的下溢标志。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 trait,以及 @lf_arith 的 checked 与 contextual trait。针对这些 trait 编写的代码可以在 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 trait 把 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 会隐藏标志。 当你需要知道收窄是否不精确或上溢时,请使用 round_ctx。
  • 已知缺陷。 在当前分支上,acos 在输入为 NaN、无穷或 ∣x∣>1|x| > 1 时会无限递归;atan2 在两种 IEEE 特殊情形下对带符号零处理有误;pow 会拒绝带有绝对值至少为 2312^{31} 的整数指数的负底数,以及底数为 −0-0 且指数为大多数非整数的情形,并把 163/4=816^{3/4} = 8 这类精确结果标记为不精确。API 参考列出了详细信息。

后续步骤