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.

Instances For
    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
                        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
                                    Instances For
                                      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 }
                                      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 < 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
                                          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
                                            @[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
                                            def PNat.xgcd (a : ℕ+) (b : ℕ+) :

                                            Extended Euclidean algorithm

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