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 : Nat

  wp is a variable which changes through the algorithm.

• x : Nat

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

• y : Nat

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

• zp : Nat

  zp is a variable which changes through the algorithm.

• ap : Nat

  ap is a variable which changes through the algorithm.

• bp : Nat

  bp is a variable which changes through the algorithm.

def PNat.instInhabitedXgcdType : Inhabited XgcdType

Equations

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

def PNat.XgcdType.instSizeOf : SizeOf XgcdType

Equations

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

def PNat.XgcdType.instRepr : Repr 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 : PNat) (x y : Nat) (z a b : PNat) : XgcdType

Another mk using ℕ and ℕ+

Equations

• Eq (PNat.XgcdType.mk' w x y z a b) { wp := w.val.pred, x := x, y := y, zp := z.val.pred, ap := a.val.pred, bp := b.val.pred }

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

w = wp + 1

Equations

• Eq u.w u.wp.succPNat

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

z = zp + 1

Equations

• Eq u.z u.zp.succPNat

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

a = ap + 1

Equations

• Eq u.a u.ap.succPNat

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

b = bp + 1

Equations

• Eq u.b u.bp.succPNat

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

r = a % b: remainder

Equations

• Eq u.r (HMod.hMod (HAdd.hAdd u.ap 1) (HAdd.hAdd u.bp 1))

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

q = ap / bp: quotient

Equations

• Eq u.q (HDiv.hDiv (HAdd.hAdd u.ap 1) (HAdd.hAdd u.bp 1))

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

qp = q - 1

Equations

• Eq u.qp (HSub.hSub u.q 1)

def PNat.XgcdType.vp (u : XgcdType) : Prod Nat Nat

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

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

def PNat.XgcdType.v (u : XgcdType) : Prod Nat Nat

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

Equations

• Eq u.v { fst := HAdd.hAdd (HMul.hMul u.w.val u.a.val) (HMul.hMul u.x u.b.val), snd := HAdd.hAdd (HMul.hMul u.y u.a.val) (HMul.hMul u.z.val u.b.val) }

def PNat.XgcdType.succ₂ (t : Prod Nat Nat) : Prod Nat Nat

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

Equations

• Eq (PNat.XgcdType.succ₂ t) { fst := t.fst.succ, snd := t.snd.succ }

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

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

IsSpecial holds if the matrix has determinant one.

Equations

• Eq u.IsSpecial (Eq (HAdd.hAdd (HAdd.hAdd u.wp u.zp) (HMul.hMul u.wp u.zp)) (HMul.hMul u.x u.y))

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

IsSpecial' is an alternative of IsSpecial.

Equations

• Eq u.IsSpecial' (Eq (HMul.hMul u.w u.z) (HMul.hMul u.x u.y).succPNat)

theorem PNat.XgcdType.isSpecial_iff (u : XgcdType) : Iff 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.

Equations

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

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

IsReduced' is an alternative of IsReduced.

Equations

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

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

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

flip flips the placement of variables during the algorithm.

Equations

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

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

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

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

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

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

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

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

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

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

theorem PNat.XgcdType.rq_eq (u : XgcdType) : Eq (HAdd.hAdd u.r (HMul.hMul (HAdd.hAdd u.bp 1) u.q)) (HAdd.hAdd u.ap 1)

Properties of division with remainder for a / b.

theorem PNat.XgcdType.qp_eq (u : XgcdType) (hr : Eq u.r 0) : Eq u.q (HAdd.hAdd 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

• Eq (PNat.XgcdType.start a b) { wp := 0, x := 0, y := 0, zp := 0, ap := HSub.hSub a.val 1, bp := HSub.hSub b.val 1 }

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

theorem PNat.XgcdType.start_v (a b : PNat) : Eq (start a b).v { fst := a.val, snd := b.val }

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

finish happens when the reducing process ends.

Equations

• Eq u.finish { wp := u.wp, x := HAdd.hAdd (HMul.hMul (HAdd.hAdd u.wp 1) u.qp) u.x, y := u.y, zp := HAdd.hAdd (HMul.hMul u.y u.qp) u.zp, ap := u.bp, bp := u.bp }

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 : Eq u.r 0) : Eq 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.

Equations

• Eq u.step { wp := HAdd.hAdd (HMul.hMul u.y u.q) u.zp, x := u.y, y := HAdd.hAdd (HMul.hMul (HAdd.hAdd u.wp 1) u.q) u.x, zp := u.wp, ap := u.bp, bp := HSub.hSub u.r 1 }

theorem PNat.XgcdType.step_wf (u : XgcdType) (hr : Ne u.r 0) : LT.lt (SizeOf.sizeOf u.step) (SizeOf.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 : Ne u.r 0) : Eq 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.

Equations

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

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

theorem PNat.XgcdType.reduce_b {u : XgcdType} (h : Ne u.r 0) : Eq 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) : Eq u.reduce.v u.v

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

Extended Euclidean algorithm

Equations

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

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

gcdD a b = gcd a b

Equations

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

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

Final value of w

Equations

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

def PNat.gcdX (a b : PNat) : Nat

Final value of x

Equations

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

def PNat.gcdY (a b : PNat) : Nat

Final value of y

Equations

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

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

Final value of z

Equations

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

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

Final value of a / d

Equations

• Eq (a.gcdA' b) (HAdd.hAdd (a.xgcd b).wp (a.xgcd b).x).succPNat

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

Final value of b / d

Equations

• Eq (a.gcdB' b) (HAdd.hAdd (a.xgcd b).y (a.xgcd b).zp).succPNat

theorem PNat.gcdA'_coe (a b : PNat) : Eq (a.gcdA' b).val (HAdd.hAdd (a.gcdW b).val (a.gcdX b))

theorem PNat.gcdB'_coe (a b : PNat) : Eq (a.gcdB' b).val (HAdd.hAdd (a.gcdY b) (a.gcdZ b).val)

theorem PNat.gcd_props (a b : PNat) : 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; And (Eq (HMul.hMul w z) (HMul.hMul x y).succPNat) (And (Eq a (HMul.hMul a' d)) (And (Eq b (HMul.hMul b' d)) (And (Eq (HMul.hMul z a') (HMul.hMul x b'.val).succPNat) (And (Eq (HMul.hMul w b') (HMul.hMul y a'.val).succPNat) (And (Eq (HMul.hMul z.val a.val) (HAdd.hAdd (HMul.hMul x b.val) d.val)) (Eq (HMul.hMul w.val b.val) (HAdd.hAdd (HMul.hMul y a.val) d.val)))))))

theorem PNat.gcd_eq (a b : PNat) : Eq (a.gcdD b) (a.gcd b)

theorem PNat.gcd_det_eq (a b : PNat) : Eq (HMul.hMul (a.gcdW b) (a.gcdZ b)) (HMul.hMul (a.gcdX b) (a.gcdY b)).succPNat

theorem PNat.gcd_a_eq (a b : PNat) : Eq a (HMul.hMul (a.gcdA' b) (a.gcd b))

theorem PNat.gcd_b_eq (a b : PNat) : Eq b (HMul.hMul (a.gcdB' b) (a.gcd b))

theorem PNat.gcd_rel_left' (a b : PNat) : Eq (HMul.hMul (a.gcdZ b) (a.gcdA' b)) (HMul.hMul (a.gcdX b) (a.gcdB' b).val).succPNat

theorem PNat.gcd_rel_right' (a b : PNat) : Eq (HMul.hMul (a.gcdW b) (a.gcdB' b)) (HMul.hMul (a.gcdY b) (a.gcdA' b).val).succPNat

theorem PNat.gcd_rel_left (a b : PNat) : Eq (HMul.hMul (a.gcdZ b).val a.val) (HAdd.hAdd (HMul.hMul (a.gcdX b) b.val) (a.gcd b).val)

theorem PNat.gcd_rel_right (a b : PNat) : Eq (HMul.hMul (a.gcdW b).val b.val) (HAdd.hAdd (HMul.hMul (a.gcdY b) a.val) (a.gcd b).val)