~/bend-docscommunity

src/math/f64.bend source

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

import Baseimport ./num.bend as Nimport ./u64.bend as Wimport ./w64.bend as Ximport ./natural.bend as M# IEEE-754 binary64 in software, on two U32 words (Bend 2 has only F32).# The algorithms are Berkeley SoftFloat 3e's (f64_add/f64_sub through# addMags/subMags, f64_mul, roundPackToF64 and normRoundPackToF64, round to# nearest with ties to even), on 64-bit significands held as two U32 limbs# (w64.bend). Division is exact long division by 32-bit quotient digits and# the square root one step of Zimmermann's Karatsuba square root with its# exact remainder, each followed by a sticky bit, so every result is correctly rounded: subnormals, overflow# to infinity, signed zeros and NaN follow IEEE 754 (x + (-x) is +0, 0/0 and# inf - inf are NaN; a NaN result is the canonical quiet NaN).##   of_bits(hi, lo), word_hi(x), word_lo(x)          the bit pattern#   is_nan, is_inf, is_finite, is_zero, signbit#   neg, abs, copysign#   add, sub, mul, div, sqrt               correctly rounded#   lt, le, eq                             IEEE comparisons (NaN unordered)#   of_nat(n)                              correctly rounded (exact below 2^53)#   f64_op, f64_is                         the num.bend instance## Exponents that may go below zero before rounding are kept as naturals# offset by OFF = 4096. Tested against the machine's doubles# (tools/check_f64.py); not proved.type F64 is Data:  Bits{lo: U32, hi: U32}def word_lo(x: F64) -> U32:  match x:    case Bits{l, h}:      ldef word_hi(x: F64) -> U32:  match x:    case Bits{l, h}:      hdef of_bits(+hi: U32, +lo: U32) -> F64:  Bits{lo, hi}def off() -> Nat:  4096ndef sgn(s: Bool) -> U32:  match s:    case True{}:      2147483648    case False{}:      0# ---- fields ----def signbit(+x: F64) -> Bool:  U32.is_le(2147483648, word_hi(x))def exp_field(+x: F64) -> Nat:  U32.to_nat(U32.and(U32.div(word_hi(x), 1048576), 2047))def frac(+x: F64) -> W.U64:  W.U64{word_lo(x), U32.and(word_hi(x), 1048575)}def mag(+x: F64) -> W.U64:  W.U64{word_lo(x), U32.and(word_hi(x), 2147483647)}def nz(+a: W.U64) -> Bool:  Bool.not(X.is_zero(a))def is_nan(+x: F64) -> Bool:  Bool.and(Nat.is_eq(exp_field(x), 2047n), nz(frac(x)))def is_inf(+x: F64) -> Bool:  Bool.and(Nat.is_eq(exp_field(x), 2047n), X.is_zero(frac(x)))def is_finite(+x: F64) -> Bool:  Nat.is_lt(exp_field(x), 2047n)def is_zero(+x: F64) -> Bool:  X.is_zero(mag(x))# ---- constants and sign ----def nan() -> F64:  Bits{0, 2146959360}def inf(s: Bool) -> F64:  Bits{0, U32.add(2146435072, sgn(s))}def zero(s: Bool) -> F64:  Bits{0, sgn(s)}def one() -> F64:  Bits{0, 1072693248}def neg(+x: F64) -> F64:  Bits{word_lo(x), U32.xor(word_hi(x), 2147483648)}def abs(+x: F64) -> F64:  Bits{word_lo(x), U32.and(word_hi(x), 2147483647)}def copysign(+x: F64, +y: F64) -> F64:  Bits{word_lo(x), U32.add(U32.and(word_hi(x), 2147483647), U32.and(word_hi(y), 2147483648))}def nan_or(+x: F64, bad: Bool) -> F64:  match bad:    case True{}:      nan()    case False{}:      x# ---- packing and rounding (SoftFloat's packToF64, roundPackToF64) ----def pack64(+w: W.U64) -> F64:  Bits{X.lo(w), X.hi(w)}# sign, exponent field e and significand m; the hidden bit of m (bit 52)# carries into the exponent field, as packToF64's addition doesdef pack(s: Bool, +e: Nat, +m: W.U64) -> F64:  pack64(X.add(W.U64{0, U32.add(sgn(s), U32.mul(U32.from_nat(e), 1048576))}, m))def zero_e(+e: Nat, z: Bool) -> Nat:  match z:    case True{}:      0n    case False{}:      edef rp_fin3(+s: Bool, +e: Nat, +r: W.U64) -> F64:  pack(s, zero_e(e, X.is_zero(r)), r)def rp_fin2(+s: Bool, +e: Nat, +r: W.U64, tie: Bool) -> F64:  rp_fin3(s, e, X.clear0(r, tie))# sig has its top bit at 62 and e is the plain exponent (0 <= e <= 0x7FD):# add half an ulp of the kept 53 bits, shift, and clear bit 0 on a tiedef rp_fin(+s: Bool, +e: Nat, +sig: W.U64) -> F64:  rp_fin2(s, e, X.shr(X.add(sig, W.U64{512, 0}), 10n), U32.is_eq(U32.and(X.lo(sig), 1023), 512))def rp_over(+s: Bool, +e: Nat, +sig: W.U64, over: Bool) -> F64:  match over:    case True{}:      inf(s)    case False{}:      rp_fin(s, Nat.sub(e, off()), sig)def rp_neg(+s: Bool, +e: Nat, +sig: W.U64, below: Bool) -> F64:  match below:    case True{}:      rp_fin(s, 0n, X.shr_jam(sig, Nat.sub(off(), e)))    case False{}:      rp_over(s, e, sig, Bool.or(Nat.is_lt(Nat.add(off(), 2045n), e), Bool.and(Nat.is_eq(e, Nat.add(off(), 2045n)), X.le(W.U64{0, 2147483648}, X.add(sig, W.U64{512, 0})))))# e is the biased exponent minus one, offset by OFF; sig has its top bit at 62def round_pack(+s: Bool, +e: Nat, +sig: W.U64) -> F64:  rp_neg(s, e, sig, Nat.is_lt(e, off()))# sig below 2^62 is shifted up oncedef rp62_pick(+s: Bool, +e: Nat, +sig: W.U64, low: Bool) -> F64:  match low:    case True{}:      round_pack(s, Nat.sub(e, 1n), X.add(sig, sig))    case False{}:      round_pack(s, e, sig)def rp62(+s: Bool, +e: Nat, +sig: W.U64) -> F64:  rp62_pick(s, e, sig, X.lt(sig, W.U64{0, 1073741824}))def nrp_exp(+e: Nat, +sd: Nat, z: Bool) -> Nat:  match z:    case True{}:      0n    case False{}:      Nat.sub(Nat.sub(e, sd), off())def nrp_pick(+s: Bool, +e: Nat, +sig: W.U64, +sd: Nat, direct: Bool) -> F64:  match direct:    case True{}:      pack(s, nrp_exp(e, sd, X.is_zero(sig)), X.shl(sig, Nat.sub(sd, 10n)))    case False{}:      round_pack(s, Nat.sub(e, sd), X.shl(sig, sd))def nrp_sd(+s: Bool, +e: Nat, +sig: W.U64, +sd: Nat) -> F64:  nrp_pick(s, e, sig, sd, Bool.and(Nat.is_le(10n, sd), Bool.and(Nat.is_le(Nat.add(off(), sd), e), Nat.is_lt(Nat.sub(e, sd), Nat.add(off(), 2045n)))))# SoftFloat's normRoundPackToF64 (e offset)def norm_round_pack(+s: Bool, +e: Nat, +sig: W.U64) -> F64:  nrp_sd(s, e, sig, Nat.sub(X.clz(sig), 1n))# the offset exponent and 53-bit significand (hidden bit at 52) of a nonzero# finite value, subnormals normalized (SoftFloat's normSubnormalF64Sig)def norm_e(+e: Nat, +f: W.U64) -> Nat:  match e:    case 0n:      Nat.sub(Nat.add(off(), 12n), X.clz(f))    case _:      Nat.add(off(), e)def norm_f(+e: Nat, +f: W.U64) -> W.U64:  match e:    case 0n:      X.shl(f, Nat.sub(X.clz(f), 11n))    case _:      X.add(f, W.U64{0, 1048576})# ---- addition (SoftFloat's addMagsF64 and subMagsF64) ----def am_small(+e: Nat, +f9: W.U64) -> W.U64:  match e:    case 0n:      X.add(f9, f9)    case _:      X.add(f9, W.U64{0, 536870912})# |x| + |y| with sign s, the larger exponent eLdef am_big(+big: F64, +s: Bool, +el: Nat, +fl: W.U64, +es: Nat, +fs: W.U64, top: Bool) -> F64:  match top:    case True{}:      nan_or(big, nz(fl))    case False{}:      rp62(s, Nat.add(off(), el), X.add(X.add(W.U64{0, 536870912}, X.shl(fl, 9n)), X.shr_jam(am_small(es, X.shl(fs, 9n)), Nat.sub(el, es))))def am_eq_n(+x: F64, +s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64, top: Bool) -> F64:  match top:    case True{}:      nan_or(x, Bool.or(nz(fa), nz(fb)))    case False{}:      round_pack(s, Nat.add(off(), e), X.shl(X.add(W.U64{0, 2097152}, X.add(fa, fb)), 9n))def am_eq(+x: F64, +s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64) -> F64:  match e:    case 0n:      pack(s, 0n, X.add(fa, fb))    case _:      am_eq_n(x, s, e, fa, fb, Nat.is_eq(e, 2047n))def am_case(+x: F64, +y: F64, +s: Bool, +ea: Nat, +eb: Nat, c: Cmp) -> F64:  match c:    case EQ{}:      am_eq(x, s, ea, frac(x), frac(y))    case LT{}:      am_big(y, s, eb, frac(y), ea, frac(x), Nat.is_eq(eb, 2047n))    case GT{}:      am_big(x, s, ea, frac(x), eb, frac(y), Nat.is_eq(ea, 2047n))def add_mags(+x: F64, +y: F64, +s: Bool) -> F64:  am_case(x, y, s, exp_field(x), exp_field(y), Nat.cmp(exp_field(x), exp_field(y)))def sm_small(+e: Nat, +f10: W.U64) -> W.U64:  match e:    case 0n:      X.add(f10, f10)    case _:      X.add(f10, W.U64{0, 1073741824})# |big| - |small| with sign s, the larger exponent eLdef sm_big(+big: F64, +s: Bool, +el: Nat, +fl: W.U64, +es: Nat, +fs: W.U64, top: Bool) -> F64:  match top:    case True{}:      nan_or(big, nz(fl))    case False{}:      norm_round_pack(s, Nat.add(off(), Nat.sub(el, 1n)), X.sub(X.add(X.shl(fl, 10n), W.U64{0, 1073741824}), X.shr_jam(sm_small(es, X.shl(fs, 10n)), Nat.sub(el, es))))def sm_exact3(+s: Bool, +e1: Nat, +d: W.U64, +sd: Nat, under: Bool) -> F64:  match under:    case True{}:      pack(s, 0n, X.shl(d, e1))    case False{}:      pack(s, Nat.sub(e1, sd), X.shl(d, sd))# an exact difference of equal exponents, normalized without roundingdef sm_exact(+s: Bool, +e1: Nat, +d: W.U64) -> F64:  sm_exact3(s, e1, d, Nat.sub(X.clz(d), 11n), Nat.is_lt(e1, Nat.sub(X.clz(d), 11n)))def sm_eq2(+s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64, c: Cmp) -> F64:  match c:    case EQ{}:      zero(False{})    case LT{}:      sm_exact(Bool.not(s), Nat.sub(e, 1n), X.sub(fb, fa))    case GT{}:      sm_exact(s, Nat.sub(e, 1n), X.sub(fa, fb))def cmp64(+a: W.U64, +b: W.U64) -> Cmp:  X.cmp(a, b)def sm_eq(+s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64, top: Bool) -> F64:  match top:    case True{}:      nan()    case False{}:      sm_eq2(s, e, fa, fb, cmp64(fa, fb))def sm_case(+x: F64, +y: F64, +s: Bool, +ea: Nat, +eb: Nat, c: Cmp) -> F64:  match c:    case EQ{}:      sm_eq(s, ea, frac(x), frac(y), Nat.is_eq(ea, 2047n))    case LT{}:      sm_big(y, Bool.not(s), eb, frac(y), ea, frac(x), Nat.is_eq(eb, 2047n))    case GT{}:      sm_big(x, s, ea, frac(x), eb, frac(y), Nat.is_eq(ea, 2047n))def sub_mags(+x: F64, +y: F64, +s: Bool) -> F64:  sm_case(x, y, s, exp_field(x), exp_field(y), Nat.cmp(exp_field(x), exp_field(y)))def add_pick(+x: F64, +y: F64, same: Bool) -> F64:  match same:    case True{}:      add_mags(x, y, signbit(x))    case False{}:      sub_mags(x, y, signbit(x))def add(+x: F64, +y: F64) -> F64:  add_pick(x, y, Bool.not(Bool.xor(signbit(x), signbit(y))))def sub(+x: F64, +y: F64) -> F64:  add(x, neg(y))# ---- multiplication (SoftFloat's f64_mul) ----def mul_n2(+s: Bool, +e: Nat, p: W.U64 & W.U64) -> F64:  (+pl, +ph) = p  rp62(s, e, X.or_bit(ph, nz(pl)))def mul_n(+s: Bool, +ea: Nat, +ma: W.U64, +eb: Nat, +mb: W.U64) -> F64:  mul_n2(s, Nat.sub(Nat.add(ea, eb), Nat.add(off(), 1023n)), X.mul128(X.shl(ma, 10n), X.shl(mb, 11n)))def fin_zero(+e: Nat, +f: W.U64) -> Bool:  Bool.and(Nat.is_eq(e, 0n), X.is_zero(f))def mul_z(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, z: Bool) -> F64:  match z:    case True{}:      zero(s)    case False{}:      mul_n(s, norm_e(ea, fa), norm_f(ea, fa), norm_e(eb, fb), norm_f(eb, fb))def inf_or_nan(+s: Bool, bad: Bool) -> F64:  match bad:    case True{}:      nan()    case False{}:      inf(s)def mul_b(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:  match top:    case True{}:      inf_or_nan(s, Bool.or(nz(fb), fin_zero(ea, fa)))    case False{}:      mul_z(s, ea, fa, eb, fb, Bool.or(fin_zero(ea, fa), fin_zero(eb, fb)))def mul_cls(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:  match top:    case True{}:      inf_or_nan(s, Bool.or(Bool.or(nz(fa), Bool.and(Nat.is_eq(eb, 2047n), nz(fb))), fin_zero(eb, fb)))    case False{}:      mul_b(s, ea, fa, eb, fb, Nat.is_eq(eb, 2047n))def mul(+x: F64, +y: F64) -> F64:  mul_cls(Bool.xor(signbit(x), signbit(y)), exp_field(x), frac(x), exp_field(y), frac(y), Nat.is_eq(exp_field(x), 2047n))# ---- division ----# the next quotient digit of r * 2^32 / b (r < b, b >= 2^52) and the remainderdef digit_fin(+r: W.U64, +b: W.U64, +d: U32) -> U32 & W.U64:  (d, X.sub(W.U64{0, X.lo(r)}, X.fst_q(X.mul_32_64(d, b))))def digit(+r: W.U64, +b: W.U64, +t: Nat) -> U32 & W.U64:  digit_fin(r, b, X.q96(W.U64{0, X.lo(r)}, X.hi(r), b, t))def dq_fin(+s: Bool, +e: Nat, +d1: U32, p: U32 & W.U64) -> F64:  (+d2, +r2) = p  round_pack(s, e, X.or_bit(X.add(W.U64{0, 1073741824}, X.shr(W.U64{d2, d1}, 2n)), Bool.or(nz(r2), Bool.not(U32.is_zero(U32.and(d2, 3))))))def dq_mid(+s: Bool, +e: Nat, +b: W.U64, +t: Nat, p: U32 & W.U64) -> F64:  (+d1, +r1) = p  dq_fin(s, e, d1, digit(r1, b, t))# a in [b, 2b): 2^62 + floor((a - b) * 2^62 / b), sticky, by two 32-bit digitsdef div_qt(+s: Bool, +e: Nat, +a: W.U64, +b: W.U64, +t: Nat) -> F64:  dq_mid(s, e, b, t, digit(X.sub(a, b), b, t))def div_q(+s: Bool, +e: Nat, +a: W.U64, +b: W.U64) -> F64:  div_qt(s, e, a, b, X.bitlen(X.hi(b)))def div_ab(+s: Bool, +ea: Nat, +a: W.U64, +eb: Nat, +b: W.U64, less: Bool) -> F64:  match less:    case True{}:      div_q(s, Nat.sub(Nat.add(ea, Nat.add(off(), 1021n)), eb), X.add(a, a), b)    case False{}:      div_q(s, Nat.sub(Nat.add(ea, Nat.add(off(), 1022n)), eb), a, b)def div_n(+s: Bool, +ea: Nat, +a: W.U64, +eb: Nat, +b: W.U64) -> F64:  div_ab(s, ea, a, eb, b, X.lt(a, b))def div_za(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, za: Bool) -> F64:  match za:    case True{}:      zero(s)    case False{}:      div_n(s, norm_e(ea, fa), norm_f(ea, fa), norm_e(eb, fb), norm_f(eb, fb))def div_z(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, zb: Bool) -> F64:  match zb:    case True{}:      inf_or_nan(s, fin_zero(ea, fa))    case False{}:      div_za(s, ea, fa, eb, fb, fin_zero(ea, fa))def div_b(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:  match top:    case True{}:      nan_or(zero(s), nz(fb))    case False{}:      div_z(s, ea, fa, eb, fb, fin_zero(eb, fb))def div_cls(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:  match top:    case True{}:      inf_or_nan(s, Bool.or(nz(fa), Nat.is_eq(eb, 2047n)))    case False{}:      div_b(s, ea, fa, eb, fb, Nat.is_eq(eb, 2047n))def div(+x: F64, +y: F64) -> F64:  div_cls(Bool.xor(signbit(x), signbit(y)), exp_field(x), frac(x), exp_field(y), frac(y), Nat.is_eq(exp_field(x), 2047n))# ---- square root ----def sq_fin(+e: Nat, +t: W.U64, sticky: Bool) -> F64:  round_pack(False{}, e, X.or_bit(t, sticky))# the root is S0 - 2: remainder (2 S0 - 3) - (D - (2 S0 - 1)) with c = 2 S0 - 1def sq_d2(+e: Nat, +s0: W.U64, +d: W.U64, +c: W.U64) -> F64:  sq_fin(e, X.sub(s0, W.U64{2, 0}), Bool.not(X.eq(X.sub(c, W.U64{2, 0}), X.sub(d, c))))# S0^2 - N = D > 0: the root is S0 - 1 when D <= 2 S0 - 1 (remainder c - D)def sq_d1(+e: Nat, +s0: W.U64, +d: W.U64, +c: W.U64, small: Bool) -> F64:  match small:    case True{}:      sq_fin(e, X.sub(s0, W.U64{1, 0}), Bool.not(X.eq(c, d)))    case False{}:      sq_d2(e, s0, d, c)def sq_neg(+e: Nat, +s0: W.U64, +d: W.U64) -> F64:  sq_d1(e, s0, d, X.sub(X.add(s0, s0), W.U64{1, 0}), X.le(d, X.sub(X.add(s0, s0), W.U64{1, 0})))# N - S0^2 = a - b for a = u 2^32, b = q^2: the root is S0 when b <= adef sq_rem(+e: Nat, +s0: W.U64, +a: W.U64, +b: W.U64, ge: Bool) -> F64:  match ge:    case True{}:      sq_fin(e, s0, Bool.not(X.eq(a, b)))    case False{}:      sq_neg(e, s0, X.sub(b, a))def sq_qu(+e: Nat, +s: U32, +q: U32, +u: U32) -> F64:  sq_rem(e, W.U64{q, s}, W.U64{0, u}, X.mul32(q, q), X.le(X.mul32(q, q), W.U64{0, u}))# q = 2^32 (r = 2 s) is taken as q = 2^32 - 1 with u = 2 sdef sq_clamp(+e: Nat, +s: U32, +q: W.U64, +u: U32, small: Bool) -> F64:  match small:    case True{}:      sq_qu(e, s, X.lo(q), u)    case False{}:      sq_qu(e, s, 4294967295, U32.add(s, s))def sq_div(+e: Nat, +s: U32, p: W.U64 & U32) -> F64:  (+q, +u) = p  sq_clamp(e, s, q, u, U32.is_zero(X.hi(q)))# one step of Zimmermann's SqrtRem (Karatsuba square root) from s = isqrt(nh):# r = nh - s^2, (q, u) = divmod(r 2^32, 2 s) give S0 = s 2^32 + q with# N - S0^2 = u 2^32 - q^2 for N = nh 2^64; S0 - 2 <= isqrt(N) <= S0 and the# exact remainder decides the root and its sticky bit without squaring S0def sq_root_s(+e: Nat, +nh: W.U64, +s: U32) -> F64:  sq_div(e, s, X.div32(W.U64{0, X.lo(X.sub(nh, X.mul32(s, s)))}, U32.add(s, s)))def sq_root(+e: Nat, +nh: W.U64) -> F64:  sq_root_s(e, nh, X.lo(X.isqrt(nh)))# m * 2^(E - OFF - 1075) with the unbiased exponent made even (E odd); the# root of m * 2^72 has its top bit at 62def sq_even(+e: Nat, +m: W.U64) -> F64:  sq_root(Nat.add(Nat.div(Nat.sub(Nat.add(e, off()), 1075n), 2n), 1048n), X.shl(m, 8n))def sq_n(+e: Nat, +m: W.U64, odd: Bool) -> F64:  match odd:    case True{}:      sq_even(e, m)    case False{}:      sq_even(Nat.sub(e, 1n), X.add(m, m))def sq_pos(+e: Nat, +f: W.U64, negative: Bool) -> F64:  match negative:    case True{}:      nan()    case False{}:      sq_n(norm_e(e, f), norm_f(e, f), Nat.is_eq(Nat.mod(norm_e(e, f), 2n), 1n))def sq_sign(+x: F64, negative: Bool, +e: Nat, +f: W.U64, zr: Bool) -> F64:  match zr:    case True{}:      x    case False{}:      sq_pos(e, f, negative)def sq_cls(+x: F64, s: Bool, +e: Nat, +f: W.U64, top: Bool) -> F64:  match top:    case True{}:      nan_or(x, Bool.or(nz(f), s))    case False{}:      sq_sign(x, s, e, f, fin_zero(e, f))def sqrt(+x: F64) -> F64:  sq_cls(x, signbit(x), exp_field(x), frac(x), Nat.is_eq(exp_field(x), 2047n))# ---- comparison and conversion ----def lt_s(+x: F64, +y: F64, sa: Bool, sb: Bool) -> Bool:  match sa sb:    case True{} False{}:      True{}    case False{} True{}:      False{}    case False{} False{}:      X.lt(mag(x), mag(y))    case True{} True{}:      X.lt(mag(y), mag(x))def le_s(+x: F64, +y: F64, sa: Bool, sb: Bool) -> Bool:  match sa sb:    case True{} False{}:      True{}    case False{} True{}:      False{}    case False{} False{}:      X.le(mag(x), mag(y))    case True{} True{}:      X.le(mag(y), mag(x))def unordered(+x: F64, +y: F64) -> Bool:  Bool.or(is_nan(x), is_nan(y))def zeros(+x: F64, +y: F64) -> Bool:  Bool.and(is_zero(x), is_zero(y))def lt_z(+x: F64, +y: F64, bad: Bool, z: Bool) -> Bool:  match bad z:    case True{} _:      False{}    case False{} True{}:      False{}    case False{} False{}:      lt_s(x, y, signbit(x), signbit(y))def lt(+x: F64, +y: F64) -> Bool:  lt_z(x, y, unordered(x, y), zeros(x, y))def le_z(+x: F64, +y: F64, bad: Bool, z: Bool) -> Bool:  match bad z:    case True{} _:      False{}    case False{} True{}:      True{}    case False{} False{}:      le_s(x, y, signbit(x), signbit(y))def le(+x: F64, +y: F64) -> Bool:  le_z(x, y, unordered(x, y), zeros(x, y))def eq_z(+x: F64, +y: F64, bad: Bool, z: Bool) -> Bool:  match bad z:    case True{} _:      False{}    case False{} True{}:      True{}    case False{} False{}:      Bool.and(U32.is_eq(word_lo(x), word_lo(y)), U32.is_eq(word_hi(x), word_hi(y)))def eq(+x: F64, +y: F64) -> Bool:  eq_z(x, y, unordered(x, y), zeros(x, y))# n mod 2^k and n div 2^k by k halvings (no 2^32 constant: the proof# checker would expand it in unary)def low_bits(k: Nat, +n: Nat) -> Nat:  match k:    case 0n:      0n    case 1n+p:      Nat.add(Nat.mod(n, 2n), Nat.double(low_bits(p, Nat.div(n, 2n))))def high_bits(k: Nat, +n: Nat) -> Nat:  match k:    case 0n:      n    case 1n+p:      high_bits(p, Nat.div(n, 2n))# the jam of n >> d (the bits shifted out OR-ed into bit 0)def jam_nat(+d: Nat, +n: Nat) -> Nat:  Nat.add(Nat.mul(2n, Nat.div(high_bits(d, n), 2n)), Nat.max(Nat.mod(high_bits(d, n), 2n), Nat.min(low_bits(d, n), 1n)))def word64(+s: Nat) -> W.U64:  W.U64{U32.from_nat(low_bits(32n, s)), U32.from_nat(high_bits(32n, s))}# n with its top bit moved to bit 62 (b = bit_length(n)): shifted up when# b <= 63 (every runtime Nat), jammed down otherwisedef sig63(+n: Nat, +b: Nat, fits: Bool) -> W.U64:  match fits:    case True{}:      X.shl(word64(n), Nat.sub(63n, b))    case False{}:      word64(jam_nat(Nat.sub(b, 63n), n))def of_nat_z(+n: Nat, z: Bool) -> F64:  match z:    case True{}:      zero(False{})    case False{}:      round_pack(False{}, Nat.add(5117n, M.bit_length(n)), sig63(n, M.bit_length(n), Nat.is_le(M.bit_length(n), 63n)))# the double nearest to n (exact below 2^53)def of_nat(+n: Nat) -> F64:  of_nat_z(n, Nat.is_eq(n, 0n))# 2^k for k <= 1023def pow2(+k: Nat) -> F64:  pack(False{}, Nat.add(1022n, k), W.U64{0, 1048576})# ---- decoding and exact rounding (the tools of everything below) ----# the significand with its hidden bit (bit 52 of a normal number) and the# scale of its unit, the exponent plus Z = 3000 (1926 for subnormals), as in# spec/math/f64.bend's mant and xexpdef dmant_z(+f: W.U64, z: Bool) -> W.U64:  match z:    case True{}:      f    case False{}:      X.add(f, W.U64{0, 1048576})def dmant(+x: F64) -> W.U64:  dmant_z(frac(x), Nat.is_eq(exp_field(x), 0n))def dexp_z(+e: Nat, z: Bool) -> Nat:  match z:    case True{}:      1926n    case False{}:      Nat.add(e, 1925n)def dexp(+x: F64) -> Nat:  dexp_z(exp_field(x), Nat.is_eq(exp_field(x), 0n))def rw_top(+s: Bool, +x: Nat, +w: W.U64, top: Bool) -> F64:  match top:    case True{}:      norm_round_pack(s, Nat.add(x, 2181n), X.shr_jam(w, 1n))    case False{}:      norm_round_pack(s, Nat.add(x, 2180n), w)def rw_z(+s: Bool, +x: Nat, +w: W.U64, z: Bool) -> F64:  match z:    case True{}:      zero(s)    case False{}:      rw_top(s, x, w, X.le(W.U64{0, 2147483648}, w))# the double nearest to (-1)^s * w * 2^(x - 3000), for x >= 63: every exact# result below is an integer of at most 64 bits at a scale, rounded oncedef round_w(+s: Bool, +x: Nat, +w: W.U64) -> F64:  rw_z(s, x, w, X.is_zero(w))# ---- rounding to an integral value (IEEE roundToIntegral) ----type RMode is Data:  Trunc{}  Floor{}  Ceil{}  Even{}def b64(b: Bool) -> W.U64:  match b:    case True{}:      W.U64{1, 0}    case False{}:      W.U64{0, 0}# whether the integer part q (remainder r, half h of the unit) moves up onedef ri_up(m: RMode, +s: Bool, +q: W.U64, +r: W.U64, +h: W.U64) -> Bool:  match m:    case Trunc{}:      False{}    case Floor{}:      Bool.and(s, Bool.not(X.is_zero(r)))    case Ceil{}:      Bool.and(Bool.not(s), Bool.not(X.is_zero(r)))    case Even{}:      Bool.or(X.lt(h, r), Bool.and(X.eq(r, h), X.odd(q)))def ri_q(m: RMode, +s: Bool, +q: W.U64, +r: W.U64, +h: W.U64) -> F64:  round_w(s, 3000n, X.add(q, b64(ri_up(m, s, q, r, h))))# k < 64 fractional bits: q = w >> k, r = w - (q << k), h = 2^(k-1)def ri_k(m: RMode, +s: Bool, +w: W.U64, +k: Nat, small: Bool) -> F64:  match small:    case True{}:      ri_q(m, s, X.shr(w, k), X.sub(w, X.shl(X.shr(w, k), k)), X.shl(W.U64{1, 0}, Nat.sub(k, 1n)))    case False{}:      ri_q(m, s, W.U64{0, 0}, w, W.U64{0, 2147483648})def ri_fin(m: RMode, +x: F64, int: Bool) -> F64:  match int:    case True{}:      x    case False{}:      ri_k(m, signbit(x), dmant(x), Nat.sub(3000n, dexp(x)), Nat.is_lt(Nat.sub(3000n, dexp(x)), 64n))def ri_cls(m: RMode, +x: F64, top: Bool) -> F64:  match top:    case True{}:      nan_or(x, nz(frac(x)))    case False{}:      ri_fin(m, x, Nat.is_le(3000n, dexp(x)))def to_integral(m: RMode, +x: F64) -> F64:  ri_cls(m, x, Nat.is_eq(exp_field(x), 2047n))def trunc(+x: F64) -> F64:  to_integral(Trunc{}, x)def floor(+x: F64) -> F64:  to_integral(Floor{}, x)def ceil(+x: F64) -> F64:  to_integral(Ceil{}, x)# ties to even (Python's round(x) and C's rint)def round(+x: F64) -> F64:  to_integral(Even{}, x)# ---- conversions to and from unsigned integers ----def tu_neg(+w: W.U64, bad: Bool) -> Result<&2, &2, N.NumError, W.U64>:  match bad:    case True{}:      Fail{N.Overflow{}}    case False{}:      Done{w}def tu_int(+s: Bool, +w: W.U64) -> Result<&2, &2, N.NumError, W.U64>:  tu_neg(w, Bool.and(s, Bool.not(X.is_zero(w))))# |x| >= 2^52: w << k with k + bit_length(w) <= 64def tu_big(+s: Bool, +w: W.U64, +k: Nat, fits: Bool) -> Result<&2, &2, N.NumError, W.U64>:  match fits:    case True{}:      tu_int(s, X.shl(w, k))    case False{}:      Fail{N.Overflow{}}def tu_fin(+x: F64, big: Bool) -> Result<&2, &2, N.NumError, W.U64>:  match big:    case True{}:      tu_big(signbit(x), dmant(x), Nat.sub(dexp(x), 3000n), Nat.is_le(Nat.add(Nat.sub(dexp(x), 3000n), Nat.sub(64n, X.clz(dmant(x)))), 64n))    case False{}:      tu_int(signbit(x), X.shr(dmant(x), Nat.sub(3000n, dexp(x))))def tu_top(isn: Bool) -> Result<&2, &2, N.NumError, W.U64>:  match isn:    case True{}:      Fail{N.BadDomain{}}    case False{}:      Fail{N.Overflow{}}def tu_cls(+x: F64, top: Bool) -> Result<&2, &2, N.NumError, W.U64>:  match top:    case True{}:      tu_top(nz(frac(x)))    case False{}:      tu_fin(x, Nat.is_le(3000n, dexp(x)))# x truncated toward zero as an unsigned 64-bit integer: NaN is a domain# error, a value outside [0, 2^64) (infinities included) an overflowdef to_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:  tu_cls(x, Nat.is_eq(exp_field(x), 2047n))def tu32_w(+w: W.U64, fits: Bool) -> Result<&2, &2, N.NumError, U32>:  match fits:    case True{}:      Done{X.lo(w)}    case False{}:      Fail{N.Overflow{}}def tu32(r: Result<&2, &2, N.NumError, W.U64>) -> Result<&2, &2, N.NumError, U32>:  match r:    case Fail{e}:      Fail{e}    case Done{+w}:      tu32_w(w, U32.is_zero(X.hi(w)))def to_u32(+x: F64) -> Result<&2, &2, N.NumError, U32>:  tu32(to_u64(x))def floor_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:  to_u64(floor(x))def ceil_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:  to_u64(ceil(x))def round_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:  to_u64(round(x))# the double nearest to an unsigned 64-bit integer (exact below 2^53)def of_u64(+w: W.U64) -> F64:  round_w(False{}, 3000n, w)def of_u32(+u: U32) -> F64:  of_u64(W.U64{u, 0})# ---- exponents: frexp, ldexp, ulp ----# a signed exponent: -mag when negtype Exp is Data:  Exp{neg: Bool, mag: Nat}def exp_pick(+t: Nat, pos: Bool) -> Exp:  match pos:    case True{}:      Exp{False{}, Nat.sub(t, 3000n)}    case False{}:      Exp{True{}, Nat.sub(3000n, t)}# t - 3000 as a signed exponentdef exp_of(+t: Nat) -> Exp:  exp_pick(t, Nat.is_le(3000n, t))def fx_fin(+x: F64, +w: W.U64, +b: Nat) -> F64 & Exp:  (round_w(signbit(x), Nat.sub(3000n, b), w), exp_of(Nat.add(dexp(x), b)))def fx_z(+x: F64, z: Bool) -> F64 & Exp:  match z:    case True{}:      (x, Exp{False{}, 0n})    case False{}:      fx_fin(x, dmant(x), Nat.sub(64n, X.clz(dmant(x))))def fx_cls(+x: F64, top: Bool) -> F64 & Exp:  match top:    case True{}:      (nan_or(x, nz(frac(x))), Exp{False{}, 0n})    case False{}:      fx_z(x, is_zero(x))# (m, e) with x = m * 2^e and 0.5 <= |m| < 1; zeros and infinities give (x, 0)def frexp(+x: F64) -> F64 & Exp:  fx_cls(x, Nat.is_eq(exp_field(x), 2047n))def ld_neg(+x: F64, +k: Nat, tiny: Bool) -> F64:  match tiny:    case True{}:      zero(signbit(x))    case False{}:      round_w(signbit(x), Nat.sub(dexp(x), k), dmant(x))def ld_fin(+x: F64, +neg: Bool, +k: Nat) -> F64:  match neg:    case True{}:      ld_neg(x, k, Nat.is_lt(dexp(x), Nat.add(k, 63n)))    case False{}:      round_w(signbit(x), Nat.add(dexp(x), k), dmant(x))def ld_e(+x: F64, +e: Exp) -> F64:  match e:    case Exp{+neg, +k}:      ld_fin(x, neg, k)def ld_z(+x: F64, +e: Exp, z: Bool) -> F64:  match z:    case True{}:      x    case False{}:      ld_e(x, e)def ld_cls(+x: F64, +e: Exp, top: Bool) -> F64:  match top:    case True{}:      nan_or(x, nz(frac(x)))    case False{}:      ld_z(x, e, is_zero(x))# x * 2^e, rounded once (overflow to infinity, gradual underflow)def ldexp(+x: F64, +e: Exp) -> F64:  ld_cls(x, e, Nat.is_eq(exp_field(x), 2047n))def ulp_cls(+x: F64, top: Bool) -> F64:  match top:    case True{}:      nan_or(inf(False{}), nz(frac(x)))    case False{}:      round_w(False{}, dexp(x), W.U64{1, 0})# the value of the least significant bit of x (Python's math.ulp)def ulp(+x: F64) -> F64:  ulp_cls(x, Nat.is_eq(exp_field(x), 2047n))# ---- neighbours and NaN-ignoring extrema ----def with_sign(+s: Bool, +w: W.U64) -> F64:  Bits{X.lo(w), U32.add(X.hi(w), sgn(s))}def na_step(+x: F64, up: Bool) -> F64:  match up:    case True{}:      with_sign(signbit(x), X.add(mag(x), W.U64{1, 0}))    case False{}:      with_sign(signbit(x), X.sub(mag(x), W.U64{1, 0}))def na_z(+x: F64, +y: F64, z: Bool) -> F64:  match z:    case True{}:      Bits{1, sgn(signbit(y))}    case False{}:      na_step(x, Bool.xor(lt(x, y), signbit(x)))def na_eq(+x: F64, +y: F64, same: Bool) -> F64:  match same:    case True{}:      y    case False{}:      na_z(x, y, is_zero(x))def na_nan(+x: F64, +y: F64, bad: Bool) -> F64:  match bad:    case True{}:      nan()    case False{}:      na_eq(x, y, eq(x, y))# the next double after x toward y (IEEE nextAfter; y when x == y)def nextafter(+x: F64, +y: F64) -> F64:  na_nan(x, y, unordered(x, y))def fpick(+x: F64, +y: F64, first: Bool) -> F64:  match first:    case True{}:      x    case False{}:      ydef fmin_z(+x: F64, +y: F64, z: Bool) -> F64:  match z:    case True{}:      zero(Bool.or(signbit(x), signbit(y)))    case False{}:      fpick(x, y, lt(x, y))def fmin_ny(+x: F64, +y: F64, ny: Bool) -> F64:  match ny:    case True{}:      x    case False{}:      fmin_z(x, y, zeros(x, y))def fmin_nx(+x: F64, +y: F64, nx: Bool) -> F64:  match nx:    case True{}:      nan_or(y, is_nan(y))    case False{}:      fmin_ny(x, y, is_nan(y))# IEEE 754-2019 minimumNumber: a NaN operand is ignored, -0 < +0def fmin(+x: F64, +y: F64) -> F64:  fmin_nx(x, y, is_nan(x))def fmax_z(+x: F64, +y: F64, z: Bool) -> F64:  match z:    case True{}:      zero(Bool.and(signbit(x), signbit(y)))    case False{}:      fpick(x, y, lt(y, x))def fmax_ny(+x: F64, +y: F64, ny: Bool) -> F64:  match ny:    case True{}:      x    case False{}:      fmax_z(x, y, zeros(x, y))def fmax_nx(+x: F64, +y: F64, nx: Bool) -> F64:  match nx:    case True{}:      nan_or(y, is_nan(y))    case False{}:      fmax_ny(x, y, is_nan(y))# IEEE 754-2019 maximumNumberdef fmax(+x: F64, +y: F64) -> F64:  fmax_nx(x, y, is_nan(x))# ---- classification, bit casts, integrality ----def is_normal(+x: F64) -> Bool:  Bool.and(Bool.not(Nat.is_eq(exp_field(x), 0n)), Nat.is_lt(exp_field(x), 2047n))def is_subnormal(+x: F64) -> Bool:  Bool.and(Nat.is_eq(exp_field(x), 0n), nz(frac(x)))def to_bits(+x: F64) -> W.U64:  W.U64{word_lo(x), word_hi(x)}def of_bits64(+w: W.U64) -> F64:  Bits{X.lo(w), X.hi(w)}def ii_k(+w: W.U64, +k: Nat, small: Bool) -> Bool:  match small:    case True{}:      X.eq(X.shl(X.shr(w, k), k), w)    case False{}:      X.is_zero(w)def ii_fin(+x: F64, int: Bool) -> Bool:  match int:    case True{}:      True{}    case False{}:      ii_k(dmant(x), Nat.sub(3000n, dexp(x)), Nat.is_lt(Nat.sub(3000n, dexp(x)), 64n))# finite with no fractional bitsdef is_integer(+x: F64) -> Bool:  Bool.and(is_finite(x), ii_fin(x, Nat.is_le(3000n, dexp(x))))# ---- the fractional part and the exact remainders ----def mf_frac(+s: Bool, +w: W.U64, +u: Nat, +k: Nat, small: Bool) -> F64:  match small:    case True{}:      round_w(s, u, X.sub(w, X.shl(X.shr(w, k), k)))    case False{}:      round_w(s, u, w)def mf_fin(+x: F64, int: Bool) -> F64 & F64:  match int:    case True{}:      (zero(signbit(x)), x)    case False{}:      (mf_frac(signbit(x), dmant(x), dexp(x), Nat.sub(3000n, dexp(x)), Nat.is_lt(Nat.sub(3000n, dexp(x)), 64n)), trunc(x))def mf_cls(+x: F64, top: Bool) -> F64 & F64:  match top:    case True{}:      (nan_or(zero(signbit(x)), nz(frac(x))), nan_or(x, nz(frac(x))))    case False{}:      mf_fin(x, Nat.is_le(3000n, dexp(x)))# (fractional part, integral part), both with the sign of x (Python's modf)def modf(+x: F64) -> F64 & F64:  mf_cls(x, Nat.is_eq(exp_field(x), 2047n))# (q, r) of mx * 2^d by b from (q0, r0) of mx by b: ten bits a step, so each# partial remainder r * 2^t stays below 2^63def fm_go(fuel: Nat, +b: W.U64, +d: Nat, st: W.U64 & W.U64) -> W.U64 & W.U64:  match fuel d:    case 0n _:      st    case 1n+f 0n:      st    case 1n+f 1n+e:      fm_go(f, b, Nat.sub(1n+e, Nat.min(1n+e, 10n)), X.divmod(X.shl(X.psnd(st), Nat.min(1n+e, 10n)), b))# the last quotient digit and remainder of mant(x) * 2^(e(x) - e(y)) by mant(y)def fm_qr(+x: F64, +y: F64) -> W.U64 & W.U64:  fm_go(Nat.sub(dexp(x), dexp(y)), dmant(y), Nat.sub(dexp(x), dexp(y)), X.divmod(dmant(x), dmant(y)))def fmod_fin(+x: F64, +y: F64, far: Bool) -> F64:  match far:    case True{}:      x    case False{}:      round_w(signbit(x), dexp(y), X.psnd(fm_qr(x, y)))def fmod_z(+x: F64, +y: F64, zx: Bool) -> F64:  match zx:    case True{}:      x    case False{}:      fmod_fin(x, y, Nat.is_lt(dexp(x), dexp(y)))def fmod_yi(+x: F64, +y: F64, yi: Bool) -> F64:  match yi:    case True{}:      x    case False{}:      fmod_z(x, y, is_zero(x))def fmod_bad(+x: F64, +y: F64, bad: Bool) -> F64:  match bad:    case True{}:      nan()    case False{}:      fmod_yi(x, y, is_inf(y))# x - n y for n = trunc(x / y), exact, with the sign of x (C's fmod)def fmod(+x: F64, +y: F64) -> F64:  fmod_bad(x, y, Bool.or(Bool.or(unordered(x, y), is_inf(x)), is_zero(y)))def rm_pick(+s: Bool, +u: Nat, +b: W.U64, +r: W.U64, flip: Bool) -> F64:  match flip:    case True{}:      round_w(Bool.not(s), u, X.sub(b, r))    case False{}:      round_w(s, u, r)# remainder r of divisor b at scale u, last quotient digit q: round the# quotient to nearest (ties to even) by flipping to r - bdef rm_fix(+s: Bool, +u: Nat, +b: W.U64, +r: W.U64, +q: W.U64) -> F64:  rm_pick(s, u, b, r, Bool.or(X.lt(b, X.add(r, r)), Bool.and(X.eq(X.add(r, r), b), X.odd(q))))def rm_qr(+x: F64, +y: F64, p: W.U64 & W.U64) -> F64:  (+q, +r) = p  rm_fix(signbit(x), dexp(y), dmant(y), r, q)def rm_near(+x: F64, +y: F64, one: Bool) -> F64:  match one:    case True{}:      rm_fix(signbit(x), dexp(x), X.add(dmant(y), dmant(y)), dmant(x), W.U64{0, 0})    case False{}:      xdef rm_far(+x: F64, +y: F64, far: Bool) -> F64:  match far:    case True{}:      rm_near(x, y, Nat.is_eq(Nat.sub(dexp(y), dexp(x)), 1n))    case False{}:      rm_qr(x, y, fm_qr(x, y))def rm_z(+x: F64, +y: F64, zx: Bool) -> F64:  match zx:    case True{}:      x    case False{}:      rm_far(x, y, Nat.is_lt(dexp(x), dexp(y)))def rm_yi(+x: F64, +y: F64, yi: Bool) -> F64:  match yi:    case True{}:      x    case False{}:      rm_z(x, y, is_zero(x))def rm_bad(+x: F64, +y: F64, bad: Bool) -> F64:  match bad:    case True{}:      nan()    case False{}:      rm_yi(x, y, is_inf(y))# IEEE remainder: x - n y for n = x / y rounded to nearest, ties to evendef remainder(+x: F64, +y: F64) -> F64:  rm_bad(x, y, Bool.or(Bool.or(unordered(x, y), is_inf(x)), is_zero(y)))# ---- closeness and exact ratios ----def ic_d(+a: F64, +b: F64, +rel: F64, +at: F64, +d: F64) -> Bool:  Bool.or(Bool.or(le(d, abs(mul(rel, b))), le(d, abs(mul(rel, a)))), le(d, at))def ic_inf(+a: F64, +b: F64, +rel: F64, +at: F64, big: Bool) -> Bool:  match big:    case True{}:      False{}    case False{}:      ic_d(a, b, rel, at, abs(sub(b, a)))def ic_eq(+a: F64, +b: F64, +rel: F64, +at: F64, same: Bool) -> Bool:  match same:    case True{}:      True{}    case False{}:      ic_inf(a, b, rel, at, Bool.or(is_inf(a), is_inf(b)))def ic_tol(+a: F64, +b: F64, +rel: F64, +at: F64, bad: Bool) -> Result<&2, &2, N.NumError, Bool>:  match bad:    case True{}:      Fail{N.BadDomain{}}    case False{}:      Done{ic_eq(a, b, rel, at, eq(a, b))}# Python's math.isclose (CPython's algorithm): a negative tolerance is a# domain errordef isclose(+a: F64, +b: F64, +rel: F64, +at: F64) -> Result<&2, &2, N.NumError, Bool>:  ic_tol(a, b, rel, at, Bool.or(lt(rel, zero(False{})), lt(at, zero(False{}))))# +-num / 2^den in lowest terms (num odd or den = 0)type Ratio is Data:  Ratio{neg: Bool, num: W.U64, den: Nat}# trailing zeros of w (at most fuel), one halving at a timedef ctz_go(fuel: Nat, +w: W.U64, odd: Bool) -> Nat:  match fuel odd:    case 0n _:      0n    case 1n+f True{}:      0n    case 1n+f False{}:      1n+ctz_go(f, X.half(w), X.odd(X.half(w)))def ctz(+w: W.U64) -> Nat:  ctz_go(64n, w, X.odd(w))def ar_small(+s: Bool, +w: W.U64, +k: Nat, +j: Nat) -> Result<&2, &2, N.NumError, Ratio>:  Done{Ratio{s, X.shr(w, j), Nat.sub(k, j)}}def ar_big(+s: Bool, +w: W.U64, +k: Nat, fits: Bool) -> Result<&2, &2, N.NumError, Ratio>:  match fits:    case True{}:      Done{Ratio{s, X.shl(w, k), 0n}}    case False{}:      Fail{N.Overflow{}}def ar_fin(+x: F64, big: Bool) -> Result<&2, &2, N.NumError, Ratio>:  match big:    case True{}:      ar_big(signbit(x), dmant(x), Nat.sub(dexp(x), 3000n), Nat.is_le(Nat.add(Nat.sub(dexp(x), 3000n), Nat.sub(64n, X.clz(dmant(x)))), 64n))    case False{}:      ar_small(signbit(x), dmant(x), Nat.sub(3000n, dexp(x)), Nat.min(ctz(dmant(x)), Nat.sub(3000n, dexp(x))))def ar_z(+x: F64, z: Bool) -> Result<&2, &2, N.NumError, Ratio>:  match z:    case True{}:      Done{Ratio{False{}, W.U64{0, 0}, 0n}}    case False{}:      ar_fin(x, Nat.is_le(3000n, dexp(x)))def ar_top(isn: Bool) -> Result<&2, &2, N.NumError, Ratio>:  match isn:    case True{}:      Fail{N.BadDomain{}}    case False{}:      Fail{N.Overflow{}}def ar_cls(+x: F64, top: Bool) -> Result<&2, &2, N.NumError, Ratio>:  match top:    case True{}:      ar_top(nz(frac(x)))    case False{}:      ar_z(x, is_zero(x))# x = +-num / 2^den exactly, in lowest terms (Python's as_integer_ratio with# the denominator's exponent); NaN is a domain error, infinities and integers# of 2^64 or more overflowdef as_integer_ratio(+x: F64) -> Result<&2, &2, N.NumError, Ratio>:  ar_cls(x, Nat.is_eq(exp_field(x), 2047n))# ---- the num.bend instance ----def f64_op(o: N.Op<F64>) -> F64:  match o:    case N.ZeroOp{}:      zero(False{})    case N.One{}:      one()    case N.Add{a, b}:      add(a, b)    case N.Sub{a, b}:      sub(a, b)    case N.Mul{a, b}:      mul(a, b)    case N.Neg{a}:      neg(a)    case N.Abs{a}:      abs(a)    case N.Quot{a, b}:      div(a, b)    case N.Rem{a, b}:      nan()    case N.Half{a}:      mul(a, Bits{0, 1071644672})    case N.MulMod{a, b, m}:      nan()    case N.Pow2{k}:      pow2(k)    case N.Sqrt{a}:      sqrt(a)    case N.PowMod{a, e, m}:      nan()    case N.GcdSmall{a, b}:      nan()def f64_is(t: N.Test<F64>) -> Bool:  match t:    case N.Lt{a, b}:      lt(a, b)    case N.AddOver{a, b}:      False{}    case N.MulOver{a, b}:      False{}    case N.Odd{a}:      False{}    case N.IsZero{a}:      is_zero(a)    case N.FastDiv{}:      True{}    case N.Mont{m}:      False{}    case N.Small{a, b}:      False{}