# lib/pythag.bend -- primitive Pythagorean triples, on top of lib/gcd.bend.
#
# THE THEOREM (Pyth, at the bottom of the file). If
#
#     A^2 + B^2 = C^2,   A, B, C > 0,   gcd(A, B) = 1,   B even
#
# then there are m > n > 0, coprime, with
#
#     A + n^2 = m^2      (i.e. A = m^2 - n^2)
#     B = 2 m n
#     C = m^2 + n^2
#
# Everything is stated additively, so no truncated subtraction is ever
# needed: "A = m^2 - n^2" is "A + n^2 = m^2", and the two halves of C - A
# and C + A arrive as witnesses.
#
# The proof: A and C are odd (A because it is coprime to the even B, C
# because C^2 = A^2 + B^2 is odd), so C - A and C + A are even; write
# C = A + 2u and v = A + u, so that u + v = C and u + A = v. Expanding
# C^2 = (A + 2u)^2 in A^2 + B^2 = C^2 leaves B^2 = 4 u v, and B = 2k gives
# k^2 = u v. A common divisor of u and v divides u + v = C and v - u = A,
# and gcd(A, C) = 1 (a common divisor of A and C divides C^2 - A^2 = B^2,
# and gcd(A^2, B^2) = 1), so u and v are coprime. The square theorem of
# lib/gcd.bend then makes u = n^2 and v = m^2, and the three identities
# above fall out.
#
# The file also proves the ONE congruence fact the descent needs:
# odd^2 + odd^2 is never a square (mod 4 it would be 2). It is used twice
# in demos/math/fermat4.bend -- once to see that x and y cannot both be
# odd, and once to see that the smaller leg of the second triple is even.

import Base
import ./nat.bend as N
import ./sqrt2.bend as S2
import ./gcd.bend as G

# ---------------------------------------------------------------------
# doubling algebra
# ---------------------------------------------------------------------

# 2x + 2y = 2(x + y)
def py.double_dist(+x: Nat, +y: Nat) -> {Nat.add(Nat.double(x), Nat.double(y)) == Nat.double(Nat.add(x, y)) : Nat}:
  %Equal.sym(Nat, Nat.double(Nat.add(x, y)), Nat.add(Nat.add(x, y), Nat.add(x, y)), N.Nat.double_add(Nat.add(x, y))) : {Nat.add(Nat.double(x), Nat.double(y)) == _ : Nat}
  %Equal.sym(Nat, Nat.double(x), Nat.add(x, x), N.Nat.double_add(x)) : {Nat.add(_, Nat.double(y)) == Nat.add(Nat.add(x, y), Nat.add(x, y)) : Nat}
  %Equal.sym(Nat, Nat.double(y), Nat.add(y, y), N.Nat.double_add(y)) : {Nat.add(Nat.add(x, x), _) == Nat.add(Nat.add(x, y), Nat.add(x, y)) : Nat}
  N.Nat.add_add_swap(x, x, y, y)

# a (2b) = 2(ab)
def py.mul_double(+a: Nat, +b: Nat) -> {Nat.mul(a, Nat.double(b)) == Nat.double(Nat.mul(a, b)) : Nat}:
  %Equal.sym(Nat, Nat.double(b), Nat.add(b, b), N.Nat.double_add(b)) : {Nat.mul(a, _) == Nat.double(Nat.mul(a, b)) : Nat}
  %Equal.sym(Nat, Nat.double(Nat.mul(a, b)), Nat.add(Nat.mul(a, b), Nat.mul(a, b)), N.Nat.double_add(Nat.mul(a, b))) : {Nat.mul(a, Nat.add(b, b)) == _ : Nat}
  N.Nat.mul_add(a, b, b)

# (2a) b = 2(ab)
def py.double_mul(+a: Nat, +b: Nat) -> {Nat.mul(Nat.double(a), b) == Nat.double(Nat.mul(a, b)) : Nat}:
  Equal.trans(Nat, Nat.mul(Nat.double(a), b), Nat.mul(b, Nat.double(a)), Nat.double(Nat.mul(a, b)),
    N.Nat.mul_comm(Nat.double(a), b),
    Equal.trans(Nat, Nat.mul(b, Nat.double(a)), Nat.double(Nat.mul(b, a)), Nat.double(Nat.mul(a, b)),
      py.mul_double(b, a),
      Equal.cong(Nat, Nat, z => Nat.double(z), Nat.mul(b, a), Nat.mul(a, b), N.Nat.mul_comm(b, a))))

# ---------------------------------------------------------------------
# parity: shapes, and "2 divides an even number"
# ---------------------------------------------------------------------

def py.not_false(b: Bool, e: {Bool.not(b) == False{} : Bool}) -> {b == True{} : Bool}:
  match b:
    case True{}:
      {==}
    case False{}:
      Empty.absurd({False{} == True{} : Bool}, N.Bool.true_ne_false(e))

def py.not_true(b: Bool, e: {Bool.not(b) == True{} : Bool}) -> {b == False{} : Bool}:
  match b:
    case True{}:
      Empty.absurd({True{} == False{} : Bool}, N.Bool.false_ne_true(e))
    case False{}:
      {==}

def py.dvd2_of_even(+n: Nat, e: {N.Nat.even(n) == True{} : Bool}) -> G.Dvd<2n, n>:
  G.Dvd{N.Nat.half(n),
    Equal.trans(Nat, Nat.mul(2n, N.Nat.half(n)), Nat.double(N.Nat.half(n)), n,
      N.Nat.two_mul(N.Nat.half(n)),
      Equal.sym(Nat, n, Nat.double(N.Nat.half(n)), N.Nat.even_half(n, e)))}

def py.even_shape(+n: Nat, e: {N.Nat.even(n) == True{} : Bool}, -P: Type, k: @a: Nat -> @h: {n == Nat.double(a) : Nat} -> P) -> P:
  k(N.Nat.half(n), N.Nat.even_half(n, e))

def py.odd_shape(n: Nat, e: {N.Nat.even(n) == False{} : Bool}, -P: Type, k: @a: Nat -> @h: {n == 1n+Nat.double(a) : Nat} -> P) -> P:
  match n:
    case 0n:
      Empty.absurd(P, N.Bool.true_ne_false(e))
    case 1n+p:
      +p2 = p
      k(N.Nat.half(p2),
        Equal.cong(Nat, Nat, z => 1n+z, p2, Nat.double(N.Nat.half(p2)),
          N.Nat.even_half(p2, py.not_false(N.Nat.even(p2), e))))

# ---------------------------------------------------------------------
# odd^2 + odd^2 is never a square: (1+2a)^2 + (1+2b)^2 = 2 + 4(...) is
# twice an odd number, while a square that is even is four times a square
# ---------------------------------------------------------------------

# (1 + 2a)^2 = 1 + 4(a + a^2)
def py.sq_odd(+a: Nat) -> {Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)) == 1n+Nat.double(Nat.double(Nat.add(a, Nat.mul(a, a)))) : Nat}:
  Equal.cong(Nat, Nat, z => 1n+z, Nat.add(Nat.double(a), Nat.mul(Nat.double(a), 1n+Nat.double(a))), Nat.double(Nat.double(Nat.add(a, Nat.mul(a, a)))),
    Equal.trans(Nat, Nat.add(Nat.double(a), Nat.mul(Nat.double(a), 1n+Nat.double(a))), Nat.add(Nat.double(a), Nat.add(Nat.double(a), Nat.mul(Nat.double(a), Nat.double(a)))), Nat.double(Nat.double(Nat.add(a, Nat.mul(a, a)))),
      Equal.cong(Nat, Nat, z => Nat.add(Nat.double(a), z), Nat.mul(Nat.double(a), 1n+Nat.double(a)), Nat.add(Nat.double(a), Nat.mul(Nat.double(a), Nat.double(a))), N.Nat.mul_succ(Nat.double(a), Nat.double(a))),
      Equal.trans(Nat, Nat.add(Nat.double(a), Nat.add(Nat.double(a), Nat.mul(Nat.double(a), Nat.double(a)))), Nat.add(Nat.double(a), Nat.add(Nat.double(a), Nat.double(Nat.double(Nat.mul(a, a))))), Nat.double(Nat.double(Nat.add(a, Nat.mul(a, a)))),
        Equal.cong(Nat, Nat, z => Nat.add(Nat.double(a), Nat.add(Nat.double(a), z)), Nat.mul(Nat.double(a), Nat.double(a)), Nat.double(Nat.double(Nat.mul(a, a))), S2.Nat.sq_double(a)),
        Equal.trans(Nat, Nat.add(Nat.double(a), Nat.add(Nat.double(a), Nat.double(Nat.double(Nat.mul(a, a))))), Nat.add(Nat.add(Nat.double(a), Nat.double(a)), Nat.double(Nat.double(Nat.mul(a, a)))), Nat.double(Nat.double(Nat.add(a, Nat.mul(a, a)))),
          N.Nat.add_assoc(Nat.double(a), Nat.double(a), Nat.double(Nat.double(Nat.mul(a, a)))),
          Equal.trans(Nat, Nat.add(Nat.add(Nat.double(a), Nat.double(a)), Nat.double(Nat.double(Nat.mul(a, a)))), Nat.add(Nat.double(Nat.double(a)), Nat.double(Nat.double(Nat.mul(a, a)))), Nat.double(Nat.double(Nat.add(a, Nat.mul(a, a)))),
            Equal.cong(Nat, Nat, z => Nat.add(z, Nat.double(Nat.double(Nat.mul(a, a)))), Nat.add(Nat.double(a), Nat.double(a)), Nat.double(Nat.double(a)),
              Equal.sym(Nat, Nat.double(Nat.double(a)), Nat.add(Nat.double(a), Nat.double(a)), N.Nat.double_add(Nat.double(a)))),
            Equal.trans(Nat, Nat.add(Nat.double(Nat.double(a)), Nat.double(Nat.double(Nat.mul(a, a)))), Nat.double(Nat.add(Nat.double(a), Nat.double(Nat.mul(a, a)))), Nat.double(Nat.double(Nat.add(a, Nat.mul(a, a)))),
              py.double_dist(Nat.double(a), Nat.double(Nat.mul(a, a))),
              Equal.cong(Nat, Nat, z => Nat.double(z), Nat.add(Nat.double(a), Nat.double(Nat.mul(a, a))), Nat.double(Nat.add(a, Nat.mul(a, a))), py.double_dist(a, Nat.mul(a, a)))))))))

def py.odd_odd.core(+a: Nat, +b: Nat, +w: Nat, e: {Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))) == Nat.mul(Nat.double(w), Nat.double(w)) : Nat}) -> Empty:
  +pp = Nat.add(a, Nat.mul(a, a))
  +qq = Nat.add(b, Nat.mul(b, b))
  +key = N.Nat.double_inj(1n+Nat.double(Nat.add(pp, qq)), Nat.double(Nat.mul(w, w)),
    Equal.trans(Nat, Nat.double(1n+Nat.double(Nat.add(pp, qq))), Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.double(Nat.double(Nat.mul(w, w))),
      Equal.sym(Nat, Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.double(1n+Nat.double(Nat.add(pp, qq))),
        Equal.trans(Nat, Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.add(1n+Nat.double(Nat.double(pp)), 1n+Nat.double(Nat.double(qq))), Nat.double(1n+Nat.double(Nat.add(pp, qq))),
          Equal.trans(Nat, Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.add(1n+Nat.double(Nat.double(pp)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.add(1n+Nat.double(Nat.double(pp)), 1n+Nat.double(Nat.double(qq))),
            Equal.cong(Nat, Nat, z => Nat.add(z, Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), 1n+Nat.double(Nat.double(pp)), py.sq_odd(a)),
            Equal.cong(Nat, Nat, z => Nat.add(1n+Nat.double(Nat.double(pp)), z), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b)), 1n+Nat.double(Nat.double(qq)), py.sq_odd(b))),
          Equal.trans(Nat, Nat.add(1n+Nat.double(Nat.double(pp)), 1n+Nat.double(Nat.double(qq))), 2n+Nat.add(Nat.double(Nat.double(pp)), Nat.double(Nat.double(qq))), Nat.double(1n+Nat.double(Nat.add(pp, qq))),
            Equal.cong(Nat, Nat, z => 1n+z, Nat.add(Nat.double(Nat.double(pp)), 1n+Nat.double(Nat.double(qq))), 1n+Nat.add(Nat.double(Nat.double(pp)), Nat.double(Nat.double(qq))), N.Nat.add_succ(Nat.double(Nat.double(pp)), Nat.double(Nat.double(qq)))),
            Equal.cong(Nat, Nat, z => 2n+z, Nat.add(Nat.double(Nat.double(pp)), Nat.double(Nat.double(qq))), Nat.double(Nat.double(Nat.add(pp, qq))),
              Equal.trans(Nat, Nat.add(Nat.double(Nat.double(pp)), Nat.double(Nat.double(qq))), Nat.double(Nat.add(Nat.double(pp), Nat.double(qq))), Nat.double(Nat.double(Nat.add(pp, qq))),
                py.double_dist(Nat.double(pp), Nat.double(qq)),
                Equal.cong(Nat, Nat, z => Nat.double(z), Nat.add(Nat.double(pp), Nat.double(qq)), Nat.double(Nat.add(pp, qq)), py.double_dist(pp, qq))))))),
      Equal.trans(Nat, Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.mul(Nat.double(w), Nat.double(w)), Nat.double(Nat.double(Nat.mul(w, w))),
        e, S2.Nat.sq_double(w))))
  N.Bool.false_ne_true(
    Equal.trans(Bool, False{}, N.Nat.even(Nat.double(Nat.mul(w, w))), True{},
      Equal.trans(Bool, False{}, N.Nat.even(1n+Nat.double(Nat.add(pp, qq))), N.Nat.even(Nat.double(Nat.mul(w, w))),
        Equal.sym(Bool, N.Nat.even(1n+Nat.double(Nat.add(pp, qq))), False{},
          Equal.cong(Bool, Bool, z => Bool.not(z), N.Nat.even(Nat.double(Nat.add(pp, qq))), True{}, N.Nat.even_double(Nat.add(pp, qq)))),
        Equal.cong(Nat, Bool, z => N.Nat.even(z), 1n+Nat.double(Nat.add(pp, qq)), Nat.double(Nat.mul(w, w)), key)),
      N.Nat.even_double(Nat.mul(w, w))))

def py.odd_odd.fin(+x: Nat, +y: Nat, +z: Nat, +a: Nat, +b: Nat, +w: Nat, hx: {x == 1n+Nat.double(a) : Nat}, hy: {y == 1n+Nat.double(b) : Nat}, hz: {z == Nat.double(w) : Nat}, e: {Nat.add(Nat.mul(x, x), Nat.mul(y, y)) == Nat.mul(z, z) : Nat}) -> Empty:
  py.odd_odd.core(a, b, w,
    Equal.trans(Nat, Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.add(Nat.mul(x, x), Nat.mul(y, y)), Nat.mul(Nat.double(w), Nat.double(w)),
      Equal.trans(Nat, Nat.add(Nat.mul(1n+Nat.double(a), 1n+Nat.double(a)), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.add(Nat.mul(x, x), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), Nat.add(Nat.mul(x, x), Nat.mul(y, y)),
        Equal.cong(Nat, Nat, t => Nat.add(Nat.mul(t, t), Nat.mul(1n+Nat.double(b), 1n+Nat.double(b))), 1n+Nat.double(a), x, Equal.sym(Nat, x, 1n+Nat.double(a), hx)),
        Equal.cong(Nat, Nat, t => Nat.add(Nat.mul(x, x), Nat.mul(t, t)), 1n+Nat.double(b), y, Equal.sym(Nat, y, 1n+Nat.double(b), hy))),
      Equal.trans(Nat, Nat.add(Nat.mul(x, x), Nat.mul(y, y)), Nat.mul(z, z), Nat.mul(Nat.double(w), Nat.double(w)),
        e, Equal.cong(Nat, Nat, t => Nat.mul(t, t), z, Nat.double(w), hz))))

# THE CONGRUENCE FACT: x and y odd makes x^2 + y^2 no square.
def py.odd_odd(+x: Nat, +y: Nat, +z: Nat, +ex: {N.Nat.even(x) == False{} : Bool}, +ey: {N.Nat.even(y) == False{} : Bool}, +e: {Nat.add(Nat.mul(x, x), Nat.mul(y, y)) == Nat.mul(z, z) : Nat}) -> Empty:
  py.odd_shape(x, ex, Empty, a => hx =>
    py.odd_shape(y, ey, Empty, b => hy =>
      py.even_shape(z,
        Equal.trans(Bool, N.Nat.even(z), N.Nat.even(Nat.mul(z, z)), True{},
          Equal.sym(Bool, N.Nat.even(Nat.mul(z, z)), N.Nat.even(z), N.Nat.even_sq(z)),
          Equal.trans(Bool, N.Nat.even(Nat.mul(z, z)), N.Nat.even(Nat.add(Nat.mul(x, x), Nat.mul(y, y))), True{},
            Equal.cong(Nat, Bool, t => N.Nat.even(t), Nat.mul(z, z), Nat.add(Nat.mul(x, x), Nat.mul(y, y)), Equal.sym(Nat, Nat.add(Nat.mul(x, x), Nat.mul(y, y)), Nat.mul(z, z), e)),
            Equal.trans(Bool, N.Nat.even(Nat.add(Nat.mul(x, x), Nat.mul(y, y))), N.Bool.same(N.Nat.even(Nat.mul(x, x)), N.Nat.even(Nat.mul(y, y))), True{},
              N.Nat.even_add(Nat.mul(x, x), Nat.mul(y, y)),
              Equal.trans(Bool, N.Bool.same(N.Nat.even(Nat.mul(x, x)), N.Nat.even(Nat.mul(y, y))), N.Bool.same(N.Nat.even(x), N.Nat.even(y)), True{},
                Equal.trans(Bool, N.Bool.same(N.Nat.even(Nat.mul(x, x)), N.Nat.even(Nat.mul(y, y))), N.Bool.same(N.Nat.even(x), N.Nat.even(Nat.mul(y, y))), N.Bool.same(N.Nat.even(x), N.Nat.even(y)),
                  Equal.cong(Bool, Bool, t => N.Bool.same(t, N.Nat.even(Nat.mul(y, y))), N.Nat.even(Nat.mul(x, x)), N.Nat.even(x), N.Nat.even_sq(x)),
                  Equal.cong(Bool, Bool, t => N.Bool.same(N.Nat.even(x), t), N.Nat.even(Nat.mul(y, y)), N.Nat.even(y), N.Nat.even_sq(y))),
                Equal.trans(Bool, N.Bool.same(N.Nat.even(x), N.Nat.even(y)), N.Bool.same(False{}, N.Nat.even(y)), True{},
                  Equal.cong(Bool, Bool, t => N.Bool.same(t, N.Nat.even(y)), N.Nat.even(x), False{}, ex),
                  Equal.cong(Bool, Bool, t => N.Bool.same(False{}, t), N.Nat.even(y), False{}, ey)))))),
        Empty, w => hz =>
          py.odd_odd.fin(x, y, z, a, b, w, hx, hy, hz, e))))

# ---------------------------------------------------------------------
# (a + 2u)^2 = a^2 + 4 u (a + u)  -- the identity that turns
# A^2 + B^2 = C^2 with C = A + 2u into B^2 = 4 u v, v = A + u
# ---------------------------------------------------------------------

def py.expand.inner(+xa: Nat, +u: Nat) -> {Nat.add(Nat.mul(xa, Nat.double(u)), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))) == Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u)))) : Nat}:
  %Equal.sym(Nat, Nat.mul(xa, Nat.double(u)), Nat.double(Nat.mul(xa, u)), py.mul_double(xa, u)) : {Nat.add(_, Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))) == Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u)))) : Nat}
  %N.Nat.mul_comm(u, xa) : {Nat.add(Nat.double(_), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))) == Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u)))) : Nat}
  %Equal.sym(Nat, Nat.mul(Nat.double(u), xa), Nat.double(Nat.mul(u, xa)), py.double_mul(u, xa)) : {Nat.add(Nat.double(Nat.mul(u, xa)), Nat.add(_, Nat.mul(Nat.double(u), Nat.double(u)))) == Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u)))) : Nat}
  %Equal.sym(Nat, Nat.mul(Nat.double(u), Nat.double(u)), Nat.double(Nat.double(Nat.mul(u, u))), S2.Nat.sq_double(u)) : {Nat.add(Nat.double(Nat.mul(u, xa)), Nat.add(Nat.double(Nat.mul(u, xa)), _)) == Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u)))) : Nat}
  Equal.trans(Nat,
    Nat.add(Nat.double(Nat.mul(u, xa)), Nat.add(Nat.double(Nat.mul(u, xa)), Nat.double(Nat.double(Nat.mul(u, u))))),
    Nat.double(Nat.add(Nat.mul(u, xa), Nat.add(Nat.mul(u, xa), Nat.double(Nat.mul(u, u))))),
    Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u)))),
    Equal.trans(Nat,
      Nat.add(Nat.double(Nat.mul(u, xa)), Nat.add(Nat.double(Nat.mul(u, xa)), Nat.double(Nat.double(Nat.mul(u, u))))),
      Nat.add(Nat.double(Nat.mul(u, xa)), Nat.double(Nat.add(Nat.mul(u, xa), Nat.double(Nat.mul(u, u))))),
      Nat.double(Nat.add(Nat.mul(u, xa), Nat.add(Nat.mul(u, xa), Nat.double(Nat.mul(u, u))))),
      Equal.cong(Nat, Nat, z => Nat.add(Nat.double(Nat.mul(u, xa)), z),
        Nat.add(Nat.double(Nat.mul(u, xa)), Nat.double(Nat.double(Nat.mul(u, u)))),
        Nat.double(Nat.add(Nat.mul(u, xa), Nat.double(Nat.mul(u, u)))),
        py.double_dist(Nat.mul(u, xa), Nat.double(Nat.mul(u, u)))),
      py.double_dist(Nat.mul(u, xa), Nat.add(Nat.mul(u, xa), Nat.double(Nat.mul(u, u))))),
    Equal.cong(Nat, Nat, z => Nat.double(z),
      Nat.add(Nat.mul(u, xa), Nat.add(Nat.mul(u, xa), Nat.double(Nat.mul(u, u)))),
      Nat.double(Nat.mul(u, Nat.add(xa, u))),
      Equal.trans(Nat,
        Nat.add(Nat.mul(u, xa), Nat.add(Nat.mul(u, xa), Nat.double(Nat.mul(u, u)))),
        Nat.add(Nat.add(Nat.mul(u, xa), Nat.mul(u, xa)), Nat.double(Nat.mul(u, u))),
        Nat.double(Nat.mul(u, Nat.add(xa, u))),
        N.Nat.add_assoc(Nat.mul(u, xa), Nat.mul(u, xa), Nat.double(Nat.mul(u, u))),
        Equal.trans(Nat,
          Nat.add(Nat.add(Nat.mul(u, xa), Nat.mul(u, xa)), Nat.double(Nat.mul(u, u))),
          Nat.add(Nat.double(Nat.mul(u, xa)), Nat.double(Nat.mul(u, u))),
          Nat.double(Nat.mul(u, Nat.add(xa, u))),
          Equal.cong(Nat, Nat, z => Nat.add(z, Nat.double(Nat.mul(u, u))),
            Nat.add(Nat.mul(u, xa), Nat.mul(u, xa)), Nat.double(Nat.mul(u, xa)),
            Equal.sym(Nat, Nat.double(Nat.mul(u, xa)), Nat.add(Nat.mul(u, xa), Nat.mul(u, xa)), N.Nat.double_add(Nat.mul(u, xa)))),
          Equal.trans(Nat,
            Nat.add(Nat.double(Nat.mul(u, xa)), Nat.double(Nat.mul(u, u))),
            Nat.double(Nat.add(Nat.mul(u, xa), Nat.mul(u, u))),
            Nat.double(Nat.mul(u, Nat.add(xa, u))),
            py.double_dist(Nat.mul(u, xa), Nat.mul(u, u)),
            Equal.cong(Nat, Nat, z => Nat.double(z), Nat.add(Nat.mul(u, xa), Nat.mul(u, u)), Nat.mul(u, Nat.add(xa, u)),
              Equal.sym(Nat, Nat.mul(u, Nat.add(xa, u)), Nat.add(Nat.mul(u, xa), Nat.mul(u, u)), N.Nat.mul_add(u, xa, u))))))))

def py.expand(+xa: Nat, +u: Nat) -> {Nat.mul(Nat.add(xa, Nat.double(u)), Nat.add(xa, Nat.double(u))) == Nat.add(Nat.mul(xa, xa), Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u))))) : Nat}:
  Equal.trans(Nat,
    Nat.mul(Nat.add(xa, Nat.double(u)), Nat.add(xa, Nat.double(u))),
    Nat.add(Nat.mul(xa, Nat.add(xa, Nat.double(u))), Nat.mul(Nat.double(u), Nat.add(xa, Nat.double(u)))),
    Nat.add(Nat.mul(xa, xa), Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u))))),
    N.Nat.add_mul(xa, Nat.double(u), Nat.add(xa, Nat.double(u))),
    Equal.trans(Nat,
      Nat.add(Nat.mul(xa, Nat.add(xa, Nat.double(u))), Nat.mul(Nat.double(u), Nat.add(xa, Nat.double(u)))),
      Nat.add(Nat.add(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u))), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))),
      Nat.add(Nat.mul(xa, xa), Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u))))),
      Equal.trans(Nat,
        Nat.add(Nat.mul(xa, Nat.add(xa, Nat.double(u))), Nat.mul(Nat.double(u), Nat.add(xa, Nat.double(u)))),
        Nat.add(Nat.add(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u))), Nat.mul(Nat.double(u), Nat.add(xa, Nat.double(u)))),
        Nat.add(Nat.add(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u))), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))),
        Equal.cong(Nat, Nat, z => Nat.add(z, Nat.mul(Nat.double(u), Nat.add(xa, Nat.double(u)))),
          Nat.mul(xa, Nat.add(xa, Nat.double(u))), Nat.add(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u))),
          N.Nat.mul_add(xa, xa, Nat.double(u))),
        Equal.cong(Nat, Nat, z => Nat.add(Nat.add(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u))), z),
          Nat.mul(Nat.double(u), Nat.add(xa, Nat.double(u))), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u))),
          N.Nat.mul_add(Nat.double(u), xa, Nat.double(u)))),
      Equal.trans(Nat,
        Nat.add(Nat.add(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u))), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))),
        Nat.add(Nat.mul(xa, xa), Nat.add(Nat.mul(xa, Nat.double(u)), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u))))),
        Nat.add(Nat.mul(xa, xa), Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u))))),
        Equal.sym(Nat,
          Nat.add(Nat.mul(xa, xa), Nat.add(Nat.mul(xa, Nat.double(u)), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u))))),
          Nat.add(Nat.add(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u))), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))),
          N.Nat.add_assoc(Nat.mul(xa, xa), Nat.mul(xa, Nat.double(u)), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u))))),
        Equal.cong(Nat, Nat, z => Nat.add(Nat.mul(xa, xa), z),
          Nat.add(Nat.mul(xa, Nat.double(u)), Nat.add(Nat.mul(Nat.double(u), xa), Nat.mul(Nat.double(u), Nat.double(u)))),
          Nat.double(Nat.double(Nat.mul(u, Nat.add(xa, u)))),
          py.expand.inner(xa, u)))))

# ---------------------------------------------------------------------
# the pieces of the parametrisation
# ---------------------------------------------------------------------

def py.same_false(x: Bool, e: {N.Bool.same(False{}, x) == False{} : Bool}) -> {x == True{} : Bool}:
  match x:
    case True{}:
      {==}
    case False{}:
      Empty.absurd({False{} == True{} : Bool}, N.Bool.true_ne_false(e))

def py.lt_add_pos(+a: Nat, +b: Nat) -> N.Nat.Lt(a, Nat.add(a, 1n+b)):
  %Equal.sym(Nat, Nat.add(a, 1n+b), 1n+Nat.add(a, b), N.Nat.add_succ(a, b)) : N.Nat.Le(1n+a, _)
  N.Nat.le_add_r(a, b)

# A^2 < C^2, because B^2 is positive
def py.sq_lt(+av: Nat, +bv: Nat, +cv: Nat, hsq: {Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)) == Nat.mul(1n+cv, 1n+cv) : Nat}) -> N.Nat.Lt(Nat.mul(1n+av, 1n+av), Nat.mul(1n+cv, 1n+cv)):
  %hsq : N.Nat.Lt(Nat.mul(1n+av, 1n+av), _)
  py.lt_add_pos(Nat.mul(1n+av, 1n+av), Nat.add(bv, Nat.mul(bv, 1n+bv)))

# and therefore A < C
def py.lt_of_sq(+a: Nat, +c: Nat, h: N.Nat.Lt(Nat.mul(a, a), Nat.mul(c, c))) -> N.Nat.Lt(a, c):
  S2.Nat.decide_lt(a, c, N.Nat.Lt(a, c),
    lt => lt,
    ge => Empty.absurd(N.Nat.Lt(a, c),
      N.Nat.lt_irrefl(Nat.mul(a, a),
        N.Nat.le_trans(1n+Nat.mul(a, a), Nat.mul(c, c), Nat.mul(a, a), h, S2.Nat.le_sq(c, a, ge)))))

# A is odd: 2 would otherwise divide both A and the even B
def py.oddA.go(+av: Nat, +bv: Nat, hcop: G.Cop<1n+av, 1n+bv>, hbe: {N.Nat.even(1n+bv) == True{} : Bool}, c: Bool, ec: {N.Nat.even(1n+av) == c : Bool}) -> {N.Nat.even(1n+av) == False{} : Bool}:
  match c:
    case False{}:
      ec
    case True{}:
      Empty.absurd({N.Nat.even(1n+av) == False{} : Bool},
        N.Nat.succ_neq_zero(0n,
          N.Nat.succ_inj(1n, 0n,
            G.gc.cop_common(1n+av, 1n+bv, 2n, hcop,
              py.dvd2_of_even(1n+av, ec), py.dvd2_of_even(1n+bv, hbe)))))

def py.oddA(+av: Nat, +bv: Nat, hcop: G.Cop<1n+av, 1n+bv>, hbe: {N.Nat.even(1n+bv) == True{} : Bool}) -> {N.Nat.even(1n+av) == False{} : Bool}:
  py.oddA.go(av, bv, hcop, hbe, N.Nat.even(1n+av), {==})

# C is odd: C^2 = A^2 + B^2 is odd + even
def py.oddC(+av: Nat, +bv: Nat, +cv: Nat, hsq: {Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)) == Nat.mul(1n+cv, 1n+cv) : Nat}, exA: {N.Nat.even(1n+av) == False{} : Bool}, hbe: {N.Nat.even(1n+bv) == True{} : Bool}) -> {N.Nat.even(1n+cv) == False{} : Bool}:
  Equal.trans(Bool, N.Nat.even(1n+cv), N.Nat.even(Nat.mul(1n+cv, 1n+cv)), False{},
    Equal.sym(Bool, N.Nat.even(Nat.mul(1n+cv, 1n+cv)), N.Nat.even(1n+cv), N.Nat.even_sq(1n+cv)),
    Equal.trans(Bool, N.Nat.even(Nat.mul(1n+cv, 1n+cv)), N.Nat.even(Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv))), False{},
      Equal.cong(Nat, Bool, t => N.Nat.even(t), Nat.mul(1n+cv, 1n+cv), Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)),
        Equal.sym(Nat, Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)), Nat.mul(1n+cv, 1n+cv), hsq)),
      Equal.trans(Bool, N.Nat.even(Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv))), N.Bool.same(N.Nat.even(Nat.mul(1n+av, 1n+av)), N.Nat.even(Nat.mul(1n+bv, 1n+bv))), False{},
        N.Nat.even_add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)),
        Equal.trans(Bool, N.Bool.same(N.Nat.even(Nat.mul(1n+av, 1n+av)), N.Nat.even(Nat.mul(1n+bv, 1n+bv))), N.Bool.same(N.Nat.even(1n+av), N.Nat.even(1n+bv)), False{},
          Equal.trans(Bool, N.Bool.same(N.Nat.even(Nat.mul(1n+av, 1n+av)), N.Nat.even(Nat.mul(1n+bv, 1n+bv))), N.Bool.same(N.Nat.even(1n+av), N.Nat.even(Nat.mul(1n+bv, 1n+bv))), N.Bool.same(N.Nat.even(1n+av), N.Nat.even(1n+bv)),
            Equal.cong(Bool, Bool, t => N.Bool.same(t, N.Nat.even(Nat.mul(1n+bv, 1n+bv))), N.Nat.even(Nat.mul(1n+av, 1n+av)), N.Nat.even(1n+av), N.Nat.even_sq(1n+av)),
            Equal.cong(Bool, Bool, t => N.Bool.same(N.Nat.even(1n+av), t), N.Nat.even(Nat.mul(1n+bv, 1n+bv)), N.Nat.even(1n+bv), N.Nat.even_sq(1n+bv))),
          Equal.trans(Bool, N.Bool.same(N.Nat.even(1n+av), N.Nat.even(1n+bv)), N.Bool.same(False{}, N.Nat.even(1n+bv)), False{},
            Equal.cong(Bool, Bool, t => N.Bool.same(t, N.Nat.even(1n+bv)), N.Nat.even(1n+av), False{}, exA),
            Equal.cong(Bool, Bool, t => N.Bool.same(False{}, t), N.Nat.even(1n+bv), True{}, hbe))))))

# C - A is even and positive: C = A + 2(1+up)
def py.get_u.go(+av: Nat, +cv: Nat, exA: {N.Nat.even(1n+av) == False{} : Bool}, exC: {N.Nat.even(1n+cv) == False{} : Bool}, +t: Nat, +ht: {Nat.add(2n+av, t) == 1n+cv : Nat}, -P: Type, k: @up: Nat -> @hu: {Nat.add(1n+av, Nat.double(1n+up)) == 1n+cv : Nat} -> P) -> P:
  +ht2 = Equal.trans(Nat, Nat.add(1n+av, 1n+t), Nat.add(2n+av, t), 1n+cv,
    Equal.cong(Nat, Nat, z => 1n+z, Nat.add(av, 1n+t), 1n+Nat.add(av, t), N.Nat.add_succ(av, t)),
    ht)
  py.odd_shape(t,
    py.not_true(N.Nat.even(t),
      py.same_false(N.Nat.even(1n+t),
        Equal.trans(Bool, N.Bool.same(False{}, N.Nat.even(1n+t)), N.Bool.same(N.Nat.even(1n+av), N.Nat.even(1n+t)), False{},
          Equal.cong(Bool, Bool, z => N.Bool.same(z, N.Nat.even(1n+t)), False{}, N.Nat.even(1n+av), Equal.sym(Bool, N.Nat.even(1n+av), False{}, exA)),
          Equal.trans(Bool, N.Bool.same(N.Nat.even(1n+av), N.Nat.even(1n+t)), N.Nat.even(Nat.add(1n+av, 1n+t)), False{},
            Equal.sym(Bool, N.Nat.even(Nat.add(1n+av, 1n+t)), N.Bool.same(N.Nat.even(1n+av), N.Nat.even(1n+t)), N.Nat.even_add(1n+av, 1n+t)),
            Equal.trans(Bool, N.Nat.even(Nat.add(1n+av, 1n+t)), N.Nat.even(1n+cv), False{},
              Equal.cong(Nat, Bool, z => N.Nat.even(z), Nat.add(1n+av, 1n+t), 1n+cv, ht2),
              exC))))),
    P, s => hs =>
      k(s,
        Equal.trans(Nat, Nat.add(1n+av, Nat.double(1n+s)), Nat.add(1n+av, 1n+t), 1n+cv,
          Equal.sym(Nat, Nat.add(1n+av, 1n+t), Nat.add(1n+av, Nat.double(1n+s)),
            Equal.cong(Nat, Nat, z => Nat.add(1n+av, 1n+z), t, 1n+Nat.double(s), hs)),
          ht2)))

def py.get_u.open(+av: Nat, +cv: Nat, exA: {N.Nat.even(1n+av) == False{} : Bool}, exC: {N.Nat.even(1n+cv) == False{} : Bool}, -P: Type, k: @up: Nat -> @hu: {Nat.add(1n+av, Nat.double(1n+up)) == 1n+cv : Nat} -> P, w: (&t: Nat -> {Nat.add(2n+av, t) == 1n+cv : Nat})) -> P:
  (t, ht) = w
  py.get_u.go(av, cv, exA, exC, t, ht, P, k)

def py.get_u(+av: Nat, +cv: Nat, hlt: N.Nat.Lt(1n+av, 1n+cv), exA: {N.Nat.even(1n+av) == False{} : Bool}, exC: {N.Nat.even(1n+cv) == False{} : Bool}, -P: Type, k: @up: Nat -> @hu: {Nat.add(1n+av, Nat.double(1n+up)) == 1n+cv : Nat} -> P) -> P:
  py.get_u.open(av, cv, exA, exC, P, k, N.Nat.le_diff(2n+av, 1n+cv, hlt))

# gcd(A, C) = 1: a common divisor of A and C divides A^2 and C^2, hence
# also B^2 = C^2 - A^2, and A^2, B^2 are coprime
def py.copAC.k(+av: Nat, +bv: Nat, +cv: Nat, +hsq: {Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)) == Nat.mul(1n+cv, 1n+cv) : Nat}, +hcop: G.Cop<1n+av, 1n+bv>, +g: Nat, d1: G.Dvd<g, 1n+av>, d2: G.Dvd<g, 1n+cv>) -> {g == 1n : Nat}:
  +dA2 = G.gc.dvd_mul_r(g, 1n+av, 1n+av, d1)
  G.gc.cop_common(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv), g,
    G.gc.cop_sq2(1n+av, 1n+bv, hcop),
    dA2,
    G.gc.dvd_sub(g, Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv), Nat.mul(1n+cv, 1n+cv), hsq, dA2,
      G.gc.dvd_mul_r(g, 1n+cv, 1n+cv, d2)))

def py.copAC(+av: Nat, +bv: Nat, +cv: Nat, +hsq: {Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)) == Nat.mul(1n+cv, 1n+cv) : Nat}, +hcop: G.Cop<1n+av, 1n+bv>) -> G.Cop<1n+av, 1n+cv>:
  G.gc.cop_intro(av, cv, g => d1 => d2 => py.copAC.k(av, bv, cv, hsq, hcop, g, d1, d2))

# gcd(u, v) = 1: a common divisor of u and v divides u + v = C and v - u = A
def py.copUV.k(+av: Nat, +cv: Nat, +up: Nat, +huv: {Nat.add(1n+up, Nat.add(1n+av, 1n+up)) == 1n+cv : Nat}, +cAC: G.Cop<1n+av, 1n+cv>, +g: Nat, d1: G.Dvd<g, 1n+up>, d2: G.Dvd<g, Nat.add(1n+av, 1n+up)>) -> {g == 1n : Nat}:
  +d1c = d1
  +d2c = d2
  G.gc.cop_common(1n+av, 1n+cv, g, cAC,
    G.gc.dvd_sub(g, 1n+up, 1n+av, Nat.add(1n+av, 1n+up), N.Nat.add_comm(1n+up, 1n+av), d1c, d2c),
    G.gc.dvd_eq(g, Nat.add(1n+up, Nat.add(1n+av, 1n+up)), 1n+cv, huv,
      G.gc.dvd_add(g, 1n+up, Nat.add(1n+av, 1n+up), d1c, d2c)))

def py.copUV(+av: Nat, +cv: Nat, +up: Nat, +huv: {Nat.add(1n+up, Nat.add(1n+av, 1n+up)) == 1n+cv : Nat}, +cAC: G.Cop<1n+av, 1n+cv>) -> G.Cop<1n+up, Nat.add(1n+av, 1n+up)>:
  G.gc.cop_intro(up, Nat.add(av, 1n+up), g => d1 => d2 => py.copUV.k(av, cv, up, huv, cAC, g, d1, d2))

# B^2 = 4 u v
def py.bsq(+av: Nat, +bv: Nat, +cv: Nat, +up: Nat, hsq: {Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)) == Nat.mul(1n+cv, 1n+cv) : Nat}, hu: {Nat.add(1n+av, Nat.double(1n+up)) == 1n+cv : Nat}) -> {Nat.mul(1n+bv, 1n+bv) == Nat.double(Nat.double(Nat.mul(1n+up, Nat.add(1n+av, 1n+up)))) : Nat}:
  N.Nat.add_cancel_l(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv), Nat.double(Nat.double(Nat.mul(1n+up, Nat.add(1n+av, 1n+up)))),
    Equal.trans(Nat, Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)), Nat.mul(1n+cv, 1n+cv), Nat.add(Nat.mul(1n+av, 1n+av), Nat.double(Nat.double(Nat.mul(1n+up, Nat.add(1n+av, 1n+up))))),
      hsq,
      Equal.trans(Nat, Nat.mul(1n+cv, 1n+cv), Nat.mul(Nat.add(1n+av, Nat.double(1n+up)), Nat.add(1n+av, Nat.double(1n+up))), Nat.add(Nat.mul(1n+av, 1n+av), Nat.double(Nat.double(Nat.mul(1n+up, Nat.add(1n+av, 1n+up))))),
        Equal.cong(Nat, Nat, z => Nat.mul(z, z), 1n+cv, Nat.add(1n+av, Nat.double(1n+up)),
          Equal.sym(Nat, Nat.add(1n+av, Nat.double(1n+up)), 1n+cv, hu)),
        py.expand(1n+av, 1n+up))))

# k^2 = u v, where B = 2k
def py.kk(+av: Nat, +bv: Nat, +up: Nat, +k: Nat, hk: {1n+bv == Nat.double(k) : Nat}, hb: {Nat.mul(1n+bv, 1n+bv) == Nat.double(Nat.double(Nat.mul(1n+up, Nat.add(1n+av, 1n+up)))) : Nat}) -> {Nat.mul(k, k) == Nat.mul(1n+up, Nat.add(1n+av, 1n+up)) : Nat}:
  N.Nat.double_inj(Nat.mul(k, k), Nat.mul(1n+up, Nat.add(1n+av, 1n+up)),
    N.Nat.double_inj(Nat.double(Nat.mul(k, k)), Nat.double(Nat.mul(1n+up, Nat.add(1n+av, 1n+up))),
      Equal.trans(Nat, Nat.double(Nat.double(Nat.mul(k, k))), Nat.mul(1n+bv, 1n+bv), Nat.double(Nat.double(Nat.mul(1n+up, Nat.add(1n+av, 1n+up)))),
        Equal.trans(Nat, Nat.double(Nat.double(Nat.mul(k, k))), Nat.mul(Nat.double(k), Nat.double(k)), Nat.mul(1n+bv, 1n+bv),
          Equal.sym(Nat, Nat.mul(Nat.double(k), Nat.double(k)), Nat.double(Nat.double(Nat.mul(k, k))), S2.Nat.sq_double(k)),
          Equal.cong(Nat, Nat, z => Nat.mul(z, z), Nat.double(k), 1n+bv, Equal.sym(Nat, 1n+bv, Nat.double(k), hk))),
        hb)))

# reading a square root off IsSq, in 1+ form
def py.sq_pos.go(+target: Nat, s: Nat, +h: {Nat.mul(s, s) == 1n+target : Nat}, -P: Type, k: @sv: Nat -> @hs: {Nat.mul(1n+sv, 1n+sv) == 1n+target : Nat} -> P) -> P:
  match s:
    case 0n:
      Empty.absurd(P, N.Nat.zero_neq_succ(target, h))
    case 1n+sv:
      k(sv, h)

def py.sq_pos(+target: Nat, sq: G.IsSq<1n+target>, -P: Type, k: @sv: Nat -> @hs: {Nat.mul(1n+sv, 1n+sv) == 1n+target : Nat} -> P) -> P:
  match sq:
    case G.IsSq{s, h}:
      py.sq_pos.go(target, s, h, P, k)

# ---------------------------------------------------------------------
# THE PARAMETRISATION
# ---------------------------------------------------------------------

type Pyth<-xa: Nat, -xb: Nat, -xc: Nat> is Data:
  Pyth{mv: Nat, nv: Nat,
       hcop: G.Cop<1n+mv, 1n+nv>,
       ha: {Nat.add(xa, Nat.mul(1n+nv, 1n+nv)) == Nat.mul(1n+mv, 1n+mv) : Nat},
       hb: {xb == Nat.double(Nat.mul(1n+mv, 1n+nv)) : Nat},
       hcc: {xc == Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)) : Nat}}

def py.pyth.fin(+av: Nat, +bv: Nat, +cv: Nat, +up: Nat, +k: Nat, +nv: Nat, +mv: Nat, +hn: {Nat.mul(1n+nv, 1n+nv) == 1n+up : Nat}, +hm: {Nat.mul(1n+mv, 1n+mv) == Nat.add(1n+av, 1n+up) : Nat}, hcuv: G.Cop<1n+up, Nat.add(1n+av, 1n+up)>, hk: {1n+bv == Nat.double(k) : Nat}, hkk: {Nat.mul(k, k) == Nat.mul(1n+up, Nat.add(1n+av, 1n+up)) : Nat}, huv: {Nat.add(1n+up, Nat.add(1n+av, 1n+up)) == 1n+cv : Nat}) -> Pyth<1n+av, 1n+bv, 1n+cv>:
  Pyth{mv, nv,
    G.gc.cop_of_sq(1n+mv, 1n+nv,
      G.gc.cop_sym(Nat.mul(1n+nv, 1n+nv), Nat.mul(1n+mv, 1n+mv),
        G.gc.cop_eq_r(Nat.mul(1n+nv, 1n+nv), Nat.add(1n+av, 1n+up), Nat.mul(1n+mv, 1n+mv),
          Equal.sym(Nat, Nat.mul(1n+mv, 1n+mv), Nat.add(1n+av, 1n+up), hm),
          G.gc.cop_eq_l(1n+up, Nat.mul(1n+nv, 1n+nv), Nat.add(1n+av, 1n+up),
            Equal.sym(Nat, Nat.mul(1n+nv, 1n+nv), 1n+up, hn),
            hcuv)))),
    Equal.trans(Nat, Nat.add(1n+av, Nat.mul(1n+nv, 1n+nv)), Nat.add(1n+av, 1n+up), Nat.mul(1n+mv, 1n+mv),
      Equal.cong(Nat, Nat, z => Nat.add(1n+av, z), Nat.mul(1n+nv, 1n+nv), 1n+up, hn),
      Equal.sym(Nat, Nat.mul(1n+mv, 1n+mv), Nat.add(1n+av, 1n+up), hm)),
    Equal.trans(Nat, 1n+bv, Nat.double(k), Nat.double(Nat.mul(1n+mv, 1n+nv)),
      hk,
      Equal.cong(Nat, Nat, z => Nat.double(z), k, Nat.mul(1n+mv, 1n+nv),
        G.gc.sq_inj(k, Nat.mul(1n+mv, 1n+nv),
          Equal.trans(Nat, Nat.mul(k, k), Nat.mul(1n+up, Nat.add(1n+av, 1n+up)), Nat.mul(Nat.mul(1n+mv, 1n+nv), Nat.mul(1n+mv, 1n+nv)),
            hkk,
            Equal.trans(Nat, Nat.mul(1n+up, Nat.add(1n+av, 1n+up)), Nat.mul(Nat.mul(1n+nv, 1n+nv), Nat.add(1n+av, 1n+up)), Nat.mul(Nat.mul(1n+mv, 1n+nv), Nat.mul(1n+mv, 1n+nv)),
              Equal.cong(Nat, Nat, z => Nat.mul(z, Nat.add(1n+av, 1n+up)), 1n+up, Nat.mul(1n+nv, 1n+nv),
                Equal.sym(Nat, Nat.mul(1n+nv, 1n+nv), 1n+up, hn)),
              Equal.trans(Nat, Nat.mul(Nat.mul(1n+nv, 1n+nv), Nat.add(1n+av, 1n+up)), Nat.mul(Nat.mul(1n+nv, 1n+nv), Nat.mul(1n+mv, 1n+mv)), Nat.mul(Nat.mul(1n+mv, 1n+nv), Nat.mul(1n+mv, 1n+nv)),
                Equal.cong(Nat, Nat, z => Nat.mul(Nat.mul(1n+nv, 1n+nv), z), Nat.add(1n+av, 1n+up), Nat.mul(1n+mv, 1n+mv),
                  Equal.sym(Nat, Nat.mul(1n+mv, 1n+mv), Nat.add(1n+av, 1n+up), hm)),
                Equal.trans(Nat, Nat.mul(Nat.mul(1n+nv, 1n+nv), Nat.mul(1n+mv, 1n+mv)), Nat.mul(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)), Nat.mul(Nat.mul(1n+mv, 1n+nv), Nat.mul(1n+mv, 1n+nv)),
                  N.Nat.mul_comm(Nat.mul(1n+nv, 1n+nv), Nat.mul(1n+mv, 1n+mv)),
                  Equal.sym(Nat, Nat.mul(Nat.mul(1n+mv, 1n+nv), Nat.mul(1n+mv, 1n+nv)), Nat.mul(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)), G.gc.nsq_prod(1n+mv, 1n+nv))))))))),
    Equal.trans(Nat, 1n+cv, Nat.add(1n+up, Nat.add(1n+av, 1n+up)), Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)),
      Equal.sym(Nat, Nat.add(1n+up, Nat.add(1n+av, 1n+up)), 1n+cv, huv),
      Equal.trans(Nat, Nat.add(1n+up, Nat.add(1n+av, 1n+up)), Nat.add(Nat.mul(1n+nv, 1n+nv), Nat.add(1n+av, 1n+up)), Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)),
        Equal.cong(Nat, Nat, z => Nat.add(z, Nat.add(1n+av, 1n+up)), 1n+up, Nat.mul(1n+nv, 1n+nv),
          Equal.sym(Nat, Nat.mul(1n+nv, 1n+nv), 1n+up, hn)),
        Equal.trans(Nat, Nat.add(Nat.mul(1n+nv, 1n+nv), Nat.add(1n+av, 1n+up)), Nat.add(Nat.mul(1n+nv, 1n+nv), Nat.mul(1n+mv, 1n+mv)), Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)),
          Equal.cong(Nat, Nat, z => Nat.add(Nat.mul(1n+nv, 1n+nv), z), Nat.add(1n+av, 1n+up), Nat.mul(1n+mv, 1n+mv),
            Equal.sym(Nat, Nat.mul(1n+mv, 1n+mv), Nat.add(1n+av, 1n+up), hm)),
          N.Nat.add_comm(Nat.mul(1n+nv, 1n+nv), Nat.mul(1n+mv, 1n+mv)))))}

# A^2 + B^2 = C^2 with A, B, C > 0, gcd(A, B) = 1 and B even
def py.pyth(+av: Nat, +bv: Nat, +cv: Nat, +hsq: {Nat.add(Nat.mul(1n+av, 1n+av), Nat.mul(1n+bv, 1n+bv)) == Nat.mul(1n+cv, 1n+cv) : Nat}, +hcop: G.Cop<1n+av, 1n+bv>, +hbe: {N.Nat.even(1n+bv) == True{} : Bool}) -> Pyth<1n+av, 1n+bv, 1n+cv>:
  +exA = py.oddA(av, bv, hcop, hbe)
  py.get_u(av, cv,
    py.lt_of_sq(1n+av, 1n+cv, py.sq_lt(av, bv, cv, hsq)),
    exA,
    py.oddC(av, bv, cv, hsq, exA, hbe),
    Pyth<1n+av, 1n+bv, 1n+cv>, up => hu =>
      +up2 = up
      +hu2 = hu
      +huv = Equal.trans(Nat, Nat.add(1n+up2, Nat.add(1n+av, 1n+up2)), Nat.add(1n+av, Nat.add(1n+up2, 1n+up2)), 1n+cv,
        N.Nat.add_swap(1n+up2, 1n+av, 1n+up2),
        Equal.trans(Nat, Nat.add(1n+av, Nat.add(1n+up2, 1n+up2)), Nat.add(1n+av, Nat.double(1n+up2)), 1n+cv,
          Equal.cong(Nat, Nat, z => Nat.add(1n+av, z), Nat.add(1n+up2, 1n+up2), Nat.double(1n+up2),
            Equal.sym(Nat, Nat.double(1n+up2), Nat.add(1n+up2, 1n+up2), N.Nat.double_add(1n+up2))),
          hu2))
      +hcuv = py.copUV(av, cv, up2, huv, py.copAC(av, bv, cv, hsq, hcop))
      py.even_shape(1n+bv, hbe, Pyth<1n+av, 1n+bv, 1n+cv>, k => hk =>
        +k2 = k
        +hk2 = hk
        +hkk = py.kk(av, bv, up2, k2, hk2, py.bsq(av, bv, cv, up2, hsq, hu2))
        py.sq_pos(up2,
          G.gc.sq_of_cop(up2, Nat.add(1n+av, 1n+up2), k2, hcuv,
            Equal.sym(Nat, Nat.mul(k2, k2), Nat.mul(1n+up2, Nat.add(1n+av, 1n+up2)), hkk)),
          Pyth<1n+av, 1n+bv, 1n+cv>, nv => hn =>
            py.sq_pos(Nat.add(av, 1n+up2),
              G.gc.sq_of_cop(Nat.add(av, 1n+up2), 1n+up2, k2,
                G.gc.cop_sym(1n+up2, Nat.add(1n+av, 1n+up2), hcuv),
                Equal.trans(Nat, Nat.mul(Nat.add(1n+av, 1n+up2), 1n+up2), Nat.mul(1n+up2, Nat.add(1n+av, 1n+up2)), Nat.mul(k2, k2),
                  N.Nat.mul_comm(Nat.add(1n+av, 1n+up2), 1n+up2),
                  Equal.sym(Nat, Nat.mul(k2, k2), Nat.mul(1n+up2, Nat.add(1n+av, 1n+up2)), hkk))),
              Pyth<1n+av, 1n+bv, 1n+cv>, mv => hm =>
                py.pyth.fin(av, bv, cv, up2, k2, nv, mv, hn, hm, hcuv, hk2, hkk, huv)))))

# the same theorem with positivity as a hypothesis rather than a shape,
# which is how the descent meets it
def py.pyth_pos(a: Nat, pa: N.Nat.Lt(0n, a), b: Nat, pb: N.Nat.Lt(0n, b), c: Nat, pc: N.Nat.Lt(0n, c), +hsq: {Nat.add(Nat.mul(a, a), Nat.mul(b, b)) == Nat.mul(c, c) : Nat}, +hcop: G.Cop<a, b>, +hbe: {N.Nat.even(b) == True{} : Bool}) -> Pyth<a, b, c>:
  match a:
    case 0n:
      match pa:
    case 1n+av:
      match b:
        case 0n:
          match pb:
        case 1n+bv:
          match c:
            case 0n:
              match pc:
            case 1n+cv:
              py.pyth(av, bv, cv, hsq, hcop, hbe)