The standard model of floating-point arithmetic with gradual underflow

Rounding a real x onto the gradual-underflow grid of fmt satisfies

round x = x * (1 + δ) + η, with |δ| ≤ u, |η| ≤ η₀, and δ * η = 0.

For nearest rounding, u = 2^(-p) with p = fracWidth + 1, and η₀ = 2^(emin - 1) is half the subnormal spacing, where emin = minSubnormalExponent. For binary32 these are 2^(-24) and 2^(-150). Any valid rounding function, including the directed ones, satisfies the same model with 2 * u and 2 * η₀. At most one of the two error terms is nonzero: η = 0 in the normal range and δ = 0 below it. The real grid has no overflow, so the grid theorems need no range hypothesis. The executable corollaries use IEEE descriptors and require the stated domain and finiteness hypotheses.

Addition and subtraction are sharper. Below the normal range, the exact sum of two grid points is itself a grid point, so the underflow term vanishes and fl(x + y) = (x + y) * (1 + δ).

References

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Section 2.1, equations (2.8) and (2.9).
• J. W. Demmel, "Underflow and the Reliability of Numerical Software," SIAM J. Sci. Stat. Comput. 5(4), 1984. https://doi.org/10.1137/0905062
• J. R. Hauser, "Handling Floating-Point Exceptions in Numeric Programs," ACM TOPLAS 18(2), 1996. https://doi.org/10.1145/227699.227701

@[reducible, inline]

noncomputable abbrev FloatLib.Floats.Formats.BinaryInterchange.Model.unitRoundoffAt (fmt : FloatFormat) : ℝ

The unit roundoff 2^(-precision) of nearest rounding to fmt.

@[reducible, inline]

noncomputable abbrev FloatLib.Floats.Formats.BinaryInterchange.Model.underflowErrorAt (fmt : FloatFormat) : ℝ

Half the subnormal spacing of fmt, 2^(minSubnormalExponent - 1). This bounds the absolute error of nearest rounding below the normal range.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.unitRoundoffAt_eq (fmt : FloatFormat) : unitRoundoffAt fmt = 2 ^ (-(↑fmt.fracWidth + 1))

The unit roundoff as a power of two.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.underflowErrorAt_eq (fmt : FloatFormat) : underflowErrorAt fmt = 2 ^ (fmt.minSubnormalExponent - 1)

The underflow allowance as a power of two.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.two_mul_unitRoundoffAt (fmt : FloatFormat) : 2 * unitRoundoffAt fmt = Flocq.bpow Numerics.binaryRadix (1 - Int.ofNat (fmt.fracWidth + 1))

Twice the unit roundoff is the relative spacing 2^(1 - precision) of the normal grid.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.two_mul_underflowErrorAt (fmt : FloatFormat) : 2 * underflowErrorAt fmt = Flocq.bpow Numerics.binaryRadix fmt.minSubnormalExponent

Twice the underflow allowance is the subnormal spacing 2^minSubnormalExponent.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.ulpAt_eq_of_abs_lt_minNormalAt (fmt : FloatFormat) {x : ℝ} (hx : |x| < minNormalAt fmt) : ulpAt fmt x = Flocq.bpow Numerics.binaryRadix fmt.minSubnormalExponent

Below the normal range, one ULP is the fixed subnormal spacing.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.ulpAt_le_of_minNormalAt_le (fmt : FloatFormat) {x : ℝ} (hx : minNormalAt fmt ≤ |x|) : ulpAt fmt x ≤ Flocq.bpow Numerics.binaryRadix (1 - Int.ofNat (fmt.fracWidth + 1)) * |x|

In the normal range, one ULP is at most 2^(1 - precision) times the magnitude.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.round_standardModel_of_nearest (fmt : FloatFormat) (rnd : ℝ → ℤ) [Flocq.ValidRndToNearest rnd] (x : ℝ) : ∃ (δ : ℝ) (η : ℝ), Flocq.round rnd x = x * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt fmt ∧ |η| ≤ underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for any nearest rounding onto the grid of fmt, with either tie-breaking rule: round x = x * (1 + δ) + η, |δ| ≤ 2^(-precision), |η| ≤ 2^(minSubnormalExponent - 1), and δ * η = 0.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.round_standardModel (fmt : FloatFormat) (rnd : ℝ → ℤ) [Flocq.ValidRnd rnd] (x : ℝ) : ∃ (δ : ℝ) (η : ℝ), Flocq.round rnd x = x * (1 + δ) + η ∧ |δ| ≤ 2 * unitRoundoffAt fmt ∧ |η| ≤ 2 * underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for any valid rounding onto the grid of fmt, including the directed roundings: round x = x * (1 + δ) + η, |δ| ≤ 2^(1 - precision), |η| ≤ 2^minSubnormalExponent, and δ * η = 0.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.roundAt_standardModel (fmt : FloatFormat) (x : ℝ) : ∃ (δ : ℝ) (η : ℝ), roundAt fmt x = x * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt fmt ∧ |η| ≤ underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for nearest-even rounding onto the grid of fmt: roundAt fmt x = x * (1 + δ) + η, |δ| ≤ 2^(-precision), |η| ≤ 2^(minSubnormalExponent - 1), and δ * η = 0. For binary32 the bounds are 2^(-24) and 2^(-150).

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.roundAt_down_standardModel (fmt : FloatFormat) (x : ℝ) : ∃ (δ : ℝ) (η : ℝ), roundAt fmt x IEEERoundingMode.down = x * (1 + δ) + η ∧ |δ| ≤ 2 * unitRoundoffAt fmt ∧ |η| ≤ 2 * underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for downward rounding onto the grid of fmt.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.roundAt_up_standardModel (fmt : FloatFormat) (x : ℝ) : ∃ (δ : ℝ) (η : ℝ), roundAt fmt x IEEERoundingMode.up = x * (1 + δ) + η ∧ |δ| ≤ 2 * unitRoundoffAt fmt ∧ |η| ≤ 2 * underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for upward rounding onto the grid of fmt.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.round_truncRound_standardModel (fmt : FloatFormat) (x : ℝ) : ∃ (δ : ℝ) (η : ℝ), Flocq.round Flocq.truncRound x = x * (1 + δ) + η ∧ |δ| ≤ 2 * unitRoundoffAt fmt ∧ |η| ≤ 2 * underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for rounding toward zero onto the grid of fmt.

Exact addition below the normal range

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.roundAt_add_eq_mul_one_add (fmt : FloatFormat) {x y : ℝ} (hx : Flocq.genericFormat Numerics.binaryRadix (fexpOf fmt) x) (hy : Flocq.genericFormat Numerics.binaryRadix (fexpOf fmt) y) : ∃ (δ : ℝ), roundAt fmt (x + y) = (x + y) * (1 + δ) ∧ |δ| ≤ unitRoundoffAt fmt

Nearest-even rounding of a sum of two grid points has a relative error of at most 2^(-precision) and no underflow term. Below the normal range the sum is itself a grid point.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_add_eq_mul_one_add {fmt : FloatFormat} (x y : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hy : y.isFinite = true) (hout : isFinite (x.add y) = true) : ∃ (δ : ℝ), toReal (x.add y) = (x.toReal + y.toReal) * (1 + δ) ∧ |δ| ≤ unitRoundoffAt fmt

For an IEEE descriptor, addition with finite operands and result has relative error at most 2^(-precision), with no underflow term: toReal (add x y) = (toReal x + toReal y) * (1 + δ).

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_sub_eq_mul_one_add {fmt : FloatFormat} (x y : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hy : y.isFinite = true) (hout : isFinite (x.sub y) = true) : ∃ (δ : ℝ), toReal (x.sub y) = (x.toReal - y.toReal) * (1 + δ) ∧ |δ| ≤ unitRoundoffAt fmt

For an IEEE descriptor, subtraction with finite operands and result has relative error at most 2^(-precision), with no underflow term: toReal (sub x y) = (toReal x - toReal y) * (1 + δ).

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.abs_toReal_add_sub_le_unitRoundoffAt {fmt : FloatFormat} (x y : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hy : y.isFinite = true) (hout : isFinite (x.add y) = true) : |toReal (x.add y) - (x.toReal + y.toReal)| ≤ unitRoundoffAt fmt * |x.toReal + y.toReal|

For an IEEE descriptor, addition with finite operands and result has absolute error at most 2^(-precision) * |x + y|.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.abs_toReal_sub_sub_le_unitRoundoffAt {fmt : FloatFormat} (x y : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hy : y.isFinite = true) (hout : isFinite (x.sub y) = true) : |toReal (x.sub y) - (x.toReal - y.toReal)| ≤ unitRoundoffAt fmt * |x.toReal - y.toReal|

For an IEEE descriptor, subtraction with finite operands and result has absolute error at most 2^(-precision) * |x - y|.

Standard model for finite IEEE executable operations

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_mul_standardModel {fmt : FloatFormat} (x y : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hy : y.isFinite = true) (hout : isFinite (x.mul y) = true) : ∃ (δ : ℝ) (η : ℝ), toReal (x.mul y) = x.toReal * y.toReal * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt fmt ∧ |η| ≤ underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for IEEE multiplication with finite operands and result.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_div_standardModel {fmt : FloatFormat} (x y : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hy : y.isFinite = true) (hy0 : y.isZero = false) (hout : isFinite (x.div y) = true) : ∃ (δ : ℝ) (η : ℝ), toReal (x.div y) = x.toReal / y.toReal * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt fmt ∧ |η| ≤ underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for IEEE division with finite operands and result and a nonzero divisor.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_sqrt_standardModel {fmt : FloatFormat} (x : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hdomain : x.isZero = true ∨ x.signBit = false) : ∃ (δ : ℝ) (η : ℝ), toReal x.sqrt = √x.toReal * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt fmt ∧ |η| ≤ underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for IEEE square root on a finite nonnegative input, including either signed zero.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_fma_standardModel {fmt : FloatFormat} (x y z : Model fmt) (hfmt : fmt.isIEEE = true) (hx : x.isFinite = true) (hy : y.isFinite = true) (hz : z.isFinite = true) (hout : isFinite (x.fma y z) = true) : ∃ (δ : ℝ) (η : ℝ), toReal (x.fma y z) = (x.toReal * y.toReal + z.toReal) * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt fmt ∧ |η| ≤ underflowErrorAt fmt ∧ δ * η = 0

Standard model with gradual underflow for IEEE fused multiply-add with finite operands and result.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_cast_standardModel_of_isIEEE_destination {src dst : FloatFormat} (x : Model src) (hdst : dst.isIEEE = true) (hx : x.isFinite = true) (hout : isFinite (cast src dst x) = true) : ∃ (δ : ℝ) (η : ℝ), toReal (cast src dst x) = x.toReal * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt dst ∧ |η| ≤ underflowErrorAt dst ∧ δ * η = 0

Standard model with gradual underflow for a finite cast to an IEEE destination, using the constants of the destination format.

theorem FloatLib.Floats.Formats.BinaryInterchange.Model.toReal_cast_standardModel {src dst : FloatFormat} (x : Model src) (hsrc : src.isIEEE = true) (hdst : dst.isIEEE = true) (hx : x.isFinite = true) (hout : isFinite (cast src dst x) = true) : ∃ (δ : ℝ) (η : ℝ), toReal (cast src dst x) = x.toReal * (1 + δ) + η ∧ |δ| ≤ unitRoundoffAt dst ∧ |η| ≤ underflowErrorAt dst ∧ δ * η = 0

IEEE-source specialization of toReal_cast_standardModel_of_isIEEE_destination.