Euclidean algorithm for ℕ

This file sets up a version of the Euclidean algorithm that only works with natural numbers. Given 0 < a, b, it computes the unique (w, x, y, z, d) such that the following identities hold:

• a = (w + x) d
• b = (y + z) d
• w * z = x * y + 1 d is then the gcd of a and b, and a' := a / d = w + x and b' := b / d = y + z are coprime.

This story is closely related to the structure of SL₂(ℕ) (as a free monoid on two generators) and the theory of continued fractions.

Main declarations

• XgcdType: Helper type in defining the gcd. Encapsulates (wp, x, y, zp, ap, bp). where wp zp, ap, bp are the variables getting changed through the algorithm.
• IsSpecial: States wp * zp = x * y + 1
• IsReduced: States ap = a ∧ bp = b

Notes

See Nat.Xgcd for a very similar algorithm allowing values in ℤ.

structure PNat.XgcdType : Type

• wp : ℕ

  wp is a variable which changes through the algorithm.

• x : ℕ

  x satisfies a / d = w + x at the final step.

• y : ℕ

  y satisfies b / d = z + y at the final step.

• zp : ℕ

  zp is a variable which changes through the algorithm.

• ap : ℕ

  ap is a variable which changes through the algorithm.

• bp : ℕ

  bp is a variable which changes through the algorithm.

A term of XgcdType is a system of six naturals. They should be thought of as representing the matrix [[w, x], [y, z]] = [[wp + 1, x], [y, zp + 1]] together with the vector [a, b] = [ap + 1, bp + 1].

instance PNat.instInhabitedXgcdType : Inhabited PNat.XgcdType

instance PNat.XgcdType.instSizeOfXgcdType : SizeOf PNat.XgcdType

instance PNat.XgcdType.instReprXgcdType : Repr PNat.XgcdType

The Repr instance converts terms to strings in a way that reflects the matrix/vector interpretation as above.

def PNat.XgcdType.mk' (w : ℕ+) (x : ℕ) (y : ℕ) (z : ℕ+) (a : ℕ+) (b : ℕ+) : PNat.XgcdType

Another mk using ℕ and ℕ+

def PNat.XgcdType.w (u : PNat.XgcdType) : ℕ+

w = wp + 1

def PNat.XgcdType.z (u : PNat.XgcdType) : ℕ+

z = zp + 1

def PNat.XgcdType.a (u : PNat.XgcdType) : ℕ+

a = ap + 1

def PNat.XgcdType.b (u : PNat.XgcdType) : ℕ+

b = bp + 1

def PNat.XgcdType.r (u : PNat.XgcdType) : ℕ

r = a % b: remainder

def PNat.XgcdType.q (u : PNat.XgcdType) : ℕ

q = ap / bp: quotient

def PNat.XgcdType.qp (u : PNat.XgcdType) : ℕ

qp = q - 1

def PNat.XgcdType.vp (u : PNat.XgcdType) : ℕ × ℕ

The map v gives the product of the matrix [[w, x], [y, z]] = [[wp + 1, x], [y, zp + 1]] and the vector [a, b] = [ap + 1, bp + 1]. The map vp gives [sp, tp] such that v = [sp + 1, tp + 1].

def PNat.XgcdType.v (u : PNat.XgcdType) : ℕ × ℕ

v = [sp + 1, tp + 1], check vp

def PNat.XgcdType.succ₂ (t : ℕ × ℕ) : ℕ × ℕ

succ₂ [t.1, t.2] = [t.1.succ, t.2.succ]

theorem PNat.XgcdType.v_eq_succ_vp (u : PNat.XgcdType) : PNat.XgcdType.v u = PNat.XgcdType.succ₂ (PNat.XgcdType.vp u)

def PNat.XgcdType.IsSpecial (u : PNat.XgcdType) : Prop

IsSpecial holds if the matrix has determinant one.

def PNat.XgcdType.IsSpecial' (u : PNat.XgcdType) : Prop

IsSpecial' is an alternative of IsSpecial.

theorem PNat.XgcdType.isSpecial_iff (u : PNat.XgcdType) : PNat.XgcdType.IsSpecial u ↔ PNat.XgcdType.IsSpecial' u

def PNat.XgcdType.IsReduced (u : PNat.XgcdType) : Prop

IsReduced holds if the two entries in the vector are the same. The reduction algorithm will produce a system with this property, whose product vector is the same as for the original system.

def PNat.XgcdType.IsReduced' (u : PNat.XgcdType) : Prop

IsReduced' is an alternative of IsReduced.

theorem PNat.XgcdType.isReduced_iff (u : PNat.XgcdType) : PNat.XgcdType.IsReduced u ↔ PNat.XgcdType.IsReduced' u

def PNat.XgcdType.flip (u : PNat.XgcdType) : PNat.XgcdType

flip flips the placement of variables during the algorithm.

@[simp]

theorem PNat.XgcdType.flip_w (u : PNat.XgcdType) : PNat.XgcdType.w (PNat.XgcdType.flip u) = PNat.XgcdType.z u

@[simp]

theorem PNat.XgcdType.flip_x (u : PNat.XgcdType) : (PNat.XgcdType.flip u).x = u.y

@[simp]

theorem PNat.XgcdType.flip_y (u : PNat.XgcdType) : (PNat.XgcdType.flip u).y = u.x

@[simp]

theorem PNat.XgcdType.flip_z (u : PNat.XgcdType) : PNat.XgcdType.z (PNat.XgcdType.flip u) = PNat.XgcdType.w u

@[simp]

theorem PNat.XgcdType.flip_a (u : PNat.XgcdType) : PNat.XgcdType.a (PNat.XgcdType.flip u) = PNat.XgcdType.b u

@[simp]

theorem PNat.XgcdType.flip_b (u : PNat.XgcdType) : PNat.XgcdType.b (PNat.XgcdType.flip u) = PNat.XgcdType.a u

theorem PNat.XgcdType.flip_isReduced (u : PNat.XgcdType) : PNat.XgcdType.IsReduced (PNat.XgcdType.flip u) ↔ PNat.XgcdType.IsReduced u

theorem PNat.XgcdType.flip_isSpecial (u : PNat.XgcdType) : PNat.XgcdType.IsSpecial (PNat.XgcdType.flip u) ↔ PNat.XgcdType.IsSpecial u

theorem PNat.XgcdType.flip_v (u : PNat.XgcdType) : PNat.XgcdType.v (PNat.XgcdType.flip u) = Prod.swap (PNat.XgcdType.v u)

theorem PNat.XgcdType.rq_eq (u : PNat.XgcdType) : PNat.XgcdType.r u + (u.bp + 1) * PNat.XgcdType.q u = u.ap + 1

Properties of division with remainder for a / b.

theorem PNat.XgcdType.qp_eq (u : PNat.XgcdType) (hr : PNat.XgcdType.r u = 0) : PNat.XgcdType.q u = PNat.XgcdType.qp u + 1

def PNat.XgcdType.start (a : ℕ+) (b : ℕ+) : PNat.XgcdType

The following function provides the starting point for our algorithm. We will apply an iterative reduction process to it, which will produce a system satisfying IsReduced. The gcd can be read off from this final system.

theorem PNat.XgcdType.start_isSpecial (a : ℕ+) (b : ℕ+) : PNat.XgcdType.IsSpecial (PNat.XgcdType.start a b)

theorem PNat.XgcdType.start_v (a : ℕ+) (b : ℕ+) : PNat.XgcdType.v (PNat.XgcdType.start a b) = (↑a, ↑b)

def PNat.XgcdType.finish (u : PNat.XgcdType) : PNat.XgcdType

finish happens when the reducing process ends.

theorem PNat.XgcdType.finish_isReduced (u : PNat.XgcdType) : PNat.XgcdType.IsReduced (PNat.XgcdType.finish u)

theorem PNat.XgcdType.finish_isSpecial (u : PNat.XgcdType) (hs : PNat.XgcdType.IsSpecial u) : PNat.XgcdType.IsSpecial (PNat.XgcdType.finish u)

theorem PNat.XgcdType.finish_v (u : PNat.XgcdType) (hr : PNat.XgcdType.r u = 0) : PNat.XgcdType.v (PNat.XgcdType.finish u) = PNat.XgcdType.v u

def PNat.XgcdType.step (u : PNat.XgcdType) : PNat.XgcdType

This is the main reduction step, which is used when u.r ≠ 0, or equivalently b does not divide a.

theorem PNat.XgcdType.step_wf (u : PNat.XgcdType) (hr : PNat.XgcdType.r u ≠ 0) : sizeOf (PNat.XgcdType.step u) < sizeOf u

We will apply the above step recursively. The following result is used to ensure that the process terminates.

theorem PNat.XgcdType.step_isSpecial (u : PNat.XgcdType) (hs : PNat.XgcdType.IsSpecial u) : PNat.XgcdType.IsSpecial (PNat.XgcdType.step u)

theorem PNat.XgcdType.step_v (u : PNat.XgcdType) (hr : PNat.XgcdType.r u ≠ 0) : PNat.XgcdType.v (PNat.XgcdType.step u) = Prod.swap (PNat.XgcdType.v u)

The reduction step does not change the product vector.

def PNat.XgcdType.reduce (u : PNat.XgcdType) : PNat.XgcdType

We can now define the full reduction function, which applies step as long as possible, and then applies finish. Note that the "have" statement puts a fact in the local context, and the equation compiler uses this fact to help construct the full definition in terms of well-founded recursion. The same fact needs to be introduced in all the inductive proofs of properties given below.

Equations

• PNat.XgcdType.reduce u = if x : PNat.XgcdType.r u = 0 then PNat.XgcdType.finish u else PNat.XgcdType.flip (PNat.XgcdType.reduce (PNat.XgcdType.step u))

theorem PNat.XgcdType.reduce_a {u : PNat.XgcdType} (h : PNat.XgcdType.r u = 0) : PNat.XgcdType.reduce u = PNat.XgcdType.finish u

theorem PNat.XgcdType.reduce_b {u : PNat.XgcdType} (h : PNat.XgcdType.r u ≠ 0) : PNat.XgcdType.reduce u = PNat.XgcdType.flip (PNat.XgcdType.reduce (PNat.XgcdType.step u))

theorem PNat.XgcdType.reduce_isReduced (u : PNat.XgcdType) : PNat.XgcdType.IsReduced (PNat.XgcdType.reduce u)

theorem PNat.XgcdType.reduce_isReduced' (u : PNat.XgcdType) : PNat.XgcdType.IsReduced' (PNat.XgcdType.reduce u)

theorem PNat.XgcdType.reduce_isSpecial (u : PNat.XgcdType) : PNat.XgcdType.IsSpecial u → PNat.XgcdType.IsSpecial (PNat.XgcdType.reduce u)

theorem PNat.XgcdType.reduce_isSpecial' (u : PNat.XgcdType) (hs : PNat.XgcdType.IsSpecial u) : PNat.XgcdType.IsSpecial' (PNat.XgcdType.reduce u)

theorem PNat.XgcdType.reduce_v (u : PNat.XgcdType) : PNat.XgcdType.v (PNat.XgcdType.reduce u) = PNat.XgcdType.v u

def PNat.xgcd (a : ℕ+) (b : ℕ+) : PNat.XgcdType

Extended Euclidean algorithm

def PNat.gcdD (a : ℕ+) (b : ℕ+) : ℕ+

gcdD a b = gcd a b

def PNat.gcdW (a : ℕ+) (b : ℕ+) : ℕ+

Final value of w

def PNat.gcdX (a : ℕ+) (b : ℕ+) : ℕ

Final value of x

def PNat.gcdY (a : ℕ+) (b : ℕ+) : ℕ

Final value of y

def PNat.gcdZ (a : ℕ+) (b : ℕ+) : ℕ+

Final value of z

def PNat.gcdA' (a : ℕ+) (b : ℕ+) : ℕ+

Final value of a / d

def PNat.gcdB' (a : ℕ+) (b : ℕ+) : ℕ+

Final value of b / d

theorem PNat.gcdA'_coe (a : ℕ+) (b : ℕ+) : ↑(PNat.gcdA' a b) = ↑(PNat.gcdW a b) + PNat.gcdX a b

theorem PNat.gcdB'_coe (a : ℕ+) (b : ℕ+) : ↑(PNat.gcdB' a b) = PNat.gcdY a b + ↑(PNat.gcdZ a b)

theorem PNat.gcd_props (a : ℕ+) (b : ℕ+) : let d := PNat.gcdD a b; let w := PNat.gcdW a b; let x := PNat.gcdX a b; let y := PNat.gcdY a b; let z := PNat.gcdZ a b; let a' := PNat.gcdA' a b; let b' := PNat.gcdB' a b; w * z = Nat.succPNat (x * y) ∧ a = a' * d ∧ b = b' * d ∧ z * a' = Nat.succPNat (x * ↑b') ∧ w * b' = Nat.succPNat (y * ↑a') ∧ ↑z * ↑a = x * ↑b + ↑d ∧ ↑w * ↑b = y * ↑a + ↑d

theorem PNat.gcd_eq (a : ℕ+) (b : ℕ+) : PNat.gcdD a b = PNat.gcd a b

theorem PNat.gcd_det_eq (a : ℕ+) (b : ℕ+) : PNat.gcdW a b * PNat.gcdZ a b = Nat.succPNat (PNat.gcdX a b * PNat.gcdY a b)

theorem PNat.gcd_a_eq (a : ℕ+) (b : ℕ+) : a = PNat.gcdA' a b * PNat.gcd a b

theorem PNat.gcd_b_eq (a : ℕ+) (b : ℕ+) : b = PNat.gcdB' a b * PNat.gcd a b

theorem PNat.gcd_rel_left' (a : ℕ+) (b : ℕ+) : PNat.gcdZ a b * PNat.gcdA' a b = Nat.succPNat (PNat.gcdX a b * ↑(PNat.gcdB' a b))

theorem PNat.gcd_rel_right' (a : ℕ+) (b : ℕ+) : PNat.gcdW a b * PNat.gcdB' a b = Nat.succPNat (PNat.gcdY a b * ↑(PNat.gcdA' a b))

theorem PNat.gcd_rel_left (a : ℕ+) (b : ℕ+) : ↑(PNat.gcdZ a b) * ↑a = PNat.gcdX a b * ↑b + ↑(PNat.gcd a b)

theorem PNat.gcd_rel_right (a : ℕ+) (b : ℕ+) : ↑(PNat.gcdW a b) * ↑b = PNat.gcdY a b * ↑a + ↑(PNat.gcd a b)