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.instSizeOf : SizeOf PNat.XgcdType

Equations

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

instance PNat.XgcdType.instRepr : 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 := (↑w).pred, x := x, y := y, zp := (↑z).pred, ap := (↑a).pred, bp := (↑b).pred }

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

w = wp + 1

Equations

• u.w = u.wp.succPNat

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

z = zp + 1

Equations

• u.z = u.zp.succPNat

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

a = ap + 1

Equations

• u.a = u.ap.succPNat

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

b = bp + 1

Equations

• u.b = u.bp.succPNat

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

r = a % b: remainder

Equations

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

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

q = ap / bp: quotient

Equations

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

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

qp = q - 1

Equations

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

Equations

• u.vp = (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

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

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

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

Equations

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

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

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

IsSpecial holds if the matrix has determinant one.

Equations

• u.IsSpecial = (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

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

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

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

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

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

IsReduced' is an alternative of IsReduced.

Equations

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

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

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

flip flips the placement of variables during the algorithm.

Equations

• u.flip = { 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) : u.flip.w = u.z

@[simp]

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

@[simp]

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

@[simp]

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

@[simp]

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

@[simp]

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

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

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

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

theorem PNat.XgcdType.rq_eq (u : PNat.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 : PNat.XgcdType) (hr : u.r = 0) : u.q = u.qp + 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.start a b).IsSpecial

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

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

finish happens when the reducing process ends.

Equations

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

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

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

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

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

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

theorem PNat.XgcdType.step_wf (u : PNat.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 : PNat.XgcdType) (hs : u.IsSpecial) : u.step.IsSpecial

theorem PNat.XgcdType.step_v (u : PNat.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 : 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

• u.reduce = if x : u.r = 0 then u.finish else u.step.reduce.flip

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

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

@[irreducible]

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

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

@[irreducible]

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

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

@[irreducible]

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

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

Extended Euclidean algorithm

Equations

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

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

gcdD a b = gcd a b

Equations

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

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

Final value of w

Equations

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

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

Final value of x

Equations

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

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

Final value of y

Equations

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

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

Final value of z

Equations

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

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

Final value of a / d

Equations

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

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

Final value of b / d

Equations

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

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 : ℕ+) : let d := a.gcdD b; let w := a.gcdW b; let x := a.gcdX b; let y := a.gcdY b; let z := a.gcdZ b; let a' := a.gcdA' b; let 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)