1	import Base
2	import ../lib/common.bend as C
3	import ../../src/math/f64.bend as F
4	import ../../src/math/natural.bend as M
5	import ../../src/math/u64.bend as WU
6	import ../../src/math/num.bend as N
7	import ./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	
70	def hi_nat(x: F.F64) -> Nat:
71	  U32.to_nat(F.word_hi(x))
72	
73	def 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
77	def sign(+x: F.F64) -> Bool:
78	  Bool.not(C.fits(31n, hi_nat(x)))
79	
80	def efield(+x: F.F64) -> Nat:
81	  C.low(11n, C.high(20n, hi_nat(x)))
82	
83	def frac(+x: F.F64) -> Nat:
84	  Nat.add(lo_nat(x), C.shift(32n, C.low(20n, hi_nat(x))))
85	
86	def 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
95	def 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	
98	def 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)
106	def 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)
114	def qnan() -> F.F64:
115	  F.Bits{0, 2146959360}
116	
117	def inf(+s: Bool) -> F.F64:
118	  pick(F.F64, s, F.Bits{0, 4293918720}, F.Bits{0, 2146435072})
119	
120	def zero(+s: Bool) -> F.F64:
121	  pick(F.F64, s, F.Bits{0, 2147483648}, F.Bits{0, 0})
122	
123	def is_nan(+x: F.F64) -> Bool:
124	  Bool.and(Nat.is_eq(efield(x), 2047n), Bool.not(Nat.is_eq(frac(x), 0n)))
125	
126	def is_inf(+x: F.F64) -> Bool:
127	  Bool.and(Nat.is_eq(efield(x), 2047n), Nat.is_eq(frac(x), 0n))
128	
129	def 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	
136	def zb() -> Nat:
137	  3000n
138	
139	# the fraction with the hidden bit 2^52 of a normal number
140	def mant(+x: F.F64) -> Nat:
141	  Nat.add(frac(x), C.shift(52n, b2n(Bool.not(Nat.is_eq(efield(x), 0n)))))
142	
143	def 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
150	def 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	
153	def 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
161	def 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
165	def 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
171	def 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
176	def 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	
179	def 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)
186	def 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	
195	def 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	
198	def 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	
201	def 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	
204	def 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)
208	def neg(+x: F.F64) -> F.F64:
209	  encode(Bool.not(sign(x)), efield(x), frac(x))
210	
211	def 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	
214	def 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
218	def kq() -> Nat:
219	  200n
220	
221	def 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	
224	def 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	
227	def 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
232	def kr() -> Nat:
233	  100n
234	
235	def 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	
238	def 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	
241	def 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	
244	def 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
250	def 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	
253	def flip(c: Cmp) -> Cmp:
254	  match c:
255	    case LT{}:
256	      GT{}
257	    case EQ{}:
258	      EQ{}
259	    case GT{}:
260	      LT{}
261	
262	def 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	
265	def ordered(+x: F.F64, +y: F.F64) -> Bool:
266	  Bool.not(Bool.or(is_nan(x), is_nan(y)))
267	
268	# ---- the contract ----
269	
270	def 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)
274	def Sub.value(+x: F.F64, +y: F.F64) -> Type:
275	  {F.sub(x, y) == add(x, neg(y)) : F.F64}
276	
277	def Mul.value(+x: F.F64, +y: F.F64) -> Type:
278	  {F.mul(x, y) == mul(x, y) : F.F64}
279	
280	def Div.value(+x: F.F64, +y: F.F64) -> Type:
281	  {F.div(x, y) == div(x, y) : F.F64}
282	
283	def Sqrt.value(+x: F.F64) -> Type:
284	  {F.sqrt(x) == sqrt(x) : F.F64}
285	
286	def 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	
289	def 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	
292	def 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	
295	def Neg.value(+x: F.F64) -> Type:
296	  {F.neg(x) == neg(x) : F.F64}
297	
298	def Abs.value(+x: F.F64) -> Type:
299	  {F.abs(x) == encode(False{}, efield(x), frac(x)) : F.F64}
300	
301	def Copysign.value(+x: F.F64, +y: F.F64) -> Type:
302	  {F.copysign(x, y) == encode(sign(y), efield(x), frac(x)) : F.F64}
303	
304	def IsNan.value(+x: F.F64) -> Type:
305	  {F.is_nan(x) == is_nan(x) : Bool}
306	
307	def IsInf.value(+x: F.F64) -> Type:
308	  {F.is_inf(x) == is_inf(x) : Bool}
309	
310	def IsFinite.value(+x: F.F64) -> Type:
311	  {F.is_finite(x) == Bool.not(Nat.is_eq(efield(x), 2047n)) : Bool}
312	
313	def IsZero.value(+x: F.F64) -> Type:
314	  {F.is_zero(x) == is_zero(x) : Bool}
315	
316	def 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)
320	def 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)
329	def 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)
343	def 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	
346	def Trunc.value(+x: F.F64) -> Type:
347	  {F.trunc(x) == to_integral(F.Trunc{}, x) : F.F64}
348	
349	def Floor.value(+x: F.F64) -> Type:
350	  {F.floor(x) == to_integral(F.Floor{}, x) : F.F64}
351	
352	def Ceil.value(+x: F.F64) -> Type:
353	  {F.ceil(x) == to_integral(F.Ceil{}, x) : F.F64}
354	
355	def 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)
361	def 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)
367	def 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	
370	def 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	
377	def 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	
384	def ToU64.value(+x: F.F64) -> Type:
385	  {rv64(F.to_u64(x)) == to_nat(x, 64n) : Result<&2, &2, N.NumError, Nat>}
386	
387	def ToU32.value(+x: F.F64) -> Type:
388	  {rv32(F.to_u32(x)) == to_nat(x, 32n) : Result<&2, &2, N.NumError, Nat>}
389	
390	def 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	
393	def 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	
396	def 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)
400	def OfU64.value(+w: WU.U64) -> Type:
401	  {F.of_u64(w) == round(False{}, SW.value(w), zb()) : F.F64}
402	
403	def 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
409	def 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)
415	def 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	
418	def 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
424	def 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	
427	def 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)
432	def 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
438	def pat(+x: F.F64) -> Nat:
439	  Nat.add(frac(x), C.shift(52n, efield(x)))
440	
441	def 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)
449	def 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	
452	def 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
457	def 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	
460	def 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	
463	def Fmin.value(+x: F.F64, +y: F.F64) -> Type:
464	  {F.fmin(x, y) == fmin(x, y) : F.F64}
465	
466	def 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	
471	def 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	
474	def 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
478	def 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	
481	def Bits.roundtrip(+x: F.F64) -> Type:
482	  {F.of_bits64(F.to_bits(x)) == x : F.F64}
483	
484	def Bits.inverse(+w: WU.U64) -> Type:
485	  {F.to_bits(F.of_bits64(w)) == w : WU.U64}
486	
487	# finite with no fractional bits
488	def 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)
495	def 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	
498	def 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
502	def rsc(+x: F.F64, +y: F.F64) -> Nat:
503	  Nat.min(xexp(x), xexp(y))
504	
505	def ra(+x: F.F64, +y: F.F64) -> Nat:
506	  C.shift(Nat.sub(xexp(x), rsc(x, y)), mant(x))
507	
508	def rb(+x: F.F64, +y: F.F64) -> Nat:
509	  C.shift(Nat.sub(xexp(y), rsc(x, y)), mant(y))
510	
511	def 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)
517	def 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	
520	def 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)
527	def 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	
530	def 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	
533	def 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	
536	def 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	
541	def lt_s(+x: F.F64, +y: F.F64) -> Bool:
542	  Bool.and(ordered(x, y), Cmp.is_lt(ord(x, y)))
543	
544	def le_s(+x: F.F64, +y: F.F64) -> Bool:
545	  Bool.and(ordered(x, y), Cmp.is_le(ord(x, y)))
546	
547	def eq_s(+x: F.F64, +y: F.F64) -> Bool:
548	  Bool.and(ordered(x, y), Cmp.is_eq(ord(x, y)))
549	
550	def 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)
556	def 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	
559	def 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	
562	def 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)
566	def 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)
574	type Rat is Data:
575	  Rat{neg: Bool, num: Nat, den: Nat}
576	
577	def 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)
588	def 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	
591	def AsIntegerRatio.value(+x: F.F64) -> Type:
592	  {rtriple(F.as_integer_ratio(x)) == as_ratio(x) : Result<&2, &2, N.NumError, Rat>}