1	import Base
2	import ../../lib/nat.bend as N
3	import ../../lib/logic.bend as L
4	import ../../lib/lemmas/proofs/nat_algebra.bend as A
5	import ../../../src/math/natural.bend as M
6	import ./arith.bend as R
7	
8	# gcd is the greatest common divisor (Lean 4 Mathlib, Mathlib/Data/Nat/GCD:
9	# Nat.gcd_dvd_left, Nat.gcd_dvd_right, Nat.dvd_gcd; the Euclid step is
10	# Nat.gcd_rec). "d divides a" is a == k d with an explicit witness k, so
11	# every divisibility fact below names its quotient.
12	
13	# ---- Euclid's step keeps the common divisors ----
14	
15	# d | a and d | n give d | a mod n: a mod n == (ka - q kn) d, q = a / n
16	def dvd_mod(+bp: Nat, +a: Nat, +d: Nat, +ka: Nat, +kn: Nat, +ea: {a == Nat.mul(ka, d) : Nat}, +en: {1n+bp == Nat.mul(kn, d) : Nat}) -> {Nat.mod(a, 1n+bp) == Nat.mul(Nat.sub(ka, Nat.mul(Nat.div(a, 1n+bp), kn)), d) : Nat}:
17	  +q = Nat.div(a, 1n+bp)
18	  %Equal.sym(Nat, Nat.mul(Nat.sub(ka, Nat.mul(q, kn)), d), Nat.sub(Nat.mul(ka, d), Nat.mul(Nat.mul(q, kn), d)), R.mul_sub(ka, Nat.mul(q, kn), d)) : {Nat.mod(a, 1n+bp) == _ : Nat}
19	  %Equal.sym(Nat, Nat.mul(Nat.mul(q, kn), d), Nat.mul(q, Nat.mul(kn, d)), A.mul_assoc(q, kn, d)) : {Nat.mod(a, 1n+bp) == Nat.sub(Nat.mul(ka, d), _) : Nat}
20	  %Equal.sym(Nat, Nat.mul(kn, d), 1n+bp, Equal.sym(Nat, 1n+bp, Nat.mul(kn, d), en)) : {Nat.mod(a, 1n+bp) == Nat.sub(Nat.mul(ka, d), Nat.mul(q, _)) : Nat}
21	  %Equal.sym(Nat, Nat.mul(ka, d), a, Equal.sym(Nat, a, Nat.mul(ka, d), ea)) : {Nat.mod(a, 1n+bp) == Nat.sub(_, Nat.mul(q, 1n+bp)) : Nat}
22	  %Equal.sym(Nat, a, Nat.add(Nat.mul(q, 1n+bp), Nat.mod(a, 1n+bp)), R.dm_eq(bp, a)) : {Nat.mod(a, 1n+bp) == Nat.sub(_, Nat.mul(q, 1n+bp)) : Nat}
23	  Equal.sym(Nat, Nat.sub(Nat.add(Nat.mul(q, 1n+bp), Nat.mod(a, 1n+bp)), Nat.mul(q, 1n+bp)), Nat.mod(a, 1n+bp), N.add_sub_cancel(Nat.mul(q, 1n+bp), Nat.mod(a, 1n+bp)))
24	
25	# d | n and d | a mod n give d | a: a == (q kn + kr) d
26	def dvd_of_mod(+bp: Nat, +a: Nat, +d: Nat, +kn: Nat, +kr: Nat, +en: {1n+bp == Nat.mul(kn, d) : Nat}, +er: {Nat.mod(a, 1n+bp) == Nat.mul(kr, d) : Nat}) -> {a == Nat.mul(Nat.add(Nat.mul(Nat.div(a, 1n+bp), kn), kr), d) : Nat}:
27	  +q = Nat.div(a, 1n+bp)
28	  %Equal.sym(Nat, Nat.mul(Nat.add(Nat.mul(q, kn), kr), d), Nat.add(Nat.mul(Nat.mul(q, kn), d), Nat.mul(kr, d)), A.mul_add_right(Nat.mul(q, kn), kr, d)) : {a == _ : Nat}
29	  %Equal.sym(Nat, Nat.mul(Nat.mul(q, kn), d), Nat.mul(q, Nat.mul(kn, d)), A.mul_assoc(q, kn, d)) : {a == Nat.add(_, Nat.mul(kr, d)) : Nat}
30	  %Equal.sym(Nat, Nat.mul(kn, d), 1n+bp, Equal.sym(Nat, 1n+bp, Nat.mul(kn, d), en)) : {a == Nat.add(Nat.mul(q, _), Nat.mul(kr, d)) : Nat}
31	  %Equal.sym(Nat, Nat.mul(kr, d), Nat.mod(a, 1n+bp), Equal.sym(Nat, Nat.mod(a, 1n+bp), Nat.mul(kr, d), er)) : {a == Nat.add(Nat.mul(q, 1n+bp), _) : Nat}
32	  R.dm_eq(bp, a)
33	
34	# ---- gcd divides both arguments ----
35	
36	# the quotients a / g and b / g, by Euclid's recursion
37	type Wit is Data:
38	  W{l: Nat, r: Nat}
39	
40	def wl(w: Wit) -> Nat:
41	  match w:
42	    case W{l, r}:
43	      l
44	
45	def wr(w: Wit) -> Nat:
46	  match w:
47	    case W{l, r}:
48	      r
49	
50	# from n == l g and a mod n == r g: a == (q l + r) g and n == l g
51	def wstep(+q: Nat, w: Wit) -> Wit:
52	  match w:
53	    case W{+l, r}:
54	      W{Nat.add(Nat.mul(q, l), r), l}
55	
56	def ws(fuel: Nat, +a: Nat, +b: Nat) -> Wit:
57	  match fuel b:
58	    case 0n _:
59	      W{1n, 0n}
60	    case 1n+f 0n:
61	      W{1n, 0n}
62	    case 1n+f 1n+ +bp:
63	      wstep(Nat.div(a, 1n+bp), ws(f, 1n+bp, Nat.mod(a, 1n+bp)))
64	
65	def divides_both(+a: Nat, +b: Nat, +g: Nat, w: Wit) -> Type:
66	  {a == Nat.mul(wl(w), g) : Nat} & {b == Nat.mul(wr(w), g) : Nat}
67	
68	def step_both(+bp: Nat, +a: Nat, +g: Nat, w: Wit, ih: divides_both(1n+bp, Nat.mod(a, 1n+bp), g, w)) -> divides_both(a, 1n+bp, g, wstep(Nat.div(a, 1n+bp), w)):
69	  match w:
70	    case W{+l, +r}:
71	      (en, er) = ih
72	      +en2 = {en : {1n+bp == Nat.mul(l, g) : Nat}}
73	      (dvd_of_mod(bp, a, g, l, r, en2, er), en2)
74	
75	def gcd_both(fuel: Nat, +a: Nat, +b: Nat, +hb: {Nat.is_le(b, fuel) == True{} : Bool}) -> divides_both(a, b, M.gcd_go(fuel, a, b), ws(fuel, a, b)):
76	  match fuel b:
77	    case 0n 0n:
78	      (Equal.sym(Nat, Nat.add(a, 0n), a, A.add_zero(a)), {==})
79	    case 0n 1n+bp:
80	      Empty.absurd(divides_both(a, 1n+bp, M.gcd_go(0n, a, 1n+bp), ws(0n, a, 1n+bp)), L.false_true(hb))
81	    case 1n+f 0n:
82	      (Equal.sym(Nat, Nat.add(a, 0n), a, A.add_zero(a)), {==})
83	    case 1n+ +f 1n+ +bp:
84	      step_both(bp, a, M.gcd_go(f, 1n+bp, Nat.mod(a, 1n+bp)), ws(f, 1n+bp, Nat.mod(a, 1n+bp)), gcd_both(f, 1n+bp, Nat.mod(a, 1n+bp), N.le_trans(Nat.mod(a, 1n+bp), bp, f, N.lt_succ_le(Nat.mod(a, 1n+bp), bp, R.dm_lt(bp, a)), hb)))
85	
86	# ---- every common divisor divides gcd ----
87	
88	# the quotient g / d, for d | a (a == ka d) and d | b (b == kb d)
89	def dw(fuel: Nat, +a: Nat, +b: Nat, +ka: Nat, +kb: Nat) -> Nat:
90	  match fuel b:
91	    case 0n _:
92	      ka
93	    case 1n+f 0n:
94	      ka
95	    case 1n+f 1n+ +bp:
96	      dw(f, 1n+bp, Nat.mod(a, 1n+bp), kb, Nat.sub(ka, Nat.mul(Nat.div(a, 1n+bp), kb)))
97	
98	def gcd_greatest(fuel: Nat, +a: Nat, +b: Nat, +hb: {Nat.is_le(b, fuel) == True{} : Bool}, +d: Nat, +ka: Nat, +kb: Nat, +ea: {a == Nat.mul(ka, d) : Nat}, +eb: {b == Nat.mul(kb, d) : Nat}) -> {M.gcd_go(fuel, a, b) == Nat.mul(dw(fuel, a, b, ka, kb), d) : Nat}:
99	  match fuel b:
100	    case 0n 0n:
101	      ea
102	    case 0n 1n+bp:
103	      Empty.absurd({M.gcd_go(0n, a, 1n+bp) == Nat.mul(dw(0n, a, 1n+bp, ka, kb), d) : Nat}, L.false_true(hb))
104	    case 1n+f 0n:
105	      ea
106	    case 1n+ +f 1n+ +bp:
107	      gcd_greatest(f, 1n+bp, Nat.mod(a, 1n+bp), N.le_trans(Nat.mod(a, 1n+bp), bp, f, N.lt_succ_le(Nat.mod(a, 1n+bp), bp, R.dm_lt(bp, a)), hb), d, kb, Nat.sub(ka, Nat.mul(Nat.div(a, 1n+bp), kb)), eb, dvd_mod(bp, a, d, ka, kb, ea, eb))
108	
109	# ---- the theorems, on gcd(a, b) = gcd_go(b, a, b) ----
110	
111	def both_l(+a: Nat, +b: Nat, +g: Nat, w: Wit, p: divides_both(a, b, g, w)) -> {a == Nat.mul(wl(w), g) : Nat}:
112	  (el, er) = p
113	  el
114	
115	def both_r(+a: Nat, +b: Nat, +g: Nat, w: Wit, p: divides_both(a, b, g, w)) -> {b == Nat.mul(wr(w), g) : Nat}:
116	  (el, er) = p
117	  er
118	
119	def gcd_dvd_left(+a: Nat, +b: Nat) -> {a == Nat.mul(wl(ws(b, a, b)), M.gcd(a, b)) : Nat}:
120	  both_l(a, b, M.gcd(a, b), ws(b, a, b), gcd_both(b, a, b, N.le_refl(b)))
121	
122	def gcd_dvd_right(+a: Nat, +b: Nat) -> {b == Nat.mul(wr(ws(b, a, b)), M.gcd(a, b)) : Nat}:
123	  both_r(a, b, M.gcd(a, b), ws(b, a, b), gcd_both(b, a, b, N.le_refl(b)))
124	
125	def dvd_gcd(+a: Nat, +b: Nat, +d: Nat, +ka: Nat, +kb: Nat, +ea: {a == Nat.mul(ka, d) : Nat}, +eb: {b == Nat.mul(kb, d) : Nat}) -> {M.gcd(a, b) == Nat.mul(dw(b, a, b, ka, kb), d) : Nat}:
126	  gcd_greatest(b, a, b, N.le_refl(b), d, ka, kb, ea, eb)