src/math/f64.bend

bend-collections@1.0.0.0 · 1,307 lines · 221 definitions

1import Base
2import bend-collections@1.0.0.0/src/math/num.bend as N
3import bend-collections@1.0.0.0/src/math/u64.bend as W
4import bend-collections@1.0.0.0/src/math/w64.bend as X
5import bend-collections@1.0.0.0/src/math/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
29type F64 is Data:
30  Bits{lo: U32, hi: U32}
31
32def word_lo(x: F64) -> U32:
33  match x:
34    case Bits{l, h}:
35      l
36
37def word_hi(x: F64) -> U32:
38  match x:
39    case Bits{l, h}:
40      h
41
42def of_bits(+hi: U32, +lo: U32) -> F64:
43  Bits{lo, hi}
44
45def off() -> Nat:
46  4096n
47
48def sgn(s: Bool) -> U32:
49  match s:
50    case True{}:
51      2147483648
52    case False{}:
53      0
54
55# ---- fields ----
56
57def signbit(+x: F64) -> Bool:
58  U32.is_le(2147483648, word_hi(x))
59
60def exp_field(+x: F64) -> Nat:
61  U32.to_nat(U32.and(U32.div(word_hi(x), 1048576), 2047))
62
63def frac(+x: F64) -> W.U64:
64  W.U64{word_lo(x), U32.and(word_hi(x), 1048575)}
65
66def mag(+x: F64) -> W.U64:
67  W.U64{word_lo(x), U32.and(word_hi(x), 2147483647)}
68
69def nz(+a: W.U64) -> Bool:
70  Bool.not(X.is_zero(a))
71
72def is_nan(+x: F64) -> Bool:
73  Bool.and(Nat.is_eq(exp_field(x), 2047n), nz(frac(x)))
74
75def is_inf(+x: F64) -> Bool:
76  Bool.and(Nat.is_eq(exp_field(x), 2047n), X.is_zero(frac(x)))
77
78def is_finite(+x: F64) -> Bool:
79  Nat.is_lt(exp_field(x), 2047n)
80
81def is_zero(+x: F64) -> Bool:
82  X.is_zero(mag(x))
83
84# ---- constants and sign ----
85
86def nan() -> F64:
87  Bits{0, 2146959360}
88
89def inf(s: Bool) -> F64:
90  Bits{0, U32.add(2146435072, sgn(s))}
91
92def zero(s: Bool) -> F64:
93  Bits{0, sgn(s)}
94
95def one() -> F64:
96  Bits{0, 1072693248}
97
98def neg(+x: F64) -> F64:
99  Bits{word_lo(x), U32.xor(word_hi(x), 2147483648)}
100
101def abs(+x: F64) -> F64:
102  Bits{word_lo(x), U32.and(word_hi(x), 2147483647)}
103
104def 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
107def 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
116def 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
121def 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
124def zero_e(+e: Nat, z: Bool) -> Nat:
125  match z:
126    case True{}:
127      0n
128    case False{}:
129      e
130
131def rp_fin3(+s: Bool, +e: Nat, +r: W.U64) -> F64:
132  pack(s, zero_e(e, X.is_zero(r)), r)
133
134def 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
139def 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
142def 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
149def 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
157def 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
161def 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
168def rp62(+s: Bool, +e: Nat, +sig: W.U64) -> F64:
169  rp62_pick(s, e, sig, X.lt(sig, W.U64{0, 1073741824}))
170
171def 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
178def 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
185def 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)
189def 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)
194def 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
201def 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
210def 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
218def 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
225def 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
232def 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
239def 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
248def 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
251def 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
259def 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
266def 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
274def 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
277def 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
286def cmp64(+a: W.U64, +b: W.U64) -> Cmp:
287  X.cmp(a, b)
288
289def 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
296def 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
305def 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
308def 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
315def add(+x: F64, +y: F64) -> F64:
316  add_pick(x, y, Bool.not(Bool.xor(signbit(x), signbit(y))))
317
318def sub(+x: F64, +y: F64) -> F64:
319  add(x, neg(y))
320
321# ---- multiplication (SoftFloat's f64_mul) ----
322
323def 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
327def 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
330def fin_zero(+e: Nat, +f: W.U64) -> Bool:
331  Bool.and(Nat.is_eq(e, 0n), X.is_zero(f))
332
333def 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
340def inf_or_nan(+s: Bool, bad: Bool) -> F64:
341  match bad:
342    case True{}:
343      nan()
344    case False{}:
345      inf(s)
346
347def 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
354def 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
361def 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
367def 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
370def 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
373def 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
377def 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
382def 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
385def 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
388def 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
395def 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
398def 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
405def 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
412def 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
419def 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
426def 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
431def 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
435def 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)
439def 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
446def 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
450def 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
457def 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
461def 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
468def 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
476def 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
479def 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
484def 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
487def 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
494def 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
501def 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
508def 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
515def 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
520def 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
531def 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
542def unordered(+x: F64, +y: F64) -> Bool:
543  Bool.or(is_nan(x), is_nan(y))
544
545def zeros(+x: F64, +y: F64) -> Bool:
546  Bool.and(is_zero(x), is_zero(y))
547
548def 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
557def lt(+x: F64, +y: F64) -> Bool:
558  lt_z(x, y, unordered(x, y), zeros(x, y))
559
560def 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
569def le(+x: F64, +y: F64) -> Bool:
570  le_z(x, y, unordered(x, y), zeros(x, y))
571
572def 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
581def 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)
586def 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
593def 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)
601def 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
604def 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
609def 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
616def 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)
624def of_nat(+n: Nat) -> F64:
625  of_nat_z(n, Nat.is_eq(n, 0n))
626
627# 2^k for k <= 1023
628def 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
636def 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
643def dmant(+x: F64) -> W.U64:
644  dmant_z(frac(x), Nat.is_eq(exp_field(x), 0n))
645
646def dexp_z(+e: Nat, z: Bool) -> Nat:
647  match z:
648    case True{}:
649      1926n
650    case False{}:
651      Nat.add(e, 1925n)
652
653def dexp(+x: F64) -> Nat:
654  dexp_z(exp_field(x), Nat.is_eq(exp_field(x), 0n))
655
656def 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
663def 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
672def 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
677type RMode is Data:
678  Trunc{}
679  Floor{}
680  Ceil{}
681  Even{}
682
683def 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
691def 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
702def 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)
706def 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
713def 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
720def 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
727def to_integral(m: RMode, +x: F64) -> F64:
728  ri_cls(m, x, Nat.is_eq(exp_field(x), 2047n))
729
730def trunc(+x: F64) -> F64:
731  to_integral(Trunc{}, x)
732
733def floor(+x: F64) -> F64:
734  to_integral(Floor{}, x)
735
736def ceil(+x: F64) -> F64:
737  to_integral(Ceil{}, x)
738
739# ties to even (Python's round(x) and C's rint)
740def round(+x: F64) -> F64:
741  to_integral(Even{}, x)
742
743# ---- conversions to and from unsigned integers ----
744
745def 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
752def 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
756def 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
763def 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
770def 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
777def 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
786def to_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:
787  tu_cls(x, Nat.is_eq(exp_field(x), 2047n))
788
789def 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
796def 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
803def to_u32(+x: F64) -> Result<&2, &2, N.NumError, U32>:
804  tu32(to_u64(x))
805
806def floor_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:
807  to_u64(floor(x))
808
809def ceil_u64(+x: F64) -> Result<&2, &2, N.NumError, W.U64>:
810  to_u64(ceil(x))
811
812def 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)
816def of_u64(+w: W.U64) -> F64:
817  round_w(False{}, 3000n, w)
818
819def 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
825type Exp is Data:
826  Exp{neg: Bool, mag: Nat}
827
828def 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
836def exp_of(+t: Nat) -> Exp:
837  exp_pick(t, Nat.is_le(3000n, t))
838
839def 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
842def 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
849def 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)
857def frexp(+x: F64) -> F64 & Exp:
858  fx_cls(x, Nat.is_eq(exp_field(x), 2047n))
859
860def 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
867def 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
874def ld_e(+x: F64, +e: Exp) -> F64:
875  match e:
876    case Exp{+neg, +k}:
877      ld_fin(x, neg, k)
878
879def 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
886def 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)
894def ldexp(+x: F64, +e: Exp) -> F64:
895  ld_cls(x, e, Nat.is_eq(exp_field(x), 2047n))
896
897def 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)
905def ulp(+x: F64) -> F64:
906  ulp_cls(x, Nat.is_eq(exp_field(x), 2047n))
907
908# ---- neighbours and NaN-ignoring extrema ----
909
910def with_sign(+s: Bool, +w: W.U64) -> F64:
911  Bits{X.lo(w), U32.add(X.hi(w), sgn(s))}
912
913def 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
920def 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
927def 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
934def 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)
942def nextafter(+x: F64, +y: F64) -> F64:
943  na_nan(x, y, unordered(x, y))
944
945def fpick(+x: F64, +y: F64, first: Bool) -> F64:
946  match first:
947    case True{}:
948      x
949    case False{}:
950      y
951
952def 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
959def 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
966def 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
974def fmin(+x: F64, +y: F64) -> F64:
975  fmin_nx(x, y, is_nan(x))
976
977def 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
984def 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
991def 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
999def fmax(+x: F64, +y: F64) -> F64:
1000  fmax_nx(x, y, is_nan(x))
1001
1002# ---- classification, bit casts, integrality ----
1003
1004def 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
1007def is_subnormal(+x: F64) -> Bool:
1008  Bool.and(Nat.is_eq(exp_field(x), 0n), nz(frac(x)))
1009
1010def to_bits(+x: F64) -> W.U64:
1011  W.U64{word_lo(x), word_hi(x)}
1012
1013def of_bits64(+w: W.U64) -> F64:
1014  Bits{X.lo(w), X.hi(w)}
1015
1016def 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
1023def 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
1031def 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
1036def 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
1043def 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
1050def 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)
1058def 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
1063def 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)
1073def 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
1076def 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
1083def 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
1090def 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
1097def 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)
1105def 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
1108def 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
1117def 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
1120def 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
1124def 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
1131def 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
1138def 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
1145def 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
1152def 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
1160def 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
1165def 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
1168def 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
1175def 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
1182def 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
1191def 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)
1195type 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
1199def 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
1208def ctz(+w: W.U64) -> Nat:
1209  ctz_go(64n, w, X.odd(w))
1210
1211def 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
1214def 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
1221def 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
1228def 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
1235def 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
1242def 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
1252def 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
1257def 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
1290def 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{}