~/jendhub

spec/math/f64.bend

bend-collections-laws-crypto@1.0.0.0 · 592 lines · 124 definitions

1import Base
2import bend-collections-laws-crypto@1.0.0.0/spec/lib/common.bend as C
3import bend-collections-laws-crypto@1.0.0.0/src/math/f64.bend as F
4import bend-collections-laws-crypto@1.0.0.0/src/math/natural.bend as M
5import bend-collections-laws-crypto@1.0.0.0/src/math/u64.bend as WU
6import bend-collections-laws-crypto@1.0.0.0/src/math/num.bend as N
7import bend-collections-laws-crypto@1.0.0.0/spec/math/w64.bend as SW
8
9# Specification of src/math/f64.bend: IEEE 754-2019 binary64 with round to
10# nearest, ties to even, written as Flocq writes it (Boldo & Melquiond,
11# Flocq's Binary / Round_NE / Bracket developments; the same shape as the
12# HOL Light and Isabelle IEEE_Floating_Point specifications): a finite
13# double denotes the exact dyadic (-1)^s * m * 2^e, every operation computes
14# its exact result on dyadics (or rationals, with a sticky bit) and rounds
15# it once (round below). Special values follow IEEE 754 section 6 and the
16# design reference python_style_math_stdlib_design.pdf (3.3, 4.4; appendix
17# A): NaN in gives the canonical quiet NaN out, x + (-x) is +0 under round
18# to nearest, inf - inf, 0 * inf, 0 / 0, inf / inf and sqrt of a negative
19# number are NaN, x / 0 is a signed infinity, sqrt(-0) is -0, comparisons
20# are false on NaN and +0 == -0.
21#
22# PROVED: every clause below is proved for every input in
23# proofs/math/typed/f64*.bend (no holes, no axioms). TESTED as well:
24# tools/check_f64.py and tools/check_f64x.py test the implementation
25# against the machine's IEEE doubles and CPython's math module on random bit
26# patterns (zeros, subnormals, infinities, NaNs, rounding ties, appendix A's
27# special values), and tools/check_f64_spec.py tests the arithmetic part of
28# this reference itself against them through a line-by-line mirror.
29#
30#   function            clauses
31#   add, sub            Add.value, Sub.value
32#   mul, div, sqrt      Mul.value, Div.value, Sqrt.value
33#   lt, le, eq          Lt.value, Le.value, Eq.value
34#   neg, abs, copysign  Neg.value, Abs.value, Copysign.value
35#   is_nan, is_inf,     IsNan.value, IsInf.value, IsFinite.value,
36#   is_finite, is_zero, IsZero.value, Signbit.value
37#   signbit
38#   of_nat              OfNat.value
39#
40# Rounding, conversions and the rest of the design reference's float
41# foundation (python_style_math_stdlib_design.pdf 3.1, 3.3, 4.3, 4.4, 10 and
42# appendix A; IEEE 754-2019 5.3.1, 5.3.3, 5.8, 9.6; C99 annex F):
43#
44#   trunc, floor, ceil, round   Trunc.value, Floor.value, Ceil.value,
45#                               Round.value (roundToIntegral, ties to even)
46#   to_u64, to_u32              ToU64.value, ToU32.value
47#   floor_u64, ceil_u64,        FloorU64.value, CeilU64.value,
48#   round_u64                   RoundU64.value
49#   of_u64, of_u32              OfU64.value, OfU32.value
50#   frexp, ldexp, ulp           Frexp.value, Ldexp.value, Ulp.value
51#   nextafter                   Nextafter.value
52#   fmin, fmax                  Fmin.value, Fmax.value (minimumNumber,
53#                               maximumNumber)
54#   is_normal, is_subnormal     IsNormal.value, IsSubnormal.value
55#   to_bits, of_bits64          Bits.value, Bits.roundtrip, Bits.inverse
56#   is_integer                  IsInteger.value
57#   modf                        Modf.value
58#   fmod, remainder             Fmod.value, Remainder.value (exact)
59#   isclose                     IsClose.value
60#   as_integer_ratio            AsIntegerRatio.value
61
62# ---- the encoding ----
63
64# The checker evaluates a closed Nat term such as 2^52 in unary, so no such
65# constant appears here: widths are stated with C.shift (x * 2^k), C.low
66# (n mod 2^k), C.high (n div 2^k) and C.fits (n < 2^k) on symbolic
67# arguments, and the special bit patterns are U32 words. Each form is the
68# plain arithmetic one (tools/check_f64_spec.py mirrors it line by line).
69
70def hi_nat(x: F.F64) -> Nat:
71  U32.to_nat(F.word_hi(x))
72
73def lo_nat(x: F.F64) -> Nat:
74  U32.to_nat(F.word_lo(x))
75
76# bit 63, bits 52..62 and bits 0..51 of the pattern
77def sign(+x: F.F64) -> Bool:
78  Bool.not(C.fits(31n, hi_nat(x)))
79
80def efield(+x: F.F64) -> Nat:
81  C.low(11n, C.high(20n, hi_nat(x)))
82
83def frac(+x: F.F64) -> Nat:
84  Nat.add(lo_nat(x), C.shift(32n, C.low(20n, hi_nat(x))))
85
86def b2n(b: Bool) -> Nat:
87  match b:
88    case True{}:
89      1n
90    case False{}:
91      0n
92
93# the double with sign s, exponent field ef < 2^11 and fraction f < 2^52:
94# (s * 2^11 + ef) * 2^52 + f
95def encode(+s: Bool, +ef: Nat, +f: Nat) -> F.F64:
96  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))))))}
97
98def pick(-T: Data, c: Bool, +a: T, +b: T) -> T:
99  match c:
100    case True{}:
101      a
102    case False{}:
103      b
104
105# pick for pairs (a Sigma is Type-sorted)
106def pickt(-T: Type, c: Bool, a: T, b: T) -> T:
107  match c:
108    case True{}:
109      a
110    case False{}:
111      b
112
113# encode(False, 2047, 2^51), encode(s, 2047, 0) and encode(s, 0, 0)
114def qnan() -> F.F64:
115  F.Bits{0, 2146959360}
116
117def inf(+s: Bool) -> F.F64:
118  pick(F.F64, s, F.Bits{0, 4293918720}, F.Bits{0, 2146435072})
119
120def zero(+s: Bool) -> F.F64:
121  pick(F.F64, s, F.Bits{0, 2147483648}, F.Bits{0, 0})
122
123def is_nan(+x: F.F64) -> Bool:
124  Bool.and(Nat.is_eq(efield(x), 2047n), Bool.not(Nat.is_eq(frac(x), 0n)))
125
126def is_inf(+x: F.F64) -> Bool:
127  Bool.and(Nat.is_eq(efield(x), 2047n), Nat.is_eq(frac(x), 0n))
128
129def is_zero(+x: F.F64) -> Bool:
130  Bool.and(Nat.is_eq(efield(x), 0n), Nat.is_eq(frac(x), 0n))
131
132# ---- the exact value of a finite double: m * 2^(xe - Z) ----
133# Z keeps every exponent a natural number: the smallest (subnormal) is
134# -1074 = 1926 - Z, and products and quotients stay above -Z.
135
136def zb() -> Nat:
137  3000n
138
139# the fraction with the hidden bit 2^52 of a normal number
140def mant(+x: F.F64) -> Nat:
141  Nat.add(frac(x), C.shift(52n, b2n(Bool.not(Nat.is_eq(efield(x), 0n)))))
142
143def xexp(+x: F.F64) -> Nat:
144  pick(Nat, Nat.is_eq(efield(x), 0n), Nat.sub(zb(), 1074n), Nat.sub(Nat.add(efield(x), zb()), 1075n))
145
146# ---- rounding (Flocq's round_NE on the binary64 format) ----
147
148# m / 2^k rounded to nearest, ties to even: the quotient q, the remainder
149# r and half the divisor h
150def rne_up(+q: Nat, +r: Nat, +h: Nat) -> Nat:
151  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)))))
152
153def rne(+m: Nat, +k: Nat) -> Nat:
154  match k:
155    case 0n:
156      m
157    case 1n+ +j:
158      rne_up(C.high(1n+j, m), C.low(1n+j, m), C.shift(j, 1n))
159
160# a normal number (ef < 2047) or an overflow to infinity
161def pack_e(+s: Bool, +ef: Nat, +f: Nat) -> F.F64:
162  pick(F.F64, Nat.is_le(2047n, ef), inf(s), encode(s, ef, f))
163
164# 2^52 <= q < 2^53 is normal with fraction q - 2^52, q < 2^52 subnormal
165def pack_n(+s: Bool, +q: Nat, +u: Nat) -> F.F64:
166  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)))
167
168# the double nearest to (-1)^s * q * 2^(u - Z) for a q of at most 53 bits at
169# the ulp u (subnormal when q < 2^52), carrying a rounded-up q = 2^53 into
170# 2^52 at the ulp 1 + u
171def pack(+s: Bool, +q: Nat, +u: Nat) -> F.F64:
172  pick(F.F64, C.fits(53n, q), pack_n(s, q, u), pack_e(s, Nat.sub(Nat.add(1n+u, 1075n), zb()), 0n))
173
174# round (-1)^s * m * 2^(x - Z): the ulp is 2^(msb - 52), but never below
175# the subnormal ulp 2^-1074
176def round_u(+s: Bool, +m: Nat, +x: Nat, +u: Nat) -> F.F64:
177  pack(s, pick(Nat, Nat.is_le(x, u), rne(m, Nat.sub(u, x)), C.shift(Nat.sub(x, u), m)), u)
178
179def round(+s: Bool, +m: Nat, +x: Nat) -> F.F64:
180  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))))
181
182# ---- the operations ----
183
184# (-1)^sa A + (-1)^sb B at a common scale x: equal magnitudes of opposite
185# signs cancel to +0 (round to nearest)
186def sub_mag(+sa: Bool, +a: Nat, +sb: Bool, +b: Nat, +x: Nat, c: Cmp) -> F.F64:
187  match c:
188    case GT{}:
189      round(sa, Nat.sub(a, b), x)
190    case LT{}:
191      round(sb, Nat.sub(b, a), x)
192    case EQ{}:
193      zero(False{})
194
195def add_mag(+sa: Bool, +a: Nat, +sb: Bool, +b: Nat, +x: Nat) -> F.F64:
196  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)))
197
198def add_fin(+x: F.F64, +y: F.F64) -> F.F64:
199  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)))
200
201def add_inf(+x: F.F64, +y: F.F64) -> F.F64:
202  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)))
203
204def add(+x: F.F64, +y: F.F64) -> F.F64:
205  pick(F.F64, Bool.or(is_nan(x), is_nan(y)), qnan(), add_inf(x, y))
206
207# the double with the sign bit flipped (NaN included)
208def neg(+x: F.F64) -> F.F64:
209  encode(Bool.not(sign(x)), efield(x), frac(x))
210
211def mul_inf(+x: F.F64, +y: F.F64, +s: Bool) -> F.F64:
212  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())))
213
214def mul(+x: F.F64, +y: F.F64) -> F.F64:
215  pick(F.F64, Bool.or(is_nan(x), is_nan(y)), qnan(), mul_inf(x, y, Bool.xor(sign(x), sign(y))))
216
217# enough quotient bits for any pair of significands, plus a sticky bit
218def kq() -> Nat:
219  200n
220
221def div_fin(+x: F.F64, +y: F.F64, +s: Bool) -> F.F64:
222  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)))
223
224def div_cls(+x: F.F64, +y: F.F64, +s: Bool) -> F.F64:
225  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))))
226
227def div(+x: F.F64, +y: F.F64) -> F.F64:
228  pick(F.F64, Bool.or(is_nan(x), is_nan(y)), qnan(), div_cls(x, y, Bool.xor(sign(x), sign(y))))
229
230# sqrt(m * 2^(x - Z)) with x - Z made even: the integer square root of the
231# significand scaled by 2^(2 K), and a sticky bit when it was not exact
232def kr() -> Nat:
233  100n
234
235def sqrt_even(+m: Nat, +x: Nat) -> F.F64:
236  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)))
237
238def sqrt_fin(+x: F.F64) -> F.F64:
239  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)))
240
241def sqrt_cls(+x: F.F64) -> F.F64:
242  pick(F.F64, is_zero(x), x, pick(F.F64, sign(x), qnan(), pick(F.F64, is_inf(x), x, sqrt_fin(x))))
243
244def sqrt(+x: F.F64) -> F.F64:
245  pick(F.F64, is_nan(x), qnan(), sqrt_cls(x))
246
247# ---- comparisons: the order of the extended reals, false on NaN ----
248
249# the order of two non-NaN doubles: infinities at the ends, zeros equal
250def mag_cmp(+x: F.F64, +y: F.F64) -> Cmp:
251  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)))))
252
253def flip(c: Cmp) -> Cmp:
254  match c:
255    case LT{}:
256      GT{}
257    case EQ{}:
258      EQ{}
259    case GT{}:
260      LT{}
261
262def ord(+x: F.F64, +y: F.F64) -> Cmp:
263  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{})))
264
265def ordered(+x: F.F64, +y: F.F64) -> Bool:
266  Bool.not(Bool.or(is_nan(x), is_nan(y)))
267
268# ---- the contract ----
269
270def Add.value(+x: F.F64, +y: F.F64) -> Type:
271  {F.add(x, y) == add(x, y) : F.F64}
272
273# x - y is x + (-y)
274def Sub.value(+x: F.F64, +y: F.F64) -> Type:
275  {F.sub(x, y) == add(x, neg(y)) : F.F64}
276
277def Mul.value(+x: F.F64, +y: F.F64) -> Type:
278  {F.mul(x, y) == mul(x, y) : F.F64}
279
280def Div.value(+x: F.F64, +y: F.F64) -> Type:
281  {F.div(x, y) == div(x, y) : F.F64}
282
283def Sqrt.value(+x: F.F64) -> Type:
284  {F.sqrt(x) == sqrt(x) : F.F64}
285
286def Lt.value(+x: F.F64, +y: F.F64) -> Type:
287  {F.lt(x, y) == Bool.and(ordered(x, y), Cmp.is_lt(ord(x, y))) : Bool}
288
289def Le.value(+x: F.F64, +y: F.F64) -> Type:
290  {F.le(x, y) == Bool.and(ordered(x, y), Cmp.is_le(ord(x, y))) : Bool}
291
292def Eq.value(+x: F.F64, +y: F.F64) -> Type:
293  {F.eq(x, y) == Bool.and(ordered(x, y), Cmp.is_eq(ord(x, y))) : Bool}
294
295def Neg.value(+x: F.F64) -> Type:
296  {F.neg(x) == neg(x) : F.F64}
297
298def Abs.value(+x: F.F64) -> Type:
299  {F.abs(x) == encode(False{}, efield(x), frac(x)) : F.F64}
300
301def Copysign.value(+x: F.F64, +y: F.F64) -> Type:
302  {F.copysign(x, y) == encode(sign(y), efield(x), frac(x)) : F.F64}
303
304def IsNan.value(+x: F.F64) -> Type:
305  {F.is_nan(x) == is_nan(x) : Bool}
306
307def IsInf.value(+x: F.F64) -> Type:
308  {F.is_inf(x) == is_inf(x) : Bool}
309
310def IsFinite.value(+x: F.F64) -> Type:
311  {F.is_finite(x) == Bool.not(Nat.is_eq(efield(x), 2047n)) : Bool}
312
313def IsZero.value(+x: F.F64) -> Type:
314  {F.is_zero(x) == is_zero(x) : Bool}
315
316def Signbit.value(+x: F.F64) -> Type:
317  {F.signbit(x) == sign(x) : Bool}
318
319# exact below 2^53, rounded above (the runtime bounds Nat by 2^48)
320def OfNat.value(+n: Nat) -> Type:
321  {F.of_nat(n) == round(False{}, n, zb()) : F.F64}
322
323# ---- rounding to an integral value (IEEE 754 roundToIntegral) ----
324
325# the integer magnitude of (-1)^s * m * 2^-k, k >= 1, in each direction: the
326# integer part high(k, m), plus one when the discarded low(k, m) is nonzero
327# and the direction is away from zero (floor of a negative, ceil of a
328# positive), or to nearest with ties to even (rne)
329def integral(m: F.RMode, +s: Bool, +n: Nat, +k: Nat) -> Nat:
330  match m:
331    case F.Trunc{}:
332      C.high(k, n)
333    case F.Floor{}:
334      Nat.add(C.high(k, n), b2n(Bool.and(s, Bool.not(Nat.is_eq(C.low(k, n), 0n)))))
335    case F.Ceil{}:
336      Nat.add(C.high(k, n), b2n(Bool.and(Bool.not(s), Bool.not(Nat.is_eq(C.low(k, n), 0n)))))
337    case F.Even{}:
338      rne(n, k)
339
340# NaN gives NaN, an infinity or a value of at least 2^52 is already integral,
341# anything else is its integer (with the sign of x, so -0.5 goes to -0 under
342# trunc, ceil and round)
343def to_integral(m: F.RMode, +x: F.F64) -> F.F64:
344  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())))
345
346def Trunc.value(+x: F.F64) -> Type:
347  {F.trunc(x) == to_integral(F.Trunc{}, x) : F.F64}
348
349def Floor.value(+x: F.F64) -> Type:
350  {F.floor(x) == to_integral(F.Floor{}, x) : F.F64}
351
352def Ceil.value(+x: F.F64) -> Type:
353  {F.ceil(x) == to_integral(F.Ceil{}, x) : F.F64}
354
355def Round.value(+x: F.F64) -> Type:
356  {F.round(x) == to_integral(F.Even{}, x) : F.F64}
357
358# ---- conversions ----
359
360# the integer part of |x| (x finite)
361def int_part(+x: F.F64) -> Nat:
362  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)))
363
364# x truncated toward zero as a w-bit unsigned integer: NaN is a domain error;
365# an infinity, a truncation below zero or one of 2^w or more an overflow
366# (3.1: int(x) truncates; 2.6: fixed widths raise instead of wrapping)
367def to_nat(+x: F.F64, +w: Nat) -> Result<&2, &2, N.NumError, Nat>:
368  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)}))
369
370def rv64(r: Result<&2, &2, N.NumError, WU.U64>) -> Result<&2, &2, N.NumError, Nat>:
371  match r:
372    case Fail{e}:
373      Fail{e}
374    case Done{w}:
375      Done{SW.value(w)}
376
377def rv32(r: Result<&2, &2, N.NumError, U32>) -> Result<&2, &2, N.NumError, Nat>:
378  match r:
379    case Fail{e}:
380      Fail{e}
381    case Done{u}:
382      Done{U32.to_nat(u)}
383
384def ToU64.value(+x: F.F64) -> Type:
385  {rv64(F.to_u64(x)) == to_nat(x, 64n) : Result<&2, &2, N.NumError, Nat>}
386
387def ToU32.value(+x: F.F64) -> Type:
388  {rv32(F.to_u32(x)) == to_nat(x, 32n) : Result<&2, &2, N.NumError, Nat>}
389
390def FloorU64.value(+x: F.F64) -> Type:
391  {rv64(F.floor_u64(x)) == to_nat(to_integral(F.Floor{}, x), 64n) : Result<&2, &2, N.NumError, Nat>}
392
393def CeilU64.value(+x: F.F64) -> Type:
394  {rv64(F.ceil_u64(x)) == to_nat(to_integral(F.Ceil{}, x), 64n) : Result<&2, &2, N.NumError, Nat>}
395
396def RoundU64.value(+x: F.F64) -> Type:
397  {rv64(F.round_u64(x)) == to_nat(to_integral(F.Even{}, x), 64n) : Result<&2, &2, N.NumError, Nat>}
398
399# the double nearest to an unsigned integer (exact below 2^53)
400def OfU64.value(+w: WU.U64) -> Type:
401  {F.of_u64(w) == round(False{}, SW.value(w), zb()) : F.F64}
402
403def OfU32.value(+u: U32) -> Type:
404  {F.of_u32(u) == round(False{}, U32.to_nat(u), zb()) : F.F64}
405
406# ---- exponents ----
407
408# t - Z as a signed exponent
409def exp_of(+t: Nat) -> F.Exp:
410  pick(F.Exp, Nat.is_le(zb(), t), F.Exp{False{}, Nat.sub(t, zb())}, F.Exp{True{}, Nat.sub(zb(), t)})
411
412# (m, e) with x = m 2^e and 1/2 <= |m| < 1: m is mant(x) scaled by
413# 2^-bit_length(mant(x)) (exact), e the rest of the scale; NaN gives (NaN, 0),
414# zeros and infinities (x, 0) (4.4)
415def frexp(+x: F.F64) -> F.F64 & F.Exp:
416  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)))))))
417
418def Frexp.value(+x: F.F64) -> Type:
419  {F.frexp(x) == frexp(x) : F.F64 & F.Exp}
420
421# x 2^e rounded once: overflow to infinity, gradual underflow (4.4); a
422# scale below 2^-Z is a product under 2^-2947, which rounds to a zero like
423# the scale 0 it is cut to
424def ldexp(+x: F.F64, +neg: Bool, +k: Nat) -> F.F64:
425  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)))))
426
427def Ldexp.value(+x: F.F64, +neg: Bool, +k: Nat) -> Type:
428  {F.ldexp(x, F.Exp{neg, k}) == ldexp(x, neg, k) : F.F64}
429
430# the weight of the last bit of x: 2^(xexp(x) - Z) (the smallest subnormal
431# for zeros), inf for infinities, NaN for NaN (4.4)
432def Ulp.value(+x: F.F64) -> Type:
433  {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}
434
435# ---- neighbours and NaN-ignoring extrema ----
436
437# the magnitude bits: exponent field and fraction as one integer below 2^63
438def pat(+x: F.F64) -> Nat:
439  Nat.add(frac(x), C.shift(52n, efield(x)))
440
441def of_pat(+s: Bool, +p: Nat) -> F.F64:
442  encode(s, C.high(52n, p), C.low(52n, p))
443
444# the ordered-set neighbour (IEEE nextAfter, C99 annex F): NaN in, NaN out;
445# y when x == y; from a zero the smallest subnormal with y's sign; otherwise
446# one step of the magnitude bits, up when moving away from zero (the step past
447# the largest finite double is infinity, the step below the smallest
448# subnormal a zero of x's sign)
449def nextafter(+x: F.F64, +y: F.F64) -> F.F64:
450  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))))))
451
452def Nextafter.value(+x: F.F64, +y: F.F64) -> Type:
453  {F.nextafter(x, y) == nextafter(x, y) : F.F64}
454
455# IEEE 754-2019 minimumNumber / maximumNumber (10: fmin/fmax fix the
456# order dependence of min/max with NaN): a NaN operand is ignored, -0 < +0
457def fmin(+x: F.F64, +y: F.F64) -> F.F64:
458  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))))
459
460def fmax(+x: F.F64, +y: F.F64) -> F.F64:
461  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))))
462
463def Fmin.value(+x: F.F64, +y: F.F64) -> Type:
464  {F.fmin(x, y) == fmin(x, y) : F.F64}
465
466def Fmax.value(+x: F.F64, +y: F.F64) -> Type:
467  {F.fmax(x, y) == fmax(x, y) : F.F64}
468
469# ---- classification, bit casts, integrality ----
470
471def IsNormal.value(+x: F.F64) -> Type:
472  {F.is_normal(x) == Bool.and(Bool.not(Nat.is_eq(efield(x), 0n)), Nat.is_lt(efield(x), 2047n)) : Bool}
473
474def IsSubnormal.value(+x: F.F64) -> Type:
475  {F.is_subnormal(x) == Bool.and(Nat.is_eq(efield(x), 0n), Bool.not(Nat.is_eq(frac(x), 0n))) : Bool}
476
477# the 64-bit pattern: sign, exponent field, fraction
478def Bits.value(+x: F.F64) -> Type:
479  {SW.value(F.to_bits(x)) == Nat.add(pat(x), C.shift(63n, b2n(sign(x)))) : Nat}
480
481def Bits.roundtrip(+x: F.F64) -> Type:
482  {F.of_bits64(F.to_bits(x)) == x : F.F64}
483
484def Bits.inverse(+w: WU.U64) -> Type:
485  {F.to_bits(F.of_bits64(w)) == w : WU.U64}
486
487# finite with no fractional bits
488def IsInteger.value(+x: F.F64) -> Type:
489  {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}
490
491# ---- the fractional part and the exact remainders ----
492
493# (fractional part, integral part), both with the sign of x, both exact:
494# modf(-2.0) = (-0.0, -2.0), modf(inf) = (0.0, inf) (4.3)
495def modf(+x: F.F64) -> F.F64 & F.F64:
496  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))))
497
498def Modf.value(+x: F.F64) -> Type:
499  {F.modf(x) == modf(x) : F.F64 & F.F64}
500
501# the significands of x and y at their common scale
502def rsc(+x: F.F64, +y: F.F64) -> Nat:
503  Nat.min(xexp(x), xexp(y))
504
505def ra(+x: F.F64, +y: F.F64) -> Nat:
506  C.shift(Nat.sub(xexp(x), rsc(x, y)), mant(x))
507
508def rb(+x: F.F64, +y: F.F64) -> Nat:
509  C.shift(Nat.sub(xexp(y), rsc(x, y)), mant(y))
510
511def rbad(+x: F.F64, +y: F.F64) -> Bool:
512  Bool.or(Bool.or(Bool.not(ordered(x, y)), is_inf(x)), is_zero(y))
513
514# x - n y with n = trunc(x / y): the exact remainder of the scaled
515# significands, with the sign of x; NaN for NaN, an infinite x or a zero y;
516# x for an infinite y (4.3: fmod(x, inf) == x)
517def fmod(+x: F.F64, +y: F.F64) -> F.F64:
518  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))))
519
520def Fmod.value(+x: F.F64, +y: F.F64) -> Type:
521  {F.fmod(x, y) == fmod(x, y) : F.F64}
522
523# IEEE remainder, x - n y with n = x / y rounded to nearest, ties to even:
524# from r = A mod B and q = A div B, n is q + 1 when r is past half of B (or
525# exactly half with q odd), which leaves B - r with the opposite sign; a zero
526# result has the sign of x (4.3)
527def rflip(+r: Nat, +b: Nat, +q: Nat) -> Bool:
528  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)))
529
530def rnear(+s: Bool, +r: Nat, +b: Nat, +q: Nat, +u: Nat) -> F.F64:
531  pick(F.F64, rflip(r, b, q), round(Bool.not(s), Nat.sub(b, r), u), round(s, r, u))
532
533def remainder(+x: F.F64, +y: F.F64) -> F.F64:
534  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))))
535
536def Remainder.value(+x: F.F64, +y: F.F64) -> Type:
537  {F.remainder(x, y) == remainder(x, y) : F.F64}
538
539# ---- closeness and exact ratios ----
540
541def lt_s(+x: F.F64, +y: F.F64) -> Bool:
542  Bool.and(ordered(x, y), Cmp.is_lt(ord(x, y)))
543
544def le_s(+x: F.F64, +y: F.F64) -> Bool:
545  Bool.and(ordered(x, y), Cmp.is_le(ord(x, y)))
546
547def eq_s(+x: F.F64, +y: F.F64) -> Bool:
548  Bool.and(ordered(x, y), Cmp.is_eq(ord(x, y)))
549
550def fabs(+x: F.F64) -> F.F64:
551  encode(False{}, efield(x), frac(x))
552
553# CPython's math.isclose, computed in binary64: a negative tolerance is a
554# domain error; a == b is close, an infinity is close only to itself, and
555# otherwise |b - a| <= max(|rel * b|, |rel * a|, abs_tol) (4.4)
556def ic_near(+a: F.F64, +b: F.F64, +rel: F.F64, +at: F.F64, +d: F.F64) -> Bool:
557  Bool.or(Bool.or(le_s(d, fabs(mul(rel, b))), le_s(d, fabs(mul(rel, a)))), le_s(d, at))
558
559def isclose(+a: F.F64, +b: F.F64, +rel: F.F64, +at: F.F64) -> Result<&2, &2, N.NumError, Bool>:
560  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))))))})
561
562def IsClose.value(+a: F.F64, +b: F.F64, +rel: F.F64, +at: F.F64) -> Type:
563  {F.isclose(a, b, rel, at) == isclose(a, b, rel, at) : Result<&2, &2, N.NumError, Bool>}
564
565# trailing zero bits of n (at most fuel)
566def tz(fuel: Nat, +n: Nat) -> Nat:
567  match fuel:
568    case 0n:
569      0n
570    case 1n+f:
571      pick(Nat, Nat.is_eq(Nat.mod(n, 2n), 1n), 0n, 1n+tz(f, Nat.div(n, 2n)))
572
573# (sign, numerator, exponent of the denominator)
574type Rat is Data:
575  Rat{neg: Bool, num: Nat, den: Nat}
576
577def rtriple(r: Result<&2, &2, N.NumError, F.Ratio>) -> Result<&2, &2, N.NumError, Rat>:
578  match r:
579    case Fail{e}:
580      Fail{e}
581    case Done{F.Ratio{neg, num, den}}:
582      Done{Rat{neg, SW.value(num), den}}
583
584# x = +-num / 2^den in lowest terms: the significand with as many of its
585# trailing zeros moved into the denominator as the scale allows; a zero is
586# 0/1 (Python's ints have no -0); NaN is a domain error, an infinity or an
587# integer of 2^64 or more an overflow (3.3)
588def as_ratio(+x: F.F64) -> Result<&2, &2, N.NumError, Rat>:
589  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))))}}))))
590
591def AsIntegerRatio.value(+x: F.F64) -> Type:
592  {rtriple(F.as_integer_ratio(x)) == as_ratio(x) : Result<&2, &2, N.NumError, Rat>}