1	# Arbitrary-precision signed integers: add, sub, mul, divmod, pow, modpow, gcd and decimal text.
2	import Base
3	
4	# A magnitude is a list of base-10^4 limbs in U32, least significant first,
5	# with no zero limb on top. Base 10^4 keeps every limb product and its carries
6	# below 2^32, and makes decimal show and read linear. Zero is never negative.
7	type BigInt is Data:
8	  BigInt{neg: Bool, mag: List<&2, U32>}
9	
10	# A match cannot scrutinize a computed Bool, so each branch is a closure.
11	def Mag.if(-A: Type, c: Bool, t: Unit -> A, e: Unit -> A) -> A:
12	  match c:
13	    case True{}:
14	      t(Unit{})
15	    case False{}:
16	      e(Unit{})
17	
18	def Mag.lone(z: Bool, x: U32) -> List<&2, U32>:
19	  match z:
20	    case True{}:
21	      Nil{}
22	    case False{}:
23	      Con{x, Nil{}}
24	
25	# x + r * 10^4, dropping a zero limb on top.
26	def Mag.cons(+x: U32, r: List<&2, U32>) -> List<&2, U32>:
27	  match r:
28	    case Nil{}:
29	      Mag.lone(U32.is_zero(x), x)
30	    case s:
31	      Con{x, s}
32	
33	def Mag.norm(xs: List<&2, U32>) -> List<&2, U32>:
34	  match xs:
35	    case Nil{}:
36	      Nil{}
37	    case Con{x, t}:
38	      Mag.cons(x, Mag.norm(t))
39	
40	def Mag.len(xs: List<&2, U32>) -> Nat:
41	  List.length(&2, U32, xs)
42	
43	def Mag.at(xs: List<&2, U32>, i: Nat) -> U32:
44	  Maybe.default(&2, U32, List.get(&2, U32, xs, i), 0)
45	
46	def Mag.top(xs: List<&2, U32>) -> U32:
47	  Maybe.default(&2, U32, List.last(&2, U32, xs), 0)
48	
49	def Mag.odd(xs: List<&2, U32>) -> Bool:
50	  match xs:
51	    case Nil{}:
52	      False{}
53	    case Con{x, t}:
54	      U32.is_eq(U32.and(x, 1), 1)
55	
56	# The limbs are 0 or 1 in z: z says the carry c is 0, so a stays as is.
57	def Mag.inc(a: List<&2, U32>, z: Bool, +c: U32) -> List<&2, U32>:
58	  match a z:
59	    case s True{}:
60	      s
61	    case Nil{} False{}:
62	      Con{c, Nil{}}
63	    case Con{x, xs} False{}:
64	      +t = (x + c : U32)
65	      +k = (t / 10000 : U32)
66	      Con{(t % 10000 : U32), Mag.inc(xs, U32.is_zero(k), k)}
67	
68	def Mag.add(a: List<&2, U32>, b: List<&2, U32>, +c: U32) -> List<&2, U32>:
69	  match a b:
70	    case Con{x, xs} Con{y, ys}:
71	      +t = (x + y + c : U32)
72	      Con{(t % 10000 : U32), Mag.add(xs, ys, (t / 10000 : U32))}
73	    case Nil{} s:
74	      Mag.inc(s, U32.is_zero(c), c)
75	    case s Nil{}:
76	      Mag.inc(s, U32.is_zero(c), c)
77	
78	# a - w for the borrow w; z says w is 0. Leaves zero limbs on top.
79	def Mag.dec(a: List<&2, U32>, z: Bool, +w: U32) -> List<&2, U32>:
80	  match a z:
81	    case s True{}:
82	      s
83	    case Nil{} False{}:
84	      Nil{}
85	    case Con{x, xs} False{}:
86	      +t = (x + 10000 - w : U32)
87	      +k = (1 - t / 10000 : U32)
88	      Con{(t % 10000 : U32), Mag.dec(xs, U32.is_zero(k), k)}
89	
90	# a - b - w for a >= b + w. Leaves zero limbs on top: callers norm.
91	def Mag.sub(a: List<&2, U32>, b: List<&2, U32>, +w: U32) -> List<&2, U32>:
92	  match a b:
93	    case Con{x, xs} Con{y, ys}:
94	      +t = (x + 10000 - y - w : U32)
95	      Con{(t % 10000 : U32), Mag.sub(xs, ys, (1 - t / 10000 : U32))}
96	    case s Nil{}:
97	      Mag.dec(s, U32.is_zero(w), w)
98	    case Nil{} s:
99	      Nil{}
100	
101	def Mag.then(hi: Cmp, lo: Cmp) -> Cmp:
102	  match hi:
103	    case EQ{}:
104	      lo
105	    case LT{}:
106	      LT{}
107	    case GT{}:
108	      GT{}
109	
110	# Both canonical, so the longer one is larger.
111	def Mag.cmp(a: List<&2, U32>, b: List<&2, U32>) -> Cmp:
112	  match a b:
113	    case Con{x, xs} Con{y, ys}:
114	      Mag.then(Mag.cmp(xs, ys), U32.cmp(x, y))
115	    case Nil{} Nil{}:
116	      EQ{}
117	    case Nil{} Con{y, ys}:
118	      LT{}
119	    case Con{x, xs} Nil{}:
120	      GT{}
121	
122	# a * m + c for m, c below 10^4.
123	def Mag.mul1(a: List<&2, U32>, +m: U32, +c: U32) -> List<&2, U32>:
124	  match a:
125	    case Nil{}:
126	      Mag.lone(U32.is_zero(c), c)
127	    case Con{x, xs}:
128	      +t = (x * m + c : U32)
129	      Con{(t % 10000 : U32), Mag.mul1(xs, m, (t / 10000 : U32))}
130	
131	def Mag.scale(+a: List<&2, U32>, +m: U32) -> List<&2, U32>:
132	  Mag.if(List<&2, U32>, U32.is_zero(m), u => Nil{}, u => Mag.mul1(a, m, 0))
133	
134	def Mag.school.put(p: List<&2, U32>, r: List<&2, U32>) -> List<&2, U32>:
135	  match p:
136	    case Nil{}:
137	      Mag.cons(0, r)
138	    case Con{h, t}:
139	      Mag.cons(h, Mag.add(t, r, 0))
140	
141	def Mag.school(a: List<&2, U32>, +b: List<&2, U32>) -> List<&2, U32>:
142	  match a:
143	    case Nil{}:
144	      Nil{}
145	    case Con{x, t}:
146	      Mag.school.put(Mag.scale(b, x), Mag.school(t, b))
147	
148	def Mag.zeros(n: Nat, xs: List<&2, U32>) -> List<&2, U32>:
149	  match n:
150	    case 0n:
151	      xs
152	    case 1n+p:
153	      Con{0, Mag.zeros(p, xs)}
154	
155	# xs * 10^(4 n).
156	def Mag.shift(xs: List<&2, U32>, n: Nat) -> List<&2, U32>:
157	  match xs:
158	    case Nil{}:
159	      Nil{}
160	    case s:
161	      Mag.zeros(n, s)
162	
163	def Mag.lo(xs: List<&2, U32>, k: Nat) -> List<&2, U32>:
164	  Mag.norm(List.take(&2, U32, xs, k))
165	
166	# z0 + (s - z0 - z2) B^k + z2 B^2k, with B = 10^(4 k).
167	def Mag.kara(+k: Nat, +z0: List<&2, U32>, +z2: List<&2, U32>,
168	  s: List<&2, U32>) -> List<&2, U32>:
169	  +mid = Mag.norm(Mag.sub(Mag.norm(Mag.sub(s, z0, 0)), z2, 0))
170	  Mag.add(z0, Mag.shift(Mag.add(mid, Mag.shift(z2, k), 0), k), 0)
171	
172	# Karatsuba while both sides have at least 32 limbs, schoolbook below. f caps
173	# the depth; 32 halvings outrun any list that fits in memory.
174	def Mag.kmul(+f: Nat, +a: List<&2, U32>, +b: List<&2, U32>) -> List<&2, U32>:
175	  match f:
176	    case 0n:
177	      Mag.school(a, b)
178	    case 1n+p:
179	      +k = Nat.div(Nat.min(Mag.len(a), Mag.len(b)), 2n)
180	      +a0 = Mag.lo(a, k)
181	      +a1 = List.drop(&2, U32, a, k)
182	      +b0 = Mag.lo(b, k)
183	      +b1 = List.drop(&2, U32, b, k)
184	      Mag.if(List<&2, U32>, Nat.is_lt(k, 16n),
185	        u => Mag.school(a, b),
186	        u => Mag.kara(k, Mag.kmul(p, a0, b0), Mag.kmul(p, a1, b1),
187	          Mag.kmul(p, Mag.add(a0, a1, 0), Mag.add(b0, b1, 0))))
188	
189	def Mag.mul(a: List<&2, U32>, b: List<&2, U32>) -> List<&2, U32>:
190	  Mag.kmul(32n, a, b)
191	
192	# Short division of a big-endian list by d below 10^4 with remainder rem:
193	# the quotient comes out little-endian, with zero limbs on top.
194	def Mag.dsmall(xs: List<&2, U32>, +d: U32, +rem: U32, acc: List<&2, U32>) ->
195	  List<&2, U32>:
196	  match xs:
197	    case Nil{}:
198	      acc
199	    case Con{x, t}:
200	      +v = (rem * 10000 + x : U32)
201	      Mag.dsmall(t, d, (v % d : U32), Con{(v / d : U32), acc})
202	
203	def Mag.half(+xs: List<&2, U32>) -> List<&2, U32>:
204	  Mag.norm(Mag.dsmall(List.reverse(&2, U32, xs), 2, 0, Nil{}))
205	
206	# Knuth's Algorithm D step D6: the estimate d is at most 2 over the true
207	# quotient digit, so f = 3 steps of p - bn always settle it.
208	def Mag.fix(+f: Nat, +d: U32, +p: List<&2, U32>, +r: List<&2, U32>,
209	  +bn: List<&2, U32>, +q: List<&2, U32>) -> List<&2, U32> & List<&2, U32>:
210	  match f:
211	    case 0n:
212	      (Con{d, q}, Mag.norm(Mag.sub(r, p, 0)))
213	    case 1n+g:
214	      Mag.if(List<&2, U32> & List<&2, U32>, Cmp.is_gt(Mag.cmp(p, r)),
215	        u => Mag.fix(g, (d - 1 : U32), Mag.norm(Mag.sub(p, bn, 0)), r, bn, q),
216	        u => (Con{d, q}, Mag.norm(Mag.sub(r, p, 0))))
217	
218	# One quotient digit: bring limb x down into r, estimate from r's top two
219	# limbs over bn's top limb vt, and correct.
220	def Mag.step(x: U32, +bn: List<&2, U32>, +n: Nat, +n1: Nat, +vt: U32,
221	  qr: List<&2, U32> & List<&2, U32>) -> List<&2, U32> & List<&2, U32>:
222	  (q, r) = qr
223	  +r2 = Mag.cons(x, r)
224	  +d = U32.min(((Mag.at(r2, n) * 10000 + Mag.at(r2, n1)) / vt : U32), 9999)
225	  Mag.fix(3n, d, Mag.scale(bn, d), r2, bn, q)
226	
227	# Long division over the dividend's limbs, most significant first. Quotient
228	# digits are consed, so the quotient comes out little-endian.
229	def Mag.long(xs: List<&2, U32>, +bn: List<&2, U32>, +n: Nat, +n1: Nat,
230	  +vt: U32, qr: List<&2, U32> & List<&2, U32>) ->
231	  List<&2, U32> & List<&2, U32>:
232	  match xs:
233	    case Nil{}:
234	      qr
235	    case Con{x, t}:
236	      Mag.long(t, bn, n, n1, vt, Mag.step(x, bn, n, n1, vt, qr))
237	
238	def Mag.divmod.fin(+f: U32, qr: List<&2, U32> & List<&2, U32>) ->
239	  List<&2, U32> & List<&2, U32>:
240	  (q, r) = qr
241	  (Mag.norm(q), Mag.norm(Mag.dsmall(List.reverse(&2, U32, r), f, 0, Nil{})))
242	
243	# a / b and a % b for b above 0. Both scale by f so that b's top limb is at
244	# least 10^4 / 2 (Knuth D1), and the remainder scales back.
245	def Mag.divmod(a: List<&2, U32>, +b: List<&2, U32>) ->
246	  List<&2, U32> & List<&2, U32>:
247	  +f = (10000 / (Mag.top(b) + 1) : U32)
248	  +bn = Mag.mul1(b, f, 0)
249	  +n = Mag.len(bn)
250	  Mag.divmod.fin(f, Mag.long(List.reverse(&2, U32, Mag.mul1(a, f, 0)), bn, n,
251	    Nat.sub(n, 1n), Mag.top(bn),
252	    {(Nil{}, Nil{}) : List<&2, U32> & List<&2, U32>}))
253	
254	def Mag.snd(qr: List<&2, U32> & List<&2, U32>) -> List<&2, U32>:
255	  (q, r) = qr
256	  r
257	
258	def Mag.rem(a: List<&2, U32>, b: List<&2, U32>) -> List<&2, U32>:
259	  Mag.snd(Mag.divmod(a, b))
260	
261	# b^e * acc by squaring; e <= f, and e halves as f counts down.
262	def Mag.pow(+f: Nat, +b: List<&2, U32>, +e: Nat, +acc: List<&2, U32>) ->
263	  List<&2, U32>:
264	  match f:
265	    case 0n:
266	      acc
267	    case 1n+g:
268	      +h = Nat.div(e, 2n)
269	      +acc2 = Mag.if(List<&2, U32>, Nat.is_eq(Nat.mod(e, 2n), 1n),
270	        u => Mag.mul(acc, b), u => acc)
271	      Mag.if(List<&2, U32>, Nat.is_eq(h, 0n), u => acc2,
272	        u => Mag.pow(g, Mag.mul(b, b), h, acc2))
273	
274	# b^e * acc mod m, with b and acc below m. f counts e's bits: 14 per limb,
275	# since 10^4 < 2^14.
276	def Mag.modpow(+f: Nat, +b: List<&2, U32>, +e: List<&2, U32>,
277	  +acc: List<&2, U32>, +m: List<&2, U32>) -> List<&2, U32>:
278	  match f:
279	    case 0n:
280	      acc
281	    case 1n+g:
282	      +h = Mag.half(e)
283	      +acc2 = Mag.if(List<&2, U32>, Mag.odd(e),
284	        u => Mag.rem(Mag.mul(acc, b), m), u => acc)
285	      Mag.if(List<&2, U32>, List.is_empty(&2, U32, h), u => acc2,
286	        u => Mag.modpow(g, Mag.rem(Mag.mul(b, b), m), h, acc2, m))
287	
288	# Euclid. Lame: under 5 steps per decimal digit of the smaller, so f from
289	# 20 per limb of both never runs out.
290	def Mag.gcd(+f: Nat, +a: List<&2, U32>, +b: List<&2, U32>) -> List<&2, U32>:
291	  match f:
292	    case 0n:
293	      a
294	    case 1n+g:
295	      Mag.if(List<&2, U32>, List.is_empty(&2, U32, b), u => a,
296	        u => Mag.gcd(g, b, Mag.rem(a, b)))
297	
298	def Mag.dig(d: U32) -> Char:
299	  Chr{(48 + d : U32)}
300	
301	def Mag.d4(+x: U32, acc: String) -> String:
302	  SCon{Mag.dig((x / 1000 : U32)), SCon{Mag.dig(((x / 100) % 10 : U32)),
303	    SCon{Mag.dig(((x / 10) % 10 : U32)), SCon{Mag.dig((x % 10 : U32)), acc}}}}
304	
305	# The top limb unpadded, every lower limb as four digits.
306	def Mag.show(xs: List<&2, U32>, acc: String) -> String:
307	  match xs:
308	    case Nil{}:
309	      acc
310	    case Con{x, Nil{}}:
311	      U32.show(x) ++ acc
312	    case Con{x, t}:
313	      Mag.show(t, Mag.d4(x, acc))
314	
315	# The digit values of s, last first, and whether every char was a digit.
316	def Mag.digits(s: String, ok: Bool, acc: List<&2, U32>) ->
317	  Bool & List<&2, U32>:
318	  match s:
319	    case SNil{}:
320	      (ok, acc)
321	    case SCon{Chr{c}, t}:
322	      +d = (c - 48 : U32)
323	      Mag.digits(t, Bool.and(ok, U32.is_lt(d, 10)), Con{d, acc})
324	
325	# Digits, least significant first, four to a limb.
326	def Mag.group(ds: List<&2, U32>) -> List<&2, U32>:
327	  match ds:
328	    case Con{a, Con{b, Con{c, Con{d, t}}}}:
329	      Con{(a + 10 * b + 100 * c + 1000 * d : U32), Mag.group(t)}
330	    case Con{a, Con{b, Con{c, Nil{}}}}:
331	      Con{(a + 10 * b + 100 * c : U32), Nil{}}
332	    case Con{a, Con{b, Nil{}}}:
333	      Con{(a + 10 * b : U32), Nil{}}
334	    case Con{a, Nil{}}:
335	      Con{a, Nil{}}
336	    case Nil{}:
337	      Nil{}
338	
339	def BigInt.mk(neg: Bool, m: List<&2, U32>) -> BigInt:
340	  match m:
341	    case Nil{}:
342	      BigInt{False{}, Nil{}}
343	    case s:
344	      BigInt{neg, s}
345	
346	def BigInt.zero() -> BigInt:
347	  BigInt{False{}, Nil{}}
348	
349	def BigInt.one() -> BigInt:
350	  BigInt{False{}, Con{1, Nil{}}}
351	
352	def BigInt.from_u32(+x: U32) -> BigInt:
353	  BigInt{False{}, Mag.norm(Con{(x % 10000 : U32),
354	    Con{((x / 10000) % 10000 : U32), Con{(x / 100000000 : U32), Nil{}}}})}
355	
356	def BigInt.neg(x: BigInt) -> BigInt:
357	  match x:
358	    case BigInt{s, m}:
359	      BigInt.mk(Bool.not(s), m)
360	
361	def BigInt.abs(x: BigInt) -> BigInt:
362	  match x:
363	    case BigInt{s, m}:
364	      BigInt{False{}, m}
365	
366	def BigInt.is_zero(x: BigInt) -> Bool:
367	  match x:
368	    case BigInt{s, m}:
369	      List.is_empty(&2, U32, m)
370	
371	def BigInt.cmp(a: BigInt, b: BigInt) -> Cmp:
372	  match a b:
373	    case BigInt{False{}, x} BigInt{False{}, y}:
374	      Mag.cmp(x, y)
375	    case BigInt{True{}, x} BigInt{True{}, y}:
376	      Mag.cmp(y, x)
377	    case BigInt{False{}, x} BigInt{True{}, y}:
378	      GT{}
379	    case BigInt{True{}, x} BigInt{False{}, y}:
380	      LT{}
381	
382	# Signs differ: the larger magnitude keeps its sign.
383	def BigInt.add.mixed(c: Cmp, +sa: Bool, x: List<&2, U32>, y: List<&2, U32>) ->
384	  BigInt:
385	  match c:
386	    case LT{}:
387	      BigInt.mk(Bool.not(sa), Mag.norm(Mag.sub(y, x, 0)))
388	    case EQ{}:
389	      BigInt.zero()
390	    case GT{}:
391	      BigInt.mk(sa, Mag.norm(Mag.sub(x, y, 0)))
392	
393	def BigInt.add.go(diff: Bool, sa: Bool, +x: List<&2, U32>,
394	  +y: List<&2, U32>) -> BigInt:
395	  match diff:
396	    case True{}:
397	      BigInt.add.mixed(Mag.cmp(x, y), sa, x, y)
398	    case False{}:
399	      BigInt{sa, Mag.add(x, y, 0)}
400	
401	def BigInt.add(a: BigInt, b: BigInt) -> BigInt:
402	  match a b:
403	    case BigInt{+sa, x} BigInt{sb, y}:
404	      BigInt.add.go(Bool.xor(sa, sb), sa, x, y)
405	
406	def BigInt.sub(a: BigInt, b: BigInt) -> BigInt:
407	  BigInt.add(a, BigInt.neg(b))
408	
409	def BigInt.mul(a: BigInt, b: BigInt) -> BigInt:
410	  match a b:
411	    case BigInt{sa, x} BigInt{sb, y}:
412	      BigInt.mk(Bool.xor(sa, sb), Mag.mul(x, y))
413	
414	def BigInt.divmod.fin(+sa: Bool, sb: Bool,
415	  qr: List<&2, U32> & List<&2, U32>) -> BigInt & BigInt:
416	  (q, r) = qr
417	  (BigInt.mk(Bool.xor(sa, sb), q), BigInt.mk(sa, r))
418	
419	# Truncates toward zero: the remainder takes the dividend's sign.
420	def BigInt.divmod(a: BigInt, b: BigInt) -> Maybe<&1, BigInt & BigInt>:
421	  match a b:
422	    case BigInt{sa, x} BigInt{sb, Nil{}}:
423	      None{}
424	    case BigInt{sa, x} BigInt{sb, y}:
425	      Some{BigInt.divmod.fin(sa, sb, Mag.divmod(x, y))}
426	
427	def BigInt.pow(b: BigInt, +e: Nat) -> BigInt:
428	  match b:
429	    case BigInt{s, x}:
430	      BigInt.mk(Bool.and(s, Nat.is_eq(Nat.mod(e, 2n), 1n)),
431	        Mag.pow(e, x, e, Con{1, Nil{}}))
432	
433	# The base's residue in [0, m).
434	def Mag.base(neg: Bool, +r: List<&2, U32>, +m: List<&2, U32>) ->
435	  List<&2, U32>:
436	  match neg:
437	    case True{}:
438	      Mag.if(List<&2, U32>, List.is_empty(&2, U32, r), u => r,
439	        u => Mag.norm(Mag.sub(m, r, 0)))
440	    case False{}:
441	      r
442	
443	def BigInt.modpow.go(neg_e: Bool, sb: Bool, x: List<&2, U32>,
444	  +y: List<&2, U32>, +m: List<&2, U32>) -> Maybe<&2, BigInt>:
445	  match neg_e:
446	    case True{}:
447	      None{}
448	    case False{}:
449	      Some{BigInt{False{}, Mag.modpow(Nat.mul(Mag.len(y), 14n),
450	        Mag.base(sb, Mag.rem(x, m), m), y, Mag.rem(Con{1, Nil{}}, m), m)}}
451	
452	# b^e mod |m| in [0, |m|). None for m = 0 or e < 0.
453	def BigInt.modpow(b: BigInt, e: BigInt, m: BigInt) -> Maybe<&2, BigInt>:
454	  match b e m:
455	    case BigInt{sb, x} BigInt{se, y} BigInt{sm, Nil{}}:
456	      None{}
457	    case BigInt{sb, x} BigInt{se, y} BigInt{sm, z}:
458	      BigInt.modpow.go(se, sb, x, y, z)
459	
460	# Never negative; gcd(0, 0) is 0.
461	def BigInt.gcd(a: BigInt, b: BigInt) -> BigInt:
462	  match a b:
463	    case BigInt{sa, +x} BigInt{sb, +y}:
464	      BigInt{False{}, Mag.gcd(
465	        Nat.add(Nat.mul(Nat.add(Mag.len(x), Mag.len(y)), 20n), 2n), x, y)}
466	
467	def BigInt.show(x: BigInt) -> String:
468	  match x:
469	    case BigInt{s, Nil{}}:
470	      "0"
471	    case BigInt{True{}, m}:
472	      SCon{'-', Mag.show(m, SNil{})}
473	    case BigInt{False{}, m}:
474	      Mag.show(m, SNil{})
475	
476	def BigInt.read.fin(neg: Bool, r: Bool & List<&2, U32>) -> Maybe<&2, BigInt>:
477	  (ok, ds) = r
478	  match ok ds:
479	    case True{} Nil{}:
480	      None{}
481	    case True{} s:
482	      Some{BigInt.mk(neg, Mag.norm(Mag.group(s)))}
483	    case False{} s:
484	      None{}
485	
486	# An optional '-' and at least one ASCII digit; leading zeros are fine.
487	def BigInt.read(s: String) -> Maybe<&2, BigInt>:
488	  match s:
489	    case SCon{'-', t}:
490	      BigInt.read.fin(True{}, Mag.digits(t, True{}, Nil{}))
491	    case t:
492	      BigInt.read.fin(False{}, Mag.digits(t, True{}, Nil{}))