Documentation

Mathlib.Data.PNat.Xgcd

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:

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 #

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.

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 }

w = wp + 1

Equations
  • u.w = u.wp.succPNat

z = zp + 1

Equations
  • u.z = u.zp.succPNat

a = ap + 1

Equations
  • u.a = u.ap.succPNat

b = bp + 1

Equations
  • u.b = u.bp.succPNat

r = a % b: remainder

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

q = ap / bp: quotient

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

qp = q - 1

Equations
  • u.qp = u.q - 1

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)

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)

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

Equations

IsSpecial holds if the matrix has determinant one.

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

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'

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)

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'

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

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
theorem PNat.XgcdType.start_v (a b : ℕ+) :
(PNat.XgcdType.start 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 }
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 }
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]

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

Extended Euclidean algorithm

Equations
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)