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.

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