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.

@[implicit_reducible]

instance PNat.instInhabitedXgcdType : Inhabited XgcdType

def PNat.instInhabitedXgcdType.default : XgcdType

@[implicit_reducible]

instance PNat.XgcdType.instSizeOf : SizeOf XgcdType

@[implicit_reducible]

instance PNat.XgcdType.instRepr : Repr 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 : ℕ+) : XgcdType

Another mk using ℕ and ℕ+

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

w = wp + 1

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

z = zp + 1

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

a = ap + 1

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

b = bp + 1

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

r = a % b: remainder

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

q = ap / bp: quotient

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

qp = q - 1

def PNat.XgcdType.vp (u : 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 : 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 : XgcdType) : u.v = succ₂ u.vp

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

IsSpecial holds if the matrix has determinant one.

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

IsSpecial' is an alternative of IsSpecial.

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

def PNat.XgcdType.IsReduced (u : 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 : XgcdType) : Prop

IsReduced' is an alternative of IsReduced.

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

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

flip flips the placement of variables during the algorithm.

@[simp]

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

@[simp]

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

@[simp]

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

@[simp]

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

@[simp]

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

@[simp]

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

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

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

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

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

Properties of division with remainder for a / b.

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

def PNat.XgcdType.start (a b : ℕ+) : 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 : ℕ+) : (start a b).IsSpecial

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

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

finish happens when the reducing process ends.

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

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

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

def PNat.XgcdType.step (u : XgcdType) : 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 : XgcdType) (hr : u.r ≠ 0) : sizeOf u.step < 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 : XgcdType) (hs : u.IsSpecial) : u.step.IsSpecial

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

The reduction step does not change the product vector.

@[irreducible]

def PNat.XgcdType.reduce (u : XgcdType) : 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.

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

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

@[irreducible]

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

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

@[irreducible]

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

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

@[irreducible]

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

def PNat.xgcd (a b : ℕ+) : 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 : ℕ+) : ↑(a.gcdA' b) = ↑(a.gcdW b) + a.gcdX b

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

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

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

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

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

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

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

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

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

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