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

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.

Instances For
Equations
Equations

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 : ℕ+) :

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 }
Instances For

w = wp + 1

Equations
• u.w = u.wp.succPNat
Instances For

z = zp + 1

Equations
• u.z = u.zp.succPNat
Instances For

a = ap + 1

Equations
• u.a = u.ap.succPNat
Instances For

b = bp + 1

Equations
• u.b = u.bp.succPNat
Instances For

r = a % b: remainder

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

q = ap / bp: quotient

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

qp = q - 1

Equations
• u.qp = u.q - 1
Instances For

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)
Instances For

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)
Instances For

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

Equations
• = (t.1.succ, t.2.succ)
Instances For

IsSpecial holds if the matrix has determinant one.

Equations
• u.IsSpecial = (u.wp + u.zp + u.wp * u.zp = u.x * u.y)
Instances For

IsSpecial' is an alternative of IsSpecial.

Equations
• u.IsSpecial' = (u.w * u.z = (u.x * u.y).succPNat)
Instances For
theorem PNat.XgcdType.isSpecial_iff (u : PNat.XgcdType) :
u.IsSpecial u.IsSpecial'

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)
Instances For

IsReduced' is an alternative of IsReduced.

Equations
• u.IsReduced' = (u.a = u.b)
Instances For
theorem PNat.XgcdType.isReduced_iff (u : PNat.XgcdType) :
u.IsReduced u.IsReduced'

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 }
Instances For
@[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

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
• = { wp := 0, x := 0, y := 0, zp := 0, ap := a - 1, bp := b - 1 }
Instances For
theorem PNat.XgcdType.start_isSpecial (a : ℕ+) (b : ℕ+) :
().IsSpecial
theorem PNat.XgcdType.start_v (a : ℕ+) (b : ℕ+) :
().v = (a, b)

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 }
Instances For
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

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 }
Instances For
theorem PNat.XgcdType.step_wf (u : PNat.XgcdType) (hr : u.r 0) :
sizeOf u.step <

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.

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
Instances For
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
theorem PNat.XgcdType.reduce_isReduced (u : PNat.XgcdType) :
u.reduce.IsReduced
theorem PNat.XgcdType.reduce_isReduced' (u : PNat.XgcdType) :
u.reduce.IsReduced'
theorem PNat.XgcdType.reduce_isSpecial (u : PNat.XgcdType) :
u.IsSpecialu.reduce.IsSpecial
theorem PNat.XgcdType.reduce_isSpecial' (u : PNat.XgcdType) (hs : u.IsSpecial) :
u.reduce.IsSpecial'
theorem PNat.XgcdType.reduce_v (u : PNat.XgcdType) :
u.reduce.v = u.v
def PNat.xgcd (a : ℕ+) (b : ℕ+) :

Extended Euclidean algorithm

Equations
• a.xgcd b = ().reduce
Instances For
def PNat.gcdD (a : ℕ+) (b : ℕ+) :

gcdD a b = gcd a b

Equations
• a.gcdD b = (a.xgcd b).a
Instances For
def PNat.gcdW (a : ℕ+) (b : ℕ+) :

Final value of w

Equations
• a.gcdW b = (a.xgcd b).w
Instances For
def PNat.gcdX (a : ℕ+) (b : ℕ+) :

Final value of x

Equations
• a.gcdX b = (a.xgcd b).x
Instances For
def PNat.gcdY (a : ℕ+) (b : ℕ+) :

Final value of y

Equations
• a.gcdY b = (a.xgcd b).y
Instances For
def PNat.gcdZ (a : ℕ+) (b : ℕ+) :

Final value of z

Equations
• a.gcdZ b = (a.xgcd b).z
Instances For
def PNat.gcdA' (a : ℕ+) (b : ℕ+) :

Final value of a / d

Equations
• a.gcdA' b = ((a.xgcd b).wp + (a.xgcd b).x).succPNat
Instances For
def PNat.gcdB' (a : ℕ+) (b : ℕ+) :

Final value of b / d

Equations
• a.gcdB' b = ((a.xgcd b).y + (a.xgcd b).zp).succPNat
Instances For
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)