~/bend-docscommunity

spec/math/f64.bend source

spec/math/f64.bend on the hub · documented module

import Baseimport ../lib/common.bend as Cimport ../../src/math/f64.bend as Fimport ../../src/math/natural.bend as Mimport ../../src/math/u64.bend as WUimport ../../src/math/num.bend as Nimport ./w64.bend as SW# Specification of src/math/f64.bend: IEEE 754-2019 binary64 with round to# nearest, ties to even, written as Flocq writes it (Boldo & Melquiond,# Flocq's Binary / Round_NE / Bracket developments; the same shape as the# HOL Light and Isabelle IEEE_Floating_Point specifications): a finite# double denotes the exact dyadic (-1)^s * m * 2^e, every operation computes# its exact result on dyadics (or rationals, with a sticky bit) and rounds# it once (round below). Special values follow IEEE 754 section 6 and the# design reference python_style_math_stdlib_design.pdf (3.3, 4.4; appendix# A): NaN in gives the canonical quiet NaN out, x + (-x) is +0 under round# to nearest, inf - inf, 0 * inf, 0 / 0, inf / inf and sqrt of a negative# number are NaN, x / 0 is a signed infinity, sqrt(-0) is -0, comparisons# are false on NaN and +0 == -0.## PROVED: every clause below is proved for every input in# proofs/math/typed/f64*.bend (no holes, no axioms). TESTED as well:# tools/check_f64.py and tools/check_f64x.py test the implementation# against the machine's IEEE doubles and CPython's math module on random bit# patterns (zeros, subnormals, infinities, NaNs, rounding ties, appendix A's# special values), and tools/check_f64_spec.py tests the arithmetic part of# this reference itself against them through a line-by-line mirror.##   function            clauses#   add, sub            Add.value, Sub.value#   mul, div, sqrt      Mul.value, Div.value, Sqrt.value#   lt, le, eq          Lt.value, Le.value, Eq.value#   neg, abs, copysign  Neg.value, Abs.value, Copysign.value#   is_nan, is_inf,     IsNan.value, IsInf.value, IsFinite.value,#   is_finite, is_zero, IsZero.value, Signbit.value#   signbit#   of_nat              OfNat.value## Rounding, conversions and the rest of the design reference's float# foundation (python_style_math_stdlib_design.pdf 3.1, 3.3, 4.3, 4.4, 10 and# appendix A; IEEE 754-2019 5.3.1, 5.3.3, 5.8, 9.6; C99 annex F):##   trunc, floor, ceil, round   Trunc.value, Floor.value, Ceil.value,#                               Round.value (roundToIntegral, ties to even)#   to_u64, to_u32              ToU64.value, ToU32.value#   floor_u64, ceil_u64,        FloorU64.value, CeilU64.value,#   round_u64                   RoundU64.value#   of_u64, of_u32              OfU64.value, OfU32.value#   frexp, ldexp, ulp           Frexp.value, Ldexp.value, Ulp.value#   nextafter                   Nextafter.value#   fmin, fmax                  Fmin.value, Fmax.value (minimumNumber,#                               maximumNumber)#   is_normal, is_subnormal     IsNormal.value, IsSubnormal.value#   to_bits, of_bits64          Bits.value, Bits.roundtrip, Bits.inverse#   is_integer                  IsInteger.value#   modf                        Modf.value#   fmod, remainder             Fmod.value, Remainder.value (exact)#   isclose                     IsClose.value#   as_integer_ratio            AsIntegerRatio.value# ---- the encoding ----# The checker evaluates a closed Nat term such as 2^52 in unary, so no such# constant appears here: widths are stated with C.shift (x * 2^k), C.low# (n mod 2^k), C.high (n div 2^k) and C.fits (n < 2^k) on symbolic# arguments, and the special bit patterns are U32 words. Each form is the# plain arithmetic one (tools/check_f64_spec.py mirrors it line by line).def hi_nat(x: F.F64) -> Nat:  U32.to_nat(F.word_hi(x))def lo_nat(x: F.F64) -> Nat:  U32.to_nat(F.word_lo(x))# bit 63, bits 52..62 and bits 0..51 of the patterndef sign(+x: F.F64) -> Bool:  Bool.not(C.fits(31n, hi_nat(x)))def efield(+x: F.F64) -> Nat:  C.low(11n, C.high(20n, hi_nat(x)))def frac(+x: F.F64) -> Nat:  Nat.add(lo_nat(x), C.shift(32n, C.low(20n, hi_nat(x))))def b2n(b: Bool) -> Nat:  match b:    case True{}:      1n    case False{}:      0n# the double with sign s, exponent field ef < 2^11 and fraction f < 2^52:# (s * 2^11 + ef) * 2^52 + fdef encode(+s: Bool, +ef: Nat, +f: Nat) -> F.F64:  F.Bits{U32.from_nat(C.low(32n, f)), U32.from_nat(Nat.add(C.high(32n, f), C.shift(20n, Nat.add(ef, C.shift(11n, b2n(s))))))}def pick(-T: Data, c: Bool, +a: T, +b: T) -> T:  match c:    case True{}:      a    case False{}:      b# pick for pairs (a Sigma is Type-sorted)def pickt(-T: Type, c: Bool, a: T, b: T) -> T:  match c:    case True{}:      a    case False{}:      b# encode(False, 2047, 2^51), encode(s, 2047, 0) and encode(s, 0, 0)def qnan() -> F.F64:  F.Bits{0, 2146959360}def inf(+s: Bool) -> F.F64:  pick(F.F64, s, F.Bits{0, 4293918720}, F.Bits{0, 2146435072})def zero(+s: Bool) -> F.F64:  pick(F.F64, s, F.Bits{0, 2147483648}, F.Bits{0, 0})def is_nan(+x: F.F64) -> Bool:  Bool.and(Nat.is_eq(efield(x), 2047n), Bool.not(Nat.is_eq(frac(x), 0n)))def is_inf(+x: F.F64) -> Bool:  Bool.and(Nat.is_eq(efield(x), 2047n), Nat.is_eq(frac(x), 0n))def is_zero(+x: F.F64) -> Bool:  Bool.and(Nat.is_eq(efield(x), 0n), Nat.is_eq(frac(x), 0n))# ---- the exact value of a finite double: m * 2^(xe - Z) ----# Z keeps every exponent a natural number: the smallest (subnormal) is# -1074 = 1926 - Z, and products and quotients stay above -Z.def zb() -> Nat:  3000n# the fraction with the hidden bit 2^52 of a normal numberdef mant(+x: F.F64) -> Nat:  Nat.add(frac(x), C.shift(52n, b2n(Bool.not(Nat.is_eq(efield(x), 0n)))))def xexp(+x: F.F64) -> Nat:  pick(Nat, Nat.is_eq(efield(x), 0n), Nat.sub(zb(), 1074n), Nat.sub(Nat.add(efield(x), zb()), 1075n))# ---- rounding (Flocq's round_NE on the binary64 format) ----# m / 2^k rounded to nearest, ties to even: the quotient q, the remainder# r and half the divisor hdef rne_up(+q: Nat, +r: Nat, +h: Nat) -> Nat:  Nat.add(q, b2n(Bool.or(Nat.is_lt(h, r), Bool.and(Nat.is_eq(r, h), Nat.is_eq(Nat.mod(q, 2n), 1n)))))def rne(+m: Nat, +k: Nat) -> Nat:  match k:    case 0n:      m    case 1n+ +j:      rne_up(C.high(1n+j, m), C.low(1n+j, m), C.shift(j, 1n))# a normal number (ef < 2047) or an overflow to infinitydef pack_e(+s: Bool, +ef: Nat, +f: Nat) -> F.F64:  pick(F.F64, Nat.is_le(2047n, ef), inf(s), encode(s, ef, f))# 2^52 <= q < 2^53 is normal with fraction q - 2^52, q < 2^52 subnormaldef pack_n(+s: Bool, +q: Nat, +u: Nat) -> F.F64:  pick(F.F64, C.fits(52n, q), encode(s, 0n, q), pack_e(s, Nat.sub(Nat.add(u, 1075n), zb()), C.low(52n, q)))# the double nearest to (-1)^s * q * 2^(u - Z) for a q of at most 53 bits at# the ulp u (subnormal when q < 2^52), carrying a rounded-up q = 2^53 into# 2^52 at the ulp 1 + udef pack(+s: Bool, +q: Nat, +u: Nat) -> F.F64:  pick(F.F64, C.fits(53n, q), pack_n(s, q, u), pack_e(s, Nat.sub(Nat.add(1n+u, 1075n), zb()), 0n))# round (-1)^s * m * 2^(x - Z): the ulp is 2^(msb - 52), but never below# the subnormal ulp 2^-1074def round_u(+s: Bool, +m: Nat, +x: Nat, +u: Nat) -> F.F64:  pack(s, pick(Nat, Nat.is_le(x, u), rne(m, Nat.sub(u, x)), C.shift(Nat.sub(x, u), m)), u)def round(+s: Bool, +m: Nat, +x: Nat) -> F.F64:  pick(F.F64, Nat.is_eq(m, 0n), zero(s), round_u(s, m, x, Nat.max(Nat.sub(Nat.add(x, M.bit_length(m)), 53n), Nat.sub(zb(), 1074n))))# ---- the operations ----# (-1)^sa A + (-1)^sb B at a common scale x: equal magnitudes of opposite# signs cancel to +0 (round to nearest)def sub_mag(+sa: Bool, +a: Nat, +sb: Bool, +b: Nat, +x: Nat, c: Cmp) -> F.F64:  match c:    case GT{}:      round(sa, Nat.sub(a, b), x)    case LT{}:      round(sb, Nat.sub(b, a), x)    case EQ{}:      zero(False{})def add_mag(+sa: Bool, +a: Nat, +sb: Bool, +b: Nat, +x: Nat) -> F.F64:  pick(F.F64, Bool.not(Bool.xor(sa, sb)), round(sa, Nat.add(a, b), x), sub_mag(sa, a, sb, b, x, Nat.cmp(a, b)))def add_fin(+x: F.F64, +y: F.F64) -> F.F64:  add_mag(sign(x), C.shift(Nat.sub(xexp(x), Nat.min(xexp(x), xexp(y))), mant(x)), sign(y), C.shift(Nat.sub(xexp(y), Nat.min(xexp(x), xexp(y))), mant(y)), Nat.min(xexp(x), xexp(y)))def add_inf(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, is_inf(x), pick(F.F64, Bool.and(is_inf(y), Bool.not(Bool.not(Bool.xor(sign(x), sign(y))))), qnan(), x), pick(F.F64, is_inf(y), y, add_fin(x, y)))def add(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, Bool.or(is_nan(x), is_nan(y)), qnan(), add_inf(x, y))# the double with the sign bit flipped (NaN included)def neg(+x: F.F64) -> F.F64:  encode(Bool.not(sign(x)), efield(x), frac(x))def mul_inf(+x: F.F64, +y: F.F64, +s: Bool) -> F.F64:  pick(F.F64, Bool.or(is_inf(x), is_inf(y)), pick(F.F64, Bool.or(is_zero(x), is_zero(y)), qnan(), inf(s)), round(s, Nat.mul(mant(x), mant(y)), Nat.sub(Nat.add(xexp(x), xexp(y)), zb())))def mul(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, Bool.or(is_nan(x), is_nan(y)), qnan(), mul_inf(x, y, Bool.xor(sign(x), sign(y))))# enough quotient bits for any pair of significands, plus a sticky bitdef kq() -> Nat:  200ndef div_fin(+x: F.F64, +y: F.F64, +s: Bool) -> F.F64:  round(s, Nat.add(Nat.mul(2n, Nat.div(C.shift(kq(), mant(x)), mant(y))), Nat.min(Nat.mod(C.shift(kq(), mant(x)), mant(y)), 1n)), Nat.sub(Nat.add(xexp(x), zb()), Nat.add(Nat.add(xexp(y), kq()), 1n)))def div_cls(+x: F.F64, +y: F.F64, +s: Bool) -> F.F64:  pick(F.F64, is_inf(x), pick(F.F64, is_inf(y), qnan(), inf(s)), pick(F.F64, is_inf(y), zero(s), pick(F.F64, is_zero(y), pick(F.F64, is_zero(x), qnan(), inf(s)), div_fin(x, y, s))))def div(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, Bool.or(is_nan(x), is_nan(y)), qnan(), div_cls(x, y, Bool.xor(sign(x), sign(y))))# sqrt(m * 2^(x - Z)) with x - Z made even: the integer square root of the# significand scaled by 2^(2 K), and a sticky bit when it was not exactdef kr() -> Nat:  100ndef sqrt_even(+m: Nat, +x: Nat) -> F.F64:  round(False{}, Nat.add(Nat.mul(2n, M.isqrt(C.shift(Nat.mul(2n, kr()), m))), b2n(Bool.not(Nat.is_eq(Nat.pow(M.isqrt(C.shift(Nat.mul(2n, kr()), m)), 2n), C.shift(Nat.mul(2n, kr()), m))))), Nat.sub(Nat.add(Nat.div(x, 2n), Nat.div(zb(), 2n)), Nat.add(kr(), 1n)))def sqrt_fin(+x: F.F64) -> F.F64:  pick(F.F64, Nat.is_eq(Nat.mod(xexp(x), 2n), 0n), sqrt_even(mant(x), xexp(x)), sqrt_even(Nat.mul(2n, mant(x)), Nat.sub(xexp(x), 1n)))def sqrt_cls(+x: F.F64) -> F.F64:  pick(F.F64, is_zero(x), x, pick(F.F64, sign(x), qnan(), pick(F.F64, is_inf(x), x, sqrt_fin(x))))def sqrt(+x: F.F64) -> F.F64:  pick(F.F64, is_nan(x), qnan(), sqrt_cls(x))# ---- comparisons: the order of the extended reals, false on NaN ----# the order of two non-NaN doubles: infinities at the ends, zeros equaldef mag_cmp(+x: F.F64, +y: F.F64) -> Cmp:  pick(Cmp, is_inf(x), pick(Cmp, is_inf(y), EQ{}, GT{}), pick(Cmp, is_inf(y), LT{}, Nat.cmp(C.shift(Nat.sub(xexp(x), Nat.min(xexp(x), xexp(y))), mant(x)), C.shift(Nat.sub(xexp(y), Nat.min(xexp(x), xexp(y))), mant(y)))))def flip(c: Cmp) -> Cmp:  match c:    case LT{}:      GT{}    case EQ{}:      EQ{}    case GT{}:      LT{}def ord(+x: F.F64, +y: F.F64) -> Cmp:  pick(Cmp, Bool.and(is_zero(x), is_zero(y)), EQ{}, pick(Cmp, Bool.not(Bool.xor(sign(x), sign(y))), pick(Cmp, sign(x), flip(mag_cmp(x, y)), mag_cmp(x, y)), pick(Cmp, sign(x), LT{}, GT{})))def ordered(+x: F.F64, +y: F.F64) -> Bool:  Bool.not(Bool.or(is_nan(x), is_nan(y)))# ---- the contract ----def Add.value(+x: F.F64, +y: F.F64) -> Type:  {F.add(x, y) == add(x, y) : F.F64}# x - y is x + (-y)def Sub.value(+x: F.F64, +y: F.F64) -> Type:  {F.sub(x, y) == add(x, neg(y)) : F.F64}def Mul.value(+x: F.F64, +y: F.F64) -> Type:  {F.mul(x, y) == mul(x, y) : F.F64}def Div.value(+x: F.F64, +y: F.F64) -> Type:  {F.div(x, y) == div(x, y) : F.F64}def Sqrt.value(+x: F.F64) -> Type:  {F.sqrt(x) == sqrt(x) : F.F64}def Lt.value(+x: F.F64, +y: F.F64) -> Type:  {F.lt(x, y) == Bool.and(ordered(x, y), Cmp.is_lt(ord(x, y))) : Bool}def Le.value(+x: F.F64, +y: F.F64) -> Type:  {F.le(x, y) == Bool.and(ordered(x, y), Cmp.is_le(ord(x, y))) : Bool}def Eq.value(+x: F.F64, +y: F.F64) -> Type:  {F.eq(x, y) == Bool.and(ordered(x, y), Cmp.is_eq(ord(x, y))) : Bool}def Neg.value(+x: F.F64) -> Type:  {F.neg(x) == neg(x) : F.F64}def Abs.value(+x: F.F64) -> Type:  {F.abs(x) == encode(False{}, efield(x), frac(x)) : F.F64}def Copysign.value(+x: F.F64, +y: F.F64) -> Type:  {F.copysign(x, y) == encode(sign(y), efield(x), frac(x)) : F.F64}def IsNan.value(+x: F.F64) -> Type:  {F.is_nan(x) == is_nan(x) : Bool}def IsInf.value(+x: F.F64) -> Type:  {F.is_inf(x) == is_inf(x) : Bool}def IsFinite.value(+x: F.F64) -> Type:  {F.is_finite(x) == Bool.not(Nat.is_eq(efield(x), 2047n)) : Bool}def IsZero.value(+x: F.F64) -> Type:  {F.is_zero(x) == is_zero(x) : Bool}def Signbit.value(+x: F.F64) -> Type:  {F.signbit(x) == sign(x) : Bool}# exact below 2^53, rounded above (the runtime bounds Nat by 2^48)def OfNat.value(+n: Nat) -> Type:  {F.of_nat(n) == round(False{}, n, zb()) : F.F64}# ---- rounding to an integral value (IEEE 754 roundToIntegral) ----# the integer magnitude of (-1)^s * m * 2^-k, k >= 1, in each direction: the# integer part high(k, m), plus one when the discarded low(k, m) is nonzero# and the direction is away from zero (floor of a negative, ceil of a# positive), or to nearest with ties to even (rne)def integral(m: F.RMode, +s: Bool, +n: Nat, +k: Nat) -> Nat:  match m:    case F.Trunc{}:      C.high(k, n)    case F.Floor{}:      Nat.add(C.high(k, n), b2n(Bool.and(s, Bool.not(Nat.is_eq(C.low(k, n), 0n)))))    case F.Ceil{}:      Nat.add(C.high(k, n), b2n(Bool.and(Bool.not(s), Bool.not(Nat.is_eq(C.low(k, n), 0n)))))    case F.Even{}:      rne(n, k)# NaN gives NaN, an infinity or a value of at least 2^52 is already integral,# anything else is its integer (with the sign of x, so -0.5 goes to -0 under# trunc, ceil and round)def to_integral(m: F.RMode, +x: F.F64) -> F.F64:  pick(F.F64, is_nan(x), qnan(), pick(F.F64, Bool.or(is_inf(x), Nat.is_le(zb(), xexp(x))), x, round(sign(x), integral(m, sign(x), mant(x), Nat.sub(zb(), xexp(x))), zb())))def Trunc.value(+x: F.F64) -> Type:  {F.trunc(x) == to_integral(F.Trunc{}, x) : F.F64}def Floor.value(+x: F.F64) -> Type:  {F.floor(x) == to_integral(F.Floor{}, x) : F.F64}def Ceil.value(+x: F.F64) -> Type:  {F.ceil(x) == to_integral(F.Ceil{}, x) : F.F64}def Round.value(+x: F.F64) -> Type:  {F.round(x) == to_integral(F.Even{}, x) : F.F64}# ---- conversions ----# the integer part of |x| (x finite)def int_part(+x: F.F64) -> Nat:  pick(Nat, Nat.is_le(zb(), xexp(x)), C.shift(Nat.sub(xexp(x), zb()), mant(x)), C.high(Nat.sub(zb(), xexp(x)), mant(x)))# x truncated toward zero as a w-bit unsigned integer: NaN is a domain error;# an infinity, a truncation below zero or one of 2^w or more an overflow# (3.1: int(x) truncates; 2.6: fixed widths raise instead of wrapping)def to_nat(+x: F.F64, +w: Nat) -> Result<&2, &2, N.NumError, Nat>:  pick(Result<&2, &2, N.NumError, Nat>, is_nan(x), Fail{N.BadDomain{}}, pick(Result<&2, &2, N.NumError, Nat>, Bool.or(is_inf(x), Bool.or(Bool.and(sign(x), Bool.not(Nat.is_eq(int_part(x), 0n))), Bool.not(C.fits(w, int_part(x))))), Fail{N.Overflow{}}, Done{int_part(x)}))def rv64(r: Result<&2, &2, N.NumError, WU.U64>) -> Result<&2, &2, N.NumError, Nat>:  match r:    case Fail{e}:      Fail{e}    case Done{w}:      Done{SW.value(w)}def rv32(r: Result<&2, &2, N.NumError, U32>) -> Result<&2, &2, N.NumError, Nat>:  match r:    case Fail{e}:      Fail{e}    case Done{u}:      Done{U32.to_nat(u)}def ToU64.value(+x: F.F64) -> Type:  {rv64(F.to_u64(x)) == to_nat(x, 64n) : Result<&2, &2, N.NumError, Nat>}def ToU32.value(+x: F.F64) -> Type:  {rv32(F.to_u32(x)) == to_nat(x, 32n) : Result<&2, &2, N.NumError, Nat>}def FloorU64.value(+x: F.F64) -> Type:  {rv64(F.floor_u64(x)) == to_nat(to_integral(F.Floor{}, x), 64n) : Result<&2, &2, N.NumError, Nat>}def CeilU64.value(+x: F.F64) -> Type:  {rv64(F.ceil_u64(x)) == to_nat(to_integral(F.Ceil{}, x), 64n) : Result<&2, &2, N.NumError, Nat>}def RoundU64.value(+x: F.F64) -> Type:  {rv64(F.round_u64(x)) == to_nat(to_integral(F.Even{}, x), 64n) : Result<&2, &2, N.NumError, Nat>}# the double nearest to an unsigned integer (exact below 2^53)def OfU64.value(+w: WU.U64) -> Type:  {F.of_u64(w) == round(False{}, SW.value(w), zb()) : F.F64}def OfU32.value(+u: U32) -> Type:  {F.of_u32(u) == round(False{}, U32.to_nat(u), zb()) : F.F64}# ---- exponents ----# t - Z as a signed exponentdef exp_of(+t: Nat) -> F.Exp:  pick(F.Exp, Nat.is_le(zb(), t), F.Exp{False{}, Nat.sub(t, zb())}, F.Exp{True{}, Nat.sub(zb(), t)})# (m, e) with x = m 2^e and 1/2 <= |m| < 1: m is mant(x) scaled by# 2^-bit_length(mant(x)) (exact), e the rest of the scale; NaN gives (NaN, 0),# zeros and infinities (x, 0) (4.4)def frexp(+x: F.F64) -> F.F64 & F.Exp:  pickt(F.F64 & F.Exp, is_nan(x), (qnan(), F.Exp{False{}, 0n}), pickt(F.F64 & F.Exp, Bool.or(is_inf(x), is_zero(x)), (x, F.Exp{False{}, 0n}), (round(sign(x), mant(x), Nat.sub(zb(), M.bit_length(mant(x)))), exp_of(Nat.add(xexp(x), M.bit_length(mant(x)))))))def Frexp.value(+x: F.F64) -> Type:  {F.frexp(x) == frexp(x) : F.F64 & F.Exp}# x 2^e rounded once: overflow to infinity, gradual underflow (4.4); a# scale below 2^-Z is a product under 2^-2947, which rounds to a zero like# the scale 0 it is cut todef ldexp(+x: F.F64, +neg: Bool, +k: Nat) -> F.F64:  pick(F.F64, is_nan(x), qnan(), pick(F.F64, Bool.or(is_inf(x), is_zero(x)), x, round(sign(x), mant(x), pick(Nat, neg, Nat.sub(xexp(x), k), Nat.add(xexp(x), k)))))def Ldexp.value(+x: F.F64, +neg: Bool, +k: Nat) -> Type:  {F.ldexp(x, F.Exp{neg, k}) == ldexp(x, neg, k) : F.F64}# the weight of the last bit of x: 2^(xexp(x) - Z) (the smallest subnormal# for zeros), inf for infinities, NaN for NaN (4.4)def Ulp.value(+x: F.F64) -> Type:  {F.ulp(x) == pick(F.F64, is_nan(x), qnan(), pick(F.F64, is_inf(x), inf(False{}), round(False{}, 1n, xexp(x)))) : F.F64}# ---- neighbours and NaN-ignoring extrema ----# the magnitude bits: exponent field and fraction as one integer below 2^63def pat(+x: F.F64) -> Nat:  Nat.add(frac(x), C.shift(52n, efield(x)))def of_pat(+s: Bool, +p: Nat) -> F.F64:  encode(s, C.high(52n, p), C.low(52n, p))# the ordered-set neighbour (IEEE nextAfter, C99 annex F): NaN in, NaN out;# y when x == y; from a zero the smallest subnormal with y's sign; otherwise# one step of the magnitude bits, up when moving away from zero (the step past# the largest finite double is infinity, the step below the smallest# subnormal a zero of x's sign)def nextafter(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, Bool.not(ordered(x, y)), qnan(), pick(F.F64, Cmp.is_eq(ord(x, y)), y, pick(F.F64, is_zero(x), encode(sign(y), 0n, 1n), of_pat(sign(x), pick(Nat, Bool.xor(Cmp.is_lt(ord(x, y)), sign(x)), Nat.add(pat(x), 1n), Nat.sub(pat(x), 1n))))))def Nextafter.value(+x: F.F64, +y: F.F64) -> Type:  {F.nextafter(x, y) == nextafter(x, y) : F.F64}# IEEE 754-2019 minimumNumber / maximumNumber (10: fmin/fmax fix the# order dependence of min/max with NaN): a NaN operand is ignored, -0 < +0def fmin(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, is_nan(x), pick(F.F64, is_nan(y), qnan(), y), pick(F.F64, is_nan(y), x, pick(F.F64, Bool.and(is_zero(x), is_zero(y)), zero(Bool.or(sign(x), sign(y))), pick(F.F64, Cmp.is_lt(ord(x, y)), x, y))))def fmax(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, is_nan(x), pick(F.F64, is_nan(y), qnan(), y), pick(F.F64, is_nan(y), x, pick(F.F64, Bool.and(is_zero(x), is_zero(y)), zero(Bool.and(sign(x), sign(y))), pick(F.F64, Cmp.is_lt(ord(y, x)), x, y))))def Fmin.value(+x: F.F64, +y: F.F64) -> Type:  {F.fmin(x, y) == fmin(x, y) : F.F64}def Fmax.value(+x: F.F64, +y: F.F64) -> Type:  {F.fmax(x, y) == fmax(x, y) : F.F64}# ---- classification, bit casts, integrality ----def IsNormal.value(+x: F.F64) -> Type:  {F.is_normal(x) == Bool.and(Bool.not(Nat.is_eq(efield(x), 0n)), Nat.is_lt(efield(x), 2047n)) : Bool}def IsSubnormal.value(+x: F.F64) -> Type:  {F.is_subnormal(x) == Bool.and(Nat.is_eq(efield(x), 0n), Bool.not(Nat.is_eq(frac(x), 0n))) : Bool}# the 64-bit pattern: sign, exponent field, fractiondef Bits.value(+x: F.F64) -> Type:  {SW.value(F.to_bits(x)) == Nat.add(pat(x), C.shift(63n, b2n(sign(x)))) : Nat}def Bits.roundtrip(+x: F.F64) -> Type:  {F.of_bits64(F.to_bits(x)) == x : F.F64}def Bits.inverse(+w: WU.U64) -> Type:  {F.to_bits(F.of_bits64(w)) == w : WU.U64}# finite with no fractional bitsdef IsInteger.value(+x: F.F64) -> Type:  {F.is_integer(x) == Bool.and(Nat.is_lt(efield(x), 2047n), Bool.or(Nat.is_le(zb(), xexp(x)), Nat.is_eq(C.low(Nat.sub(zb(), xexp(x)), mant(x)), 0n))) : Bool}# ---- the fractional part and the exact remainders ----# (fractional part, integral part), both with the sign of x, both exact:# modf(-2.0) = (-0.0, -2.0), modf(inf) = (0.0, inf) (4.3)def modf(+x: F.F64) -> F.F64 & F.F64:  pickt(F.F64 & F.F64, is_nan(x), (qnan(), qnan()), pickt(F.F64 & F.F64, Bool.or(is_inf(x), Nat.is_le(zb(), xexp(x))), (zero(sign(x)), x), (round(sign(x), C.low(Nat.sub(zb(), xexp(x)), mant(x)), xexp(x)), to_integral(F.Trunc{}, x))))def Modf.value(+x: F.F64) -> Type:  {F.modf(x) == modf(x) : F.F64 & F.F64}# the significands of x and y at their common scaledef rsc(+x: F.F64, +y: F.F64) -> Nat:  Nat.min(xexp(x), xexp(y))def ra(+x: F.F64, +y: F.F64) -> Nat:  C.shift(Nat.sub(xexp(x), rsc(x, y)), mant(x))def rb(+x: F.F64, +y: F.F64) -> Nat:  C.shift(Nat.sub(xexp(y), rsc(x, y)), mant(y))def rbad(+x: F.F64, +y: F.F64) -> Bool:  Bool.or(Bool.or(Bool.not(ordered(x, y)), is_inf(x)), is_zero(y))# x - n y with n = trunc(x / y): the exact remainder of the scaled# significands, with the sign of x; NaN for NaN, an infinite x or a zero y;# x for an infinite y (4.3: fmod(x, inf) == x)def fmod(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, rbad(x, y), qnan(), pick(F.F64, is_inf(y), x, round(sign(x), Nat.mod(ra(x, y), rb(x, y)), rsc(x, y))))def Fmod.value(+x: F.F64, +y: F.F64) -> Type:  {F.fmod(x, y) == fmod(x, y) : F.F64}# IEEE remainder, x - n y with n = x / y rounded to nearest, ties to even:# from r = A mod B and q = A div B, n is q + 1 when r is past half of B (or# exactly half with q odd), which leaves B - r with the opposite sign; a zero# result has the sign of x (4.3)def rflip(+r: Nat, +b: Nat, +q: Nat) -> Bool:  Bool.or(Nat.is_lt(b, Nat.add(r, r)), Bool.and(Nat.is_eq(Nat.add(r, r), b), Nat.is_eq(Nat.mod(q, 2n), 1n)))def rnear(+s: Bool, +r: Nat, +b: Nat, +q: Nat, +u: Nat) -> F.F64:  pick(F.F64, rflip(r, b, q), round(Bool.not(s), Nat.sub(b, r), u), round(s, r, u))def remainder(+x: F.F64, +y: F.F64) -> F.F64:  pick(F.F64, rbad(x, y), qnan(), pick(F.F64, is_inf(y), x, rnear(sign(x), Nat.mod(ra(x, y), rb(x, y)), rb(x, y), Nat.div(ra(x, y), rb(x, y)), rsc(x, y))))def Remainder.value(+x: F.F64, +y: F.F64) -> Type:  {F.remainder(x, y) == remainder(x, y) : F.F64}# ---- closeness and exact ratios ----def lt_s(+x: F.F64, +y: F.F64) -> Bool:  Bool.and(ordered(x, y), Cmp.is_lt(ord(x, y)))def le_s(+x: F.F64, +y: F.F64) -> Bool:  Bool.and(ordered(x, y), Cmp.is_le(ord(x, y)))def eq_s(+x: F.F64, +y: F.F64) -> Bool:  Bool.and(ordered(x, y), Cmp.is_eq(ord(x, y)))def fabs(+x: F.F64) -> F.F64:  encode(False{}, efield(x), frac(x))# CPython's math.isclose, computed in binary64: a negative tolerance is a# domain error; a == b is close, an infinity is close only to itself, and# otherwise |b - a| <= max(|rel * b|, |rel * a|, abs_tol) (4.4)def ic_near(+a: F.F64, +b: F.F64, +rel: F.F64, +at: F.F64, +d: F.F64) -> Bool:  Bool.or(Bool.or(le_s(d, fabs(mul(rel, b))), le_s(d, fabs(mul(rel, a)))), le_s(d, at))def isclose(+a: F.F64, +b: F.F64, +rel: F.F64, +at: F.F64) -> Result<&2, &2, N.NumError, Bool>:  pick(Result<&2, &2, N.NumError, Bool>, Bool.or(lt_s(rel, zero(False{})), lt_s(at, zero(False{}))), Fail{N.BadDomain{}}, Done{pick(Bool, eq_s(a, b), True{}, pick(Bool, Bool.or(is_inf(a), is_inf(b)), False{}, ic_near(a, b, rel, at, fabs(add(b, neg(a))))))})def IsClose.value(+a: F.F64, +b: F.F64, +rel: F.F64, +at: F.F64) -> Type:  {F.isclose(a, b, rel, at) == isclose(a, b, rel, at) : Result<&2, &2, N.NumError, Bool>}# trailing zero bits of n (at most fuel)def tz(fuel: Nat, +n: Nat) -> Nat:  match fuel:    case 0n:      0n    case 1n+f:      pick(Nat, Nat.is_eq(Nat.mod(n, 2n), 1n), 0n, 1n+tz(f, Nat.div(n, 2n)))# (sign, numerator, exponent of the denominator)type Rat is Data:  Rat{neg: Bool, num: Nat, den: Nat}def rtriple(r: Result<&2, &2, N.NumError, F.Ratio>) -> Result<&2, &2, N.NumError, Rat>:  match r:    case Fail{e}:      Fail{e}    case Done{F.Ratio{neg, num, den}}:      Done{Rat{neg, SW.value(num), den}}# x = +-num / 2^den in lowest terms: the significand with as many of its# trailing zeros moved into the denominator as the scale allows; a zero is# 0/1 (Python's ints have no -0); NaN is a domain error, an infinity or an# integer of 2^64 or more an overflow (3.3)def as_ratio(+x: F.F64) -> Result<&2, &2, N.NumError, Rat>:  pick(Result<&2, &2, N.NumError, Rat>, is_nan(x), Fail{N.BadDomain{}}, pick(Result<&2, &2, N.NumError, Rat>, is_inf(x), Fail{N.Overflow{}}, pick(Result<&2, &2, N.NumError, Rat>, is_zero(x), Done{Rat{False{}, 0n, 0n}}, pick(Result<&2, &2, N.NumError, Rat>, Nat.is_le(zb(), xexp(x)), pick(Result<&2, &2, N.NumError, Rat>, C.fits(64n, int_part(x)), Done{Rat{sign(x), int_part(x), 0n}}, Fail{N.Overflow{}}), Done{Rat{sign(x), C.high(Nat.min(tz(64n, mant(x)), Nat.sub(zb(), xexp(x))), mant(x)), Nat.sub(Nat.sub(zb(), xexp(x)), Nat.min(tz(64n, mant(x)), Nat.sub(zb(), xexp(x))))}}))))def AsIntegerRatio.value(+x: F.F64) -> Type:  {rtriple(F.as_integer_ratio(x)) == as_ratio(x) : Result<&2, &2, N.NumError, Rat>}