src/math/generic.bend

bend-collections-laws-crypto@1.0.0.0 · 610 lines · 97 definitions

1import Base
2import bend-collections-laws-crypto@1.0.0.0/src/math/num.bend as N
3
4# The math of natural.bend written once for any numeric type. Every function
5# takes the type and its interface as templates (~T, ~op, ~test; num.bend), so
6# each instance (instances.bend: U32, U64, F32; f64.bend: F64) compiles to its
7# own copy with the operations inlined:
8#
9#   G.gcd(~U32, ~I.u32_op, ~I.u32_is, a, b)
10#
11#   ordered-ring functions (every instance, floats included):
12#     min, max, clamp, abs, sign, sum, prod, pow (exponent a Nat)
13#   integer functions (the unsigned instances):
14#     gcd, lcm, gcd_all, lcm_all, bit_length, isqrt, iroot, ilog,
15#     factorial, perm, comb, pow_mod, mod_inverse, divmod
16#
17# Semantics are Python's, with checked results for fixed widths (design
18# reference section 2.6): a result that does not fit is Fail{Overflow}, never
19# a wrapped value. comb reduces by a gcd at every step, so it overflows only
20# when the result does; pow_mod and mod_inverse use the instance's MulMod, so
21# no intermediate leaves [0, m). Loops are bounded by fuel for 64-bit widths.
22# The Nat versions in natural.bend are the proved reference; these are tested
23# against Python (tools/check_generic.py). sum is a left fold (Python's float
24# sum is compensated; use it for exact float sums once fsum exists).
25
26type QuotRem<-T: Data> is Data:
27  TQR{quot: T, rem: T}
28
29def ok(-T: Data, m: Maybe<&2, T>) -> Result<&2, &2, N.NumError, T>:
30  match m:
31    case None{}:
32      Fail{N.Overflow{}}
33    case Some{x}:
34      Done{x}
35
36def or_zero(-T: Data, +z: T, m: Maybe<&2, T>) -> T:
37  match m:
38    case None{}:
39      z
40    case Some{x}:
41      x
42
43# enough steps for every loop at 64 bits (Euclid < 1.5 w + 2, bit loops w)
44def fuel() -> Nat:
45  140n
46
47# x when it fits, else None
48def fits(-T: Data, over: Bool, +x: T) -> Maybe<&2, T>:
49  match over:
50    case True{}:
51      None{}
52    case False{}:
53      Some{x}
54
55def cadd(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Maybe<&2, T>:
56  fits(T, test(N.AddOver{a, b}), op(N.Add{a, b}))
57
58def cmul(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Maybe<&2, T>:
59  fits(T, test(N.MulOver{a, b}), op(N.Mul{a, b}))
60
61def lt(~T: Data, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Bool:
62  test(N.Lt{a, b})
63
64def le(~T: Data, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Bool:
65  Bool.not(test(N.Lt{b, a}))
66
67def pos(~T: Data, ~test: N.Test<T> -> Bool, +a: T) -> Bool:
68  Bool.not(test(N.IsZero{a}))
69
70# ---- ordered-ring functions ----
71
72def pick(-T: Data, c: Bool, +a: T, +b: T) -> T:
73  match c:
74    case True{}:
75      a
76    case False{}:
77      b
78
79# Python's min(a, b) keeps a unless b < a; max keeps a unless a < b.
80def min(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> T:
81  pick(T, test(N.Lt{b, a}), b, a)
82
83def max(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> T:
84  pick(T, test(N.Lt{a, b}), b, a)
85
86def clamp_ok(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T, +lo: T, +hi: T, bad: Bool) -> Result<&2, &2, N.NumError, T>:
87  match bad:
88    case True{}:
89      Fail{N.BadDomain{}}
90    case False{}:
91      Done{min(~T, ~op, ~test, max(~T, ~op, ~test, x, lo), hi)}
92
93# min(max(x, lo), hi); hi < lo is a Domain error
94def clamp(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T, +lo: T, +hi: T) -> Result<&2, &2, N.NumError, T>:
95  clamp_ok(~T, ~op, ~test, x, lo, hi, test(N.Lt{hi, lo}))
96
97def abs(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T) -> T:
98  op(N.Abs{x})
99
100def sign_below(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T, below: Bool) -> T:
101  match below:
102    case True{}:
103      op(N.Neg{op(N.One{})})
104    case False{}:
105      x
106
107def sign_above(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T, above: Bool) -> T:
108  match above:
109    case True{}:
110      op(N.One{})
111    case False{}:
112      sign_below(~T, ~op, ~test, x, test(N.Lt{x, op(N.ZeroOp{})}))
113
114# 1 above zero, -1 below, x itself otherwise (0, -0.0, NaN)
115def sign(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T) -> T:
116  sign_above(~T, ~op, ~test, x, test(N.Lt{op(N.ZeroOp{}), x}))
117
118def sum_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>, acc: Maybe<&2, T>) -> Maybe<&2, T>:
119  match xs acc:
120    case Nil{} _:
121      acc
122    case Con{x, rest} None{}:
123      None{}
124    case Con{+x, rest} Some{+a}:
125      sum_go(~T, ~op, ~test, rest, cadd(~T, ~op, ~test, a, x))
126
127def sum(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>) -> Result<&2, &2, N.NumError, T>:
128  ok(T, sum_go(~T, ~op, ~test, xs, Some{op(N.ZeroOp{})}))
129
130def prod_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>, acc: Maybe<&2, T>) -> Maybe<&2, T>:
131  match xs acc:
132    case Nil{} _:
133      acc
134    case Con{x, rest} None{}:
135      None{}
136    case Con{+x, rest} Some{+a}:
137      prod_go(~T, ~op, ~test, rest, cmul(~T, ~op, ~test, a, x))
138
139def prod(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>) -> Result<&2, &2, N.NumError, T>:
140  ok(T, prod_go(~T, ~op, ~test, xs, Some{op(N.One{})}))
141
142# (value, still fits) after multiplying by x when bit is set
143def mul_st(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, bit: Bool, +x: T, +xok: Bool, st: T & Bool) -> T & Bool:
144  match bit st:
145    case False{} _:
146      st
147    case True{} Tuple{+a, +aok}:
148      (op(N.Mul{a, x}), Bool.and(Bool.and(aok, xok), Bool.not(test(N.MulOver{a, x}))))
149
150# the base squared, only while higher exponent bits remain
151def sq_val(~T: Data, ~op: N.Op<T> -> T, more: Bool, +b: T) -> T:
152  match more:
153    case False{}:
154      b
155    case True{}:
156      op(N.Mul{b, b})
157
158def sq_ok(~T: Data, ~test: N.Test<T> -> Bool, more: Bool, +b: T, +bok: Bool) -> Bool:
159  match more:
160    case False{}:
161      bok
162    case True{}:
163      Bool.and(bok, Bool.not(test(N.MulOver{b, b})))
164
165def st_fin(-T: Data, st: T & Bool) -> Maybe<&2, T>:
166  match st:
167    case Tuple{+a, +aok}:
168      fits(T, Bool.not(aok), a)
169
170# square-and-multiply, low bits first; an overflowing square means the
171# result overflows too (a later bit multiplies it in). Callers pass fuel
172# 1 + k: every step halves k, so the loop always ends on k == 0.
173def pow_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +k: Nat, +b: T, +bok: Bool, st: T & Bool) -> Maybe<&2, T>:
174  match fuel k:
175    case 0n _:
176      st_fin(T, st)
177    case 1n+f 0n:
178      st_fin(T, st)
179    case 1n+f 1n+j:
180      pow_go(~T, ~op, ~test, f, Nat.div(1n+j, 2n), sq_val(~T, ~op, Nat.is_lt(1n, 1n+j), b), sq_ok(~T, ~test, Nat.is_lt(1n, 1n+j), b, bok), mul_st(~T, ~op, ~test, Nat.is_eq(Nat.mod(1n+j, 2n), 1n), b, bok, st))
181
182def pow(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T, +k: Nat) -> Result<&2, &2, N.NumError, T>:
183  ok(T, pow_go(~T, ~op, ~test, 1n+k, k, x, True{}, (op(N.One{}), True{})))
184
185# ---- integer functions (unsigned) ----
186
187def one(~T: Data, ~op: N.Op<T> -> T) -> T:
188  op(N.One{})
189
190def inc(~T: Data, ~op: N.Op<T> -> T, +x: T) -> T:
191  op(N.Add{x, op(N.One{})})
192
193def dec(~T: Data, ~op: N.Op<T> -> T, +x: T) -> T:
194  op(N.Sub{x, op(N.One{})})
195
196# (x, x == 0)
197def zst(~T: Data, ~test: N.Test<T> -> Bool, +x: T) -> T & Bool:
198  (x, test(N.IsZero{x}))
199
200# Euclid's algorithm on (a, (b, b == 0))
201def gcd_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +a: T, st: T & Bool) -> T:
202  match fuel st:
203    case 0n _:
204      a
205    case 1n+f Tuple{b, True{}}:
206      a
207    case 1n+f Tuple{+b, False{}}:
208      gcd_go(~T, ~op, ~test, f, b, zst(~T, ~test, op(N.Rem{a, b})))
209
210# ---- binary gcd (Stein): halve, subtract, compare; no division ----
211#
212# Knuth TAOCP vol. 2, 4.5.2 Algorithm B; proved in proofs/math/typed/bgcd*.bend.
213
214def hst(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T) -> T & Bool:
215  (op(N.Half{x}), test(N.Odd{op(N.Half{x})}))
216
217# x with its factors of two removed (x != 0)
218def strip_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, st: T & Bool) -> T:
219  match fuel st:
220    case 0n Tuple{x, o}:
221      x
222    case 1n+f Tuple{x, True{}}:
223      x
224    case 1n+f Tuple{+x, False{}}:
225      strip_go(~T, ~op, ~test, f, hst(~T, ~op, ~test, x))
226
227def strip(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T) -> T:
228  strip_go(~T, ~op, ~test, fuel(), (x, test(N.Odd{x})))
229
230# (d, d == 0, a and d small)
231def zsm(~T: Data, ~test: N.Test<T> -> Bool, +a: T, +d: T) -> T & (Bool & Bool):
232  (d, (test(N.IsZero{d}), test(N.Small{a, d})))
233
234def bnx2(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T, less: Bool) -> T & (T & (Bool & Bool)):
235  match less:
236    case True{}:
237      (b, zsm(~T, ~test, b, op(N.Sub{a, b})))
238    case False{}:
239      (a, zsm(~T, ~test, a, op(N.Sub{b, a})))
240
241def bnx(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> T & (T & (Bool & Bool)):
242  bnx2(~T, ~op, ~test, a, b, test(N.Lt{b, a}))
243
244# (a odd, (d, d == 0, small)): once a and d fit the instance's small type
245# (Small), their gcd is taken there (GcdSmall: U64 runs Euclid on U32)
246def bloop(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, st: T & (T & (Bool & Bool))) -> T:
247  match fuel st:
248    case 0n Tuple{a, x}:
249      a
250    case 1n+f Tuple{a, Tuple{d, Tuple{True{}, sm}}}:
251      a
252    case 1n+f Tuple{+a, Tuple{+d, Tuple{False{}, True{}}}}:
253      op(N.GcdSmall{a, d})
254    case 1n+f Tuple{+a, Tuple{+d, Tuple{False{}, False{}}}}:
255      bloop(~T, ~op, ~test, f, bnx(~T, ~op, ~test, a, strip(~T, ~op, ~test, d)))
256
257def bl_start(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> T:
258  bloop(~T, ~op, ~test, fuel(), (a, zsm(~T, ~test, a, b)))
259
260def dbl(~T: Data, ~op: N.Op<T> -> T, c: Nat, +x: T) -> T:
261  match c:
262    case 0n:
263      x
264    case 1n+k:
265      dbl(~T, ~op, k, op(N.Add{x, x}))
266
267def ev2(~T: Data, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Bool:
268  Bool.not(Bool.or(test(N.Odd{a}), test(N.Odd{b})))
269
270def tw(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fl: Nat, +a: T, +b: T, +c: Nat, both: Bool) -> T:
271  match fl both:
272    case 0n _:
273      dbl(~T, ~op, c, bl_start(~T, ~op, ~test, strip(~T, ~op, ~test, a), b))
274    case 1n+f False{}:
275      dbl(~T, ~op, c, bl_start(~T, ~op, ~test, strip(~T, ~op, ~test, a), b))
276    case 1n+f True{}:
277      tw(~T, ~op, ~test, f, op(N.Half{a}), op(N.Half{b}), 1n+c, ev2(~T, ~test, op(N.Half{a}), op(N.Half{b})))
278
279def bg_b(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T, bz: Bool) -> T:
280  match bz:
281    case True{}:
282      a
283    case False{}:
284      tw(~T, ~op, ~test, fuel(), a, b, 0n, ev2(~T, ~test, a, b))
285
286def bg_a(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T, az: Bool) -> T:
287  match az:
288    case True{}:
289      b
290    case False{}:
291      bg_b(~T, ~op, ~test, a, b, test(N.IsZero{b}))
292
293def gcd_rb(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T, bz: Bool) -> T:
294  match bz:
295    case True{}:
296      a
297    case False{}:
298      bg_a(~T, ~op, ~test, b, op(N.Rem{a, b}), test(N.IsZero{b}))
299
300# one Euclid step first (so a tiny b costs one division, not a bit per
301# subtraction), then the binary loop
302def gcd_bin(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> T:
303  gcd_rb(~T, ~op, ~test, a, b, test(N.IsZero{b}))
304
305def gcd_pick(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T, fast: Bool) -> T:
306  match fast:
307    case True{}:
308      gcd_go(~T, ~op, ~test, fuel(), a, zst(~T, ~test, b))
309    case False{}:
310      gcd_bin(~T, ~op, ~test, a, b)
311
312# Euclid when division is native, the binary gcd otherwise
313def gcd(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> T:
314  gcd_pick(~T, ~op, ~test, a, b, test(N.FastDiv{}))
315
316def lcm_ok(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T, zero: Bool) -> Result<&2, &2, N.NumError, T>:
317  match zero:
318    case True{}:
319      Done{op(N.ZeroOp{})}
320    case False{}:
321      ok(T, cmul(~T, ~op, ~test, op(N.Quot{a, gcd(~T, ~op, ~test, a, b)}), b))
322
323def lcm(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Result<&2, &2, N.NumError, T>:
324  lcm_ok(~T, ~op, ~test, a, b, Bool.or(test(N.IsZero{a}), test(N.IsZero{b})))
325
326def gcd_all_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>, +acc: T) -> T:
327  match xs:
328    case Nil{}:
329      acc
330    case Con{+x, rest}:
331      gcd_all_go(~T, ~op, ~test, rest, gcd(~T, ~op, ~test, acc, x))
332
333def gcd_all(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>) -> T:
334  gcd_all_go(~T, ~op, ~test, xs, op(N.ZeroOp{}))
335
336def lcm_all_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>, acc: Result<&2, &2, N.NumError, T>) -> Result<&2, &2, N.NumError, T>:
337  match xs acc:
338    case Nil{} _:
339      acc
340    case Con{x, rest} Fail{e}:
341      Fail{e}
342    case Con{+x, rest} Done{+a}:
343      lcm_all_go(~T, ~op, ~test, rest, lcm(~T, ~op, ~test, a, x))
344
345def lcm_all(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, xs: List<&2, T>) -> Result<&2, &2, N.NumError, T>:
346  lcm_all_go(~T, ~op, ~test, xs, Done{op(N.One{})})
347
348# the number of halvings to reach zero
349def bit_length_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +k: Nat, st: T & Bool) -> Nat:
350  match fuel st:
351    case 0n _:
352      k
353    case 1n+f Tuple{n, True{}}:
354      k
355    case 1n+f Tuple{+n, False{}}:
356      bit_length_go(~T, ~op, ~test, f, 1n+k, zst(~T, ~test, op(N.Half{n})))
357
358def bit_length(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T) -> Nat:
359  bit_length_go(~T, ~op, ~test, fuel(), 0n, zst(~T, ~test, n))
360
361# the instance's integer square root (a hardware estimate and an exact
362# integer correction for U32 and U64)
363def isqrt(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T) -> T:
364  op(N.Sqrt{n})
365
366# ---- roots and logarithms ----
367
368# r^k <= n, with an overflowing power counted as > n
369def root_le_fin(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, p: Maybe<&2, T>) -> Bool:
370  match p:
371    case None{}:
372      False{}
373    case Some{+x}:
374      le(~T, ~test, x, n)
375
376def root_le(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: Nat, +r: T) -> Bool:
377  root_le_fin(~T, ~op, ~test, n, pow_go(~T, ~op, ~test, 1n+k, k, r, True{}, (op(N.One{}), True{})))
378
379# lo + (hi - lo) / 2, never above hi
380def imid(~T: Data, ~op: N.Op<T> -> T, +lo: T, +hi: T) -> T:
381  op(N.Add{lo, op(N.Half{op(N.Sub{hi, lo})})})
382
383# the last r in [lo, hi) with r^k <= n, given it holds at lo and fails at
384# hi: bisection while hi - lo > 1
385def search_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fu: Nat, +n: T, +k: Nat, +lo: T, +hi: T, +m: T, +more: Bool, +hit: Bool) -> T:
386  match fu more hit:
387    case 0n _ _:
388      lo
389    case 1n+f False{} _:
390      lo
391    case 1n+f True{} True{}:
392      search_go(~T, ~op, ~test, f, n, k, m, hi, imid(~T, ~op, m, hi), test(N.Lt{inc(~T, ~op, m), hi}), root_le(~T, ~op, ~test, n, k, imid(~T, ~op, m, hi)))
393    case 1n+f True{} False{}:
394      search_go(~T, ~op, ~test, f, n, k, lo, m, imid(~T, ~op, lo, m), test(N.Lt{inc(~T, ~op, lo), m}), root_le(~T, ~op, ~test, n, k, imid(~T, ~op, lo, m)))
395
396def search(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: Nat, +hi: T) -> T:
397  search_go(~T, ~op, ~test, fuel(), n, k, op(N.ZeroOp{}), hi, imid(~T, ~op, op(N.ZeroOp{}), hi), test(N.Lt{op(N.One{}), hi}), root_le(~T, ~op, ~test, n, k, imid(~T, ~op, op(N.ZeroOp{}), hi)))
398
399# 2^(bits(n) / k + 1) > the k-th root of n
400def iroot_k(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: Nat) -> T:
401  match k:
402    case 0n:
403      op(N.ZeroOp{})
404    case 1n:
405      n
406    case 2n+ +kp:
407      search(~T, ~op, ~test, n, 2n+kp, op(N.Pow2{1n+Nat.div(bit_length(~T, ~op, ~test, n), 2n+kp)}))
408
409def iroot(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: Nat) -> Result<&2, &2, N.NumError, T>:
410  match k:
411    case 0n:
412      Fail{N.BadDomain{}}
413    case 1n+ +kp:
414      Done{iroot_k(~T, ~op, ~test, n, 1n+kp)}
415
416# the largest k with b^k <= n; p <= n / b is asked before p * b is formed
417def ilog_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +q: T, +b: T, +k: Nat, +p: T, up: Bool) -> Nat:
418  match fuel up:
419    case 0n _:
420      k
421    case 1n+f False{}:
422      k
423    case 1n+f True{}:
424      ilog_go(~T, ~op, ~test, f, q, b, 1n+k, op(N.Mul{p, b}), le(~T, ~test, op(N.Mul{p, b}), q))
425
426def ilog_ok(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +b: T, bad: Bool) -> Result<&2, &2, N.NumError, Nat>:
427  match bad:
428    case True{}:
429      Fail{N.BadDomain{}}
430    case False{}:
431      Done{ilog_go(~T, ~op, ~test, fuel(), op(N.Quot{n, b}), b, 0n, op(N.One{}), le(~T, ~test, op(N.One{}), op(N.Quot{n, b})))}
432
433def ilog(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +b: T) -> Result<&2, &2, N.NumError, Nat>:
434  ilog_ok(~T, ~op, ~test, n, b, Bool.or(test(N.IsZero{n}), test(N.Lt{b, op(N.Add{op(N.One{}), op(N.One{})})})))
435
436# ---- factorial, permutations, combinations ----
437
438# (acc * x, still fits)
439def mul_by(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T, st: T & Bool) -> T & Bool:
440  match st:
441    case Tuple{+a, +aok}:
442      (op(N.Mul{a, x}), Bool.and(aok, Bool.not(test(N.MulOver{a, x}))))
443
444# 2 * 3 * ... * n: i runs up to n, stopping at the first overflow
445def fact_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +i: T, +n: T, st: T & Bool, more: Bool) -> Maybe<&2, T>:
446  match fuel st more:
447    case 0n _ _:
448      None{}
449    case 1n+f Tuple{a, False{}} _:
450      None{}
451    case 1n+f Tuple{+a, True{}} False{}:
452      Some{a}
453    case 1n+f Tuple{+a, True{}} True{}:
454      fact_go(~T, ~op, ~test, f, inc(~T, ~op, i), n, mul_by(~T, ~op, ~test, inc(~T, ~op, i), (a, True{})), test(N.Lt{inc(~T, ~op, i), n}))
455
456def factorial(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T) -> Result<&2, &2, N.NumError, T>:
457  ok(T, fact_go(~T, ~op, ~test, fuel(), op(N.One{}), n, (op(N.One{}), True{}), test(N.Lt{op(N.One{}), n})))
458
459# n (n - 1) ... (n - k + 1): m runs down while k counts the factors left
460def perm_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +k: T, +m: T, st: T & Bool, more: Bool) -> Maybe<&2, T>:
461  match fuel st more:
462    case 0n _ _:
463      None{}
464    case 1n+f Tuple{a, False{}} _:
465      None{}
466    case 1n+f Tuple{+a, True{}} False{}:
467      Some{a}
468    case 1n+f Tuple{+a, True{}} True{}:
469      perm_go(~T, ~op, ~test, f, dec(~T, ~op, k), dec(~T, ~op, m), mul_by(~T, ~op, ~test, m, (a, True{})), Bool.not(test(N.IsZero{dec(~T, ~op, k)})))
470
471def perm_big(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: T, big: Bool) -> Result<&2, &2, N.NumError, T>:
472  match big:
473    case True{}:
474      Done{op(N.ZeroOp{})}
475    case False{}:
476      ok(T, perm_go(~T, ~op, ~test, fuel(), k, n, (op(N.One{}), True{}), Bool.not(test(N.IsZero{k}))))
477
478def perm(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: T) -> Result<&2, &2, N.NumError, T>:
479  perm_big(~T, ~op, ~test, n, k, test(N.Lt{n, k}))
480
481# C(n, i+1) = (r / g) * ((n - i) / ((i + 1) / g)) with r = C(n, i) and
482# g = gcd(r, i + 1): exact, and it overflows only if C(n, i+1) does (C(n, i)
483# grows for i < n / 2)
484def comb_mul(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +r: T, +t: T, +j: T, +g: T) -> T & Bool:
485  mul_by(~T, ~op, ~test, op(N.Quot{t, op(N.Quot{j, g})}), (op(N.Quot{r, g}), True{}))
486
487def comb_next(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +r: T, +t: T, +j: T) -> T & Bool:
488  comb_mul(~T, ~op, ~test, r, t, j, gcd(~T, ~op, ~test, r, j))
489
490# i runs up to k (k <= n / 2), r = C(n, i)
491def comb_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +n: T, +i: T, +k: T, st: T & Bool, more: Bool) -> Maybe<&2, T>:
492  match fuel st more:
493    case 0n _ _:
494      None{}
495    case 1n+f Tuple{r, False{}} _:
496      None{}
497    case 1n+f Tuple{+r, True{}} False{}:
498      Some{r}
499    case 1n+f Tuple{+r, True{}} True{}:
500      comb_go(~T, ~op, ~test, f, n, inc(~T, ~op, i), k, comb_next(~T, ~op, ~test, r, op(N.Sub{n, i}), inc(~T, ~op, i)), test(N.Lt{inc(~T, ~op, i), k}))
501
502def comb_k(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: T) -> Result<&2, &2, N.NumError, T>:
503  ok(T, comb_go(~T, ~op, ~test, fuel(), n, op(N.ZeroOp{}), k, (op(N.One{}), True{}), test(N.Lt{op(N.ZeroOp{}), k})))
504
505def comb_big(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: T, big: Bool) -> Result<&2, &2, N.NumError, T>:
506  match big:
507    case True{}:
508      Done{op(N.ZeroOp{})}
509    case False{}:
510      comb_k(~T, ~op, ~test, n, min(~T, ~op, ~test, k, op(N.Sub{n, k})))
511
512def comb(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +n: T, +k: T) -> Result<&2, &2, N.NumError, T>:
513  comb_big(~T, ~op, ~test, n, k, test(N.Lt{n, k}))
514
515# ---- modular arithmetic ----
516
517def mm_bit(~T: Data, ~op: N.Op<T> -> T, odd: Bool, +m: T, +b: T, +acc: T) -> T:
518  match odd:
519    case False{}:
520      acc
521    case True{}:
522      op(N.MulMod{acc, b, m})
523
524# right-to-left binary exponentiation on (e, e == 0), every product reduced
525# mod m
526def pow_mod_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +m: T, +b: T, +acc: T, st: T & Bool) -> T:
527  match fuel st:
528    case 0n _:
529      acc
530    case 1n+f Tuple{e, True{}}:
531      acc
532    case 1n+f Tuple{+e, False{}}:
533      pow_mod_go(~T, ~op, ~test, f, m, op(N.MulMod{b, b, m}), mm_bit(~T, ~op, test(N.Odd{e}), m, b, acc), zst(~T, ~test, op(N.Half{e})))
534
535def pow_mod_pick(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +b: T, +e: T, +m: T, mont: Bool) -> T:
536  match mont:
537    case True{}:
538      op(N.PowMod{op(N.Rem{b, m}), e, m})
539    case False{}:
540      pow_mod_go(~T, ~op, ~test, fuel(), m, op(N.Rem{b, m}), op(N.Rem{op(N.One{}), m}), zst(~T, ~test, e))
541
542def pow_mod_ok(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +b: T, +e: T, +m: T, zero: Bool) -> Result<&2, &2, N.NumError, T>:
543  match zero:
544    case True{}:
545      Fail{N.DivByZero{}}
546    case False{}:
547      Done{pow_mod_pick(~T, ~op, ~test, b, e, m, test(N.Mont{m}))}
548
549# Python's pow(b, e, m) for m > 0
550def pow_mod(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +b: T, +e: T, +m: T) -> Result<&2, &2, N.NumError, T>:
551  pow_mod_ok(~T, ~op, ~test, b, e, m, test(N.IsZero{m}))
552
553# s0 - x mod m for s0, x < m, without leaving [0, m)
554def submod(~T: Data, ~op: N.Op<T> -> T, +m: T, +s0: T, +x: T, below: Bool) -> T:
555  match below:
556    case False{}:
557      op(N.Sub{s0, x})
558    case True{}:
559      op(N.Add{s0, op(N.Sub{m, x})})
560
561def inv_sub(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +m: T, +s0: T, +x: T) -> T:
562  submod(~T, ~op, m, s0, x, test(N.Lt{s0, x}))
563
564def inv_step(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +m: T, +q: T, +s0: T, +s1: T) -> T:
565  inv_sub(~T, ~op, ~test, m, s0, op(N.MulMod{op(N.Rem{q, m}), s1, m}))
566
567# extended Euclid on (m, a mod m), keeping the coefficient of a mod m;
568# st = (r1, r1 == 0)
569def inv_go(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, fuel: Nat, +m: T, +r0: T, +s0: T, +s1: T, st: T & Bool) -> T & T:
570  match fuel st:
571    case 0n _:
572      (r0, s0)
573    case 1n+f Tuple{r1, True{}}:
574      (r0, s0)
575    case 1n+f Tuple{+r1, False{}}:
576      inv_go(~T, ~op, ~test, f, m, r1, s1, inv_step(~T, ~op, ~test, m, op(N.Quot{r0, r1}), s0, s1), zst(~T, ~test, op(N.Rem{r0, r1})))
577
578def inv_one(~T: Data, ~op: N.Op<T> -> T, +m: T, +s: T, one: Bool) -> Result<&2, &2, N.NumError, T>:
579  match one:
580    case True{}:
581      Done{op(N.Rem{s, m})}
582    case False{}:
583      Fail{N.NoInverse{}}
584
585def inv_fin(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +m: T, r: T & T) -> Result<&2, &2, N.NumError, T>:
586  (+g, +s) = r
587  inv_one(~T, ~op, m, s, Bool.and(Bool.not(test(N.Lt{g, op(N.One{})})), Bool.not(test(N.Lt{op(N.One{}), g}))))
588
589def inv_ok(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +m: T, zero: Bool) -> Result<&2, &2, N.NumError, T>:
590  match zero:
591    case True{}:
592      Fail{N.DivByZero{}}
593    case False{}:
594      inv_fin(~T, ~op, ~test, m, inv_go(~T, ~op, ~test, fuel(), m, m, op(N.ZeroOp{}), op(N.One{}), zst(~T, ~test, op(N.Rem{a, m}))))
595
596# pow(a, -1, m): the x < m with a * x == 1 (mod m)
597def mod_inverse(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +m: T) -> Result<&2, &2, N.NumError, T>:
598  inv_ok(~T, ~op, ~test, a, m, test(N.IsZero{m}))
599
600# ---- division ----
601
602def divmod_ok(~T: Data, ~op: N.Op<T> -> T, +a: T, +b: T, zero: Bool) -> Result<&2, &2, N.NumError, QuotRem<T>>:
603  match zero:
604    case True{}:
605      Fail{N.DivByZero{}}
606    case False{}:
607      Done{TQR{op(N.Quot{a, b}), op(N.Rem{a, b})}}
608
609def divmod(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Result<&2, &2, N.NumError, QuotRem<T>>:
610  divmod_ok(~T, ~op, a, b, test(N.IsZero{b}))