1	import Base
2	import ./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	
26	type QuotRem<-T: Data> is Data:
27	  TQR{quot: T, rem: T}
28	
29	def 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	
36	def 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)
44	def fuel() -> Nat:
45	  140n
46	
47	# x when it fits, else None
48	def 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	
55	def 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	
58	def 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	
61	def lt(~T: Data, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Bool:
62	  test(N.Lt{a, b})
63	
64	def le(~T: Data, ~test: N.Test<T> -> Bool, +a: T, +b: T) -> Bool:
65	  Bool.not(test(N.Lt{b, a}))
66	
67	def pos(~T: Data, ~test: N.Test<T> -> Bool, +a: T) -> Bool:
68	  Bool.not(test(N.IsZero{a}))
69	
70	# ---- ordered-ring functions ----
71	
72	def 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.
80	def 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	
83	def 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	
86	def 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
94	def 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	
97	def abs(~T: Data, ~op: N.Op<T> -> T, ~test: N.Test<T> -> Bool, +x: T) -> T:
98	  op(N.Abs{x})
99	
100	def 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	
107	def 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)
115	def 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	
118	def 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	
127	def 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	
130	def 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	
139	def 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
143	def 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
151	def 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	
158	def 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	
165	def 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.
173	def 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	
182	def 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	
187	def one(~T: Data, ~op: N.Op<T> -> T) -> T:
188	  op(N.One{})
189	
190	def inc(~T: Data, ~op: N.Op<T> -> T, +x: T) -> T:
191	  op(N.Add{x, op(N.One{})})
192	
193	def dec(~T: Data, ~op: N.Op<T> -> T, +x: T) -> T:
194	  op(N.Sub{x, op(N.One{})})
195	
196	# (x, x == 0)
197	def 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))
201	def 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	
214	def 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)
218	def 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	
227	def 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)
231	def 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	
234	def 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	
241	def 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)
246	def 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	
257	def 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	
260	def 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	
267	def 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	
270	def 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	
279	def 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	
286	def 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	
293	def 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
302	def 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	
305	def 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
313	def 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	
316	def 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	
323	def 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	
326	def 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	
333	def 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	
336	def 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	
345	def 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
349	def 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	
358	def 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)
363	def 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
369	def 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	
376	def 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
380	def 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
385	def 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	
396	def 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
400	def 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	
409	def 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
417	def 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	
426	def 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	
433	def 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)
439	def 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
445	def 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	
456	def 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
460	def 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	
471	def 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	
478	def 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)
484	def 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	
487	def 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)
491	def 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	
502	def 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	
505	def 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	
512	def 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	
517	def 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
526	def 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	
535	def 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	
542	def 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
550	def 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)
554	def 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	
561	def 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	
564	def 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)
569	def 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	
578	def 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	
585	def 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	
589	def 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)
597	def 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	
602	def 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	
609	def 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}))