1	import Base
2	import ./num.bend as N
3	import ./u64.bend as W
4	import ./w64.bend as X
5	import ./natural.bend as M
6	
7	# IEEE-754 binary64 in software, on two U32 words (Bend 2 has only F32).
8	# The algorithms are Berkeley SoftFloat 3e's (f64_add/f64_sub through
9	# addMags/subMags, f64_mul, roundPackToF64 and normRoundPackToF64, round to
10	# nearest with ties to even), on 64-bit significands held as two U32 limbs
11	# (w64.bend). Division is exact long division by 32-bit quotient digits and
12	# the square root one step of Zimmermann's Karatsuba square root with its
13	# exact remainder, each followed by a sticky bit, so every result is correctly rounded: subnormals, overflow
14	# to infinity, signed zeros and NaN follow IEEE 754 (x + (-x) is +0, 0/0 and
15	# inf - inf are NaN; a NaN result is the canonical quiet NaN).
16	#
17	#   of_bits(hi, lo), word_hi(x), word_lo(x)          the bit pattern
18	#   is_nan, is_inf, is_finite, is_zero, signbit
19	#   neg, abs, copysign
20	#   add, sub, mul, div, sqrt               correctly rounded
21	#   lt, le, eq                             IEEE comparisons (NaN unordered)
22	#   of_nat(n)                              correctly rounded (exact below 2^53)
23	#   f64_op, f64_is                         the num.bend instance
24	#
25	# Exponents that may go below zero before rounding are kept as naturals
26	# offset by OFF = 4096. Tested against the machine's doubles
27	# (tools/check_f64.py); not proved.
28	
29	type F64 is Data:
30	  Bits{lo: U32, hi: U32}
31	
32	def word_lo(x: F64) -> U32:
33	  match x:
34	    case Bits{l, h}:
35	      l
36	
37	def word_hi(x: F64) -> U32:
38	  match x:
39	    case Bits{l, h}:
40	      h
41	
42	def of_bits(+hi: U32, +lo: U32) -> F64:
43	  Bits{lo, hi}
44	
45	def off() -> Nat:
46	  4096n
47	
48	def sgn(s: Bool) -> U32:
49	  match s:
50	    case True{}:
51	      2147483648
52	    case False{}:
53	      0
54	
55	# ---- fields ----
56	
57	def signbit(+x: F64) -> Bool:
58	  U32.is_le(2147483648, word_hi(x))
59	
60	def exp_field(+x: F64) -> Nat:
61	  U32.to_nat(U32.and(U32.div(word_hi(x), 1048576), 2047))
62	
63	def frac(+x: F64) -> W.U64:
64	  W.U64{word_lo(x), U32.and(word_hi(x), 1048575)}
65	
66	def mag(+x: F64) -> W.U64:
67	  W.U64{word_lo(x), U32.and(word_hi(x), 2147483647)}
68	
69	def nz(+a: W.U64) -> Bool:
70	  Bool.not(X.is_zero(a))
71	
72	def is_nan(+x: F64) -> Bool:
73	  Bool.and(Nat.is_eq(exp_field(x), 2047n), nz(frac(x)))
74	
75	def is_inf(+x: F64) -> Bool:
76	  Bool.and(Nat.is_eq(exp_field(x), 2047n), X.is_zero(frac(x)))
77	
78	def is_finite(+x: F64) -> Bool:
79	  Nat.is_lt(exp_field(x), 2047n)
80	
81	def is_zero(+x: F64) -> Bool:
82	  X.is_zero(mag(x))
83	
84	# ---- constants and sign ----
85	
86	def nan() -> F64:
87	  Bits{0, 2146959360}
88	
89	def inf(s: Bool) -> F64:
90	  Bits{0, U32.add(2146435072, sgn(s))}
91	
92	def zero(s: Bool) -> F64:
93	  Bits{0, sgn(s)}
94	
95	def one() -> F64:
96	  Bits{0, 1072693248}
97	
98	def neg(+x: F64) -> F64:
99	  Bits{word_lo(x), U32.xor(word_hi(x), 2147483648)}
100	
101	def abs(+x: F64) -> F64:
102	  Bits{word_lo(x), U32.and(word_hi(x), 2147483647)}
103	
104	def copysign(+x: F64, +y: F64) -> F64:
105	  Bits{word_lo(x), U32.add(U32.and(word_hi(x), 2147483647), U32.and(word_hi(y), 2147483648))}
106	
107	def nan_or(+x: F64, bad: Bool) -> F64:
108	  match bad:
109	    case True{}:
110	      nan()
111	    case False{}:
112	      x
113	
114	# ---- packing and rounding (SoftFloat's packToF64, roundPackToF64) ----
115	
116	def pack64(+w: W.U64) -> F64:
117	  Bits{X.lo(w), X.hi(w)}
118	
119	# sign, exponent field e and significand m; the hidden bit of m (bit 52)
120	# carries into the exponent field, as packToF64's addition does
121	def pack(s: Bool, +e: Nat, +m: W.U64) -> F64:
122	  pack64(X.add(W.U64{0, U32.add(sgn(s), U32.mul(U32.from_nat(e), 1048576))}, m))
123	
124	def zero_e(+e: Nat, z: Bool) -> Nat:
125	  match z:
126	    case True{}:
127	      0n
128	    case False{}:
129	      e
130	
131	def rp_fin3(+s: Bool, +e: Nat, +r: W.U64) -> F64:
132	  pack(s, zero_e(e, X.is_zero(r)), r)
133	
134	def rp_fin2(+s: Bool, +e: Nat, +r: W.U64, tie: Bool) -> F64:
135	  rp_fin3(s, e, X.clear0(r, tie))
136	
137	# sig has its top bit at 62 and e is the plain exponent (0 <= e <= 0x7FD):
138	# add half an ulp of the kept 53 bits, shift, and clear bit 0 on a tie
139	def rp_fin(+s: Bool, +e: Nat, +sig: W.U64) -> F64:
140	  rp_fin2(s, e, X.shr(X.add(sig, W.U64{512, 0}), 10n), U32.is_eq(U32.and(X.lo(sig), 1023), 512))
141	
142	def rp_over(+s: Bool, +e: Nat, +sig: W.U64, over: Bool) -> F64:
143	  match over:
144	    case True{}:
145	      inf(s)
146	    case False{}:
147	      rp_fin(s, Nat.sub(e, off()), sig)
148	
149	def rp_neg(+s: Bool, +e: Nat, +sig: W.U64, below: Bool) -> F64:
150	  match below:
151	    case True{}:
152	      rp_fin(s, 0n, X.shr_jam(sig, Nat.sub(off(), e)))
153	    case False{}:
154	      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})))))
155	
156	# e is the biased exponent minus one, offset by OFF; sig has its top bit at 62
157	def round_pack(+s: Bool, +e: Nat, +sig: W.U64) -> F64:
158	  rp_neg(s, e, sig, Nat.is_lt(e, off()))
159	
160	# sig below 2^62 is shifted up once
161	def rp62_pick(+s: Bool, +e: Nat, +sig: W.U64, low: Bool) -> F64:
162	  match low:
163	    case True{}:
164	      round_pack(s, Nat.sub(e, 1n), X.add(sig, sig))
165	    case False{}:
166	      round_pack(s, e, sig)
167	
168	def rp62(+s: Bool, +e: Nat, +sig: W.U64) -> F64:
169	  rp62_pick(s, e, sig, X.lt(sig, W.U64{0, 1073741824}))
170	
171	def nrp_exp(+e: Nat, +sd: Nat, z: Bool) -> Nat:
172	  match z:
173	    case True{}:
174	      0n
175	    case False{}:
176	      Nat.sub(Nat.sub(e, sd), off())
177	
178	def nrp_pick(+s: Bool, +e: Nat, +sig: W.U64, +sd: Nat, direct: Bool) -> F64:
179	  match direct:
180	    case True{}:
181	      pack(s, nrp_exp(e, sd, X.is_zero(sig)), X.shl(sig, Nat.sub(sd, 10n)))
182	    case False{}:
183	      round_pack(s, Nat.sub(e, sd), X.shl(sig, sd))
184	
185	def nrp_sd(+s: Bool, +e: Nat, +sig: W.U64, +sd: Nat) -> F64:
186	  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)))))
187	
188	# SoftFloat's normRoundPackToF64 (e offset)
189	def norm_round_pack(+s: Bool, +e: Nat, +sig: W.U64) -> F64:
190	  nrp_sd(s, e, sig, Nat.sub(X.clz(sig), 1n))
191	
192	# the offset exponent and 53-bit significand (hidden bit at 52) of a nonzero
193	# finite value, subnormals normalized (SoftFloat's normSubnormalF64Sig)
194	def norm_e(+e: Nat, +f: W.U64) -> Nat:
195	  match e:
196	    case 0n:
197	      Nat.sub(Nat.add(off(), 12n), X.clz(f))
198	    case _:
199	      Nat.add(off(), e)
200	
201	def norm_f(+e: Nat, +f: W.U64) -> W.U64:
202	  match e:
203	    case 0n:
204	      X.shl(f, Nat.sub(X.clz(f), 11n))
205	    case _:
206	      X.add(f, W.U64{0, 1048576})
207	
208	# ---- addition (SoftFloat's addMagsF64 and subMagsF64) ----
209	
210	def am_small(+e: Nat, +f9: W.U64) -> W.U64:
211	  match e:
212	    case 0n:
213	      X.add(f9, f9)
214	    case _:
215	      X.add(f9, W.U64{0, 536870912})
216	
217	# |x| + |y| with sign s, the larger exponent eL
218	def am_big(+big: F64, +s: Bool, +el: Nat, +fl: W.U64, +es: Nat, +fs: W.U64, top: Bool) -> F64:
219	  match top:
220	    case True{}:
221	      nan_or(big, nz(fl))
222	    case False{}:
223	      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))))
224	
225	def am_eq_n(+x: F64, +s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64, top: Bool) -> F64:
226	  match top:
227	    case True{}:
228	      nan_or(x, Bool.or(nz(fa), nz(fb)))
229	    case False{}:
230	      round_pack(s, Nat.add(off(), e), X.shl(X.add(W.U64{0, 2097152}, X.add(fa, fb)), 9n))
231	
232	def am_eq(+x: F64, +s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64) -> F64:
233	  match e:
234	    case 0n:
235	      pack(s, 0n, X.add(fa, fb))
236	    case _:
237	      am_eq_n(x, s, e, fa, fb, Nat.is_eq(e, 2047n))
238	
239	def am_case(+x: F64, +y: F64, +s: Bool, +ea: Nat, +eb: Nat, c: Cmp) -> F64:
240	  match c:
241	    case EQ{}:
242	      am_eq(x, s, ea, frac(x), frac(y))
243	    case LT{}:
244	      am_big(y, s, eb, frac(y), ea, frac(x), Nat.is_eq(eb, 2047n))
245	    case GT{}:
246	      am_big(x, s, ea, frac(x), eb, frac(y), Nat.is_eq(ea, 2047n))
247	
248	def add_mags(+x: F64, +y: F64, +s: Bool) -> F64:
249	  am_case(x, y, s, exp_field(x), exp_field(y), Nat.cmp(exp_field(x), exp_field(y)))
250	
251	def sm_small(+e: Nat, +f10: W.U64) -> W.U64:
252	  match e:
253	    case 0n:
254	      X.add(f10, f10)
255	    case _:
256	      X.add(f10, W.U64{0, 1073741824})
257	
258	# |big| - |small| with sign s, the larger exponent eL
259	def sm_big(+big: F64, +s: Bool, +el: Nat, +fl: W.U64, +es: Nat, +fs: W.U64, top: Bool) -> F64:
260	  match top:
261	    case True{}:
262	      nan_or(big, nz(fl))
263	    case False{}:
264	      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))))
265	
266	def sm_exact3(+s: Bool, +e1: Nat, +d: W.U64, +sd: Nat, under: Bool) -> F64:
267	  match under:
268	    case True{}:
269	      pack(s, 0n, X.shl(d, e1))
270	    case False{}:
271	      pack(s, Nat.sub(e1, sd), X.shl(d, sd))
272	
273	# an exact difference of equal exponents, normalized without rounding
274	def sm_exact(+s: Bool, +e1: Nat, +d: W.U64) -> F64:
275	  sm_exact3(s, e1, d, Nat.sub(X.clz(d), 11n), Nat.is_lt(e1, Nat.sub(X.clz(d), 11n)))
276	
277	def sm_eq2(+s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64, c: Cmp) -> F64:
278	  match c:
279	    case EQ{}:
280	      zero(False{})
281	    case LT{}:
282	      sm_exact(Bool.not(s), Nat.sub(e, 1n), X.sub(fb, fa))
283	    case GT{}:
284	      sm_exact(s, Nat.sub(e, 1n), X.sub(fa, fb))
285	
286	def cmp64(+a: W.U64, +b: W.U64) -> Cmp:
287	  X.cmp(a, b)
288	
289	def sm_eq(+s: Bool, +e: Nat, +fa: W.U64, +fb: W.U64, top: Bool) -> F64:
290	  match top:
291	    case True{}:
292	      nan()
293	    case False{}:
294	      sm_eq2(s, e, fa, fb, cmp64(fa, fb))
295	
296	def sm_case(+x: F64, +y: F64, +s: Bool, +ea: Nat, +eb: Nat, c: Cmp) -> F64:
297	  match c:
298	    case EQ{}:
299	      sm_eq(s, ea, frac(x), frac(y), Nat.is_eq(ea, 2047n))
300	    case LT{}:
301	      sm_big(y, Bool.not(s), eb, frac(y), ea, frac(x), Nat.is_eq(eb, 2047n))
302	    case GT{}:
303	      sm_big(x, s, ea, frac(x), eb, frac(y), Nat.is_eq(ea, 2047n))
304	
305	def sub_mags(+x: F64, +y: F64, +s: Bool) -> F64:
306	  sm_case(x, y, s, exp_field(x), exp_field(y), Nat.cmp(exp_field(x), exp_field(y)))
307	
308	def add_pick(+x: F64, +y: F64, same: Bool) -> F64:
309	  match same:
310	    case True{}:
311	      add_mags(x, y, signbit(x))
312	    case False{}:
313	      sub_mags(x, y, signbit(x))
314	
315	def add(+x: F64, +y: F64) -> F64:
316	  add_pick(x, y, Bool.not(Bool.xor(signbit(x), signbit(y))))
317	
318	def sub(+x: F64, +y: F64) -> F64:
319	  add(x, neg(y))
320	
321	# ---- multiplication (SoftFloat's f64_mul) ----
322	
323	def mul_n2(+s: Bool, +e: Nat, p: W.U64 & W.U64) -> F64:
324	  (+pl, +ph) = p
325	  rp62(s, e, X.or_bit(ph, nz(pl)))
326	
327	def mul_n(+s: Bool, +ea: Nat, +ma: W.U64, +eb: Nat, +mb: W.U64) -> F64:
328	  mul_n2(s, Nat.sub(Nat.add(ea, eb), Nat.add(off(), 1023n)), X.mul128(X.shl(ma, 10n), X.shl(mb, 11n)))
329	
330	def fin_zero(+e: Nat, +f: W.U64) -> Bool:
331	  Bool.and(Nat.is_eq(e, 0n), X.is_zero(f))
332	
333	def mul_z(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, z: Bool) -> F64:
334	  match z:
335	    case True{}:
336	      zero(s)
337	    case False{}:
338	      mul_n(s, norm_e(ea, fa), norm_f(ea, fa), norm_e(eb, fb), norm_f(eb, fb))
339	
340	def inf_or_nan(+s: Bool, bad: Bool) -> F64:
341	  match bad:
342	    case True{}:
343	      nan()
344	    case False{}:
345	      inf(s)
346	
347	def mul_b(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:
348	  match top:
349	    case True{}:
350	      inf_or_nan(s, Bool.or(nz(fb), fin_zero(ea, fa)))
351	    case False{}:
352	      mul_z(s, ea, fa, eb, fb, Bool.or(fin_zero(ea, fa), fin_zero(eb, fb)))
353	
354	def mul_cls(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:
355	  match top:
356	    case True{}:
357	      inf_or_nan(s, Bool.or(Bool.or(nz(fa), Bool.and(Nat.is_eq(eb, 2047n), nz(fb))), fin_zero(eb, fb)))
358	    case False{}:
359	      mul_b(s, ea, fa, eb, fb, Nat.is_eq(eb, 2047n))
360	
361	def mul(+x: F64, +y: F64) -> F64:
362	  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))
363	
364	# ---- division ----
365	
366	# the next quotient digit of r * 2^32 / b (r < b, b >= 2^52) and the remainder
367	def digit_fin(+r: W.U64, +b: W.U64, +d: U32) -> U32 & W.U64:
368	  (d, X.sub(W.U64{0, X.lo(r)}, X.fst_q(X.mul_32_64(d, b))))
369	
370	def digit(+r: W.U64, +b: W.U64, +t: Nat) -> U32 & W.U64:
371	  digit_fin(r, b, X.q96(W.U64{0, X.lo(r)}, X.hi(r), b, t))
372	
373	def dq_fin(+s: Bool, +e: Nat, +d1: U32, p: U32 & W.U64) -> F64:
374	  (+d2, +r2) = p
375	  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))))))
376	
377	def dq_mid(+s: Bool, +e: Nat, +b: W.U64, +t: Nat, p: U32 & W.U64) -> F64:
378	  (+d1, +r1) = p
379	  dq_fin(s, e, d1, digit(r1, b, t))
380	
381	# a in [b, 2b): 2^62 + floor((a - b) * 2^62 / b), sticky, by two 32-bit digits
382	def div_qt(+s: Bool, +e: Nat, +a: W.U64, +b: W.U64, +t: Nat) -> F64:
383	  dq_mid(s, e, b, t, digit(X.sub(a, b), b, t))
384	
385	def div_q(+s: Bool, +e: Nat, +a: W.U64, +b: W.U64) -> F64:
386	  div_qt(s, e, a, b, X.bitlen(X.hi(b)))
387	
388	def div_ab(+s: Bool, +ea: Nat, +a: W.U64, +eb: Nat, +b: W.U64, less: Bool) -> F64:
389	  match less:
390	    case True{}:
391	      div_q(s, Nat.sub(Nat.add(ea, Nat.add(off(), 1021n)), eb), X.add(a, a), b)
392	    case False{}:
393	      div_q(s, Nat.sub(Nat.add(ea, Nat.add(off(), 1022n)), eb), a, b)
394	
395	def div_n(+s: Bool, +ea: Nat, +a: W.U64, +eb: Nat, +b: W.U64) -> F64:
396	  div_ab(s, ea, a, eb, b, X.lt(a, b))
397	
398	def div_za(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, za: Bool) -> F64:
399	  match za:
400	    case True{}:
401	      zero(s)
402	    case False{}:
403	      div_n(s, norm_e(ea, fa), norm_f(ea, fa), norm_e(eb, fb), norm_f(eb, fb))
404	
405	def div_z(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, zb: Bool) -> F64:
406	  match zb:
407	    case True{}:
408	      inf_or_nan(s, fin_zero(ea, fa))
409	    case False{}:
410	      div_za(s, ea, fa, eb, fb, fin_zero(ea, fa))
411	
412	def div_b(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:
413	  match top:
414	    case True{}:
415	      nan_or(zero(s), nz(fb))
416	    case False{}:
417	      div_z(s, ea, fa, eb, fb, fin_zero(eb, fb))
418	
419	def div_cls(+s: Bool, +ea: Nat, +fa: W.U64, +eb: Nat, +fb: W.U64, top: Bool) -> F64:
420	  match top:
421	    case True{}:
422	      inf_or_nan(s, Bool.or(nz(fa), Nat.is_eq(eb, 2047n)))
423	    case False{}:
424	      div_b(s, ea, fa, eb, fb, Nat.is_eq(eb, 2047n))
425	
426	def div(+x: F64, +y: F64) -> F64:
427	  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))
428	
429	# ---- square root ----
430	
431	def sq_fin(+e: Nat, +t: W.U64, sticky: Bool) -> F64:
432	  round_pack(False{}, e, X.or_bit(t, sticky))
433	
434	# the root is S0 - 2: remainder (2 S0 - 3) - (D - (2 S0 - 1)) with c = 2 S0 - 1
435	def sq_d2(+e: Nat, +s0: W.U64, +d: W.U64, +c: W.U64) -> F64:
436	  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))))
437	
438	# S0^2 - N = D > 0: the root is S0 - 1 when D <= 2 S0 - 1 (remainder c - D)
439	def sq_d1(+e: Nat, +s0: W.U64, +d: W.U64, +c: W.U64, small: Bool) -> F64:
440	  match small:
441	    case True{}:
442	      sq_fin(e, X.sub(s0, W.U64{1, 0}), Bool.not(X.eq(c, d)))
443	    case False{}:
444	      sq_d2(e, s0, d, c)
445	
446	def sq_neg(+e: Nat, +s0: W.U64, +d: W.U64) -> F64:
447	  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})))
448	
449	# N - S0^2 = a - b for a = u 2^32, b = q^2: the root is S0 when b <= a
450	def sq_rem(+e: Nat, +s0: W.U64, +a: W.U64, +b: W.U64, ge: Bool) -> F64:
451	  match ge:
452	    case True{}:
453	      sq_fin(e, s0, Bool.not(X.eq(a, b)))
454	    case False{}:
455	      sq_neg(e, s0, X.sub(b, a))
456	
457	def sq_qu(+e: Nat, +s: U32, +q: U32, +u: U32) -> F64:
458	  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}))
459	
460	# q = 2^32 (r = 2 s) is taken as q = 2^32 - 1 with u = 2 s
461	def sq_clamp(+e: Nat, +s: U32, +q: W.U64, +u: U32, small: Bool) -> F64:
462	  match small:
463	    case True{}:
464	      sq_qu(e, s, X.lo(q), u)
465	    case False{}:
466	      sq_qu(e, s, 4294967295, U32.add(s, s))
467	
468	def sq_div(+e: Nat, +s: U32, p: W.U64 & U32) -> F64:
469	  (+q, +u) = p
470	  sq_clamp(e, s, q, u, U32.is_zero(X.hi(q)))
471	
472	# one step of Zimmermann's SqrtRem (Karatsuba square root) from s = isqrt(nh):
473	# r = nh - s^2, (q, u) = divmod(r 2^32, 2 s) give S0 = s 2^32 + q with
474	# N - S0^2 = u 2^32 - q^2 for N = nh 2^64; S0 - 2 <= isqrt(N) <= S0 and the
475	# exact remainder decides the root and its sticky bit without squaring S0
476	def sq_root_s(+e: Nat, +nh: W.U64, +s: U32) -> F64:
477	  sq_div(e, s, X.div32(W.U64{0, X.lo(X.sub(nh, X.mul32(s, s)))}, U32.add(s, s)))
478	
479	def sq_root(+e: Nat, +nh: W.U64) -> F64:
480	  sq_root_s(e, nh, X.lo(X.isqrt(nh)))
481	
482	# m * 2^(E - OFF - 1075) with the unbiased exponent made even (E odd); the
483	# root of m * 2^72 has its top bit at 62
484	def sq_even(+e: Nat, +m: W.U64) -> F64:
485	  sq_root(Nat.add(Nat.div(Nat.sub(Nat.add(e, off()), 1075n), 2n), 1048n), X.shl(m, 8n))
486	
487	def sq_n(+e: Nat, +m: W.U64, odd: Bool) -> F64:
488	  match odd:
489	    case True{}:
490	      sq_even(e, m)
491	    case False{}:
492	      sq_even(Nat.sub(e, 1n), X.add(m, m))
493	
494	def sq_pos(+e: Nat, +f: W.U64, negative: Bool) -> F64:
495	  match negative:
496	    case True{}:
497	      nan()
498	    case False{}:
499	      sq_n(norm_e(e, f), norm_f(e, f), Nat.is_eq(Nat.mod(norm_e(e, f), 2n), 1n))
500	
501	def sq_sign(+x: F64, negative: Bool, +e: Nat, +f: W.U64, zr: Bool) -> F64:
502	  match zr:
503	    case True{}:
504	      x
505	    case False{}:
506	      sq_pos(e, f, negative)
507	
508	def sq_cls(+x: F64, s: Bool, +e: Nat, +f: W.U64, top: Bool) -> F64:
509	  match top:
510	    case True{}:
511	      nan_or(x, Bool.or(nz(f), s))
512	    case False{}:
513	      sq_sign(x, s, e, f, fin_zero(e, f))
514	
515	def sqrt(+x: F64) -> F64:
516	  sq_cls(x, signbit(x), exp_field(x), frac(x), Nat.is_eq(exp_field(x), 2047n))
517	
518	# ---- comparison and conversion ----
519	
520	def lt_s(+x: F64, +y: F64, sa: Bool, sb: Bool) -> Bool:
521	  match sa sb:
522	    case True{} False{}:
523	      True{}
524	    case False{} True{}:
525	      False{}
526	    case False{} False{}:
527	      X.lt(mag(x), mag(y))
528	    case True{} True{}:
529	      X.lt(mag(y), mag(x))
530	
531	def le_s(+x: F64, +y: F64, sa: Bool, sb: Bool) -> Bool:
532	  match sa sb:
533	    case True{} False{}:
534	      True{}
535	    case False{} True{}:
536	      False{}
537	    case False{} False{}:
538	      X.le(mag(x), mag(y))
539	    case True{} True{}:
540	      X.le(mag(y), mag(x))
541	
542	def unordered(+x: F64, +y: F64) -> Bool:
543	  Bool.or(is_nan(x), is_nan(y))
544	
545	def zeros(+x: F64, +y: F64) -> Bool:
546	  Bool.and(is_zero(x), is_zero(y))
547	
548	def lt_z(+x: F64, +y: F64, bad: Bool, z: Bool) -> Bool:
549	  match bad z:
550	    case True{} _:
551	      False{}
552	    case False{} True{}:
553	      False{}
554	    case False{} False{}:
555	      lt_s(x, y, signbit(x), signbit(y))
556	
557	def lt(+x: F64, +y: F64) -> Bool:
558	  lt_z(x, y, unordered(x, y), zeros(x, y))
559	
560	def le_z(+x: F64, +y: F64, bad: Bool, z: Bool) -> Bool:
561	  match bad z:
562	    case True{} _:
563	      False{}
564	    case False{} True{}:
565	      True{}
566	    case False{} False{}:
567	      le_s(x, y, signbit(x), signbit(y))
568	
569	def le(+x: F64, +y: F64) -> Bool:
570	  le_z(x, y, unordered(x, y), zeros(x, y))
571	
572	def eq_z(+x: F64, +y: F64, bad: Bool, z: Bool) -> Bool:
573	  match bad z:
574	    case True{} _:
575	      False{}
576	    case False{} True{}:
577	      True{}
578	    case False{} False{}:
579	      Bool.and(U32.is_eq(word_lo(x), word_lo(y)), U32.is_eq(word_hi(x), word_hi(y)))
580	
581	def eq(+x: F64, +y: F64) -> Bool:
582	  eq_z(x, y, unordered(x, y), zeros(x, y))
583	
584	# n mod 2^k and n div 2^k by k halvings (no 2^32 constant: the proof
585	# checker would expand it in unary)
586	def low_bits(k: Nat, +n: Nat) -> Nat:
587	  match k:
588	    case 0n:
589	      0n
590	    case 1n+p:
591	      Nat.add(Nat.mod(n, 2n), Nat.double(low_bits(p, Nat.div(n, 2n))))
592	
593	def high_bits(k: Nat, +n: Nat) -> Nat:
594	  match k:
595	    case 0n:
596	      n
597	    case 1n+p:
598	      high_bits(p, Nat.div(n, 2n))
599	
600	# the jam of n >> d (the bits shifted out OR-ed into bit 0)
601	def jam_nat(+d: Nat, +n: Nat) -> Nat:
602	  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)))
603	
604	def word64(+s: Nat) -> W.U64:
605	  W.U64{U32.from_nat(low_bits(32n, s)), U32.from_nat(high_bits(32n, s))}
606	
607	# n with its top bit moved to bit 62 (b = bit_length(n)): shifted up when
608	# b <= 63 (every runtime Nat), jammed down otherwise
609	def sig63(+n: Nat, +b: Nat, fits: Bool) -> W.U64:
610	  match fits:
611	    case True{}:
612	      X.shl(word64(n), Nat.sub(63n, b))
613	    case False{}:
614	      word64(jam_nat(Nat.sub(b, 63n), n))
615	
616	def of_nat_z(+n: Nat, z: Bool) -> F64:
617	  match z:
618	    case True{}:
619	      zero(False{})
620	    case False{}:
621	      round_pack(False{}, Nat.add(5117n, M.bit_length(n)), sig63(n, M.bit_length(n), Nat.is_le(M.bit_length(n), 63n)))
622	
623	# the double nearest to n (exact below 2^53)
624	def of_nat(+n: Nat) -> F64:
625	  of_nat_z(n, Nat.is_eq(n, 0n))
626	
627	# 2^k for k <= 1023
628	def pow2(+k: Nat) -> F64:
629	  pack(False{}, Nat.add(1022n, k), W.U64{0, 1048576})
630	
631	# ---- decoding and exact rounding (the tools of everything below) ----
632	
633	# the significand with its hidden bit (bit 52 of a normal number) and the
634	# scale of its unit, the exponent plus Z = 3000 (1926 for subnormals), as in
635	# spec/math/f64.bend's mant and xexp
636	def dmant_z(+f: W.U64, z: Bool) -> W.U64:
637	  match z:
638	    case True{}:
639	      f
640	    case False{}:
641	      X.add(f, W.U64{0, 1048576})
642	
643	def dmant(+x: F64) -> W.U64:
644	  dmant_z(frac(x), Nat.is_eq(exp_field(x), 0n))
645	
646	def dexp_z(+e: Nat, z: Bool) -> Nat:
647	  match z:
648	    case True{}:
649	      1926n
650	    case False{}:
651	      Nat.add(e, 1925n)
652	
653	def dexp(+x: F64) -> Nat:
654	  dexp_z(exp_field(x), Nat.is_eq(exp_field(x), 0n))
655	
656	def rw_top(+s: Bool, +x: Nat, +w: W.U64, top: Bool) -> F64:
657	  match top:
658	    case True{}:
659	      norm_round_pack(s, Nat.add(x, 2181n), X.shr_jam(w, 1n))
660	    case False{}:
661	      norm_round_pack(s, Nat.add(x, 2180n), w)
662	
663	def rw_z(+s: Bool, +x: Nat, +w: W.U64, z: Bool) -> F64:
664	  match z:
665	    case True{}:
666	      zero(s)
667	    case False{}:
668	      rw_top(s, x, w, X.le(W.U64{0, 2147483648}, w))
669	
670	# the double nearest to (-1)^s * w * 2^(x - 3000), for x >= 63: every exact
671	# result below is an integer of at most 64 bits at a scale, rounded once
672	def round_w(+s: Bool, +x: Nat, +w: W.U64) -> F64:
673	  rw_z(s, x, w, X.is_zero(w))
674	
675	# ---- rounding to an integral value (IEEE roundToIntegral) ----
676	
677	type RMode is Data:
678	  Trunc{}
679	  Floor{}
680	  Ceil{}
681	  Even{}
682	
683	def b64(b: Bool) -> W.U64:
684	  match b:
685	    case True{}:
686	      W.U64{1, 0}
687	    case False{}:
688	      W.U64{0, 0}
689	
690	# whether the integer part q (remainder r, half h of the unit) moves up one
691	def ri_up(m: RMode, +s: Bool, +q: W.U64, +r: W.U64, +h: W.U64) -> Bool:
692	  match m:
693	    case Trunc{}:
694	      False{}
695	    case Floor{}:
696	      Bool.and(s, Bool.not(X.is_zero(r)))
697	    case Ceil{}:
698	      Bool.and(Bool.not(s), Bool.not(X.is_zero(r)))
699	    case Even{}:
700	      Bool.or(X.lt(h, r), Bool.and(X.eq(r, h), X.odd(q)))
701	
702	def ri_q(m: RMode, +s: Bool, +q: W.U64, +r: W.U64, +h: W.U64) -> F64:
703	  round_w(s, 3000n, X.add(q, b64(ri_up(m, s, q, r, h))))
704	
705	# k < 64 fractional bits: q = w >> k, r = w - (q << k), h = 2^(k-1)
706	def ri_k(m: RMode, +s: Bool, +w: W.U64, +k: Nat, small: Bool) -> F64:
707	  match small:
708	    case True{}:
709	      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)))
710	    case False{}:
711	      ri_q(m, s, W.U64{0, 0}, w, W.U64{0, 2147483648})
712	
713	def ri_fin(m: RMode, +x: F64, int: Bool) -> F64:
714	  match int:
715	    case True{}:
716	      x
717	    case False{}:
718	      ri_k(m, signbit(x), dmant(x), Nat.sub(3000n, dexp(x)), Nat.is_lt(Nat.sub(3000n, dexp(x)), 64n))
719	
720	def ri_cls(m: RMode, +x: F64, top: Bool) -> F64:
721	  match top:
722	    case True{}:
723	      nan_or(x, nz(frac(x)))
724	    case False{}:
725	      ri_fin(m, x, Nat.is_le(3000n, dexp(x)))
726	
727	def to_integral(m: RMode, +x: F64) -> F64:
728	  ri_cls(m, x, Nat.is_eq(exp_field(x), 2047n))
729	
730	def trunc(+x: F64) -> F64:
731	  to_integral(Trunc{}, x)
732	
733	def floor(+x: F64) -> F64:
734	  to_integral(Floor{}, x)
735	
736	def ceil(+x: F64) -> F64:
737	  to_integral(Ceil{}, x)
738	
739	# ties to even (Python's round(x) and C's rint)
740	def round(+x: F64) -> F64:
741	  to_integral(Even{}, x)
742	
743	# ---- conversions to and from unsigned integers ----
744	
745	def tu_neg(+w: W.U64, bad: Bool) -> Result<&2, &2, N.NumError, W.U64>:
746	  match bad:
747	    case True{}:
748	      Fail{N.Overflow{}}
749	    case False{}:
750	      Done{w}
751	
752	def tu_int(+s: Bool, +w: W.U64) -> Result<&2, &2, N.NumError, W.U64>:
753	  tu_neg(w, Bool.and(s, Bool.not(X.is_zero(w))))
754	
755	# |x| >= 2^52: w << k with k + bit_length(w) <= 64
756	def tu_big(+s: Bool, +w: W.U64, +k: Nat, fits: Bool) -> Result<&2, &2, N.NumError, W.U64>:
757	  match fits:
758	    case True{}:
759	      tu_int(s, X.shl(w, k))
760	    case False{}:
761	      Fail{N.Overflow{}}
762	
763	def tu_fin(+x: F64, big: Bool) -> Result<&2, &2, N.NumError, W.U64>:
764	  match big:
765	    case True{}:
766	      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))
767	    case False{}:
768	      tu_int(signbit(x), X.shr(dmant(x), Nat.sub(3000n, dexp(x))))
769	
770	def tu_top(isn: Bool) -> Result<&2, &2, N.NumError, W.U64>:
771	  match isn:
772	    case True{}:
773	      Fail{N.BadDomain{}}
774	    case False{}:
775	      Fail{N.Overflow{}}
776	
777	def tu_cls(+x: F64, top: Bool) -> Result<&2, &2, N.NumError, W.U64>:
778	  match top:
779	    case True{}:
780	      tu_top(nz(frac(x)))
781	    case False{}:
782	      tu_fin(x, Nat.is_le(3000n, dexp(x)))
783	
784	# x truncated toward zero as an unsigned 64-bit integer: NaN is a domain
785	# error, a value outside [0, 2^64) (infinities included) an overflow
786	def to_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:
787	  tu_cls(x, Nat.is_eq(exp_field(x), 2047n))
788	
789	def tu32_w(+w: W.U64, fits: Bool) -> Result<&2, &2, N.NumError, U32>:
790	  match fits:
791	    case True{}:
792	      Done{X.lo(w)}
793	    case False{}:
794	      Fail{N.Overflow{}}
795	
796	def tu32(r: Result<&2, &2, N.NumError, W.U64>) -> Result<&2, &2, N.NumError, U32>:
797	  match r:
798	    case Fail{e}:
799	      Fail{e}
800	    case Done{+w}:
801	      tu32_w(w, U32.is_zero(X.hi(w)))
802	
803	def to_u32(+x: F64) -> Result<&2, &2, N.NumError, U32>:
804	  tu32(to_u64(x))
805	
806	def floor_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:
807	  to_u64(floor(x))
808	
809	def ceil_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:
810	  to_u64(ceil(x))
811	
812	def round_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:
813	  to_u64(round(x))
814	
815	# the double nearest to an unsigned 64-bit integer (exact below 2^53)
816	def of_u64(+w: W.U64) -> F64:
817	  round_w(False{}, 3000n, w)
818	
819	def of_u32(+u: U32) -> F64:
820	  of_u64(W.U64{u, 0})
821	
822	# ---- exponents: frexp, ldexp, ulp ----
823	
824	# a signed exponent: -mag when neg
825	type Exp is Data:
826	  Exp{neg: Bool, mag: Nat}
827	
828	def exp_pick(+t: Nat, pos: Bool) -> Exp:
829	  match pos:
830	    case True{}:
831	      Exp{False{}, Nat.sub(t, 3000n)}
832	    case False{}:
833	      Exp{True{}, Nat.sub(3000n, t)}
834	
835	# t - 3000 as a signed exponent
836	def exp_of(+t: Nat) -> Exp:
837	  exp_pick(t, Nat.is_le(3000n, t))
838	
839	def fx_fin(+x: F64, +w: W.U64, +b: Nat) -> F64 & Exp:
840	  (round_w(signbit(x), Nat.sub(3000n, b), w), exp_of(Nat.add(dexp(x), b)))
841	
842	def fx_z(+x: F64, z: Bool) -> F64 & Exp:
843	  match z:
844	    case True{}:
845	      (x, Exp{False{}, 0n})
846	    case False{}:
847	      fx_fin(x, dmant(x), Nat.sub(64n, X.clz(dmant(x))))
848	
849	def fx_cls(+x: F64, top: Bool) -> F64 & Exp:
850	  match top:
851	    case True{}:
852	      (nan_or(x, nz(frac(x))), Exp{False{}, 0n})
853	    case False{}:
854	      fx_z(x, is_zero(x))
855	
856	# (m, e) with x = m * 2^e and 0.5 <= |m| < 1; zeros and infinities give (x, 0)
857	def frexp(+x: F64) -> F64 & Exp:
858	  fx_cls(x, Nat.is_eq(exp_field(x), 2047n))
859	
860	def ld_neg(+x: F64, +k: Nat, tiny: Bool) -> F64:
861	  match tiny:
862	    case True{}:
863	      zero(signbit(x))
864	    case False{}:
865	      round_w(signbit(x), Nat.sub(dexp(x), k), dmant(x))
866	
867	def ld_fin(+x: F64, +neg: Bool, +k: Nat) -> F64:
868	  match neg:
869	    case True{}:
870	      ld_neg(x, k, Nat.is_lt(dexp(x), Nat.add(k, 63n)))
871	    case False{}:
872	      round_w(signbit(x), Nat.add(dexp(x), k), dmant(x))
873	
874	def ld_e(+x: F64, +e: Exp) -> F64:
875	  match e:
876	    case Exp{+neg, +k}:
877	      ld_fin(x, neg, k)
878	
879	def ld_z(+x: F64, +e: Exp, z: Bool) -> F64:
880	  match z:
881	    case True{}:
882	      x
883	    case False{}:
884	      ld_e(x, e)
885	
886	def ld_cls(+x: F64, +e: Exp, top: Bool) -> F64:
887	  match top:
888	    case True{}:
889	      nan_or(x, nz(frac(x)))
890	    case False{}:
891	      ld_z(x, e, is_zero(x))
892	
893	# x * 2^e, rounded once (overflow to infinity, gradual underflow)
894	def ldexp(+x: F64, +e: Exp) -> F64:
895	  ld_cls(x, e, Nat.is_eq(exp_field(x), 2047n))
896	
897	def ulp_cls(+x: F64, top: Bool) -> F64:
898	  match top:
899	    case True{}:
900	      nan_or(inf(False{}), nz(frac(x)))
901	    case False{}:
902	      round_w(False{}, dexp(x), W.U64{1, 0})
903	
904	# the value of the least significant bit of x (Python's math.ulp)
905	def ulp(+x: F64) -> F64:
906	  ulp_cls(x, Nat.is_eq(exp_field(x), 2047n))
907	
908	# ---- neighbours and NaN-ignoring extrema ----
909	
910	def with_sign(+s: Bool, +w: W.U64) -> F64:
911	  Bits{X.lo(w), U32.add(X.hi(w), sgn(s))}
912	
913	def na_step(+x: F64, up: Bool) -> F64:
914	  match up:
915	    case True{}:
916	      with_sign(signbit(x), X.add(mag(x), W.U64{1, 0}))
917	    case False{}:
918	      with_sign(signbit(x), X.sub(mag(x), W.U64{1, 0}))
919	
920	def na_z(+x: F64, +y: F64, z: Bool) -> F64:
921	  match z:
922	    case True{}:
923	      Bits{1, sgn(signbit(y))}
924	    case False{}:
925	      na_step(x, Bool.xor(lt(x, y), signbit(x)))
926	
927	def na_eq(+x: F64, +y: F64, same: Bool) -> F64:
928	  match same:
929	    case True{}:
930	      y
931	    case False{}:
932	      na_z(x, y, is_zero(x))
933	
934	def na_nan(+x: F64, +y: F64, bad: Bool) -> F64:
935	  match bad:
936	    case True{}:
937	      nan()
938	    case False{}:
939	      na_eq(x, y, eq(x, y))
940	
941	# the next double after x toward y (IEEE nextAfter; y when x == y)
942	def nextafter(+x: F64, +y: F64) -> F64:
943	  na_nan(x, y, unordered(x, y))
944	
945	def fpick(+x: F64, +y: F64, first: Bool) -> F64:
946	  match first:
947	    case True{}:
948	      x
949	    case False{}:
950	      y
951	
952	def fmin_z(+x: F64, +y: F64, z: Bool) -> F64:
953	  match z:
954	    case True{}:
955	      zero(Bool.or(signbit(x), signbit(y)))
956	    case False{}:
957	      fpick(x, y, lt(x, y))
958	
959	def fmin_ny(+x: F64, +y: F64, ny: Bool) -> F64:
960	  match ny:
961	    case True{}:
962	      x
963	    case False{}:
964	      fmin_z(x, y, zeros(x, y))
965	
966	def fmin_nx(+x: F64, +y: F64, nx: Bool) -> F64:
967	  match nx:
968	    case True{}:
969	      nan_or(y, is_nan(y))
970	    case False{}:
971	      fmin_ny(x, y, is_nan(y))
972	
973	# IEEE 754-2019 minimumNumber: a NaN operand is ignored, -0 < +0
974	def fmin(+x: F64, +y: F64) -> F64:
975	  fmin_nx(x, y, is_nan(x))
976	
977	def fmax_z(+x: F64, +y: F64, z: Bool) -> F64:
978	  match z:
979	    case True{}:
980	      zero(Bool.and(signbit(x), signbit(y)))
981	    case False{}:
982	      fpick(x, y, lt(y, x))
983	
984	def fmax_ny(+x: F64, +y: F64, ny: Bool) -> F64:
985	  match ny:
986	    case True{}:
987	      x
988	    case False{}:
989	      fmax_z(x, y, zeros(x, y))
990	
991	def fmax_nx(+x: F64, +y: F64, nx: Bool) -> F64:
992	  match nx:
993	    case True{}:
994	      nan_or(y, is_nan(y))
995	    case False{}:
996	      fmax_ny(x, y, is_nan(y))
997	
998	# IEEE 754-2019 maximumNumber
999	def fmax(+x: F64, +y: F64) -> F64:
1000	  fmax_nx(x, y, is_nan(x))
1001	
1002	# ---- classification, bit casts, integrality ----
1003	
1004	def is_normal(+x: F64) -> Bool:
1005	  Bool.and(Bool.not(Nat.is_eq(exp_field(x), 0n)), Nat.is_lt(exp_field(x), 2047n))
1006	
1007	def is_subnormal(+x: F64) -> Bool:
1008	  Bool.and(Nat.is_eq(exp_field(x), 0n), nz(frac(x)))
1009	
1010	def to_bits(+x: F64) -> W.U64:
1011	  W.U64{word_lo(x), word_hi(x)}
1012	
1013	def of_bits64(+w: W.U64) -> F64:
1014	  Bits{X.lo(w), X.hi(w)}
1015	
1016	def ii_k(+w: W.U64, +k: Nat, small: Bool) -> Bool:
1017	  match small:
1018	    case True{}:
1019	      X.eq(X.shl(X.shr(w, k), k), w)
1020	    case False{}:
1021	      X.is_zero(w)
1022	
1023	def ii_fin(+x: F64, int: Bool) -> Bool:
1024	  match int:
1025	    case True{}:
1026	      True{}
1027	    case False{}:
1028	      ii_k(dmant(x), Nat.sub(3000n, dexp(x)), Nat.is_lt(Nat.sub(3000n, dexp(x)), 64n))
1029	
1030	# finite with no fractional bits
1031	def is_integer(+x: F64) -> Bool:
1032	  Bool.and(is_finite(x), ii_fin(x, Nat.is_le(3000n, dexp(x))))
1033	
1034	# ---- the fractional part and the exact remainders ----
1035	
1036	def mf_frac(+s: Bool, +w: W.U64, +u: Nat, +k: Nat, small: Bool) -> F64:
1037	  match small:
1038	    case True{}:
1039	      round_w(s, u, X.sub(w, X.shl(X.shr(w, k), k)))
1040	    case False{}:
1041	      round_w(s, u, w)
1042	
1043	def mf_fin(+x: F64, int: Bool) -> F64 & F64:
1044	  match int:
1045	    case True{}:
1046	      (zero(signbit(x)), x)
1047	    case False{}:
1048	      (mf_frac(signbit(x), dmant(x), dexp(x), Nat.sub(3000n, dexp(x)), Nat.is_lt(Nat.sub(3000n, dexp(x)), 64n)), trunc(x))
1049	
1050	def mf_cls(+x: F64, top: Bool) -> F64 & F64:
1051	  match top:
1052	    case True{}:
1053	      (nan_or(zero(signbit(x)), nz(frac(x))), nan_or(x, nz(frac(x))))
1054	    case False{}:
1055	      mf_fin(x, Nat.is_le(3000n, dexp(x)))
1056	
1057	# (fractional part, integral part), both with the sign of x (Python's modf)
1058	def modf(+x: F64) -> F64 & F64:
1059	  mf_cls(x, Nat.is_eq(exp_field(x), 2047n))
1060	
1061	# (q, r) of mx * 2^d by b from (q0, r0) of mx by b: ten bits a step, so each
1062	# partial remainder r * 2^t stays below 2^63
1063	def fm_go(fuel: Nat, +b: W.U64, +d: Nat, st: W.U64 & W.U64) -> W.U64 & W.U64:
1064	  match fuel d:
1065	    case 0n _:
1066	      st
1067	    case 1n+f 0n:
1068	      st
1069	    case 1n+f 1n+e:
1070	      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))
1071	
1072	# the last quotient digit and remainder of mant(x) * 2^(e(x) - e(y)) by mant(y)
1073	def fm_qr(+x: F64, +y: F64) -> W.U64 & W.U64:
1074	  fm_go(Nat.sub(dexp(x), dexp(y)), dmant(y), Nat.sub(dexp(x), dexp(y)), X.divmod(dmant(x), dmant(y)))
1075	
1076	def fmod_fin(+x: F64, +y: F64, far: Bool) -> F64:
1077	  match far:
1078	    case True{}:
1079	      x
1080	    case False{}:
1081	      round_w(signbit(x), dexp(y), X.psnd(fm_qr(x, y)))
1082	
1083	def fmod_z(+x: F64, +y: F64, zx: Bool) -> F64:
1084	  match zx:
1085	    case True{}:
1086	      x
1087	    case False{}:
1088	      fmod_fin(x, y, Nat.is_lt(dexp(x), dexp(y)))
1089	
1090	def fmod_yi(+x: F64, +y: F64, yi: Bool) -> F64:
1091	  match yi:
1092	    case True{}:
1093	      x
1094	    case False{}:
1095	      fmod_z(x, y, is_zero(x))
1096	
1097	def fmod_bad(+x: F64, +y: F64, bad: Bool) -> F64:
1098	  match bad:
1099	    case True{}:
1100	      nan()
1101	    case False{}:
1102	      fmod_yi(x, y, is_inf(y))
1103	
1104	# x - n y for n = trunc(x / y), exact, with the sign of x (C's fmod)
1105	def fmod(+x: F64, +y: F64) -> F64:
1106	  fmod_bad(x, y, Bool.or(Bool.or(unordered(x, y), is_inf(x)), is_zero(y)))
1107	
1108	def rm_pick(+s: Bool, +u: Nat, +b: W.U64, +r: W.U64, flip: Bool) -> F64:
1109	  match flip:
1110	    case True{}:
1111	      round_w(Bool.not(s), u, X.sub(b, r))
1112	    case False{}:
1113	      round_w(s, u, r)
1114	
1115	# remainder r of divisor b at scale u, last quotient digit q: round the
1116	# quotient to nearest (ties to even) by flipping to r - b
1117	def rm_fix(+s: Bool, +u: Nat, +b: W.U64, +r: W.U64, +q: W.U64) -> F64:
1118	  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))))
1119	
1120	def rm_qr(+x: F64, +y: F64, p: W.U64 & W.U64) -> F64:
1121	  (+q, +r) = p
1122	  rm_fix(signbit(x), dexp(y), dmant(y), r, q)
1123	
1124	def rm_near(+x: F64, +y: F64, one: Bool) -> F64:
1125	  match one:
1126	    case True{}:
1127	      rm_fix(signbit(x), dexp(x), X.add(dmant(y), dmant(y)), dmant(x), W.U64{0, 0})
1128	    case False{}:
1129	      x
1130	
1131	def rm_far(+x: F64, +y: F64, far: Bool) -> F64:
1132	  match far:
1133	    case True{}:
1134	      rm_near(x, y, Nat.is_eq(Nat.sub(dexp(y), dexp(x)), 1n))
1135	    case False{}:
1136	      rm_qr(x, y, fm_qr(x, y))
1137	
1138	def rm_z(+x: F64, +y: F64, zx: Bool) -> F64:
1139	  match zx:
1140	    case True{}:
1141	      x
1142	    case False{}:
1143	      rm_far(x, y, Nat.is_lt(dexp(x), dexp(y)))
1144	
1145	def rm_yi(+x: F64, +y: F64, yi: Bool) -> F64:
1146	  match yi:
1147	    case True{}:
1148	      x
1149	    case False{}:
1150	      rm_z(x, y, is_zero(x))
1151	
1152	def rm_bad(+x: F64, +y: F64, bad: Bool) -> F64:
1153	  match bad:
1154	    case True{}:
1155	      nan()
1156	    case False{}:
1157	      rm_yi(x, y, is_inf(y))
1158	
1159	# IEEE remainder: x - n y for n = x / y rounded to nearest, ties to even
1160	def remainder(+x: F64, +y: F64) -> F64:
1161	  rm_bad(x, y, Bool.or(Bool.or(unordered(x, y), is_inf(x)), is_zero(y)))
1162	
1163	# ---- closeness and exact ratios ----
1164	
1165	def ic_d(+a: F64, +b: F64, +rel: F64, +at: F64, +d: F64) -> Bool:
1166	  Bool.or(Bool.or(le(d, abs(mul(rel, b))), le(d, abs(mul(rel, a)))), le(d, at))
1167	
1168	def ic_inf(+a: F64, +b: F64, +rel: F64, +at: F64, big: Bool) -> Bool:
1169	  match big:
1170	    case True{}:
1171	      False{}
1172	    case False{}:
1173	      ic_d(a, b, rel, at, abs(sub(b, a)))
1174	
1175	def ic_eq(+a: F64, +b: F64, +rel: F64, +at: F64, same: Bool) -> Bool:
1176	  match same:
1177	    case True{}:
1178	      True{}
1179	    case False{}:
1180	      ic_inf(a, b, rel, at, Bool.or(is_inf(a), is_inf(b)))
1181	
1182	def ic_tol(+a: F64, +b: F64, +rel: F64, +at: F64, bad: Bool) -> Result<&2, &2, N.NumError, Bool>:
1183	  match bad:
1184	    case True{}:
1185	      Fail{N.BadDomain{}}
1186	    case False{}:
1187	      Done{ic_eq(a, b, rel, at, eq(a, b))}
1188	
1189	# Python's math.isclose (CPython's algorithm): a negative tolerance is a
1190	# domain error
1191	def isclose(+a: F64, +b: F64, +rel: F64, +at: F64) -> Result<&2, &2, N.NumError, Bool>:
1192	  ic_tol(a, b, rel, at, Bool.or(lt(rel, zero(False{})), lt(at, zero(False{}))))
1193	
1194	# +-num / 2^den in lowest terms (num odd or den = 0)
1195	type Ratio is Data:
1196	  Ratio{neg: Bool, num: W.U64, den: Nat}
1197	
1198	# trailing zeros of w (at most fuel), one halving at a time
1199	def ctz_go(fuel: Nat, +w: W.U64, odd: Bool) -> Nat:
1200	  match fuel odd:
1201	    case 0n _:
1202	      0n
1203	    case 1n+f True{}:
1204	      0n
1205	    case 1n+f False{}:
1206	      1n+ctz_go(f, X.half(w), X.odd(X.half(w)))
1207	
1208	def ctz(+w: W.U64) -> Nat:
1209	  ctz_go(64n, w, X.odd(w))
1210	
1211	def ar_small(+s: Bool, +w: W.U64, +k: Nat, +j: Nat) -> Result<&2, &2, N.NumError, Ratio>:
1212	  Done{Ratio{s, X.shr(w, j), Nat.sub(k, j)}}
1213	
1214	def ar_big(+s: Bool, +w: W.U64, +k: Nat, fits: Bool) -> Result<&2, &2, N.NumError, Ratio>:
1215	  match fits:
1216	    case True{}:
1217	      Done{Ratio{s, X.shl(w, k), 0n}}
1218	    case False{}:
1219	      Fail{N.Overflow{}}
1220	
1221	def ar_fin(+x: F64, big: Bool) -> Result<&2, &2, N.NumError, Ratio>:
1222	  match big:
1223	    case True{}:
1224	      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))
1225	    case False{}:
1226	      ar_small(signbit(x), dmant(x), Nat.sub(3000n, dexp(x)), Nat.min(ctz(dmant(x)), Nat.sub(3000n, dexp(x))))
1227	
1228	def ar_z(+x: F64, z: Bool) -> Result<&2, &2, N.NumError, Ratio>:
1229	  match z:
1230	    case True{}:
1231	      Done{Ratio{False{}, W.U64{0, 0}, 0n}}
1232	    case False{}:
1233	      ar_fin(x, Nat.is_le(3000n, dexp(x)))
1234	
1235	def ar_top(isn: Bool) -> Result<&2, &2, N.NumError, Ratio>:
1236	  match isn:
1237	    case True{}:
1238	      Fail{N.BadDomain{}}
1239	    case False{}:
1240	      Fail{N.Overflow{}}
1241	
1242	def ar_cls(+x: F64, top: Bool) -> Result<&2, &2, N.NumError, Ratio>:
1243	  match top:
1244	    case True{}:
1245	      ar_top(nz(frac(x)))
1246	    case False{}:
1247	      ar_z(x, is_zero(x))
1248	
1249	# x = +-num / 2^den exactly, in lowest terms (Python's as_integer_ratio with
1250	# the denominator's exponent); NaN is a domain error, infinities and integers
1251	# of 2^64 or more overflow
1252	def as_integer_ratio(+x: F64) -> Result<&2, &2, N.NumError, Ratio>:
1253	  ar_cls(x, Nat.is_eq(exp_field(x), 2047n))
1254	
1255	# ---- the num.bend instance ----
1256	
1257	def f64_op(o: N.Op<F64>) -> F64:
1258	  match o:
1259	    case N.ZeroOp{}:
1260	      zero(False{})
1261	    case N.One{}:
1262	      one()
1263	    case N.Add{a, b}:
1264	      add(a, b)
1265	    case N.Sub{a, b}:
1266	      sub(a, b)
1267	    case N.Mul{a, b}:
1268	      mul(a, b)
1269	    case N.Neg{a}:
1270	      neg(a)
1271	    case N.Abs{a}:
1272	      abs(a)
1273	    case N.Quot{a, b}:
1274	      div(a, b)
1275	    case N.Rem{a, b}:
1276	      nan()
1277	    case N.Half{a}:
1278	      mul(a, Bits{0, 1071644672})
1279	    case N.MulMod{a, b, m}:
1280	      nan()
1281	    case N.Pow2{k}:
1282	      pow2(k)
1283	    case N.Sqrt{a}:
1284	      sqrt(a)
1285	    case N.PowMod{a, e, m}:
1286	      nan()
1287	    case N.GcdSmall{a, b}:
1288	      nan()
1289	
1290	def f64_is(t: N.Test<F64>) -> Bool:
1291	  match t:
1292	    case N.Lt{a, b}:
1293	      lt(a, b)
1294	    case N.AddOver{a, b}:
1295	      False{}
1296	    case N.MulOver{a, b}:
1297	      False{}
1298	    case N.Odd{a}:
1299	      False{}
1300	    case N.IsZero{a}:
1301	      is_zero(a)
1302	    case N.FastDiv{}:
1303	      True{}
1304	    case N.Mont{m}:
1305	      False{}
1306	    case N.Small{a, b}:
1307	      False{}