~/jendhub

spec/math/natural.bend

bend-collections-laws-math@1.0.0.0 · 300 lines · 58 definitions

1import Base
2import bend-collections-laws-math@1.0.0.0/spec/lib/common.bend as C
3import bend-collections-laws-math@1.0.0.0/src/math/natural.bend as M
4
5# Specification of src/math/natural.bend: Python's exact integer functions
6# (math.gcd, lcm, isqrt, factorial, perm, comb, prod, the three-argument
7# pow, pow(a, -1, m), divmod, int.bit_length) on the natural numbers.
8#
9# The reference definitions (factorial, desc, choose, lsum, lprod) are
10# structural recursions in the shape of Lean 4 Mathlib; the clauses mirror
11# Mathlib's statements, the Why3 gallery's (isqrt, fast exponentiation) and
12# HACL*'s (Lib.Exponentiation). Python semantics follow
13# python_style_math_stdlib_design.pdf: the error model (section 2.3: an
14# invalid argument is an error value, never a wrapped or partial result),
15# the exact integer functions (4.2), isqrt (8.5), factorial/comb/perm (8.6),
16# gcd, modular exponentiation and inverse (8.7), divmod and pow (3.1).
17# Divisibility "d divides a" is dvd(d, a): a quotient k with a == k d.
18#
19#   function     clauses                                         mirrors
20#   gcd          Gcd.divides_left, Gcd.divides_right,           Mathlib Nat.gcd_dvd_left/right,
21#                Gcd.greatest, Gcd.mul_left                      Nat.dvd_gcd, Nat.gcd_mul_left
22#   lcm          Lcm.divides_left, Lcm.divides_right,           Nat.dvd_lcm_left/right,
23#                Lcm.least, Lcm.gcd_mul_lcm                      Nat.lcm_dvd, Nat.gcd_mul_lcm
24#   gcd_all      GcdAll.divides, GcdAll.greatest                 Finset.gcd_dvd, Finset.dvd_gcd
25#   lcm_all      LcmAll.divides, LcmAll.least                    Finset.dvd_lcm, Finset.lcm_dvd
26#   sum, prod    Sum.value, Prod.value                           List.sum, List.prod
27#   factorial    Factorial.value                                 Nat.factorial_succ
28#   perm         Perm.value, Perm.self                           Nat.descFactorial, _self
29#   comb         Comb.value, Comb.pascal, Comb.factorials,       Nat.choose_succ_succ,
30#                Comb.zero                                       Nat.choose_mul_factorial_mul_factorial,
31#                                                                Nat.choose_eq_zero_of_lt
32#   isqrt        Isqrt.le, Isqrt.lt_succ                         Nat.sqrt_le', Nat.lt_succ_sqrt'
33#   iroot        Iroot.done, Iroot.le, Iroot.lt_succ,            Why3 isqrt generalised;
34#                Iroot.zero_degree                               PDF 4.2
35#   ilog         Ilog.done, Ilog.pow_le, Ilog.lt_pow_succ,       Nat.pow_log_le_self,
36#                Ilog.zero, Ilog.small_base                      Nat.lt_pow_succ_log_self; PDF 4.2
37#   bit_length   BitLength.lt, BitLength.le                      Nat.lt_size_self, Nat.size_le
38#   pow_mod      PowMod.value, PowMod.zero_modulus               Nat.pow_mod; HACL* Lib.Exponentiation
39#   mod_inverse  ModInverse.inverse, ModInverse.reduced,         Nat.gcdA, ZMod.mul_inv_of_unit,
40#                ModInverse.not_coprime, ModInverse.zero_modulus Nat.Coprime; PDF 8.7
41#   divmod       DivMod.value, DivMod.euclid, DivMod.rem_lt,     Nat.div_add_mod, Nat.mod_lt;
42#                DivMod.zero_divisor                             PDF 3.1
43#   clamp        Clamp.value, Clamp.ge, Clamp.le, Clamp.id,      PDF 10 (clamp is a gap Python
44#                Clamp.domain                                    leaves; lo > hi is a Domain error)
45#
46# Proved in proofs/math/natural/ (proof.bend, one definition per clause
47# under the clause's name).
48
49# ---- reference definitions ----
50
51# Mathlib/Data/Nat/Factorial/Basic (Nat.factorial_succ)
52def factorial(n: Nat) -> Nat:
53  match n:
54    case 0n:
55      1n
56    case 1n+ +p:
57      Nat.mul(1n+p, factorial(p))
58
59# n.descFactorial k: (n - k) * n.descFactorial k at k + 1
60def desc(+n: Nat, k: Nat) -> Nat:
61  match k:
62    case 0n:
63      1n
64    case 1n+ +j:
65      Nat.mul(Nat.sub(n, j), desc(n, j))
66
67# Mathlib/Data/Nat/Choose/Basic (Pascal's rule, Nat.choose_succ_succ)
68def choose(n: Nat, k: Nat) -> Nat:
69  match n k:
70    case 0n 0n:
71      1n
72    case 0n 1n+j:
73      0n
74    case 1n+m 0n:
75      1n
76    case 1n+ +m 1n+ +j:
77      Nat.add(choose(m, j), choose(m, 1n+j))
78
79# List.sum and List.prod as right folds
80def lsum(xs: List<&2, Nat>) -> Nat:
81  match xs:
82    case Nil{}:
83      0n
84    case Con{x, t}:
85      Nat.add(x, lsum(t))
86
87def lprod(xs: List<&2, Nat>) -> Nat:
88  match xs:
89    case Nil{}:
90      1n
91    case Con{x, t}:
92      Nat.mul(x, lprod(t))
93
94# d divides a: a quotient k with a == k d (Mathlib's Dvd on Nat)
95def dvd(+d: Nat, +a: Nat) -> Type:
96  Sigma<&1, &1, Nat, k => {a == Nat.mul(k, d) : Nat}>
97
98# xs[i], or dflt past the end
99def nth0(+xs: List<&2, Nat>, +i: Nat, +dflt: Nat) -> Nat:
100  match xs i:
101    case Nil{} i0:
102      dflt
103    case Con{x, t} 0n:
104      x
105    case Con{x, t} 1n+j:
106      nth0(t, j, dflt)
107
108# d divides every element, with the quotients given: xs[i] == ks[i] d
109def dvd_with(+d: Nat, +xs: List<&2, Nat>, +ks: List<&2, Nat>) -> Type:
110  match xs ks:
111    case Nil{} Nil{}:
112      Unit
113    case Nil{} Con{k, kt}:
114      Unit
115    case Con{x, t} Nil{}:
116      Empty
117    case Con{+x, +t} Con{+k, +kt}:
118      {x == Nat.mul(k, d) : Nat} & dvd_with(d, t, kt)
119
120# m is a multiple of every element, with the quotients given: m == ks[i] xs[i]
121def mul_with(+m: Nat, +xs: List<&2, Nat>, +ks: List<&2, Nat>) -> Type:
122  match xs ks:
123    case Nil{} Nil{}:
124      Unit
125    case Nil{} Con{k, kt}:
126      Unit
127    case Con{x, t} Nil{}:
128      Empty
129    case Con{+x, +t} Con{+k, +kt}:
130      {m == Nat.mul(k, x) : Nat} & mul_with(m, t, kt)
131
132# ---- gcd, lcm ----
133
134def Gcd.divides_left(+a: Nat, +b: Nat) -> Type:
135  dvd(M.gcd(a, b), a)
136
137def Gcd.divides_right(+a: Nat, +b: Nat) -> Type:
138  dvd(M.gcd(a, b), b)
139
140# every common divisor divides the gcd
141def Gcd.greatest(+a: Nat, +b: Nat, +d: Nat, +ka: Nat, +kb: Nat, +ea: {a == Nat.mul(ka, d) : Nat}, +eb: {b == Nat.mul(kb, d) : Nat}) -> Type:
142  dvd(d, M.gcd(a, b))
143
144def Gcd.mul_left(+c: Nat, +a: Nat, +b: Nat) -> Type:
145  {M.gcd(Nat.mul(c, a), Nat.mul(c, b)) == Nat.mul(c, M.gcd(a, b)) : Nat}
146
147def Lcm.divides_left(+a: Nat, +b: Nat) -> Type:
148  dvd(a, M.lcm(a, b))
149
150def Lcm.divides_right(+a: Nat, +b: Nat) -> Type:
151  dvd(b, M.lcm(a, b))
152
153# the lcm divides every common multiple
154def Lcm.least(+a: Nat, +b: Nat, +m: Nat, +x: Nat, +y: Nat, +ex: {m == Nat.mul(x, a) : Nat}, +ey: {m == Nat.mul(y, b) : Nat}) -> Type:
155  dvd(M.lcm(a, b), m)
156
157def Lcm.gcd_mul_lcm(+a: Nat, +b: Nat) -> Type:
158  {Nat.mul(M.gcd(a, b), M.lcm(a, b)) == Nat.mul(a, b) : Nat}
159
160# ---- lists ----
161
162def GcdAll.divides(+xs: List<&2, Nat>, +i: Nat) -> Type:
163  dvd(M.gcd_all(xs), nth0(xs, i, 0n))
164
165def GcdAll.greatest(+d: Nat, +xs: List<&2, Nat>, +ks: List<&2, Nat>, h: dvd_with(d, xs, ks)) -> Type:
166  dvd(d, M.gcd_all(xs))
167
168def LcmAll.divides(+xs: List<&2, Nat>, +i: Nat) -> Type:
169  dvd(nth0(xs, i, 1n), M.lcm_all(xs))
170
171def LcmAll.least(+m: Nat, +xs: List<&2, Nat>, +ks: List<&2, Nat>, h: mul_with(m, xs, ks)) -> Type:
172  dvd(M.lcm_all(xs), m)
173
174def Sum.value(+xs: List<&2, Nat>) -> Type:
175  {M.sum(xs) == lsum(xs) : Nat}
176
177def Prod.value(+xs: List<&2, Nat>) -> Type:
178  {M.prod(xs) == lprod(xs) : Nat}
179
180# ---- factorial, perm, comb ----
181
182def Factorial.value(+n: Nat) -> Type:
183  {M.factorial(n) == factorial(n) : Nat}
184
185def Perm.value(+n: Nat, +k: Nat) -> Type:
186  {M.perm(n, k) == desc(n, k) : Nat}
187
188def Perm.self(+n: Nat) -> Type:
189  {M.perm(n, n) == factorial(n) : Nat}
190
191def Comb.value(+n: Nat, +k: Nat) -> Type:
192  {M.comb(n, k) == choose(n, k) : Nat}
193
194# Pascal's rule on the implementation itself
195def Comb.pascal(+m: Nat, +j: Nat) -> Type:
196  {M.comb(1n+m, 1n+j) == Nat.add(M.comb(m, j), M.comb(m, 1n+j)) : Nat}
197
198def Comb.factorials(+n: Nat, +k: Nat, +h: {Nat.is_le(k, n) == True{} : Bool}) -> Type:
199  {Nat.mul(M.comb(n, k), Nat.mul(factorial(k), factorial(Nat.sub(n, k)))) == factorial(n) : Nat}
200
201# PDF 8.6: comb(n, k) is 0 for k > n, not an error
202def Comb.zero(+n: Nat, +k: Nat, +h: {Nat.is_lt(n, k) == True{} : Bool}) -> Type:
203  {M.comb(n, k) == 0n : Nat}
204
205# ---- roots, logarithms, bits ----
206
207def Isqrt.le(+n: Nat) -> Type:
208  {Nat.is_le(Nat.pow(M.isqrt(n), 2n), n) == True{} : Bool}
209
210def Isqrt.lt_succ(+n: Nat) -> Type:
211  {Nat.is_lt(n, Nat.pow(1n+M.isqrt(n), 2n)) == True{} : Bool}
212
213# a positive degree always has a root
214def Iroot.done(+n: Nat, +kp: Nat) -> Type:
215  Sigma<&1, &1, Nat, r => {M.iroot(n, 1n+kp) == Done{r} : Result<&2, &2, M.MathError, Nat>}>
216
217def Iroot.le(+n: Nat, +kp: Nat, +r: Nat, +h: {M.iroot(n, 1n+kp) == Done{r} : Result<&2, &2, M.MathError, Nat>}) -> Type:
218  {Nat.is_le(Nat.pow(r, 1n+kp), n) == True{} : Bool}
219
220def Iroot.lt_succ(+n: Nat, +kp: Nat, +r: Nat, +h: {M.iroot(n, 1n+kp) == Done{r} : Result<&2, &2, M.MathError, Nat>}) -> Type:
221  {Nat.is_lt(n, Nat.pow(1n+r, 1n+kp)) == True{} : Bool}
222
223# PDF 2.3: the zeroth root is a domain error
224def Iroot.zero_degree(+n: Nat) -> Type:
225  {M.iroot(n, 0n) == Fail{M.Domain{}} : Result<&2, &2, M.MathError, Nat>}
226
227# a positive argument and a base of at least 2 always have a logarithm
228def Ilog.done(+np: Nat, +bq: Nat) -> Type:
229  Sigma<&1, &1, Nat, r => {M.ilog(1n+np, 2n+bq) == Done{r} : Result<&2, &2, M.MathError, Nat>}>
230
231def Ilog.pow_le(+np: Nat, +bq: Nat, +r: Nat, +h: {M.ilog(1n+np, 2n+bq) == Done{r} : Result<&2, &2, M.MathError, Nat>}) -> Type:
232  {Nat.is_le(Nat.pow(2n+bq, r), 1n+np) == True{} : Bool}
233
234def Ilog.lt_pow_succ(+np: Nat, +bq: Nat, +r: Nat, +h: {M.ilog(1n+np, 2n+bq) == Done{r} : Result<&2, &2, M.MathError, Nat>}) -> Type:
235  {Nat.is_lt(1n+np, Nat.pow(2n+bq, 1n+r)) == True{} : Bool}
236
237# PDF 2.3 / 4.5: log(0) is a domain error, as is a base below 2
238def Ilog.zero(+b: Nat) -> Type:
239  {M.ilog(0n, b) == Fail{M.Domain{}} : Result<&2, &2, M.MathError, Nat>}
240
241def Ilog.small_base(+n: Nat, +b: Nat, +hb: {Nat.is_lt(b, 2n) == True{} : Bool}) -> Type:
242  {M.ilog(n, b) == Fail{M.Domain{}} : Result<&2, &2, M.MathError, Nat>}
243
244def BitLength.lt(+n: Nat) -> Type:
245  {Nat.is_lt(n, C.pow2(M.bit_length(n))) == True{} : Bool}
246
247def BitLength.le(+np: Nat) -> Type:
248  {Nat.is_le(C.pow2(M.bit_length(1n+np)), Nat.double(1n+np)) == True{} : Bool}
249
250# ---- modular arithmetic ----
251
252def PowMod.value(+b: Nat, +e: Nat, +mp: Nat) -> Type:
253  {M.pow_mod(b, e, 1n+mp) == Done{Nat.mod(Nat.pow(b, e), 1n+mp)} : Result<&2, &2, M.MathError, Nat>}
254
255# PDF 8.7: pow(b, e, 0) raises ValueError; here ZeroDivision
256def PowMod.zero_modulus(+b: Nat, +e: Nat) -> Type:
257  {M.pow_mod(b, e, 0n) == Fail{M.ZeroDivision{}} : Result<&2, &2, M.MathError, Nat>}
258
259def ModInverse.inverse(+a: Nat, +mp: Nat, +x: Nat, +h: {M.mod_inverse(a, 1n+mp) == Done{x} : Result<&2, &2, M.MathError, Nat>}) -> Type:
260  {Nat.mod(Nat.mul(a, x), 1n+mp) == Nat.mod(1n, 1n+mp) : Nat}
261
262def ModInverse.reduced(+a: Nat, +mp: Nat, +x: Nat, +h: {M.mod_inverse(a, 1n+mp) == Done{x} : Result<&2, &2, M.MathError, Nat>}) -> Type:
263  {Nat.is_lt(x, 1n+mp) == True{} : Bool}
264
265# PDF 8.7: no inverse exactly when a and m share a factor g >= 2
266def ModInverse.not_coprime(+a: Nat, +mp: Nat, +h: {M.mod_inverse(a, 1n+mp) == Fail{M.NotInvertible{}} : Result<&2, &2, M.MathError, Nat>}) -> Type:
267  Sigma<&1, &1, Nat, g => {Nat.is_le(2n, g) == True{} : Bool} & (dvd(g, a) & dvd(g, 1n+mp))>
268
269def ModInverse.zero_modulus(+a: Nat) -> Type:
270  {M.mod_inverse(a, 0n) == Fail{M.ZeroDivision{}} : Result<&2, &2, M.MathError, Nat>}
271
272# ---- divmod and clamp ----
273
274def DivMod.value(+a: Nat, +bp: Nat) -> Type:
275  {M.divmod(a, 1n+bp) == Done{M.QR{Nat.div(a, 1n+bp), Nat.mod(a, 1n+bp)}} : Result<&2, &2, M.MathError, M.QuotRem>}
276
277def DivMod.euclid(+a: Nat, +bp: Nat) -> Type:
278  {a == Nat.add(Nat.mul(Nat.div(a, 1n+bp), 1n+bp), Nat.mod(a, 1n+bp)) : Nat}
279
280def DivMod.rem_lt(+a: Nat, +bp: Nat) -> Type:
281  {Nat.is_lt(Nat.mod(a, 1n+bp), 1n+bp) == True{} : Bool}
282
283# PDF 3.1: divmod(a, 0) raises ZeroDivisionError
284def DivMod.zero_divisor(+a: Nat) -> Type:
285  {M.divmod(a, 0n) == Fail{M.ZeroDivision{}} : Result<&2, &2, M.MathError, M.QuotRem>}
286
287def Clamp.value(+x: Nat, +lo: Nat, +hi: Nat, +h: {Nat.is_le(lo, hi) == True{} : Bool}) -> Type:
288  {M.clamp(x, lo, hi) == Done{Nat.min(Nat.max(x, lo), hi)} : Result<&2, &2, M.MathError, Nat>}
289
290def Clamp.ge(+x: Nat, +lo: Nat, +hi: Nat, +h: {Nat.is_le(lo, hi) == True{} : Bool}) -> Type:
291  {Nat.is_le(lo, Nat.min(Nat.max(x, lo), hi)) == True{} : Bool}
292
293def Clamp.le(+x: Nat, +lo: Nat, +hi: Nat) -> Type:
294  {Nat.is_le(Nat.min(Nat.max(x, lo), hi), hi) == True{} : Bool}
295
296def Clamp.id(+x: Nat, +lo: Nat, +hi: Nat, +h1: {Nat.is_le(lo, x) == True{} : Bool}, +h2: {Nat.is_le(x, hi) == True{} : Bool}) -> Type:
297  {Nat.min(Nat.max(x, lo), hi) == x : Nat}
298
299def Clamp.domain(+x: Nat, +lo: Nat, +hi: Nat, +h: {Nat.is_lt(hi, lo) == True{} : Bool}) -> Type:
300  {M.clamp(x, lo, hi) == Fail{M.Domain{}} : Result<&2, &2, M.MathError, Nat>}