# lib/gcd.bend -- divisibility, Bezout coprimality, Euclid's algorithm, and
# the theorem "if a b is a square and a, b are coprime then a is a square".
#
# This is the bottom of the ladder that lib/pythag.bend and
# fermat4.bend stand on. Two propositions carry everything:
#
#   Dvd<d, n>   d divides n:  a Nat k with d k = n
#   Cop<a, b>   a and b are coprime: INTEGERS u, v with u a + v b = 1
#
# Both are `Data` (a record of a witness and an equation), so a proof of
# either can be copied with `+` -- unlike a computed proposition such as
# Nat.Le, which has to be copied in continuation-passing style.
#
# Coprimality is Bezout's identity rather than "gcd = 1" because every use
# of it below is an algebraic one: Gauss's lemma (a | b c and a coprime to
# b give a | c) and the closure of coprimality under products are both one
# ring identity in Int, where "gcd = 1" would need a fresh induction each
# time. The integers come from int.bend + int_filled.bend, whose
# ring laws are already proven.
#
# Euclid's algorithm is run by SUBTRACTION, not by division: gcd(a, b) with
# b = a + c recurses on (a, c). That keeps the whole file free of divmod
# (so it does not need a divmod library, and inherits no
# TODOs), and the Bezout coefficients transport through a subtraction step
# with one rewrite -- u a + v c = g becomes (u - v) a + v b = g.

import Base
import ./nat.bend as N
import ./sqrt2.bend as S2
import ./int.bend as I
import ./int_filled.bend as IF

# ---------------------------------------------------------------------
# Nat odds and ends that lib/nat.bend does not have
# ---------------------------------------------------------------------

# a (1+b) = 0 forces a = 0
def gc.mul_zero_l(a: Nat, +b: Nat, e: {Nat.mul(a, 1n+b) == 0n : Nat}) -> {a == 0n : Nat}:
  match a:
    case 0n:
      {==}
    case 1n+p:
      Empty.absurd({1n+p == 0n : Nat}, N.Nat.succ_neq_zero(Nat.add(b, Nat.mul(p, 1n+b)), e))

# a b = 1 forces a = 1
def gc.mul_eq_one(a: Nat, b: Nat, e: {Nat.mul(a, b) == 1n : Nat}) -> {a == 1n : Nat}:
  match a:
    case 0n:
      Empty.absurd({0n == 1n : Nat}, N.Nat.zero_neq_succ(0n, e))
    case 1n+ap:
      match b:
        case 0n:
          Empty.absurd({1n+ap == 1n : Nat}, N.Nat.zero_neq_succ(0n,
            Equal.trans(Nat, 0n, Nat.mul(1n+ap, 0n), 1n,
              Equal.sym(Nat, Nat.mul(1n+ap, 0n), 0n, N.Nat.mul_zero(1n+ap)), e)))
        case 1n+bp:
          +ap2 = ap
          +bp2 = bp
          Equal.cong(Nat, Nat, z => 1n+z, ap2, 0n,
            gc.mul_zero_l(ap2, bp2,
              N.Nat.add_eq_zero_r(bp2, Nat.mul(ap2, 1n+bp2),
                N.Nat.succ_inj(Nat.add(bp2, Nat.mul(ap2, 1n+bp2)), 0n, e))))

# a positive number is not zero
def gc.succ_pos(+a: Nat) -> N.Nat.Lt(0n, 1n+a):
  Unit{}

# ---------------------------------------------------------------------
# the Int ring toolkit: the five identities the Bezout algebra needs
# ---------------------------------------------------------------------

# (x + y) z = x z + y z
def gc.iadd_mul(+x: I.Int, +y: I.Int, +z: I.Int) -> {I.Int.mul(I.Int.add(x, y), z) == I.Int.add(I.Int.mul(x, z), I.Int.mul(y, z)) : I.Int}:
  %I.mul_comm(z, x) : {I.Int.mul(I.Int.add(x, y), z) == I.Int.add(_, I.Int.mul(y, z)) : I.Int}
  %I.mul_comm(z, y) : {I.Int.mul(I.Int.add(x, y), z) == I.Int.add(I.Int.mul(z, x), _) : I.Int}
  %I.mul_add(z, x, y) : {I.Int.mul(I.Int.add(x, y), z) == _ : I.Int}
  I.mul_comm(I.Int.add(x, y), z)

# x (y z) = y (x z)
def gc.imul_swap(+x: I.Int, +y: I.Int, +z: I.Int) -> {I.Int.mul(x, I.Int.mul(y, z)) == I.Int.mul(y, I.Int.mul(x, z)) : I.Int}:
  %Equal.sym(I.Int, I.Int.mul(x, I.Int.mul(y, z)), I.Int.mul(I.Int.mul(x, y), z), I.mul_assoc(x, y, z)) : {_ == I.Int.mul(y, I.Int.mul(x, z)) : I.Int}
  %Equal.sym(I.Int, I.Int.mul(y, I.Int.mul(x, z)), I.Int.mul(I.Int.mul(y, x), z), I.mul_assoc(y, x, z)) : {I.Int.mul(I.Int.mul(x, y), z) == _ : I.Int}
  Equal.cong(I.Int, I.Int, w => I.Int.mul(w, z), I.Int.mul(x, y), I.Int.mul(y, x), I.mul_comm(x, y))

# (x y) z = (x z) y
def gc.imul_shift(+x: I.Int, +y: I.Int, +z: I.Int) -> {I.Int.mul(I.Int.mul(x, y), z) == I.Int.mul(I.Int.mul(x, z), y) : I.Int}:
  %I.mul_assoc(x, z, y) : {I.Int.mul(I.Int.mul(x, y), z) == _ : I.Int}
  %I.mul_assoc(x, y, z) : {_ == I.Int.mul(x, I.Int.mul(z, y)) : I.Int}
  Equal.cong(I.Int, I.Int, w => I.Int.mul(x, w), I.Int.mul(y, z), I.Int.mul(z, y), I.mul_comm(y, z))

# x + (y + z) = y + (x + z)
def gc.iadd_swap(+x: I.Int, +y: I.Int, +z: I.Int) -> {I.Int.add(x, I.Int.add(y, z)) == I.Int.add(y, I.Int.add(x, z)) : I.Int}:
  %Equal.sym(I.Int, I.Int.add(x, I.Int.add(y, z)), I.Int.add(I.Int.add(x, y), z), I.add_assoc(x, y, z)) : {_ == I.Int.add(y, I.Int.add(x, z)) : I.Int}
  %Equal.sym(I.Int, I.Int.add(y, I.Int.add(x, z)), I.Int.add(I.Int.add(y, x), z), I.add_assoc(y, x, z)) : {I.Int.add(I.Int.add(x, y), z) == _ : I.Int}
  Equal.cong(I.Int, I.Int, w => I.Int.add(w, z), I.Int.add(x, y), I.Int.add(y, x), I.add_comm(x, y))

# 1 x = x
def gc.ione_mul(+x: I.Int) -> {I.Int.mul(I.Pos{1n}, x) == x : I.Int}:
  Equal.trans(I.Int, I.Int.mul(I.Pos{1n}, x), I.Int.mul(x, I.Pos{1n}), x,
    I.mul_comm(I.Pos{1n}, x), I.mul_one(x))

# the absolute value of Pos{n} is n, which is what carries an Int equation
# back into Nat: abs(Pos{n}) = abs(Pos{d} * x) = d * abs(x)
def gc.abs_pos(n: Nat) -> {I.Int.abs(I.Pos{n}) == n : Nat}:
  {==}

# ---------------------------------------------------------------------
# divisibility: d | n is a Nat k with d k = n
# ---------------------------------------------------------------------

type Dvd<-d: Nat, -n: Nat> is Data:
  Dvd{k: Nat, h: {Nat.mul(d, k) == n : Nat}}

def gc.dvd_refl(+d: Nat) -> Dvd<d, d>:
  Dvd{1n, N.Nat.mul_one(d)}

def gc.dvd_one(+n: Nat) -> Dvd<1n, n>:
  Dvd{n, N.Nat.one_mul(n)}

def gc.dvd_zero(d: Nat) -> Dvd<d, 0n>:
  Dvd{0n, N.Nat.mul_zero(d)}

# n = m carries d | n to d | m
def gc.dvd_eq(-d: Nat, -n: Nat, -m: Nat, e: {n == m : Nat}, h: Dvd<d, n>) -> Dvd<d, m>:
  %e : Dvd<d, _>
  h

# d | n and n | c give d | c
def gc.dvd_trans(+a: Nat, +b: Nat, +c: Nat, h1: Dvd<a, b>, h2: Dvd<b, c>) -> Dvd<a, c>:
  match h1:
    case Dvd{k, e1}:
      match h2:
        case Dvd{j, e2}:
          +j2 = j
          +k2 = k
          Dvd{Nat.mul(k2, j2),
            Equal.trans(Nat, Nat.mul(a, Nat.mul(k2, j2)), Nat.mul(b, j2), c,
              Equal.trans(Nat, Nat.mul(a, Nat.mul(k2, j2)), Nat.mul(Nat.mul(a, k2), j2), Nat.mul(b, j2),
                N.Nat.mul_assoc(a, k2, j2),
                Equal.cong(Nat, Nat, z => Nat.mul(z, j2), Nat.mul(a, k2), b, e1)),
              e2)}

# d | n gives d | n m
def gc.dvd_mul_r(+d: Nat, +n: Nat, +m: Nat, h: Dvd<d, n>) -> Dvd<d, Nat.mul(n, m)>:
  match h:
    case Dvd{k, e}:
      +m2 = m
      +k2 = k
      Dvd{Nat.mul(k2, m2),
        Equal.trans(Nat, Nat.mul(d, Nat.mul(k2, m2)), Nat.mul(Nat.mul(d, k2), m2), Nat.mul(n, m2),
          N.Nat.mul_assoc(d, k2, m2),
          Equal.cong(Nat, Nat, z => Nat.mul(z, m2), Nat.mul(d, k2), n, e))}

# d | n gives d | m n
def gc.dvd_mul_l(+d: Nat, +n: Nat, +m: Nat, h: Dvd<d, n>) -> Dvd<d, Nat.mul(m, n)>:
  gc.dvd_eq(d, Nat.mul(n, m), Nat.mul(m, n), N.Nat.mul_comm(n, m), gc.dvd_mul_r(d, n, m, h))

# d | x and d | y give d | x + y
def gc.dvd_add(+d: Nat, +x: Nat, +y: Nat, h1: Dvd<d, x>, h2: Dvd<d, y>) -> Dvd<d, Nat.add(x, y)>:
  match h1:
    case Dvd{k, e1}:
      match h2:
        case Dvd{j, e2}:
          +k2 = k
          +j2 = j
          Dvd{Nat.add(k2, j2),
            Equal.trans(Nat, Nat.mul(d, Nat.add(k2, j2)), Nat.add(Nat.mul(d, k2), Nat.mul(d, j2)), Nat.add(x, y),
              N.Nat.mul_add(d, k2, j2),
              Equal.trans(Nat, Nat.add(Nat.mul(d, k2), Nat.mul(d, j2)), Nat.add(x, Nat.mul(d, j2)), Nat.add(x, y),
                Equal.cong(Nat, Nat, z => Nat.add(z, Nat.mul(d, j2)), Nat.mul(d, k2), x, e1),
                Equal.cong(Nat, Nat, z => Nat.add(x, z), Nat.mul(d, j2), y, e2)))}

# a divisor of a positive number is positive
def gc.dvd_pos(d: Nat, +n: Nat, h: Dvd<d, 1n+n>) -> N.Nat.Lt(0n, d):
  match d:
    case 0n:
      match h:
        case Dvd{k, e}:
          Empty.absurd(N.Nat.Lt(0n, 0n), N.Nat.zero_neq_succ(n, e))
    case 1n+dp:
      Unit{}

# a divisor of a positive number is at most it
def gc.dvd_le.go(+d: Nat, +n: Nat, k: Nat, e: {Nat.mul(d, k) == 1n+n : Nat}) -> N.Nat.Le(d, 1n+n):
  match k:
    case 0n:
      Empty.absurd(N.Nat.Le(d, 1n+n), N.Nat.zero_neq_succ(n,
        Equal.trans(Nat, 0n, Nat.mul(d, 0n), 1n+n,
          Equal.sym(Nat, Nat.mul(d, 0n), 0n, N.Nat.mul_zero(d)), e)))
    case 1n+kp:
      %e : N.Nat.Le(d, _)
      %Equal.sym(Nat, Nat.mul(d, 1n+kp), Nat.add(d, Nat.mul(d, kp)), N.Nat.mul_succ(d, kp)) : N.Nat.Le(d, _)
      N.Nat.le_add_r(d, Nat.mul(d, kp))

def gc.dvd_le(+d: Nat, +n: Nat, h: Dvd<d, 1n+n>) -> N.Nat.Le(d, 1n+n):
  match h:
    case Dvd{k, e}:
      gc.dvd_le.go(d, n, k, e)

# ---------------------------------------------------------------------
# coprimality: Bezout's identity over the integers
# ---------------------------------------------------------------------

type Cop<-a: Nat, -b: Nat> is Data:
  Cop{u: I.Int, v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}}

def gc.cop_sym(+a: Nat, +b: Nat, c: Cop<a, b>) -> Cop<b, a>:
  match c:
    case Cop{u, v, h}:
      +u2 = u
      +v2 = v
      Cop{v2, u2,
        %I.add_comm(I.Int.mul(u2, I.Pos{a}), I.Int.mul(v2, I.Pos{b})) : {_ == I.Pos{1n} : I.Int}
        h}

# 1 is coprime to everything
def gc.cop_one(b: Nat) -> Cop<1n, b>:
  Cop{I.Pos{1n}, I.Pos{0n}, {==}}

def gc.cop_eq_r(-a: Nat, -b: Nat, -c: Nat, e: {b == c : Nat}, h: Cop<a, b>) -> Cop<a, c>:
  %e : Cop<a, _>
  h

def gc.cop_eq_l(-a: Nat, -b: Nat, -c: Nat, e: {a == b : Nat}, h: Cop<a, c>) -> Cop<b, c>:
  %e : Cop<_, c>
  h

# a divisor of a is coprime to whatever a is coprime to:  a = c t and
# u a + v b = 1 give (u t) c + v b = 1
def gc.cop_dvd_l.go(+a: Nat, +b: Nat, +c: Nat, +t: Nat, et: {Nat.mul(c, t) == a : Nat}, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}) -> Cop<c, b>:
  Cop{I.Int.mul(u, I.Pos{t}), v,
    %I.mul_assoc(u, I.Pos{t}, I.Pos{c}) : {I.Int.add(_, I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}
    %I.mul_comm(I.Pos{c}, I.Pos{t}) : {I.Int.add(I.Int.mul(u, _), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}
    %Equal.sym(I.Int, I.Int.mul(I.Pos{c}, I.Pos{t}), I.Pos{Nat.mul(c, t)}, I.of_mul(c, t)) : {I.Int.add(I.Int.mul(u, _), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}
    %Equal.sym(I.Int, I.Pos{Nat.mul(c, t)}, I.Pos{a}, Equal.cong(Nat, I.Int, z => I.Pos{z}, Nat.mul(c, t), a, et)) : {I.Int.add(I.Int.mul(u, _), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}
    h}

def gc.cop_dvd_l(+a: Nat, +b: Nat, +c: Nat, hc: Cop<a, b>, hd: Dvd<c, a>) -> Cop<c, b>:
  match hc:
    case Cop{u, v, h}:
      match hd:
        case Dvd{t, et}:
          gc.cop_dvd_l.go(a, b, c, t, et, u, v, h)

# a common divisor of a coprime pair is 1:  a = w s, b = w t and
# u a + v b = 1 give w (u s + v t) = 1, so w = 1
def gc.cop_common.go(+a: Nat, +b: Nat, +w: Nat, +s: Nat, +t: Nat, es: {Nat.mul(w, s) == a : Nat}, et: {Nat.mul(w, t) == b : Nat}, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}) -> {w == 1n : Nat}:
  +x = I.Int.add(I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t}))
  +ee = Equal.trans(I.Int, I.Int.mul(I.Pos{w}, x), I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})), I.Pos{1n},
    %Equal.trans(I.Int, I.Int.mul(I.Pos{w}, I.Pos{s}), I.Pos{Nat.mul(w, s)}, I.Pos{a}, I.of_mul(w, s), Equal.cong(Nat, I.Int, z => I.Pos{z}, Nat.mul(w, s), a, es)) : {I.Int.mul(I.Pos{w}, x) == I.Int.add(I.Int.mul(u, _), I.Int.mul(v, I.Pos{b})) : I.Int}
    %Equal.trans(I.Int, I.Int.mul(I.Pos{w}, I.Pos{t}), I.Pos{Nat.mul(w, t)}, I.Pos{b}, I.of_mul(w, t), Equal.cong(Nat, I.Int, z => I.Pos{z}, Nat.mul(w, t), b, et)) : {I.Int.mul(I.Pos{w}, x) == I.Int.add(I.Int.mul(u, I.Int.mul(I.Pos{w}, I.Pos{s})), I.Int.mul(v, _)) : I.Int}
    %gc.imul_swap(I.Pos{w}, u, I.Pos{s}) : {I.Int.mul(I.Pos{w}, x) == I.Int.add(_, I.Int.mul(v, I.Int.mul(I.Pos{w}, I.Pos{t}))) : I.Int}
    %gc.imul_swap(I.Pos{w}, v, I.Pos{t}) : {I.Int.mul(I.Pos{w}, x) == I.Int.add(I.Int.mul(I.Pos{w}, I.Int.mul(u, I.Pos{s})), _) : I.Int}
    I.mul_add(I.Pos{w}, I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t})),
    h)
  gc.mul_eq_one(w, I.Int.abs(x),
    Equal.trans(Nat, Nat.mul(w, I.Int.abs(x)), I.Int.abs(I.Int.mul(I.Pos{w}, x)), 1n,
      Equal.sym(Nat, I.Int.abs(I.Int.mul(I.Pos{w}, x)), Nat.mul(w, I.Int.abs(x)), IF.Int.abs_mul(I.Pos{w}, x)),
      Equal.cong(I.Int, Nat, z => I.Int.abs(z), I.Int.mul(I.Pos{w}, x), I.Pos{1n}, ee)))

def gc.cop_common(+a: Nat, +b: Nat, +w: Nat, hc: Cop<a, b>, d1: Dvd<w, a>, d2: Dvd<w, b>) -> {w == 1n : Nat}:
  match hc:
    case Cop{u, v, h}:
      match d1:
        case Dvd{s, es}:
          match d2:
            case Dvd{t, et}:
              gc.cop_common.go(a, b, w, s, t, es, et, u, v, h)

# x (y z) = y (z x)
def gc.imul_rot(+x: I.Int, +y: I.Int, +z: I.Int) -> {I.Int.mul(x, I.Int.mul(y, z)) == I.Int.mul(y, I.Int.mul(z, x)) : I.Int}:
  Equal.trans(I.Int, I.Int.mul(x, I.Int.mul(y, z)), I.Int.mul(y, I.Int.mul(x, z)), I.Int.mul(y, I.Int.mul(z, x)),
    gc.imul_swap(x, y, z),
    Equal.cong(I.Int, I.Int, w => I.Int.mul(y, w), I.Int.mul(x, z), I.Int.mul(z, x), I.mul_comm(x, z)))

# ---------------------------------------------------------------------
# Gauss's lemma: a | b c and a coprime to b give a | c.
# c = c (u a + v b) = (c u) a + v (b c) = (c u) a + v (a k) = a (c u + v k),
# and the absolute value of that Int identity IS the Nat divisibility.
# ---------------------------------------------------------------------

# a k = b c turns a (v k) into c (v b)
def gc.gauss.mid(+a: Nat, +b: Nat, +c: Nat, +k: Nat, ek: {Nat.mul(a, k) == Nat.mul(b, c) : Nat}, +v: I.Int) -> {I.Int.mul(I.Pos{a}, I.Int.mul(v, I.Pos{k})) == I.Int.mul(I.Pos{c}, I.Int.mul(v, I.Pos{b})) : I.Int}:
  %gc.imul_swap(v, I.Pos{c}, I.Pos{b}) : {I.Int.mul(I.Pos{a}, I.Int.mul(v, I.Pos{k})) == _ : I.Int}
  %I.mul_comm(I.Pos{b}, I.Pos{c}) : {I.Int.mul(I.Pos{a}, I.Int.mul(v, I.Pos{k})) == I.Int.mul(v, _) : I.Int}
  %Equal.sym(I.Int, I.Int.mul(I.Pos{b}, I.Pos{c}), I.Pos{Nat.mul(b, c)}, I.of_mul(b, c)) : {I.Int.mul(I.Pos{a}, I.Int.mul(v, I.Pos{k})) == I.Int.mul(v, _) : I.Int}
  %Equal.cong(Nat, I.Int, z => I.Pos{z}, Nat.mul(a, k), Nat.mul(b, c), ek) : {I.Int.mul(I.Pos{a}, I.Int.mul(v, I.Pos{k})) == I.Int.mul(v, _) : I.Int}
  %I.of_mul(a, k) : {I.Int.mul(I.Pos{a}, I.Int.mul(v, I.Pos{k})) == I.Int.mul(v, _) : I.Int}
  gc.imul_swap(I.Pos{a}, v, I.Pos{k})

# the Int identity behind Gauss's lemma
def gc.gauss.key(+a: Nat, +b: Nat, +c: Nat, +k: Nat, ek: {Nat.mul(a, k) == Nat.mul(b, c) : Nat}, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}) -> {I.Int.mul(I.Pos{a}, I.Int.add(I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))) == I.Pos{c} : I.Int}:
  %I.mul_one(I.Pos{c}) : {I.Int.mul(I.Pos{a}, I.Int.add(I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))) == _ : I.Int}
  %h : {I.Int.mul(I.Pos{a}, I.Int.add(I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))) == I.Int.mul(I.Pos{c}, _) : I.Int}
  %Equal.sym(I.Int, I.Int.mul(I.Pos{c}, I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b}))), I.Int.add(I.Int.mul(I.Pos{c}, I.Int.mul(u, I.Pos{a})), I.Int.mul(I.Pos{c}, I.Int.mul(v, I.Pos{b}))), I.mul_add(I.Pos{c}, I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b}))) : {I.Int.mul(I.Pos{a}, I.Int.add(I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))) == _ : I.Int}
  %gc.imul_rot(I.Pos{a}, I.Pos{c}, u) : {I.Int.mul(I.Pos{a}, I.Int.add(I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))) == I.Int.add(_, I.Int.mul(I.Pos{c}, I.Int.mul(v, I.Pos{b}))) : I.Int}
  %gc.gauss.mid(a, b, c, k, ek, v) : {I.Int.mul(I.Pos{a}, I.Int.add(I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))) == I.Int.add(I.Int.mul(I.Pos{a}, I.Int.mul(I.Pos{c}, u)), _) : I.Int}
  I.mul_add(I.Pos{a}, I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))

def gc.gauss.go(+a: Nat, +b: Nat, +c: Nat, +k: Nat, ek: {Nat.mul(a, k) == Nat.mul(b, c) : Nat}, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}) -> Dvd<a, c>:
  +w = I.Int.add(I.Int.mul(I.Pos{c}, u), I.Int.mul(v, I.Pos{k}))
  +ee = gc.gauss.key(a, b, c, k, ek, u, v, h)
  Dvd{I.Int.abs(w),
    Equal.trans(Nat, Nat.mul(a, I.Int.abs(w)), I.Int.abs(I.Int.mul(I.Pos{a}, w)), c,
      Equal.sym(Nat, I.Int.abs(I.Int.mul(I.Pos{a}, w)), Nat.mul(a, I.Int.abs(w)), IF.Int.abs_mul(I.Pos{a}, w)),
      Equal.cong(I.Int, Nat, z => I.Int.abs(z), I.Int.mul(I.Pos{a}, w), I.Pos{c}, ee))}

def gc.gauss(+a: Nat, +b: Nat, +c: Nat, hc: Cop<a, b>, hd: Dvd<a, Nat.mul(b, c)>) -> Dvd<a, c>:
  match hc:
    case Cop{u, v, h}:
      match hd:
        case Dvd{k, ek}:
          gc.gauss.go(a, b, c, k, ek, u, v, h)

# ---------------------------------------------------------------------
# coprimality is closed under products: (u a + v b)(u' a + v' c) = 1
# regrouped so that everything but v v' b c carries a factor a
# ---------------------------------------------------------------------

def gc.bez_l1(+v: I.Int, +b: I.Int, +u2: I.Int, +a: I.Int) -> {I.Int.mul(I.Int.mul(v, b), I.Int.mul(u2, a)) == I.Int.mul(I.Int.mul(I.Int.mul(v, u2), b), a) : I.Int}:
  %gc.imul_shift(v, b, u2) : {I.Int.mul(I.Int.mul(v, b), I.Int.mul(u2, a)) == I.Int.mul(_, a) : I.Int}
  I.mul_assoc(I.Int.mul(v, b), u2, a)

def gc.bez_l2(+v: I.Int, +b: I.Int, +v2: I.Int, +c: I.Int) -> {I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c)) == I.Int.mul(I.Int.mul(v, v2), I.Int.mul(b, c)) : I.Int}:
  %I.mul_assoc(v, v2, I.Int.mul(b, c)) : {I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c)) == _ : I.Int}
  %gc.imul_swap(b, v2, c) : {I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c)) == I.Int.mul(v, _) : I.Int}
  Equal.sym(I.Int, I.Int.mul(v, I.Int.mul(b, I.Int.mul(v2, c))), I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c)), I.mul_assoc(v, b, I.Int.mul(v2, c)))

def gc.bez_prod(+u: I.Int, +v: I.Int, +u2: I.Int, +v2: I.Int, +a: I.Int, +b: I.Int, +c: I.Int) -> {I.Int.mul(I.Int.add(I.Int.mul(u, a), I.Int.mul(v, b)), I.Int.add(I.Int.mul(u2, a), I.Int.mul(v2, c))) == I.Int.add(I.Int.mul(I.Int.add(I.Int.mul(u, I.Int.add(I.Int.mul(u2, a), I.Int.mul(v2, c))), I.Int.mul(I.Int.mul(v, u2), b)), a), I.Int.mul(I.Int.mul(v, v2), I.Int.mul(b, c))) : I.Int}:
  +y = I.Int.add(I.Int.mul(u2, a), I.Int.mul(v2, c))
  +p1 = I.Int.mul(I.Int.mul(u, y), a)
  +p2 = I.Int.mul(I.Int.mul(I.Int.mul(v, u2), b), a)
  +p3 = I.Int.mul(I.Int.mul(v, v2), I.Int.mul(b, c))
  Equal.trans(I.Int, I.Int.mul(I.Int.add(I.Int.mul(u, a), I.Int.mul(v, b)), y), I.Int.add(I.Int.mul(I.Int.mul(u, a), y), I.Int.mul(I.Int.mul(v, b), y)), I.Int.add(I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), b)), a), p3),
    gc.iadd_mul(I.Int.mul(u, a), I.Int.mul(v, b), y),
    Equal.trans(I.Int, I.Int.add(I.Int.mul(I.Int.mul(u, a), y), I.Int.mul(I.Int.mul(v, b), y)), I.Int.add(p1, I.Int.add(p2, p3)), I.Int.add(I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), b)), a), p3),
      Equal.trans(I.Int, I.Int.add(I.Int.mul(I.Int.mul(u, a), y), I.Int.mul(I.Int.mul(v, b), y)), I.Int.add(p1, I.Int.mul(I.Int.mul(v, b), y)), I.Int.add(p1, I.Int.add(p2, p3)),
        Equal.cong(I.Int, I.Int, z => I.Int.add(z, I.Int.mul(I.Int.mul(v, b), y)), I.Int.mul(I.Int.mul(u, a), y), p1, gc.imul_shift(u, a, y)),
        Equal.cong(I.Int, I.Int, z => I.Int.add(p1, z), I.Int.mul(I.Int.mul(v, b), y), I.Int.add(p2, p3),
          Equal.trans(I.Int, I.Int.mul(I.Int.mul(v, b), y), I.Int.add(I.Int.mul(I.Int.mul(v, b), I.Int.mul(u2, a)), I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c))), I.Int.add(p2, p3),
            I.mul_add(I.Int.mul(v, b), I.Int.mul(u2, a), I.Int.mul(v2, c)),
            Equal.trans(I.Int, I.Int.add(I.Int.mul(I.Int.mul(v, b), I.Int.mul(u2, a)), I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c))), I.Int.add(p2, I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c))), I.Int.add(p2, p3),
              Equal.cong(I.Int, I.Int, z => I.Int.add(z, I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c))), I.Int.mul(I.Int.mul(v, b), I.Int.mul(u2, a)), p2, gc.bez_l1(v, b, u2, a)),
              Equal.cong(I.Int, I.Int, z => I.Int.add(p2, z), I.Int.mul(I.Int.mul(v, b), I.Int.mul(v2, c)), p3, gc.bez_l2(v, b, v2, c)))))),
      Equal.trans(I.Int, I.Int.add(p1, I.Int.add(p2, p3)), I.Int.add(I.Int.add(p1, p2), p3), I.Int.add(I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), b)), a), p3),
        I.add_assoc(p1, p2, p3),
        Equal.cong(I.Int, I.Int, z => I.Int.add(z, p3), I.Int.add(p1, p2), I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), b)), a),
          Equal.sym(I.Int, I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), b)), a), I.Int.add(p1, p2), gc.iadd_mul(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), b), a))))))

def gc.cop_mul.go(+a: Nat, +b: Nat, +c: Nat, +u: I.Int, +v: I.Int, +u2: I.Int, +v2: I.Int, h1: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{1n} : I.Int}, h2: {I.Int.add(I.Int.mul(u2, I.Pos{a}), I.Int.mul(v2, I.Pos{c})) == I.Pos{1n} : I.Int}) -> Cop<a, Nat.mul(b, c)>:
  +x = I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b}))
  +y = I.Int.add(I.Int.mul(u2, I.Pos{a}), I.Int.mul(v2, I.Pos{c}))
  Cop{I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), I.Pos{b})), I.Int.mul(v, v2),
    %I.of_mul(b, c) : {I.Int.add(I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), I.Pos{b})), I.Pos{a}), I.Int.mul(I.Int.mul(v, v2), _)) == I.Pos{1n} : I.Int}
    Equal.trans(I.Int, I.Int.add(I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), I.Pos{b})), I.Pos{a}), I.Int.mul(I.Int.mul(v, v2), I.Int.mul(I.Pos{b}, I.Pos{c}))), I.Int.mul(x, y), I.Pos{1n},
      Equal.sym(I.Int, I.Int.mul(x, y), I.Int.add(I.Int.mul(I.Int.add(I.Int.mul(u, y), I.Int.mul(I.Int.mul(v, u2), I.Pos{b})), I.Pos{a}), I.Int.mul(I.Int.mul(v, v2), I.Int.mul(I.Pos{b}, I.Pos{c}))),
        gc.bez_prod(u, v, u2, v2, I.Pos{a}, I.Pos{b}, I.Pos{c})),
      Equal.trans(I.Int, I.Int.mul(x, y), I.Int.mul(I.Pos{1n}, y), I.Pos{1n},
        Equal.cong(I.Int, I.Int, z => I.Int.mul(z, y), x, I.Pos{1n}, h1),
        Equal.trans(I.Int, I.Int.mul(I.Pos{1n}, y), I.Int.mul(I.Pos{1n}, I.Pos{1n}), I.Pos{1n},
          Equal.cong(I.Int, I.Int, z => I.Int.mul(I.Pos{1n}, z), y, I.Pos{1n}, h2),
          {==})))}

def gc.cop_mul(+a: Nat, +b: Nat, +c: Nat, hb: Cop<a, b>, hc: Cop<a, c>) -> Cop<a, Nat.mul(b, c)>:
  match hb:
    case Cop{u, v, h1}:
      match hc:
        case Cop{u2, v2, h2}:
          gc.cop_mul.go(a, b, c, u, v, u2, v2, h1, h2)

# a coprime to b is coprime to b squared
def gc.cop_sq(+a: Nat, +b: Nat, h: Cop<a, b>) -> Cop<a, Nat.mul(b, b)>:
  +h2 = h
  gc.cop_mul(a, b, b, h2, h2)

# ---------------------------------------------------------------------
# Euclid's algorithm, by subtraction, carrying Bezout's identity.
#
# Bez<a, b> is everything the algorithm produces at once: a common
# divisor g of a and b together with integers u, v making u a + v b = g.
# (That g is the GREATEST common divisor is never needed below -- what is
# needed is that the cofactors are coprime, and Bezout gives that.)
# ---------------------------------------------------------------------

type Bez<-a: Nat, -b: Nat> is Data:
  Bez{g: Nat, ga: Dvd<g, a>, gb: Dvd<g, b>, u: I.Int, v: I.Int,
      h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{g} : I.Int}}

def gc.bez_eq_r(-a: Nat, -b: Nat, -c: Nat, e: {b == c : Nat}, z: Bez<a, b>) -> Bez<a, c>:
  %e : Bez<a, _>
  z

def gc.bez_sym(+a: Nat, +b: Nat, z: Bez<a, b>) -> Bez<b, a>:
  match z:
    case Bez{g, ga, gb, u, v, h}:
      +u2 = u
      +v2 = v
      Bez{g, gb, ga, v2, u2,
        %I.add_comm(I.Int.mul(u2, I.Pos{a}), I.Int.mul(v2, I.Pos{b})) : {_ == I.Pos{g} : I.Int}
        h}

# (u - v) A + v (A + C) = u A + v C : the whole content of a subtraction step
def gc.bez_step(+u: I.Int, +v: I.Int, +xa: I.Int, +xc: I.Int) -> {I.Int.add(I.Int.mul(I.Int.sub(u, v), xa), I.Int.mul(v, I.Int.add(xa, xc))) == I.Int.add(I.Int.mul(u, xa), I.Int.mul(v, xc)) : I.Int}:
  %Equal.sym(I.Int, I.Int.sub(u, v), I.Int.add(u, I.Int.neg(v)), I.sub_def(u, v)) : {I.Int.add(I.Int.mul(_, xa), I.Int.mul(v, I.Int.add(xa, xc))) == I.Int.add(I.Int.mul(u, xa), I.Int.mul(v, xc)) : I.Int}
  %Equal.sym(I.Int, I.Int.mul(I.Int.add(u, I.Int.neg(v)), xa), I.Int.add(I.Int.mul(u, xa), I.Int.mul(I.Int.neg(v), xa)), gc.iadd_mul(u, I.Int.neg(v), xa)) : {I.Int.add(_, I.Int.mul(v, I.Int.add(xa, xc))) == I.Int.add(I.Int.mul(u, xa), I.Int.mul(v, xc)) : I.Int}
  %Equal.sym(I.Int, I.Int.mul(v, I.Int.add(xa, xc)), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc)), I.mul_add(v, xa, xc)) : {I.Int.add(I.Int.add(I.Int.mul(u, xa), I.Int.mul(I.Int.neg(v), xa)), _) == I.Int.add(I.Int.mul(u, xa), I.Int.mul(v, xc)) : I.Int}
  Equal.trans(I.Int,
    I.Int.add(I.Int.add(I.Int.mul(u, xa), I.Int.mul(I.Int.neg(v), xa)), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc))),
    I.Int.add(I.Int.mul(u, xa), I.Int.add(I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), I.Int.mul(v, xc))),
    I.Int.add(I.Int.mul(u, xa), I.Int.mul(v, xc)),
    Equal.trans(I.Int,
      I.Int.add(I.Int.add(I.Int.mul(u, xa), I.Int.mul(I.Int.neg(v), xa)), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc))),
      I.Int.add(I.Int.mul(u, xa), I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc)))),
      I.Int.add(I.Int.mul(u, xa), I.Int.add(I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), I.Int.mul(v, xc))),
      Equal.sym(I.Int,
        I.Int.add(I.Int.mul(u, xa), I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc)))),
        I.Int.add(I.Int.add(I.Int.mul(u, xa), I.Int.mul(I.Int.neg(v), xa)), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc))),
        I.add_assoc(I.Int.mul(u, xa), I.Int.mul(I.Int.neg(v), xa), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc)))),
      Equal.cong(I.Int, I.Int, z => I.Int.add(I.Int.mul(u, xa), z),
        I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.add(I.Int.mul(v, xa), I.Int.mul(v, xc))),
        I.Int.add(I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), I.Int.mul(v, xc)),
        I.add_assoc(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa), I.Int.mul(v, xc)))),
    Equal.cong(I.Int, I.Int, z => I.Int.add(I.Int.mul(u, xa), z),
      I.Int.add(I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), I.Int.mul(v, xc)),
      I.Int.mul(v, xc),
      Equal.trans(I.Int,
        I.Int.add(I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), I.Int.mul(v, xc)),
        I.Int.add(I.Pos{0n}, I.Int.mul(v, xc)),
        I.Int.mul(v, xc),
        Equal.cong(I.Int, I.Int, z => I.Int.add(z, I.Int.mul(v, xc)),
          I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), I.Pos{0n},
          Equal.trans(I.Int, I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), I.Int.mul(I.Int.add(I.Int.neg(v), v), xa), I.Pos{0n},
            Equal.sym(I.Int, I.Int.mul(I.Int.add(I.Int.neg(v), v), xa), I.Int.add(I.Int.mul(I.Int.neg(v), xa), I.Int.mul(v, xa)), gc.iadd_mul(I.Int.neg(v), v, xa)),
            Equal.trans(I.Int, I.Int.mul(I.Int.add(I.Int.neg(v), v), xa), I.Int.mul(I.Pos{0n}, xa), I.Pos{0n},
              Equal.cong(I.Int, I.Int, z => I.Int.mul(z, xa), I.Int.add(I.Int.neg(v), v), I.Pos{0n},
                Equal.trans(I.Int, I.Int.add(I.Int.neg(v), v), I.Int.add(v, I.Int.neg(v)), I.Pos{0n},
                  I.add_comm(I.Int.neg(v), v), I.add_neg(v))),
              Equal.trans(I.Int, I.Int.mul(I.Pos{0n}, xa), I.Int.mul(xa, I.Pos{0n}), I.Pos{0n},
                I.mul_comm(I.Pos{0n}, xa), I.mul_zero(xa))))),
        Equal.trans(I.Int, I.Int.add(I.Pos{0n}, I.Int.mul(v, xc)), I.Int.add(I.Int.mul(v, xc), I.Pos{0n}), I.Int.mul(v, xc),
          I.add_comm(I.Pos{0n}, I.Int.mul(v, xc)), I.add_zero(I.Int.mul(v, xc))))))

# from a common divisor of (a, c), one of (a, a + c)
def gc.bez_up(+a: Nat, +c: Nat, z: Bez<a, c>) -> Bez<a, Nat.add(a, c)>:
  match z:
    case Bez{g, ga, gb, u, v, h}:
      +ga2 = ga
      +u2 = u
      +v2 = v
      +g2 = g
      Bez{g2, ga2, gc.dvd_add(g2, a, c, ga2, gb), I.Int.sub(u2, v2), v2,
        %I.of_add(a, c) : {I.Int.add(I.Int.mul(I.Int.sub(u2, v2), I.Pos{a}), I.Int.mul(v2, _)) == I.Pos{g2} : I.Int}
        Equal.trans(I.Int,
          I.Int.add(I.Int.mul(I.Int.sub(u2, v2), I.Pos{a}), I.Int.mul(v2, I.Int.add(I.Pos{a}, I.Pos{c}))),
          I.Int.add(I.Int.mul(u2, I.Pos{a}), I.Int.mul(v2, I.Pos{c})), I.Pos{g2},
          gc.bez_step(u2, v2, I.Pos{a}, I.Pos{c}), h)}

# 1 a + 0 a = a
def gc.bez_refl(+a: Nat) -> Bez<a, a>:
  Bez{a, gc.dvd_refl(a), gc.dvd_refl(a), I.Pos{1n}, I.Pos{0n},
    Equal.cong(Nat, I.Int, z => I.Pos{z}, Nat.add(Nat.add(a, 0n), 0n), a,
      Equal.trans(Nat, Nat.add(Nat.add(a, 0n), 0n), Nat.add(a, 0n), a,
        N.Nat.add_zero(Nat.add(a, 0n)), N.Nat.add_zero(a)))}

# a < b gives b = (1+a) + d
def gc.split_lt.open(-ap: Nat, -bp: Nat, -P: Type, k: @d: Nat -> @e: {Nat.add(1n+ap, d) == bp : Nat} -> P, w: (&d: Nat -> {Nat.add(1n+ap, d) == bp : Nat})) -> P:
  (d, e) = w
  k(d, e)

def gc.split_lt(+ap: Nat, +bp: Nat, lt: N.Nat.Lt(ap, bp), -P: Type, k: @d: Nat -> @e: {Nat.add(1n+ap, d) == bp : Nat} -> P) -> P:
  gc.split_lt.open(ap, bp, P, k, N.Nat.le_diff(1n+ap, bp, lt))

# THE ALGORITHM. The fuel is ap + bp and every step strictly lowers it:
# with a < b the pair (a, b) becomes (a, b - a), and with b < a the pair is
# turned round first. Both non-trivial branches are one `bez_up`, so the
# Bezout bookkeeping is written once.
def gc.euclid.go(f: Nat, ap: Nat, bp: Nat, bound: N.Nat.Le(Nat.add(ap, bp), f)) -> Bez<1n+ap, 1n+bp>:
  match f:
    case 0n:
      match ap:
        case 0n:
          match bp:
            case 0n:
              gc.bez_refl(1n)
            case 1n+y:
              match bound:
        case 1n+x:
          match bound:
    case 1n+g:
      +g2 = g
      +ap2 = ap
      +bp2 = bp
      S2.Nat.le_copy(Nat.add(ap2, bp2), 1n+g2, bound, Bez<1n+ap2, 1n+bp2>, hb1 => hb2 =>
        S2.Nat.decide_lt(ap2, bp2, Bez<1n+ap2, 1n+bp2>,
          lt => gc.split_lt(ap2, bp2, lt, Bez<1n+ap2, 1n+bp2>, d => ed =>
            +d2 = d
            +ed2 = ed
            gc.bez_eq_r(1n+ap2, Nat.add(1n+ap2, 1n+d2), 1n+bp2,
              Equal.cong(Nat, Nat, z => 1n+z, Nat.add(ap2, 1n+d2), bp2,
                Equal.trans(Nat, Nat.add(ap2, 1n+d2), 1n+Nat.add(ap2, d2), bp2, N.Nat.add_succ(ap2, d2), ed2)),
              gc.bez_up(1n+ap2, 1n+d2,
                gc.euclid.go(g2, ap2, d2,
                  N.Nat.le_trans(1n+Nat.add(ap2, d2), Nat.add(ap2, bp2), 1n+g2,
                    %ed2 : N.Nat.Le(1n+Nat.add(ap2, d2), Nat.add(ap2, _))
                    %Equal.sym(Nat, Nat.add(ap2, 1n+Nat.add(ap2, d2)), 1n+Nat.add(ap2, Nat.add(ap2, d2)), N.Nat.add_succ(ap2, Nat.add(ap2, d2))) : N.Nat.Le(1n+Nat.add(ap2, d2), _)
                    N.Nat.le_add_l(ap2, Nat.add(ap2, d2)),
                    hb1))))),
          ge => S2.Nat.decide_lt(bp2, ap2, Bez<1n+ap2, 1n+bp2>,
            lt2 => gc.split_lt(bp2, ap2, lt2, Bez<1n+ap2, 1n+bp2>, d => ed =>
              +d2 = d
              +ed2 = ed
              gc.bez_sym(1n+bp2, 1n+ap2,
                gc.bez_eq_r(1n+bp2, Nat.add(1n+bp2, 1n+d2), 1n+ap2,
                  Equal.cong(Nat, Nat, z => 1n+z, Nat.add(bp2, 1n+d2), ap2,
                    Equal.trans(Nat, Nat.add(bp2, 1n+d2), 1n+Nat.add(bp2, d2), ap2, N.Nat.add_succ(bp2, d2), ed2)),
                  gc.bez_up(1n+bp2, 1n+d2,
                    gc.bez_sym(1n+d2, 1n+bp2,
                      gc.euclid.go(g2, d2, bp2,
                        N.Nat.le_trans(1n+Nat.add(d2, bp2), Nat.add(ap2, bp2), 1n+g2,
                          %ed2 : N.Nat.Le(1n+Nat.add(d2, bp2), Nat.add(_, bp2))
                          N.Nat.le_add_both(d2, Nat.add(bp2, d2), bp2, bp2, N.Nat.le_add_l(bp2, d2), N.Nat.le_refl(bp2)),
                          hb2))))))),
            ge2 => gc.bez_eq_r(1n+ap2, 1n+ap2, 1n+bp2,
              Equal.cong(Nat, Nat, z => 1n+z, ap2, bp2, N.Nat.le_antisym(ap2, bp2, ge2, ge)),
              gc.bez_refl(1n+ap2)))))

def gc.euclid(ap: Nat, bp: Nat) -> Bez<1n+ap, 1n+bp>:
  +ap2 = ap
  +bp2 = bp
  gc.euclid.go(Nat.add(ap2, bp2), ap2, bp2, N.Nat.le_refl(Nat.add(ap2, bp2)))

# ---------------------------------------------------------------------
# Nat algebra the square theorem needs (the same five identities as the
# Int toolkit above, on Nat, plus cancellation of a positive factor)
# ---------------------------------------------------------------------

def gc.nmul_swap(+x: Nat, +y: Nat, +z: Nat) -> {Nat.mul(x, Nat.mul(y, z)) == Nat.mul(y, Nat.mul(x, z)) : Nat}:
  %Equal.sym(Nat, Nat.mul(x, Nat.mul(y, z)), Nat.mul(Nat.mul(x, y), z), N.Nat.mul_assoc(x, y, z)) : {_ == Nat.mul(y, Nat.mul(x, z)) : Nat}
  %Equal.sym(Nat, Nat.mul(y, Nat.mul(x, z)), Nat.mul(Nat.mul(y, x), z), N.Nat.mul_assoc(y, x, z)) : {Nat.mul(Nat.mul(x, y), z) == _ : Nat}
  Equal.cong(Nat, Nat, w => Nat.mul(w, z), Nat.mul(x, y), Nat.mul(y, x), N.Nat.mul_comm(x, y))

def gc.nmul_shift(+x: Nat, +y: Nat, +z: Nat) -> {Nat.mul(Nat.mul(x, y), z) == Nat.mul(Nat.mul(x, z), y) : Nat}:
  %N.Nat.mul_assoc(x, z, y) : {Nat.mul(Nat.mul(x, y), z) == _ : Nat}
  %N.Nat.mul_assoc(x, y, z) : {_ == Nat.mul(x, Nat.mul(z, y)) : Nat}
  Equal.cong(Nat, Nat, w => Nat.mul(x, w), Nat.mul(y, z), Nat.mul(z, y), N.Nat.mul_comm(y, z))

def gc.nsq_prod(+a: Nat, +b: Nat) -> {Nat.mul(Nat.mul(a, b), Nat.mul(a, b)) == Nat.mul(Nat.mul(a, a), Nat.mul(b, b)) : Nat}:
  %N.Nat.mul_assoc(a, b, Nat.mul(a, b)) : {_ == Nat.mul(Nat.mul(a, a), Nat.mul(b, b)) : Nat}
  %N.Nat.mul_assoc(a, a, Nat.mul(b, b)) : {Nat.mul(a, Nat.mul(b, Nat.mul(a, b))) == _ : Nat}
  Equal.cong(Nat, Nat, z => Nat.mul(a, z), Nat.mul(b, Nat.mul(a, b)), Nat.mul(a, Nat.mul(b, b)), gc.nmul_swap(b, a, b))

def gc.ncancel_r(+m: Nat, x: Nat, y: Nat, e: {Nat.mul(x, 1n+m) == Nat.mul(y, 1n+m) : Nat}) -> {x == y : Nat}:
  match x:
    case 0n:
      match y:
        case 0n:
          {==}
        case 1n+yp:
          Empty.absurd({0n == 1n+yp : Nat}, N.Nat.zero_neq_succ(Nat.add(m, Nat.mul(yp, 1n+m)), e))
    case 1n+xp:
      match y:
        case 0n:
          Empty.absurd({1n+xp == 0n : Nat}, N.Nat.succ_neq_zero(Nat.add(m, Nat.mul(xp, 1n+m)), e))
        case 1n+yp:
          Equal.cong(Nat, Nat, z => 1n+z, xp, yp,
            gc.ncancel_r(m, xp, yp,
              N.Nat.add_cancel_l(m, Nat.mul(xp, 1n+m), Nat.mul(yp, 1n+m),
                N.Nat.succ_inj(Nat.add(m, Nat.mul(xp, 1n+m)), Nat.add(m, Nat.mul(yp, 1n+m)), e))))

def gc.ncancel(+m: Nat, +x: Nat, +y: Nat, e: {Nat.mul(1n+m, x) == Nat.mul(1n+m, y) : Nat}) -> {x == y : Nat}:
  gc.ncancel_r(m, x, y,
    Equal.trans(Nat, Nat.mul(x, 1n+m), Nat.mul(1n+m, x), Nat.mul(y, 1n+m),
      N.Nat.mul_comm(x, 1n+m),
      Equal.trans(Nat, Nat.mul(1n+m, x), Nat.mul(1n+m, y), Nat.mul(y, 1n+m), e, N.Nat.mul_comm(1n+m, y))))

# a factor of a positive product is positive
def gc.pos_of_mul(g: Nat, +k: Nat, +ap: Nat, e: {Nat.mul(g, k) == 1n+ap : Nat}) -> N.Nat.Lt(0n, g):
  match g:
    case 0n:
      Empty.absurd(N.Nat.Lt(0n, 0n), N.Nat.zero_neq_succ(ap, e))
    case 1n+x:
      Unit{}

def gc.pos_of_mul_r(+g: Nat, +k: Nat, +ap: Nat, e: {Nat.mul(g, k) == 1n+ap : Nat}) -> N.Nat.Lt(0n, k):
  gc.pos_of_mul(k, g, ap, Equal.trans(Nat, Nat.mul(k, g), Nat.mul(g, k), 1n+ap, N.Nat.mul_comm(k, g), e))

def gc.pos_shape(g: Nat, p: N.Nat.Lt(0n, g), -P: Type, k: @gp: Nat -> @e: {g == 1n+gp : Nat} -> P) -> P:
  match g:
    case 0n:
      match p:
    case 1n+gp:
      k(gp, {==})

# ---------------------------------------------------------------------
# cancelling a positive factor in Int (no zero divisors), and the
# coprime cofactors of a gcd
# ---------------------------------------------------------------------

def gc.pnat(z: I.Int) -> Nat:
  match z:
    case I.Pos{n}:
      n
    case I.NegS{n}:
      0n

def gc.pos_ne_zero(+gp: Nat, e: {I.Pos{1n+gp} == I.Pos{0n} : I.Int}) -> Empty:
  N.Nat.succ_neq_zero(gp, Equal.cong(I.Int, Nat, z => gc.pnat(z), I.Pos{1n+gp}, I.Pos{0n}, e))

def gc.ineg_add(+y: I.Int) -> {I.Int.add(I.Int.neg(y), y) == I.Pos{0n} : I.Int}:
  Equal.trans(I.Int, I.Int.add(I.Int.neg(y), y), I.Int.add(y, I.Int.neg(y)), I.Pos{0n},
    I.add_comm(I.Int.neg(y), y), I.add_neg(y))

def gc.isub_zero(+x: I.Int, +y: I.Int, e: {I.Int.add(x, I.Int.neg(y)) == I.Pos{0n} : I.Int}) -> {x == y : I.Int}:
  Equal.trans(I.Int, x, I.Int.add(x, I.Pos{0n}), y,
    Equal.sym(I.Int, I.Int.add(x, I.Pos{0n}), x, I.add_zero(x)),
    Equal.trans(I.Int, I.Int.add(x, I.Pos{0n}), I.Int.add(x, I.Int.add(I.Int.neg(y), y)), y,
      Equal.cong(I.Int, I.Int, z => I.Int.add(x, z), I.Pos{0n}, I.Int.add(I.Int.neg(y), y),
        Equal.sym(I.Int, I.Int.add(I.Int.neg(y), y), I.Pos{0n}, gc.ineg_add(y))),
      Equal.trans(I.Int, I.Int.add(x, I.Int.add(I.Int.neg(y), y)), I.Int.add(I.Int.add(x, I.Int.neg(y)), y), y,
        I.add_assoc(x, I.Int.neg(y), y),
        Equal.trans(I.Int, I.Int.add(I.Int.add(x, I.Int.neg(y)), y), I.Int.add(I.Pos{0n}, y), y,
          Equal.cong(I.Int, I.Int, z => I.Int.add(z, y), I.Int.add(x, I.Int.neg(y)), I.Pos{0n}, e),
          Equal.trans(I.Int, I.Int.add(I.Pos{0n}, y), I.Int.add(y, I.Pos{0n}), y,
            I.add_comm(I.Pos{0n}, y), I.add_zero(y))))))

def gc.icancel(+gp: Nat, +x: I.Int, +y: I.Int, e: {I.Int.mul(I.Pos{1n+gp}, x) == I.Int.mul(I.Pos{1n+gp}, y) : I.Int}) -> {x == y : I.Int}:
  gc.isub_zero(x, y,
    I.mul_eq_zero(I.Pos{1n+gp}, I.Int.add(x, I.Int.neg(y)),
      Equal.trans(I.Int, I.Int.mul(I.Pos{1n+gp}, I.Int.add(x, I.Int.neg(y))), I.Int.add(I.Int.mul(I.Pos{1n+gp}, x), I.Int.mul(I.Pos{1n+gp}, I.Int.neg(y))), I.Pos{0n},
        I.mul_add(I.Pos{1n+gp}, x, I.Int.neg(y)),
        Equal.trans(I.Int, I.Int.add(I.Int.mul(I.Pos{1n+gp}, x), I.Int.mul(I.Pos{1n+gp}, I.Int.neg(y))), I.Int.add(I.Int.mul(I.Pos{1n+gp}, y), I.Int.mul(I.Pos{1n+gp}, I.Int.neg(y))), I.Pos{0n},
          Equal.cong(I.Int, I.Int, z => I.Int.add(z, I.Int.mul(I.Pos{1n+gp}, I.Int.neg(y))), I.Int.mul(I.Pos{1n+gp}, x), I.Int.mul(I.Pos{1n+gp}, y), e),
          Equal.trans(I.Int, I.Int.add(I.Int.mul(I.Pos{1n+gp}, y), I.Int.mul(I.Pos{1n+gp}, I.Int.neg(y))), I.Int.mul(I.Pos{1n+gp}, I.Int.add(y, I.Int.neg(y))), I.Pos{0n},
            Equal.sym(I.Int, I.Int.mul(I.Pos{1n+gp}, I.Int.add(y, I.Int.neg(y))), I.Int.add(I.Int.mul(I.Pos{1n+gp}, y), I.Int.mul(I.Pos{1n+gp}, I.Int.neg(y))), I.mul_add(I.Pos{1n+gp}, y, I.Int.neg(y))),
            Equal.trans(I.Int, I.Int.mul(I.Pos{1n+gp}, I.Int.add(y, I.Int.neg(y))), I.Int.mul(I.Pos{1n+gp}, I.Pos{0n}), I.Pos{0n},
              Equal.cong(I.Int, I.Int, z => I.Int.mul(I.Pos{1n+gp}, z), I.Int.add(y, I.Int.neg(y)), I.Pos{0n}, I.add_neg(y)),
              I.mul_zero(I.Pos{1n+gp}))))),
      ne => gc.pos_ne_zero(gp, ne)))

# w (u s + v t) = u (w s) + v (w t)
def gc.dist2(+a: Nat, +b: Nat, +w: Nat, +s: Nat, +t: Nat, es: {Nat.mul(w, s) == a : Nat}, et: {Nat.mul(w, t) == b : Nat}, +u: I.Int, +v: I.Int) -> {I.Int.mul(I.Pos{w}, I.Int.add(I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t}))) == I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) : I.Int}:
  %Equal.trans(I.Int, I.Int.mul(I.Pos{w}, I.Pos{s}), I.Pos{Nat.mul(w, s)}, I.Pos{a}, I.of_mul(w, s), Equal.cong(Nat, I.Int, z => I.Pos{z}, Nat.mul(w, s), a, es)) : {I.Int.mul(I.Pos{w}, I.Int.add(I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t}))) == I.Int.add(I.Int.mul(u, _), I.Int.mul(v, I.Pos{b})) : I.Int}
  %Equal.trans(I.Int, I.Int.mul(I.Pos{w}, I.Pos{t}), I.Pos{Nat.mul(w, t)}, I.Pos{b}, I.of_mul(w, t), Equal.cong(Nat, I.Int, z => I.Pos{z}, Nat.mul(w, t), b, et)) : {I.Int.mul(I.Pos{w}, I.Int.add(I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t}))) == I.Int.add(I.Int.mul(u, I.Int.mul(I.Pos{w}, I.Pos{s})), I.Int.mul(v, _)) : I.Int}
  %gc.imul_swap(I.Pos{w}, u, I.Pos{s}) : {I.Int.mul(I.Pos{w}, I.Int.add(I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t}))) == I.Int.add(_, I.Int.mul(v, I.Int.mul(I.Pos{w}, I.Pos{t}))) : I.Int}
  %gc.imul_swap(I.Pos{w}, v, I.Pos{t}) : {I.Int.mul(I.Pos{w}, I.Int.add(I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t}))) == I.Int.add(I.Int.mul(I.Pos{w}, I.Int.mul(u, I.Pos{s})), _) : I.Int}
  I.mul_add(I.Pos{w}, I.Int.mul(u, I.Pos{s}), I.Int.mul(v, I.Pos{t}))

# Split<a, b>: a common divisor 1+gp of a and b whose cofactors are COPRIME.
# That is the only property of a greatest common divisor this file uses.
type Split<-a: Nat, -b: Nat> is Data:
  Split{gp: Nat, a1: Nat, b1: Nat,
        ea: {Nat.mul(1n+gp, a1) == a : Nat},
        eb: {Nat.mul(1n+gp, b1) == b : Nat},
        hc: Cop<a1, b1>}

def gc.split.go(+a: Nat, +b: Nat, +gp: Nat, +k1: Nat, +k2: Nat, e1: {Nat.mul(1n+gp, k1) == a : Nat}, e2: {Nat.mul(1n+gp, k2) == b : Nat}, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{1n+gp} : I.Int}) -> Split<a, b>:
  +e1c = e1
  +e2c = e2
  Split{gp, k1, k2, e1c, e2c,
    Cop{u, v,
      gc.icancel(gp, I.Int.add(I.Int.mul(u, I.Pos{k1}), I.Int.mul(v, I.Pos{k2})), I.Pos{1n},
        Equal.trans(I.Int, I.Int.mul(I.Pos{1n+gp}, I.Int.add(I.Int.mul(u, I.Pos{k1}), I.Int.mul(v, I.Pos{k2}))), I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})), I.Int.mul(I.Pos{1n+gp}, I.Pos{1n}),
          gc.dist2(a, b, 1n+gp, k1, k2, e1c, e2c, u, v),
          Equal.trans(I.Int, I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})), I.Pos{1n+gp}, I.Int.mul(I.Pos{1n+gp}, I.Pos{1n}),
            h, Equal.sym(I.Int, I.Int.mul(I.Pos{1n+gp}, I.Pos{1n}), I.Pos{1n+gp}, I.mul_one(I.Pos{1n+gp})))))}}

def gc.split.mk(+a: Nat, +b: Nat, g: Nat, k1: Nat, k2: Nat, e1: {Nat.mul(g, k1) == a : Nat}, e2: {Nat.mul(g, k2) == b : Nat}, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b})) == I.Pos{g} : I.Int}, pos: N.Nat.Lt(0n, g)) -> Split<a, b>:
  match g:
    case 0n:
      match pos:
    case 1n+gp:
      gc.split.go(a, b, gp, k1, k2, e1, e2, u, v, h)

def gc.split_of.fin(+ap: Nat, +b: Nat, +g: Nat, ga: Dvd<g, 1n+ap>, gb: Dvd<g, b>, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{1n+ap}), I.Int.mul(v, I.Pos{b})) == I.Pos{g} : I.Int}) -> Split<1n+ap, b>:
  match ga:
    case Dvd{k1, e1}:
      match gb:
        case Dvd{k2, e2}:
          +e1c = e1
          +k1c = k1
          gc.split.mk(1n+ap, b, g, k1c, k2, e1c, e2, u, v, h, gc.pos_of_mul(g, k1c, ap, e1c))

def gc.split_of(+ap: Nat, +b: Nat, z: Bez<1n+ap, b>) -> Split<1n+ap, b>:
  match z:
    case Bez{g, ga, gb, u, v, h}:
      gc.split_of.fin(ap, b, g, ga, gb, u, v, h)

def gc.split(ap: Nat, bp: Nat) -> Split<1n+ap, 1n+bp>:
  +ap2 = ap
  +bp2 = bp
  gc.split_of(ap2, 1n+bp2, gc.euclid(ap2, bp2))

# ---------------------------------------------------------------------
# THE SQUARE THEOREM: coprime factors of a square are squares.
#
#   a b = c^2 with gcd(a, b) = 1  implies  a is a square
#
# Proof, with G = gcd(a, c), a = G a1, c = G c1 and gcd(a1, c1) = 1:
# cancelling G in a b = c^2 gives a1 b = G c1^2, so a1 divides G c1^2 and,
# being coprime to c1 hence to c1^2, divides G (Gauss). Write G = a1 w.
# Then a = G a1 = a1^2 w and, cancelling a1 in a1 b = G c1^2, b = w c1^2.
# So w divides both a and b, which are coprime: w = 1 and a = a1^2.
# ---------------------------------------------------------------------

type IsSq<-n: Nat> is Data:
  IsSq{s: Nat, h: {Nat.mul(s, s) == n : Nat}}

def gc.ncancel_pos(+x: Nat, +y: Nat, m: Nat, pos: N.Nat.Lt(0n, m), e: {Nat.mul(m, x) == Nat.mul(m, y) : Nat}) -> {x == y : Nat}:
  match m:
    case 0n:
      match pos:
    case 1n+mp:
      gc.ncancel(mp, x, y, e)

def gc.sq_fin(+av: Nat, +b: Nat, +a1: Nat, +c1: Nat, +gp: Nat, +w: Nat, +ea: {Nat.mul(1n+gp, a1) == 1n+av : Nat}, ew: {Nat.mul(a1, w) == 1n+gp : Nat}, key1: {Nat.mul(a1, b) == Nat.mul(1n+gp, Nat.mul(c1, c1)) : Nat}, hc: Cop<1n+av, b>) -> IsSq<1n+av>:
  +ew2 = ew
  +key2 = Equal.trans(Nat, Nat.mul(Nat.mul(a1, a1), w), Nat.mul(1n+gp, a1), 1n+av,
    Equal.trans(Nat, Nat.mul(Nat.mul(a1, a1), w), Nat.mul(Nat.mul(a1, w), a1), Nat.mul(1n+gp, a1),
      Equal.sym(Nat, Nat.mul(Nat.mul(a1, w), a1), Nat.mul(Nat.mul(a1, a1), w), gc.nmul_shift(a1, w, a1)),
      Equal.cong(Nat, Nat, z => Nat.mul(z, a1), Nat.mul(a1, w), 1n+gp, ew2)),
    ea)
  +key3 = gc.ncancel_pos(b, Nat.mul(w, Nat.mul(c1, c1)), a1, gc.pos_of_mul_r(1n+gp, a1, av, ea),
    Equal.trans(Nat, Nat.mul(a1, b), Nat.mul(1n+gp, Nat.mul(c1, c1)), Nat.mul(a1, Nat.mul(w, Nat.mul(c1, c1))),
      key1,
      Equal.trans(Nat, Nat.mul(1n+gp, Nat.mul(c1, c1)), Nat.mul(Nat.mul(a1, w), Nat.mul(c1, c1)), Nat.mul(a1, Nat.mul(w, Nat.mul(c1, c1))),
        Equal.cong(Nat, Nat, z => Nat.mul(z, Nat.mul(c1, c1)), 1n+gp, Nat.mul(a1, w), Equal.sym(Nat, Nat.mul(a1, w), 1n+gp, ew2)),
        Equal.sym(Nat, Nat.mul(a1, Nat.mul(w, Nat.mul(c1, c1))), Nat.mul(Nat.mul(a1, w), Nat.mul(c1, c1)), N.Nat.mul_assoc(a1, w, Nat.mul(c1, c1))))))
  +w1 = gc.cop_common(1n+av, b, w, hc,
    Dvd{Nat.mul(a1, a1), Equal.trans(Nat, Nat.mul(w, Nat.mul(a1, a1)), Nat.mul(Nat.mul(a1, a1), w), 1n+av, N.Nat.mul_comm(w, Nat.mul(a1, a1)), key2)},
    Dvd{Nat.mul(c1, c1), Equal.sym(Nat, b, Nat.mul(w, Nat.mul(c1, c1)), key3)})
  IsSq{a1,
    Equal.trans(Nat, Nat.mul(a1, a1), Nat.mul(Nat.mul(a1, a1), w), 1n+av,
      Equal.trans(Nat, Nat.mul(a1, a1), Nat.mul(Nat.mul(a1, a1), 1n), Nat.mul(Nat.mul(a1, a1), w),
        Equal.sym(Nat, Nat.mul(Nat.mul(a1, a1), 1n), Nat.mul(a1, a1), N.Nat.mul_one(Nat.mul(a1, a1))),
        Equal.cong(Nat, Nat, z => Nat.mul(Nat.mul(a1, a1), z), 1n, w, Equal.sym(Nat, w, 1n, w1))),
      key2)}

def gc.sq_open(+av: Nat, +b: Nat, +a1: Nat, +c1: Nat, +gp: Nat, ea: {Nat.mul(1n+gp, a1) == 1n+av : Nat}, key1: {Nat.mul(a1, b) == Nat.mul(1n+gp, Nat.mul(c1, c1)) : Nat}, hc: Cop<1n+av, b>, dw: Dvd<a1, 1n+gp>) -> IsSq<1n+av>:
  match dw:
    case Dvd{w, ew}:
      gc.sq_fin(av, b, a1, c1, gp, w, ea, ew, key1, hc)

def gc.sq_main2(+av: Nat, +b: Nat, +cv: Nat, +gp: Nat, +a1: Nat, +c1: Nat, +ea: {Nat.mul(1n+gp, a1) == 1n+av : Nat}, ec: {Nat.mul(1n+gp, c1) == 1n+cv : Nat}, hcop: Cop<a1, c1>, hc: Cop<1n+av, b>, e: {Nat.mul(1n+av, b) == Nat.mul(1n+cv, 1n+cv) : Nat}) -> IsSq<1n+av>:
  +key1 = gc.ncancel(gp, Nat.mul(a1, b), Nat.mul(1n+gp, Nat.mul(c1, c1)),
    Equal.trans(Nat, Nat.mul(1n+gp, Nat.mul(a1, b)), Nat.mul(Nat.mul(1n+gp, a1), b), Nat.mul(1n+gp, Nat.mul(1n+gp, Nat.mul(c1, c1))),
      N.Nat.mul_assoc(1n+gp, a1, b),
      Equal.trans(Nat, Nat.mul(Nat.mul(1n+gp, a1), b), Nat.mul(1n+av, b), Nat.mul(1n+gp, Nat.mul(1n+gp, Nat.mul(c1, c1))),
        Equal.cong(Nat, Nat, z => Nat.mul(z, b), Nat.mul(1n+gp, a1), 1n+av, ea),
        Equal.trans(Nat, Nat.mul(1n+av, b), Nat.mul(1n+cv, 1n+cv), Nat.mul(1n+gp, Nat.mul(1n+gp, Nat.mul(c1, c1))),
          e,
          Equal.trans(Nat, Nat.mul(1n+cv, 1n+cv), Nat.mul(Nat.mul(1n+gp, c1), Nat.mul(1n+gp, c1)), Nat.mul(1n+gp, Nat.mul(1n+gp, Nat.mul(c1, c1))),
            Equal.cong(Nat, Nat, z => Nat.mul(z, z), 1n+cv, Nat.mul(1n+gp, c1), Equal.sym(Nat, Nat.mul(1n+gp, c1), 1n+cv, ec)),
            Equal.trans(Nat, Nat.mul(Nat.mul(1n+gp, c1), Nat.mul(1n+gp, c1)), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(c1, c1)), Nat.mul(1n+gp, Nat.mul(1n+gp, Nat.mul(c1, c1))),
              gc.nsq_prod(1n+gp, c1),
              Equal.sym(Nat, Nat.mul(1n+gp, Nat.mul(1n+gp, Nat.mul(c1, c1))), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(c1, c1)), N.Nat.mul_assoc(1n+gp, 1n+gp, Nat.mul(c1, c1)))))))))
  gc.sq_open(av, b, a1, c1, gp, ea, key1, hc,
    gc.gauss(a1, Nat.mul(c1, c1), 1n+gp, gc.cop_sq(a1, c1, hcop),
      Dvd{b, Equal.trans(Nat, Nat.mul(a1, b), Nat.mul(1n+gp, Nat.mul(c1, c1)), Nat.mul(Nat.mul(c1, c1), 1n+gp),
        key1, N.Nat.mul_comm(1n+gp, Nat.mul(c1, c1)))}))

def gc.sq_main(+av: Nat, +b: Nat, +cv: Nat, hc: Cop<1n+av, b>, e: {Nat.mul(1n+av, b) == Nat.mul(1n+cv, 1n+cv) : Nat}, sp: Split<1n+av, 1n+cv>) -> IsSq<1n+av>:
  match sp:
    case Split{gp, a1, c1, ea, ec, hcop}:
      gc.sq_main2(av, b, cv, gp, a1, c1, ea, ec, hcop, hc, e)

# a b = c^2 with a, b coprime and a positive makes a a square
def gc.sq_of_cop(+av: Nat, +b: Nat, c: Nat, hc: Cop<1n+av, b>, e: {Nat.mul(1n+av, b) == Nat.mul(c, c) : Nat}) -> IsSq<1n+av>:
  match c:
    case 0n:
      IsSq{1n,
        Equal.sym(Nat, 1n+av, 1n,
          gc.cop_common(1n+av, 0n, 1n+av,
            gc.cop_eq_r(1n+av, b, 0n,
              gc.mul_zero_l(b, av, Equal.trans(Nat, Nat.mul(b, 1n+av), Nat.mul(1n+av, b), 0n, N.Nat.mul_comm(b, 1n+av), e)),
              hc),
            gc.dvd_refl(1n+av), gc.dvd_zero(1n+av)))}
    case 1n+cv:
      +av2 = av
      +cv2 = cv
      gc.sq_main(av2, b, cv2, hc, e, gc.split(av2, cv2))

# ---------------------------------------------------------------------
# order and subtraction of divisibility
# ---------------------------------------------------------------------

def gc.le_mul_l(+c: Nat, +x: Nat, +y: Nat, h: N.Nat.Le(x, y)) -> N.Nat.Le(Nat.mul(c, x), Nat.mul(c, y)):
  match c:
    case 0n:
      Unit{}
    case 1n+p:
      S2.Nat.le_copy(x, y, h, N.Nat.Le(Nat.mul(1n+p, x), Nat.mul(1n+p, y)), h1 => h2 =>
        N.Nat.le_add_both(x, y, Nat.mul(p, x), Nat.mul(p, y), h1, gc.le_mul_l(p, x, y, h2)))

def gc.le_mul_cancel_l.no(+dp: Nat, +a: Nat, +b: Nat, h: N.Nat.Le(Nat.mul(1n+dp, a), Nat.mul(1n+dp, b)), lt: N.Nat.Lt(b, a)) -> Empty:
  N.Nat.lt_irrefl(Nat.mul(1n+dp, b),
    N.Nat.le_trans(1n+Nat.mul(1n+dp, b), Nat.mul(1n+dp, a), Nat.mul(1n+dp, b),
      N.Nat.le_trans(1n+Nat.mul(1n+dp, b), Nat.mul(1n+dp, 1n+b), Nat.mul(1n+dp, a),
        %Equal.sym(Nat, Nat.mul(1n+dp, 1n+b), Nat.add(1n+dp, Nat.mul(1n+dp, b)), N.Nat.mul_succ(1n+dp, b)) : N.Nat.Le(1n+Nat.mul(1n+dp, b), _)
        N.Nat.le_add_l(dp, Nat.mul(1n+dp, b)),
        gc.le_mul_l(1n+dp, 1n+b, a, lt)),
      h))

def gc.le_mul_cancel_l(+dp: Nat, +a: Nat, +b: Nat, h: N.Nat.Le(Nat.mul(1n+dp, a), Nat.mul(1n+dp, b))) -> N.Nat.Le(a, b):
  S2.Nat.decide_lt(b, a, N.Nat.Le(a, b),
    lt => Empty.absurd(N.Nat.Le(a, b), gc.le_mul_cancel_l.no(dp, a, b, h, lt)),
    ge => ge)

# d | x and d | x + y give d | y
def gc.dvd_sub.fin(+dp: Nat, +x: Nat, +y: Nat, +s: Nat, +bb: Nat, +aa: Nat, eB: {Nat.mul(1n+dp, bb) == x : Nat}, eA: {Nat.mul(1n+dp, aa) == s : Nat}, e: {Nat.add(x, y) == s : Nat}, +t: Nat, et: {Nat.add(bb, t) == aa : Nat}) -> Dvd<1n+dp, y>:
  Dvd{t,
    N.Nat.add_cancel_l(Nat.mul(1n+dp, bb), Nat.mul(1n+dp, t), y,
      Equal.trans(Nat, Nat.add(Nat.mul(1n+dp, bb), Nat.mul(1n+dp, t)), Nat.mul(1n+dp, Nat.add(bb, t)), Nat.add(Nat.mul(1n+dp, bb), y),
        Equal.sym(Nat, Nat.mul(1n+dp, Nat.add(bb, t)), Nat.add(Nat.mul(1n+dp, bb), Nat.mul(1n+dp, t)), N.Nat.mul_add(1n+dp, bb, t)),
        Equal.trans(Nat, Nat.mul(1n+dp, Nat.add(bb, t)), Nat.mul(1n+dp, aa), Nat.add(Nat.mul(1n+dp, bb), y),
          Equal.cong(Nat, Nat, z => Nat.mul(1n+dp, z), Nat.add(bb, t), aa, et),
          Equal.trans(Nat, Nat.mul(1n+dp, aa), Nat.add(x, y), Nat.add(Nat.mul(1n+dp, bb), y),
            Equal.trans(Nat, Nat.mul(1n+dp, aa), s, Nat.add(x, y), eA, Equal.sym(Nat, Nat.add(x, y), s, e)),
            Equal.cong(Nat, Nat, z => Nat.add(z, y), x, Nat.mul(1n+dp, bb), Equal.sym(Nat, Nat.mul(1n+dp, bb), x, eB))))))}

def gc.dvd_sub.open(+dp: Nat, +x: Nat, +y: Nat, +s: Nat, +bb: Nat, +aa: Nat, eB: {Nat.mul(1n+dp, bb) == x : Nat}, eA: {Nat.mul(1n+dp, aa) == s : Nat}, e: {Nat.add(x, y) == s : Nat}, w: (&t: Nat -> {Nat.add(bb, t) == aa : Nat})) -> Dvd<1n+dp, y>:
  (t, et) = w
  gc.dvd_sub.fin(dp, x, y, s, bb, aa, eB, eA, e, t, et)

def gc.dvd_sub.pos(+dp: Nat, +x: Nat, +y: Nat, +s: Nat, +bb: Nat, +aa: Nat, +eB: {Nat.mul(1n+dp, bb) == x : Nat}, +eA: {Nat.mul(1n+dp, aa) == s : Nat}, +e: {Nat.add(x, y) == s : Nat}) -> Dvd<1n+dp, y>:
  gc.dvd_sub.open(dp, x, y, s, bb, aa, eB, eA, e,
    N.Nat.le_diff(bb, aa,
      gc.le_mul_cancel_l(dp, bb, aa,
        %Equal.trans(Nat, Nat.add(Nat.mul(1n+dp, bb), y), Nat.add(x, y), Nat.mul(1n+dp, aa),
          Equal.cong(Nat, Nat, z => Nat.add(z, y), Nat.mul(1n+dp, bb), x, eB),
          Equal.trans(Nat, Nat.add(x, y), s, Nat.mul(1n+dp, aa), e, Equal.sym(Nat, Nat.mul(1n+dp, aa), s, eA))) : N.Nat.Le(Nat.mul(1n+dp, bb), _)
        N.Nat.le_add_r(Nat.mul(1n+dp, bb), y))))

def gc.dvd_sub.mk(d: Nat, +x: Nat, +y: Nat, +s: Nat, +bb: Nat, +aa: Nat, +eB: {Nat.mul(d, bb) == x : Nat}, +eA: {Nat.mul(d, aa) == s : Nat}, +e: {Nat.add(x, y) == s : Nat}) -> Dvd<d, y>:
  match d:
    case 0n:
      Dvd{0n,
        Equal.sym(Nat, y, 0n,
          N.Nat.add_eq_zero_r(x, y,
            Equal.trans(Nat, Nat.add(x, y), s, 0n, e, Equal.sym(Nat, Nat.mul(0n, aa), s, eA))))}
    case 1n+dp:
      gc.dvd_sub.pos(dp, x, y, s, bb, aa, eB, eA, e)

def gc.dvd_sub(+d: Nat, +x: Nat, +y: Nat, +s: Nat, +e: {Nat.add(x, y) == s : Nat}, d1: Dvd<d, x>, d2: Dvd<d, s>) -> Dvd<d, y>:
  match d1:
    case Dvd{bb, eB}:
      match d2:
        case Dvd{aa, eA}:
          gc.dvd_sub.mk(d, x, y, s, bb, aa, eB, eA, e)

# ---------------------------------------------------------------------
# introducing coprimality from "every common divisor is 1" -- the bridge
# that lets every coprimality fact below be a divisibility argument
# ---------------------------------------------------------------------

def gc.cop_intro.go(+up: Nat, +vp: Nat, k: @g: Nat -> @d1: Dvd<g, 1n+up> -> @d2: Dvd<g, 1n+vp> -> {g == 1n : Nat}, gp: Nat, u1: Nat, v1: Nat, ea: {Nat.mul(1n+gp, u1) == 1n+up : Nat}, eb: {Nat.mul(1n+gp, v1) == 1n+vp : Nat}, hcop: Cop<u1, v1>) -> Cop<1n+up, 1n+vp>:
  +ea2 = ea
  +eb2 = eb
  +u12 = u1
  +v12 = v1
  +e0 = N.Nat.succ_inj(gp, 0n, k(1n+gp, Dvd{u12, ea2}, Dvd{v12, eb2}))
  gc.cop_eq_l(u12, 1n+up, 1n+vp,
    Equal.trans(Nat, u12, Nat.mul(1n, u12), 1n+up,
      Equal.sym(Nat, Nat.mul(1n, u12), u12, N.Nat.one_mul(u12)),
      Equal.trans(Nat, Nat.mul(1n, u12), Nat.mul(1n+gp, u12), 1n+up,
        Equal.cong(Nat, Nat, z => Nat.mul(1n+z, u12), 0n, gp, Equal.sym(Nat, gp, 0n, e0)),
        ea2)),
    gc.cop_eq_r(u12, v12, 1n+vp,
      Equal.trans(Nat, v12, Nat.mul(1n, v12), 1n+vp,
        Equal.sym(Nat, Nat.mul(1n, v12), v12, N.Nat.one_mul(v12)),
        Equal.trans(Nat, Nat.mul(1n, v12), Nat.mul(1n+gp, v12), 1n+vp,
          Equal.cong(Nat, Nat, z => Nat.mul(1n+z, v12), 0n, gp, Equal.sym(Nat, gp, 0n, e0)),
          eb2)),
      hcop))

def gc.cop_intro.op(+up: Nat, +vp: Nat, k: @g: Nat -> @d1: Dvd<g, 1n+up> -> @d2: Dvd<g, 1n+vp> -> {g == 1n : Nat}, sp: Split<1n+up, 1n+vp>) -> Cop<1n+up, 1n+vp>:
  match sp:
    case Split{gp, u1, v1, ea, eb, hcop}:
      gc.cop_intro.go(up, vp, k, gp, u1, v1, ea, eb, hcop)

def gc.cop_intro(+up: Nat, +vp: Nat, k: @g: Nat -> @d1: Dvd<g, 1n+up> -> @d2: Dvd<g, 1n+vp> -> {g == 1n : Nat}) -> Cop<1n+up, 1n+vp>:
  gc.cop_intro.op(up, vp, k, gc.split(up, vp))

# ---------------------------------------------------------------------
# squares and coprimality
# ---------------------------------------------------------------------

def gc.cop_sq2(+a: Nat, +b: Nat, h: Cop<a, b>) -> Cop<Nat.mul(a, a), Nat.mul(b, b)>:
  gc.cop_sym(Nat.mul(b, b), Nat.mul(a, a),
    gc.cop_sq(Nat.mul(b, b), a,
      gc.cop_sym(a, Nat.mul(b, b), gc.cop_sq(a, b, h))))

def gc.cop_of_sq.go(+a: Nat, +b: Nat, +u: I.Int, +v: I.Int, h: {I.Int.add(I.Int.mul(u, I.Pos{Nat.mul(a, a)}), I.Int.mul(v, I.Pos{Nat.mul(b, b)})) == I.Pos{1n} : I.Int}) -> Cop<a, b>:
  Cop{I.Int.mul(u, I.Pos{a}), I.Int.mul(v, I.Pos{b}),
    %I.mul_assoc(u, I.Pos{a}, I.Pos{a}) : {I.Int.add(_, I.Int.mul(I.Int.mul(v, I.Pos{b}), I.Pos{b})) == I.Pos{1n} : I.Int}
    %I.mul_assoc(v, I.Pos{b}, I.Pos{b}) : {I.Int.add(I.Int.mul(u, I.Int.mul(I.Pos{a}, I.Pos{a})), _) == I.Pos{1n} : I.Int}
    %Equal.sym(I.Int, I.Int.mul(I.Pos{a}, I.Pos{a}), I.Pos{Nat.mul(a, a)}, I.of_mul(a, a)) : {I.Int.add(I.Int.mul(u, _), I.Int.mul(v, I.Int.mul(I.Pos{b}, I.Pos{b}))) == I.Pos{1n} : I.Int}
    %Equal.sym(I.Int, I.Int.mul(I.Pos{b}, I.Pos{b}), I.Pos{Nat.mul(b, b)}, I.of_mul(b, b)) : {I.Int.add(I.Int.mul(u, I.Pos{Nat.mul(a, a)}), I.Int.mul(v, _)) == I.Pos{1n} : I.Int}
    h}

def gc.cop_of_sq(+a: Nat, +b: Nat, h: Cop<Nat.mul(a, a), Nat.mul(b, b)>) -> Cop<a, b>:
  match h:
    case Cop{u, v, e}:
      gc.cop_of_sq.go(a, b, u, v, e)

# squares are injective: a^2 = b^2 forces a = b
def gc.sq_inj.no(+a: Nat, +b: Nat, e: {Nat.mul(a, a) == Nat.mul(b, b) : Nat}, lt: N.Nat.Lt(a, b)) -> Empty:
  N.Nat.lt_irrefl(Nat.mul(a, a),
    N.Nat.le_trans(1n+Nat.mul(a, a), Nat.mul(b, b), Nat.mul(a, a),
      N.Nat.le_trans(1n+Nat.mul(a, a), Nat.mul(1n+a, 1n+a), Nat.mul(b, b),
        %Equal.sym(Nat, Nat.mul(1n+a, 1n+a), Nat.add(1n+a, Nat.mul(1n+a, a)), N.Nat.mul_succ(1n+a, a)) : N.Nat.Le(1n+Nat.mul(a, a), _)
        N.Nat.le_trans(Nat.mul(a, a), Nat.add(a, Nat.mul(a, a)), Nat.add(a, Nat.add(a, Nat.mul(a, a))),
          N.Nat.le_add_l(a, Nat.mul(a, a)), N.Nat.le_add_l(a, Nat.add(a, Nat.mul(a, a)))),
        S2.Nat.le_sq(1n+a, b, lt)),
      N.Nat.eq_le(Nat.mul(b, b), Nat.mul(a, a), Equal.sym(Nat, Nat.mul(a, a), Nat.mul(b, b), e))))

def gc.sq_inj(+a: Nat, +b: Nat, +e: {Nat.mul(a, a) == Nat.mul(b, b) : Nat}) -> {a == b : Nat}:
  S2.Nat.decide_lt(a, b, {a == b : Nat},
    lt => Empty.absurd({a == b : Nat}, gc.sq_inj.no(a, b, e, lt)),
    ge => S2.Nat.decide_lt(b, a, {a == b : Nat},
      lt2 => Empty.absurd({a == b : Nat}, gc.sq_inj.no(b, a, Equal.sym(Nat, Nat.mul(a, a), Nat.mul(b, b), e), lt2)),
      ge2 => N.Nat.le_antisym(a, b, ge2, ge)))

# a^2 | b^2 forces a | b: with G = gcd(a, b), a = G a1, b = G b1 and
# gcd(a1, b1) = 1, cancelling G^2 leaves a1^2 | b1^2, and a1^2 is a common
# divisor of the coprime pair a1^2, b1^2, so a1 = 1 and a = G divides b.
def gc.sq_dvd.fin(+ap: Nat, +bp: Nat, +gp: Nat, +a1: Nat, +b1: Nat, +ea: {Nat.mul(1n+gp, a1) == 1n+ap : Nat}, +eb: {Nat.mul(1n+gp, b1) == 1n+bp : Nat}, hcop: Cop<a1, b1>, +k: Nat, ek: {Nat.mul(Nat.mul(1n+ap, 1n+ap), k) == Nat.mul(1n+bp, 1n+bp) : Nat}) -> Dvd<1n+ap, 1n+bp>:
  +one = gc.mul_eq_one(a1, a1,
    gc.cop_common(Nat.mul(a1, a1), Nat.mul(b1, b1), Nat.mul(a1, a1),
      gc.cop_sq2(a1, b1, hcop),
      gc.dvd_refl(Nat.mul(a1, a1)),
      Dvd{k,
        gc.ncancel_pos(Nat.mul(Nat.mul(a1, a1), k), Nat.mul(b1, b1), Nat.mul(1n+gp, 1n+gp), Unit{},
          Equal.trans(Nat, Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(Nat.mul(a1, a1), k)), Nat.mul(Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(a1, a1)), k), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(b1, b1)),
            N.Nat.mul_assoc(Nat.mul(1n+gp, 1n+gp), Nat.mul(a1, a1), k),
            Equal.trans(Nat, Nat.mul(Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(a1, a1)), k), Nat.mul(Nat.mul(1n+ap, 1n+ap), k), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(b1, b1)),
              Equal.cong(Nat, Nat, z => Nat.mul(z, k), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(a1, a1)), Nat.mul(1n+ap, 1n+ap),
                Equal.trans(Nat, Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(a1, a1)), Nat.mul(Nat.mul(1n+gp, a1), Nat.mul(1n+gp, a1)), Nat.mul(1n+ap, 1n+ap),
                  Equal.sym(Nat, Nat.mul(Nat.mul(1n+gp, a1), Nat.mul(1n+gp, a1)), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(a1, a1)), gc.nsq_prod(1n+gp, a1)),
                  Equal.cong(Nat, Nat, z => Nat.mul(z, z), Nat.mul(1n+gp, a1), 1n+ap, ea))),
              Equal.trans(Nat, Nat.mul(Nat.mul(1n+ap, 1n+ap), k), Nat.mul(1n+bp, 1n+bp), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(b1, b1)),
                ek,
                Equal.trans(Nat, Nat.mul(1n+bp, 1n+bp), Nat.mul(Nat.mul(1n+gp, b1), Nat.mul(1n+gp, b1)), Nat.mul(Nat.mul(1n+gp, 1n+gp), Nat.mul(b1, b1)),
                  Equal.cong(Nat, Nat, z => Nat.mul(z, z), 1n+bp, Nat.mul(1n+gp, b1), Equal.sym(Nat, Nat.mul(1n+gp, b1), 1n+bp, eb)),
                  gc.nsq_prod(1n+gp, b1))))))}))
  Dvd{b1,
    Equal.trans(Nat, Nat.mul(1n+ap, b1), Nat.mul(1n+gp, b1), 1n+bp,
      Equal.cong(Nat, Nat, z => Nat.mul(z, b1), 1n+ap, 1n+gp,
        Equal.sym(Nat, 1n+gp, 1n+ap,
          Equal.trans(Nat, 1n+gp, Nat.mul(1n+gp, 1n), 1n+ap,
            Equal.sym(Nat, Nat.mul(1n+gp, 1n), 1n+gp, N.Nat.mul_one(1n+gp)),
            Equal.trans(Nat, Nat.mul(1n+gp, 1n), Nat.mul(1n+gp, a1), 1n+ap,
              Equal.cong(Nat, Nat, z => Nat.mul(1n+gp, z), 1n, a1, Equal.sym(Nat, a1, 1n, one)),
              ea)))),
      eb)}

def gc.sq_dvd.op(+ap: Nat, +bp: Nat, sp: Split<1n+ap, 1n+bp>, h: Dvd<Nat.mul(1n+ap, 1n+ap), Nat.mul(1n+bp, 1n+bp)>) -> Dvd<1n+ap, 1n+bp>:
  match sp:
    case Split{gp, a1, b1, ea, eb, hcop}:
      match h:
        case Dvd{k, ek}:
          gc.sq_dvd.fin(ap, bp, gp, a1, b1, ea, eb, hcop, k, ek)

def gc.sq_dvd(+ap: Nat, +bp: Nat, h: Dvd<Nat.mul(1n+ap, 1n+ap), Nat.mul(1n+bp, 1n+bp)>) -> Dvd<1n+ap, 1n+bp>:
  gc.sq_dvd.op(ap, bp, gc.split(ap, bp), h)

# the same with positivity as a hypothesis rather than a shape
def gc.sq_dvd_pos(a: Nat, pa: N.Nat.Lt(0n, a), b: Nat, pb: N.Nat.Lt(0n, b), h: Dvd<Nat.mul(a, a), Nat.mul(b, b)>) -> Dvd<a, b>:
  match a:
    case 0n:
      match pa:
    case 1n+ap:
      match b:
        case 0n:
          match pb:
        case 1n+bp:
          gc.sq_dvd(ap, bp, h)