bigint.bend

bend-kit-bignum@0.1.0.0 · 492 lines · 65 definitions

1# Arbitrary-precision signed integers: add, sub, mul, divmod, pow, modpow, gcd and decimal text.
2import 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.
7type 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.
11def 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
18def 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.
26def 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
33def 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
40def Mag.len(xs: List<&2, U32>) -> Nat:
41  List.length(&2, U32, xs)
42
43def Mag.at(xs: List<&2, U32>, i: Nat) -> U32:
44  Maybe.default(&2, U32, List.get(&2, U32, xs, i), 0)
45
46def Mag.top(xs: List<&2, U32>) -> U32:
47  Maybe.default(&2, U32, List.last(&2, U32, xs), 0)
48
49def 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.
57def 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
68def 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.
79def 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.
91def 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
101def 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.
111def 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.
123def 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
131def 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
134def 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
141def 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
148def 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).
156def 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
163def 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).
167def 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.
174def 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
189def 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.
194def 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
203def 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.
208def 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.
220def 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.
229def 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
238def 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.
245def 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
254def Mag.snd(qr: List<&2, U32> & List<&2, U32>) -> List<&2, U32>:
255  (q, r) = qr
256  r
257
258def 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.
262def 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.
276def 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.
290def 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
298def Mag.dig(d: U32) -> Char:
299  Chr{(48 + d : U32)}
300
301def 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.
306def 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.
316def 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.
326def 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
339def 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
346def BigInt.zero() -> BigInt:
347  BigInt{False{}, Nil{}}
348
349def BigInt.one() -> BigInt:
350  BigInt{False{}, Con{1, Nil{}}}
351
352def 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
356def BigInt.neg(x: BigInt) -> BigInt:
357  match x:
358    case BigInt{s, m}:
359      BigInt.mk(Bool.not(s), m)
360
361def BigInt.abs(x: BigInt) -> BigInt:
362  match x:
363    case BigInt{s, m}:
364      BigInt{False{}, m}
365
366def BigInt.is_zero(x: BigInt) -> Bool:
367  match x:
368    case BigInt{s, m}:
369      List.is_empty(&2, U32, m)
370
371def 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.
383def 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
393def 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
401def 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
406def BigInt.sub(a: BigInt, b: BigInt) -> BigInt:
407  BigInt.add(a, BigInt.neg(b))
408
409def 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
414def 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.
420def 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
427def 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).
434def 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
443def 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.
453def 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.
461def 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
467def 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
476def 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.
487def 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{}))