1 import Base 2 import ../lib/common.bend as C 3 import ../../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) 52 def 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 60 def 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) 68 def 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 80 def 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 87 def 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) 95 def 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 99 def 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 109 def 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] 121 def 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 134 def Gcd.divides_left(+a: Nat, +b: Nat) -> Type: 135 dvd(M.gcd(a, b), a) 136 137 def Gcd.divides_right(+a: Nat, +b: Nat) -> Type: 138 dvd(M.gcd(a, b), b) 139 140 # every common divisor divides the gcd 141 def 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 144 def 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 147 def Lcm.divides_left(+a: Nat, +b: Nat) -> Type: 148 dvd(a, M.lcm(a, b)) 149 150 def Lcm.divides_right(+a: Nat, +b: Nat) -> Type: 151 dvd(b, M.lcm(a, b)) 152 153 # the lcm divides every common multiple 154 def 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 157 def 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 162 def GcdAll.divides(+xs: List<&2, Nat>, +i: Nat) -> Type: 163 dvd(M.gcd_all(xs), nth0(xs, i, 0n)) 164 165 def 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 168 def LcmAll.divides(+xs: List<&2, Nat>, +i: Nat) -> Type: 169 dvd(nth0(xs, i, 1n), M.lcm_all(xs)) 170 171 def 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 174 def Sum.value(+xs: List<&2, Nat>) -> Type: 175 {M.sum(xs) == lsum(xs) : Nat} 176 177 def Prod.value(+xs: List<&2, Nat>) -> Type: 178 {M.prod(xs) == lprod(xs) : Nat} 179 180 # ---- factorial, perm, comb ---- 181 182 def Factorial.value(+n: Nat) -> Type: 183 {M.factorial(n) == factorial(n) : Nat} 184 185 def Perm.value(+n: Nat, +k: Nat) -> Type: 186 {M.perm(n, k) == desc(n, k) : Nat} 187 188 def Perm.self(+n: Nat) -> Type: 189 {M.perm(n, n) == factorial(n) : Nat} 190 191 def 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 195 def 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 198 def 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 202 def 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 207 def Isqrt.le(+n: Nat) -> Type: 208 {Nat.is_le(Nat.pow(M.isqrt(n), 2n), n) == True{} : Bool} 209 210 def 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 214 def 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 217 def 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 220 def 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 224 def 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 228 def 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 231 def 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 234 def 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 238 def Ilog.zero(+b: Nat) -> Type: 239 {M.ilog(0n, b) == Fail{M.Domain{}} : Result<&2, &2, M.MathError, Nat>} 240 241 def 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 244 def BitLength.lt(+n: Nat) -> Type: 245 {Nat.is_lt(n, C.pow2(M.bit_length(n))) == True{} : Bool} 246 247 def 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 252 def 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 256 def 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 259 def 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 262 def 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 266 def 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 269 def 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 274 def 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 277 def 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 280 def 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 284 def DivMod.zero_divisor(+a: Nat) -> Type: 285 {M.divmod(a, 0n) == Fail{M.ZeroDivision{}} : Result<&2, &2, M.MathError, M.QuotRem>} 286 287 def 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 290 def 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 293 def Clamp.le(+x: Nat, +lo: Nat, +hi: Nat) -> Type: 294 {Nat.is_le(Nat.min(Nat.max(x, lo), hi), hi) == True{} : Bool} 295 296 def 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 299 def 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>}