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

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].

• 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.

instance PNat.instInhabitedXgcdType : Inhabited PNat.XgcdType

Equations

• PNat.instInhabitedXgcdType = { default := { wp := default, x := default, y := default, zp := default, ap := default, bp := default } }

instance PNat.XgcdType.instSizeOfXgcdType : SizeOf PNat.XgcdType

Equations

• PNat.XgcdType.instSizeOfXgcdType = { sizeOf := fun (u : PNat.XgcdType) => u.bp }

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.

Equations

• One or more equations did not get rendered due to their size.

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

Another mk using ℕ and ℕ+

Equations

• PNat.XgcdType.mk' w x y z a b = { wp := Nat.pred ↑w, x := x, y := y, zp := Nat.pred ↑z, ap := Nat.pred ↑a, bp := Nat.pred ↑b }

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

w = wp + 1

Equations

• PNat.XgcdType.w u = Nat.succPNat u.wp

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

z = zp + 1

Equations

• PNat.XgcdType.z u = Nat.succPNat u.zp

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

a = ap + 1

Equations

• PNat.XgcdType.a u = Nat.succPNat u.ap

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

b = bp + 1

Equations

• PNat.XgcdType.b u = Nat.succPNat u.bp

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

r = a % b: remainder

Equations

• PNat.XgcdType.r u = (u.ap + 1) % (u.bp + 1)

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

q = ap / bp: quotient

Equations

• PNat.XgcdType.q u = (u.ap + 1) / (u.bp + 1)

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

qp = q - 1

Equations

• PNat.XgcdType.qp u = PNat.XgcdType.q u - 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].

Equations

• PNat.XgcdType.vp u = (u.wp + u.x + u.ap + u.wp * u.ap + u.x * u.bp, u.y + u.zp + u.bp + u.y * u.ap + u.zp * u.bp)

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

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

Equations

• PNat.XgcdType.v u = (↑(PNat.XgcdType.w u) * ↑(PNat.XgcdType.a u) + u.x * ↑(PNat.XgcdType.b u), u.y * ↑(PNat.XgcdType.a u) + ↑(PNat.XgcdType.z u) * ↑(PNat.XgcdType.b u))

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

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

Equations

• PNat.XgcdType.succ₂ t = (Nat.succ t.1, Nat.succ t.2)

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.

Equations

• PNat.XgcdType.IsSpecial u = (u.wp + u.zp + u.wp * u.zp = u.x * u.y)

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

IsSpecial' is an alternative of IsSpecial.

Equations

• PNat.XgcdType.IsSpecial' u = (PNat.XgcdType.w u * PNat.XgcdType.z u = Nat.succPNat (u.x * u.y))

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.

Equations

• PNat.XgcdType.IsReduced u = (u.ap = u.bp)

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

IsReduced' is an alternative of IsReduced.

Equations

• PNat.XgcdType.IsReduced' u = (PNat.XgcdType.a u = PNat.XgcdType.b u)

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.

Equations

• PNat.XgcdType.flip u = { wp := u.zp, x := u.y, y := u.x, zp := u.wp, ap := u.bp, bp := u.ap }

@[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.

Equations

• PNat.XgcdType.start a b = { wp := 0, x := 0, y := 0, zp := 0, ap := ↑a - 1, bp := ↑b - 1 }

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.

Equations

• PNat.XgcdType.finish u = { wp := u.wp, x := (u.wp + 1) * PNat.XgcdType.qp u + u.x, y := u.y, zp := u.y * PNat.XgcdType.qp u + u.zp, ap := u.bp, bp := u.bp }

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.

Equations

• PNat.XgcdType.step u = { wp := u.y * PNat.XgcdType.q u + u.zp, x := u.y, y := (u.wp + 1) * PNat.XgcdType.q u + u.x, zp := u.wp, ap := u.bp, bp := PNat.XgcdType.r u - 1 }

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

Equations

• PNat.xgcd a b = PNat.XgcdType.reduce (PNat.XgcdType.start a b)

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

gcdD a b = gcd a b

Equations

• PNat.gcdD a b = PNat.XgcdType.a (PNat.xgcd a b)

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

Final value of w

Equations

• PNat.gcdW a b = PNat.XgcdType.w (PNat.xgcd a b)

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

Final value of x

Equations

• PNat.gcdX a b = (PNat.xgcd a b).x

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

Final value of y

Equations

• PNat.gcdY a b = (PNat.xgcd a b).y

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

Final value of z

Equations

• PNat.gcdZ a b = PNat.XgcdType.z (PNat.xgcd a b)

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

Final value of a / d

Equations

• PNat.gcdA' a b = Nat.succPNat ((PNat.xgcd a b).wp + (PNat.xgcd a b).x)

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

Final value of b / d

Equations

• PNat.gcdB' a b = Nat.succPNat ((PNat.xgcd a b).y + (PNat.xgcd a b).zp)

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)