# demos/math/fermat4.bend -- Fermat's Last Theorem for exponent 4.
#
#     x^4 + y^4 = z^4 has no solution in positive naturals
#
# This is Fermat's own case, and the only one he left a proof of. It is
# proved here as he proved it, by infinite descent, on the stronger
# statement that x^4 + y^4 is never a perfect SQUARE.
#
# THE ARGUMENT, in words.
#
# Suppose x^4 + y^4 = z^2 with x, y, z > 0 and gcd(x, y) = 1, and take a
# solution with z as small as possible (here: run a fuel down from z).
# Then (x^2)^2 + (y^2)^2 = z^2 is a PRIMITIVE Pythagorean triple. One of
# x, y is even -- not both, since they are coprime, and not neither, since
# odd^2 + odd^2 is 2 mod 4 and never a square (lib/pythag.bend's
# py.odd_odd) -- so name them with y even. The parametrisation of
# primitive triples (lib/pythag.bend's py.pyth) gives coprime m > n > 0
# with
#
#     x^2 = m^2 - n^2,    y^2 = 2 m n,    z = m^2 + n^2.
#
# Now x^2 + n^2 = m^2 is a second primitive triple. Its odd leg is x, so
# n is the even one -- again because odd^2 + odd^2 is no square -- and
# x is coprime to n because a common divisor would divide m^2 and n^2.
# Parametrising it gives coprime r, s > 0 with
#
#     x = r^2 - s^2,      n = 2 r s,      m = r^2 + s^2.
#
# Then y^2 = 2 m n = 4 r s m, and y is even, say y = 2Y, so
#
#     Y^2 = r * s * m
#
# with r, s, m pairwise coprime (r and s by construction, and each of them
# with m = r^2 + s^2 because a common divisor would divide the other
# square). A product of pairwise coprime numbers is a square only if each
# factor is (lib/gcd.bend's gc.sq_of_cop), so
#
#     r = a^2,   s = b^2,   m = c^2,
#
# and therefore a^4 + b^4 = r^2 + s^2 = m = c^2: a NEW solution. It is
# smaller: c <= c^2 = m <= m^2 < m^2 + n^2 = z. Descent, contradiction.
#
# The general case then follows: divide x and y by their gcd G. Since
# G^4 divides z^4 and (G^2)^2 | (z^2)^2 forces G^2 | z^2 and then G | z
# (lib/gcd.bend's gc.sq_dvd), the reduced triple is again a solution.
#
# The statement at the bottom is in the vocabulary of demos/math/mathlib.bend,
# so the website's math card can quote it directly.

import Base
import ./mathlib.bend as M
import ./nat.bend as N
import ./sqrt2.bend as S2
import ./gcd.bend as G
import ./pythag.bend as PY

# ---------------------------------------------------------------------
# mathlib's add / mul / pow are Base's, definition for definition
# ---------------------------------------------------------------------

def f4.add_bridge(a: Nat, -b: Nat) -> {M.add(a, b) == Nat.add(a, b) : Nat}:
  match a:
    case 0n:
      {==}
    case 1n+p:
      Equal.cong(Nat, Nat, x => 1n+x, M.add(p, b), Nat.add(p, b), f4.add_bridge(p, b))

def f4.mul_bridge(a: Nat, +b: Nat) -> {M.mul(a, b) == Nat.mul(a, b) : Nat}:
  match a:
    case 0n:
      {==}
    case 1n+p:
      +p2 = p
      %Equal.sym(Nat, M.mul(p2, b), Nat.mul(p2, b), f4.mul_bridge(p2, b)) : {M.add(b, _) == Nat.add(b, Nat.mul(p2, b)) : Nat}
      f4.add_bridge(b, Nat.mul(p2, b))

# a^4, mathlib's way and Base's, are the same number
def f4.pow4_bridge(+a: Nat) -> {M.pow(a, 4n) == Nat.mul(Nat.mul(a, a), Nat.mul(a, a)) : Nat}:
  Equal.trans(Nat, M.pow(a, 4n), Nat.mul(a, Nat.mul(a, Nat.mul(a, Nat.mul(a, 1n)))), Nat.mul(Nat.mul(a, a), Nat.mul(a, a)),
    Equal.trans(Nat, M.pow(a, 4n), Nat.mul(a, M.pow(a, 3n)), Nat.mul(a, Nat.mul(a, Nat.mul(a, Nat.mul(a, 1n)))),
      f4.mul_bridge(a, M.pow(a, 3n)),
      Equal.cong(Nat, Nat, z => Nat.mul(a, z), M.pow(a, 3n), Nat.mul(a, Nat.mul(a, Nat.mul(a, 1n))),
        Equal.trans(Nat, M.pow(a, 3n), Nat.mul(a, M.pow(a, 2n)), Nat.mul(a, Nat.mul(a, Nat.mul(a, 1n))),
          f4.mul_bridge(a, M.pow(a, 2n)),
          Equal.cong(Nat, Nat, z => Nat.mul(a, z), M.pow(a, 2n), Nat.mul(a, Nat.mul(a, 1n)),
            Equal.trans(Nat, M.pow(a, 2n), Nat.mul(a, M.pow(a, 1n)), Nat.mul(a, Nat.mul(a, 1n)),
              f4.mul_bridge(a, M.pow(a, 1n)),
              Equal.cong(Nat, Nat, z => Nat.mul(a, z), M.pow(a, 1n), Nat.mul(a, 1n),
                f4.mul_bridge(a, 1n))))))),
    Equal.trans(Nat, Nat.mul(a, Nat.mul(a, Nat.mul(a, Nat.mul(a, 1n)))), Nat.mul(a, Nat.mul(a, Nat.mul(a, a))), Nat.mul(Nat.mul(a, a), Nat.mul(a, a)),
      Equal.cong(Nat, Nat, z => Nat.mul(a, Nat.mul(a, Nat.mul(a, z))), Nat.mul(a, 1n), a, N.Nat.mul_one(a)),
      N.Nat.mul_assoc(a, a, Nat.mul(a, a))))

# the Bool test is faithful: a proof that the equation is impossible makes
# nat_eq return False
def f4.nat_eq_false(a: Nat, b: Nat, ne: {a == b : Nat} -> Empty) -> {M.nat_eq(a, b) == False{} : Bool}:
  match a b:
    case 0n 0n:
      Empty.absurd({True{} == False{} : Bool}, ne({==}))
    case 0n 1n+q:
      {==}
    case 1n+p 0n:
      {==}
    case 1n+p 1n+q:
      f4.nat_eq_false(p, q, e => ne(Equal.cong(Nat, Nat, x => 1n+x, p, q, e)))

# ---------------------------------------------------------------------
# the descent
#
# f4.Rec(g) is the induction hypothesis: any coprime solution whose z is
# at most g is impossible. The step below produces one from a solution
# whose z is at most 1 + g, so the recursion is on the fuel and every
# step is one application of the parametrisation, twice.
# ---------------------------------------------------------------------

def f4.Rec(g: Nat) -> Type:
  @xv2: Nat -> @yv2: Nat -> @zv2: Nat -> @b2: N.Nat.Le(zv2, g) -> @c2: G.Cop<1n+xv2, 1n+yv2> ->
  @e2: {Nat.add(Nat.mul(Nat.mul(1n+xv2, 1n+xv2), Nat.mul(1n+xv2, 1n+xv2)), Nat.mul(Nat.mul(1n+yv2, 1n+yv2), Nat.mul(1n+yv2, 1n+yv2))) == Nat.mul(1n+zv2, 1n+zv2) : Nat} -> Empty

# p (q t) = q (t p)
def f4.rot3(+p: Nat, +q: Nat, +t: Nat) -> {Nat.mul(p, Nat.mul(q, t)) == Nat.mul(q, Nat.mul(t, p)) : Nat}:
  Equal.trans(Nat, Nat.mul(p, Nat.mul(q, t)), Nat.mul(q, Nat.mul(p, t)), Nat.mul(q, Nat.mul(t, p)),
    G.gc.nmul_swap(p, q, t),
    Equal.cong(Nat, Nat, z => Nat.mul(q, z), Nat.mul(p, t), Nat.mul(t, p), N.Nat.mul_comm(p, t)))

# THE LAST STEP: a^4 + b^4 = c^2 with c < z, handed to the induction
# hypothesis. c <= c^2 = m <= m^2 < m^2 + n^2 = z.
def f4.s4(+g: Nat, +zv: Nat, +mv: Nat, +nv: Nat, +rv: Nat, +sv: Nat, +av: Nat, +bv: Nat, +cv: Nat, bound: N.Nat.Le(zv, 1n+g), hc1: {1n+zv == Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)) : Nat}, hc2: {1n+mv == Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)) : Nat}, ha: {Nat.mul(1n+av, 1n+av) == 1n+rv : Nat}, hb: {Nat.mul(1n+bv, 1n+bv) == 1n+sv : Nat}, +hcc: {Nat.mul(1n+cv, 1n+cv) == 1n+mv : Nat}, hcopab: G.Cop<1n+av, 1n+bv>, rec: f4.Rec(g)) -> Empty:
  rec(av, bv, cv,
    N.Nat.le_trans(1n+cv, zv, 1n+g,
      N.Nat.le_trans(2n+cv, 1n+Nat.mul(1n+mv, 1n+mv), 1n+zv,
        N.Nat.le_trans(1n+cv, 1n+mv, Nat.mul(1n+mv, 1n+mv),
          N.Nat.le_trans(1n+cv, Nat.mul(1n+cv, 1n+cv), 1n+mv,
            N.Nat.le_add_r(cv, Nat.mul(cv, 1n+cv)),
            N.Nat.eq_le(Nat.mul(1n+cv, 1n+cv), 1n+mv, hcc)),
          N.Nat.le_add_r(mv, Nat.mul(mv, 1n+mv))),
        %Equal.sym(Nat, 1n+zv, Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)), hc1) : N.Nat.Lt(Nat.mul(1n+mv, 1n+mv), _)
        PY.py.lt_add_pos(Nat.mul(1n+mv, 1n+mv), Nat.add(nv, Nat.mul(nv, 1n+nv)))),
      bound),
    hcopab,
    Equal.trans(Nat, Nat.add(Nat.mul(Nat.mul(1n+av, 1n+av), Nat.mul(1n+av, 1n+av)), Nat.mul(Nat.mul(1n+bv, 1n+bv), Nat.mul(1n+bv, 1n+bv))), Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(Nat.mul(1n+bv, 1n+bv), Nat.mul(1n+bv, 1n+bv))), Nat.mul(1n+cv, 1n+cv),
      Equal.cong(Nat, Nat, z => Nat.add(Nat.mul(z, z), Nat.mul(Nat.mul(1n+bv, 1n+bv), Nat.mul(1n+bv, 1n+bv))), Nat.mul(1n+av, 1n+av), 1n+rv, ha),
      Equal.trans(Nat, Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(Nat.mul(1n+bv, 1n+bv), Nat.mul(1n+bv, 1n+bv))), Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)), Nat.mul(1n+cv, 1n+cv),
        Equal.cong(Nat, Nat, z => Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(z, z)), Nat.mul(1n+bv, 1n+bv), 1n+sv, hb),
        Equal.trans(Nat, Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)), 1n+mv, Nat.mul(1n+cv, 1n+cv),
          Equal.sym(Nat, 1n+mv, Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)), hc2),
          Equal.sym(Nat, Nat.mul(1n+cv, 1n+cv), 1n+mv, hcc)))))

# r and m = r^2 + s^2 are coprime: a common divisor divides r^2, hence s^2
def f4.copRM.k(+rv: Nat, +sv: Nat, +mv: Nat, +hc2: {1n+mv == Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)) : Nat}, +hcop_rs: G.Cop<1n+rv, 1n+sv>, +gg: Nat, d1: G.Dvd<gg, 1n+rv>, d2: G.Dvd<gg, 1n+mv>) -> {gg == 1n : Nat}:
  +dr2 = G.gc.dvd_mul_r(gg, 1n+rv, 1n+rv, d1)
  G.gc.cop_common(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv), gg,
    G.gc.cop_sq2(1n+rv, 1n+sv, hcop_rs),
    dr2,
    G.gc.dvd_sub(gg, Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv), 1n+mv,
      Equal.sym(Nat, 1n+mv, Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)), hc2),
      dr2, d2))

def f4.copRM(+rv: Nat, +sv: Nat, +mv: Nat, +hc2: {1n+mv == Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)) : Nat}, +hcop_rs: G.Cop<1n+rv, 1n+sv>) -> G.Cop<1n+rv, 1n+mv>:
  G.gc.cop_intro(rv, mv, gg => d1 => d2 => f4.copRM.k(rv, sv, mv, hc2, hcop_rs, gg, d1, d2))

def f4.copSM.k(+rv: Nat, +sv: Nat, +mv: Nat, +hc2: {1n+mv == Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)) : Nat}, +hcop_rs: G.Cop<1n+rv, 1n+sv>, +gg: Nat, d1: G.Dvd<gg, 1n+sv>, d2: G.Dvd<gg, 1n+mv>) -> {gg == 1n : Nat}:
  +ds2 = G.gc.dvd_mul_r(gg, 1n+sv, 1n+sv, d1)
  G.gc.cop_common(Nat.mul(1n+sv, 1n+sv), Nat.mul(1n+rv, 1n+rv), gg,
    G.gc.cop_sq2(1n+sv, 1n+rv, G.gc.cop_sym(1n+rv, 1n+sv, hcop_rs)),
    ds2,
    G.gc.dvd_sub(gg, Nat.mul(1n+sv, 1n+sv), Nat.mul(1n+rv, 1n+rv), 1n+mv,
      Equal.trans(Nat, Nat.add(Nat.mul(1n+sv, 1n+sv), Nat.mul(1n+rv, 1n+rv)), Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)), 1n+mv,
        N.Nat.add_comm(Nat.mul(1n+sv, 1n+sv), Nat.mul(1n+rv, 1n+rv)),
        Equal.sym(Nat, 1n+mv, Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)), hc2)),
      ds2, d2))

def f4.copSM(+rv: Nat, +sv: Nat, +mv: Nat, +hc2: {1n+mv == Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)) : Nat}, +hcop_rs: G.Cop<1n+rv, 1n+sv>) -> G.Cop<1n+sv, 1n+mv>:
  G.gc.cop_intro(sv, mv, gg => d1 => d2 => f4.copSM.k(rv, sv, mv, hc2, hcop_rs, gg, d1, d2))

# Y^2 = r s m with r, s, m pairwise coprime: each is a square
def f4.s3b(+g: Nat, +zv: Nat, +mv: Nat, +nv: Nat, +rv: Nat, +sv: Nat, +yy: Nat, bound: N.Nat.Le(zv, 1n+g), hc1: {1n+zv == Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)) : Nat}, +hc2: {1n+mv == Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)) : Nat}, +hcop_rs: G.Cop<1n+rv, 1n+sv>, +key: {Nat.mul(yy, yy) == Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)) : Nat}, rec: f4.Rec(g)) -> Empty:
  +crm = f4.copRM(rv, sv, mv, hc2, hcop_rs)
  +csm = f4.copSM(rv, sv, mv, hc2, hcop_rs)
  PY.py.sq_pos(rv,
    G.gc.sq_of_cop(rv, Nat.mul(1n+sv, 1n+mv), yy,
      G.gc.cop_mul(1n+rv, 1n+sv, 1n+mv, hcop_rs, crm),
      Equal.sym(Nat, Nat.mul(yy, yy), Nat.mul(1n+rv, Nat.mul(1n+sv, 1n+mv)),
        Equal.trans(Nat, Nat.mul(yy, yy), Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)), Nat.mul(1n+rv, Nat.mul(1n+sv, 1n+mv)),
          key, f4.rot3(1n+mv, 1n+rv, 1n+sv)))),
    Empty, av => ha =>
      PY.py.sq_pos(sv,
        G.gc.sq_of_cop(sv, Nat.mul(1n+rv, 1n+mv), yy,
          G.gc.cop_mul(1n+sv, 1n+rv, 1n+mv, G.gc.cop_sym(1n+rv, 1n+sv, hcop_rs), csm),
          Equal.sym(Nat, Nat.mul(yy, yy), Nat.mul(1n+sv, Nat.mul(1n+rv, 1n+mv)),
            Equal.trans(Nat, Nat.mul(yy, yy), Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)), Nat.mul(1n+sv, Nat.mul(1n+rv, 1n+mv)),
              key,
              Equal.trans(Nat, Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)), Nat.mul(1n+mv, Nat.mul(1n+sv, 1n+rv)), Nat.mul(1n+sv, Nat.mul(1n+rv, 1n+mv)),
                Equal.cong(Nat, Nat, z => Nat.mul(1n+mv, z), Nat.mul(1n+rv, 1n+sv), Nat.mul(1n+sv, 1n+rv), N.Nat.mul_comm(1n+rv, 1n+sv)),
                f4.rot3(1n+mv, 1n+sv, 1n+rv))))),
        Empty, bv => hb =>
          PY.py.sq_pos(mv,
            G.gc.sq_of_cop(mv, Nat.mul(1n+rv, 1n+sv), yy,
              G.gc.cop_mul(1n+mv, 1n+rv, 1n+sv, G.gc.cop_sym(1n+rv, 1n+mv, crm), G.gc.cop_sym(1n+sv, 1n+mv, csm)),
              Equal.sym(Nat, Nat.mul(yy, yy), Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)), key)),
            Empty, cv => hcc =>
              +ha2 = ha
              +hb2 = hb
              +av2 = av
              +bv2 = bv
              f4.s4(g, zv, mv, nv, rv, sv, av2, bv2, cv, bound, hc1, hc2, ha2, hb2, hcc,
                G.gc.cop_of_sq(1n+av2, 1n+bv2,
                  G.gc.cop_eq_r(Nat.mul(1n+av2, 1n+av2), 1n+sv, Nat.mul(1n+bv2, 1n+bv2),
                    Equal.sym(Nat, Nat.mul(1n+bv2, 1n+bv2), 1n+sv, hb2),
                    G.gc.cop_eq_l(1n+rv, Nat.mul(1n+av2, 1n+av2), 1n+sv,
                      Equal.sym(Nat, Nat.mul(1n+av2, 1n+av2), 1n+rv, ha2),
                      hcop_rs))),
                rec))))

# y = 2Y, and Y^2 = m (r s)
def f4.s3a(+g: Nat, +yv: Nat, +zv: Nat, +mv: Nat, +nv: Nat, +rv: Nat, +sv: Nat, bound: N.Nat.Le(zv, 1n+g), hc1: {1n+zv == Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)) : Nat}, hc2: {1n+mv == Nat.add(Nat.mul(1n+rv, 1n+rv), Nat.mul(1n+sv, 1n+sv)) : Nat}, hb1: {Nat.mul(1n+yv, 1n+yv) == Nat.double(Nat.mul(1n+mv, 1n+nv)) : Nat}, hb2: {1n+nv == Nat.double(Nat.mul(1n+rv, 1n+sv)) : Nat}, hcop_rs: G.Cop<1n+rv, 1n+sv>, hyeven: {N.Nat.even(1n+yv) == True{} : Bool}, rec: f4.Rec(g)) -> Empty:
  PY.py.even_shape(1n+yv, hyeven, Empty, yy => hY =>
    +yy2 = yy
    f4.s3b(g, zv, mv, nv, rv, sv, yy2, bound, hc1, hc2, hcop_rs,
      N.Nat.double_inj(Nat.mul(yy2, yy2), Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)),
        N.Nat.double_inj(Nat.double(Nat.mul(yy2, yy2)), Nat.double(Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv))),
          Equal.trans(Nat, Nat.double(Nat.double(Nat.mul(yy2, yy2))), Nat.mul(Nat.double(yy2), Nat.double(yy2)), Nat.double(Nat.double(Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)))),
            Equal.sym(Nat, Nat.mul(Nat.double(yy2), Nat.double(yy2)), Nat.double(Nat.double(Nat.mul(yy2, yy2))), S2.Nat.sq_double(yy2)),
            Equal.trans(Nat, Nat.mul(Nat.double(yy2), Nat.double(yy2)), Nat.mul(1n+yv, 1n+yv), Nat.double(Nat.double(Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)))),
              Equal.cong(Nat, Nat, z => Nat.mul(z, z), Nat.double(yy2), 1n+yv, Equal.sym(Nat, 1n+yv, Nat.double(yy2), hY)),
              Equal.trans(Nat, Nat.mul(1n+yv, 1n+yv), Nat.double(Nat.mul(1n+mv, 1n+nv)), Nat.double(Nat.double(Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)))),
                hb1,
                Equal.trans(Nat, Nat.double(Nat.mul(1n+mv, 1n+nv)), Nat.double(Nat.mul(1n+mv, Nat.double(Nat.mul(1n+rv, 1n+sv)))), Nat.double(Nat.double(Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv)))),
                  Equal.cong(Nat, Nat, z => Nat.double(Nat.mul(1n+mv, z)), 1n+nv, Nat.double(Nat.mul(1n+rv, 1n+sv)), hb2),
                  Equal.cong(Nat, Nat, z => Nat.double(z), Nat.mul(1n+mv, Nat.double(Nat.mul(1n+rv, 1n+sv))), Nat.double(Nat.mul(1n+mv, Nat.mul(1n+rv, 1n+sv))),
                    PY.py.mul_double(1n+mv, Nat.mul(1n+rv, 1n+sv))))))))),
      rec))

# the second triple x^2 + n^2 = m^2: n is even, x is coprime to n
def f4.n1even(+xv: Nat, +nv: Nat, +mv: Nat, hxodd: {N.Nat.even(1n+xv) == False{} : Bool}, ha1: {Nat.add(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+nv, 1n+nv)) == Nat.mul(1n+mv, 1n+mv) : Nat}, c: Bool, ec: {N.Nat.even(1n+nv) == c : Bool}) -> {N.Nat.even(1n+nv) == True{} : Bool}:
  match c:
    case True{}:
      ec
    case False{}:
      Empty.absurd({N.Nat.even(1n+nv) == True{} : Bool},
        PY.py.odd_odd(1n+xv, 1n+nv, 1n+mv, hxodd, ec, ha1))

def f4.copXN.k(+xv: Nat, +nv: Nat, +mv: Nat, +ha1: {Nat.add(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+nv, 1n+nv)) == Nat.mul(1n+mv, 1n+mv) : Nat}, +hcop_mn: G.Cop<1n+mv, 1n+nv>, +gg: Nat, d1: G.Dvd<gg, 1n+xv>, d2: G.Dvd<gg, 1n+nv>) -> {gg == 1n : Nat}:
  +dn2 = G.gc.dvd_mul_r(gg, 1n+nv, 1n+nv, d2)
  G.gc.cop_common(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv), gg,
    G.gc.cop_sq2(1n+mv, 1n+nv, hcop_mn),
    G.gc.dvd_eq(gg, Nat.add(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+nv, 1n+nv)), Nat.mul(1n+mv, 1n+mv), ha1,
      G.gc.dvd_add(gg, Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+nv, 1n+nv), G.gc.dvd_mul_r(gg, 1n+xv, 1n+xv, d1), dn2)),
    dn2)

def f4.copXN(+xv: Nat, +nv: Nat, +mv: Nat, +ha1: {Nat.add(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+nv, 1n+nv)) == Nat.mul(1n+mv, 1n+mv) : Nat}, +hcop_mn: G.Cop<1n+mv, 1n+nv>) -> G.Cop<1n+xv, 1n+nv>:
  G.gc.cop_intro(xv, nv, gg => d1 => d2 => f4.copXN.k(xv, nv, mv, ha1, hcop_mn, gg, d1, d2))

def f4.s2b(+g: Nat, +xv: Nat, +yv: Nat, +zv: Nat, +mv: Nat, +nv: Nat, bound: N.Nat.Le(zv, 1n+g), hb1: {Nat.mul(1n+yv, 1n+yv) == Nat.double(Nat.mul(1n+mv, 1n+nv)) : Nat}, hc1: {1n+zv == Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)) : Nat}, hyeven: {N.Nat.even(1n+yv) == True{} : Bool}, rec: f4.Rec(g), pt: PY.Pyth<1n+xv, 1n+nv, 1n+mv>) -> Empty:
  match pt:
    case PY.Pyth{rv, sv, hcop_rs, ha2, hb2, hc2}:
      f4.s3a(g, yv, zv, mv, nv, rv, sv, bound, hc1, hc2, hb1, hb2, hcop_rs, hyeven, rec)

def f4.s2(+g: Nat, +xv: Nat, +yv: Nat, +zv: Nat, +mv: Nat, +nv: Nat, bound: N.Nat.Le(zv, 1n+g), hcop_mn: G.Cop<1n+mv, 1n+nv>, +ha1: {Nat.add(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+nv, 1n+nv)) == Nat.mul(1n+mv, 1n+mv) : Nat}, hb1: {Nat.mul(1n+yv, 1n+yv) == Nat.double(Nat.mul(1n+mv, 1n+nv)) : Nat}, hc1: {1n+zv == Nat.add(Nat.mul(1n+mv, 1n+mv), Nat.mul(1n+nv, 1n+nv)) : Nat}, hxodd: {N.Nat.even(1n+xv) == False{} : Bool}, hyeven: {N.Nat.even(1n+yv) == True{} : Bool}, rec: f4.Rec(g)) -> Empty:
  f4.s2b(g, xv, yv, zv, mv, nv, bound, hb1, hc1, hyeven, rec,
    PY.py.pyth_pos(1n+xv, Unit{}, 1n+nv, Unit{}, 1n+mv, Unit{}, ha1,
      f4.copXN(xv, nv, mv, ha1, hcop_mn),
      f4.n1even(xv, nv, mv, hxodd, ha1, N.Nat.even(1n+nv), {==})))

def f4.s1b(+g: Nat, +xv: Nat, +yv: Nat, +zv: Nat, bound: N.Nat.Le(zv, 1n+g), hxodd: {N.Nat.even(1n+xv) == False{} : Bool}, hyeven: {N.Nat.even(1n+yv) == True{} : Bool}, rec: f4.Rec(g), pt: PY.Pyth<Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+yv, 1n+yv), 1n+zv>) -> Empty:
  match pt:
    case PY.Pyth{mv, nv, hcop_mn, ha1, hb1, hc1}:
      f4.s2(g, xv, yv, zv, mv, nv, bound, hcop_mn, ha1, hb1, hc1, hxodd, hyeven, rec)

# THE STEP: the first parametrisation
def f4.s1(+g: Nat, +xv: Nat, +yv: Nat, +zv: Nat, bound: N.Nat.Le(zv, 1n+g), hcop: G.Cop<1n+xv, 1n+yv>, hxodd: {N.Nat.even(1n+xv) == False{} : Bool}, +hyeven: {N.Nat.even(1n+yv) == True{} : Bool}, e: {Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))) == Nat.mul(1n+zv, 1n+zv) : Nat}, rec: f4.Rec(g)) -> Empty:
  f4.s1b(g, xv, yv, zv, bound, hxodd, hyeven, rec,
    PY.py.pyth_pos(Nat.mul(1n+xv, 1n+xv), Unit{}, Nat.mul(1n+yv, 1n+yv), Unit{}, 1n+zv, Unit{}, e,
      G.gc.cop_sq2(1n+xv, 1n+yv, hcop),
      Equal.trans(Bool, N.Nat.even(Nat.mul(1n+yv, 1n+yv)), N.Nat.even(1n+yv), True{}, N.Nat.even_sq(1n+yv), hyeven)))

# exactly one of x, y is even: both even contradicts coprimality, and
# neither even contradicts py.odd_odd applied to (x^2)^2 + (y^2)^2 = z^2
def f4.yeven(+xv: Nat, +yv: Nat, +zv: Nat, +hxodd: {N.Nat.even(1n+xv) == False{} : Bool}, e: {Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))) == Nat.mul(1n+zv, 1n+zv) : Nat}, c: Bool, ec: {N.Nat.even(1n+yv) == c : Bool}) -> {N.Nat.even(1n+yv) == True{} : Bool}:
  match c:
    case True{}:
      ec
    case False{}:
      Empty.absurd({N.Nat.even(1n+yv) == True{} : Bool},
        PY.py.odd_odd(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+yv, 1n+yv), 1n+zv,
          Equal.trans(Bool, N.Nat.even(Nat.mul(1n+xv, 1n+xv)), N.Nat.even(1n+xv), False{}, N.Nat.even_sq(1n+xv), hxodd),
          Equal.trans(Bool, N.Nat.even(Nat.mul(1n+yv, 1n+yv)), N.Nat.even(1n+yv), False{}, N.Nat.even_sq(1n+yv), ec),
          e))

def f4.yodd(+xv: Nat, +yv: Nat, hcop: G.Cop<1n+xv, 1n+yv>, +hxe: {N.Nat.even(1n+xv) == True{} : Bool}, c: Bool, ec: {N.Nat.even(1n+yv) == c : Bool}) -> {N.Nat.even(1n+yv) == False{} : Bool}:
  match c:
    case False{}:
      ec
    case True{}:
      Empty.absurd({N.Nat.even(1n+yv) == False{} : Bool},
        N.Nat.succ_neq_zero(0n,
          N.Nat.succ_inj(1n, 0n,
            G.gc.cop_common(1n+xv, 1n+yv, 2n, hcop,
              PY.py.dvd2_of_even(1n+xv, hxe), PY.py.dvd2_of_even(1n+yv, ec)))))

# name the odd leg first, swapping x and y if necessary
def f4.branch(+g: Nat, +xv: Nat, +yv: Nat, +zv: Nat, bound: N.Nat.Le(zv, 1n+g), +hcop: G.Cop<1n+xv, 1n+yv>, +e: {Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))) == Nat.mul(1n+zv, 1n+zv) : Nat}, rec: f4.Rec(g), c: Bool, ec: {N.Nat.even(1n+xv) == c : Bool}) -> Empty:
  match c:
    case False{}:
      +ec2 = ec
      f4.s1(g, xv, yv, zv, bound, hcop, ec2,
        f4.yeven(xv, yv, zv, ec2, e, N.Nat.even(1n+yv), {==}), e, rec)
    case True{}:
      +ec2 = ec
      f4.s1(g, yv, xv, zv, bound,
        G.gc.cop_sym(1n+xv, 1n+yv, hcop),
        f4.yodd(xv, yv, hcop, ec2, N.Nat.even(1n+yv), {==}),
        ec2,
        Equal.trans(Nat, Nat.add(Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv)), Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv))), Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))), Nat.mul(1n+zv, 1n+zv),
          N.Nat.add_comm(Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv)), Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv))),
          e),
        rec)

def f4.absurd_le21(h: N.Nat.Le(2n, 1n)) -> Empty:
  match h:

# THE DESCENT. The fuel bounds z and every step produces a strictly
# smaller z, so the fuel really runs down.
def f4.go(f: Nat, +xv: Nat, +yv: Nat, zv: Nat, bound: N.Nat.Le(zv, f), hcop: G.Cop<1n+xv, 1n+yv>, e: {Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))) == Nat.mul(1n+zv, 1n+zv) : Nat}) -> Empty:
  match f:
    case 0n:
      match zv:
        case 0n:
          f4.absurd_le21(
            N.Nat.le_trans(2n, Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))), 1n,
              N.Nat.le_add_both(1n, Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), 1n, Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv)), Unit{}, Unit{}),
              N.Nat.eq_le(Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))), 1n, e)))
        case 1n+zv1:
          match bound:
    case 1n+g:
      +g2 = g
      +xv2 = xv
      f4.branch(g2, xv2, yv, zv, bound, hcop, e,
        xv3 => yv3 => zv3 => b3 => c3 => e3 => f4.go(g2, xv3, yv3, zv3, b3, c3, e3),
        N.Nat.even(1n+xv2), {==})

# no coprime solution
def f4.no_coprime(+xv: Nat, +yv: Nat, +zv: Nat, hcop: G.Cop<1n+xv, 1n+yv>, e: {Nat.add(Nat.mul(Nat.mul(1n+xv, 1n+xv), Nat.mul(1n+xv, 1n+xv)), Nat.mul(Nat.mul(1n+yv, 1n+yv), Nat.mul(1n+yv, 1n+yv))) == Nat.mul(1n+zv, 1n+zv) : Nat}) -> Empty:
  f4.go(zv, xv, yv, zv, N.Nat.le_refl(zv), hcop, e)

# ---------------------------------------------------------------------
# from any solution to a coprime one: divide by the gcd
# ---------------------------------------------------------------------

def f4.p4(+t: Nat) -> Nat:
  Nat.mul(Nat.mul(t, t), Nat.mul(t, t))

# (a b)^4 = a^4 b^4
def f4.p4_mul(+a: Nat, +b: Nat) -> {f4.p4(Nat.mul(a, b)) == Nat.mul(f4.p4(a), f4.p4(b)) : Nat}:
  Equal.trans(Nat, Nat.mul(Nat.mul(Nat.mul(a, b), Nat.mul(a, b)), Nat.mul(Nat.mul(a, b), Nat.mul(a, b))), Nat.mul(Nat.mul(Nat.mul(a, a), Nat.mul(b, b)), Nat.mul(Nat.mul(a, a), Nat.mul(b, b))), Nat.mul(Nat.mul(Nat.mul(a, a), Nat.mul(a, a)), Nat.mul(Nat.mul(b, b), Nat.mul(b, b))),
    Equal.cong(Nat, Nat, z => Nat.mul(z, z), Nat.mul(Nat.mul(a, b), Nat.mul(a, b)), Nat.mul(Nat.mul(a, a), Nat.mul(b, b)), G.gc.nsq_prod(a, b)),
    G.gc.nsq_prod(Nat.mul(a, a), Nat.mul(b, b)))

def f4.pos_mul(a: Nat, pa: N.Nat.Lt(0n, a), b: Nat, pb: N.Nat.Lt(0n, b)) -> N.Nat.Lt(0n, Nat.mul(a, b)):
  match a:
    case 0n:
      match pa:
    case 1n+ap:
      match b:
        case 0n:
          match pb:
        case 1n+bp:
          Unit{}

# G^4 (x1^4 + y1^4) = x^4 + y^4, where x = G x1 and y = G y1
def f4.gsum(+xv: Nat, +yv: Nat, +gp: Nat, +x1: Nat, +y1: Nat, ex: {Nat.mul(1n+gp, x1) == 1n+xv : Nat}, ey: {Nat.mul(1n+gp, y1) == 1n+yv : Nat}) -> {Nat.mul(f4.p4(1n+gp), Nat.add(f4.p4(x1), f4.p4(y1))) == Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)) : Nat}:
  Equal.trans(Nat, Nat.mul(f4.p4(1n+gp), Nat.add(f4.p4(x1), f4.p4(y1))), Nat.add(Nat.mul(f4.p4(1n+gp), f4.p4(x1)), Nat.mul(f4.p4(1n+gp), f4.p4(y1))), Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)),
    N.Nat.mul_add(f4.p4(1n+gp), f4.p4(x1), f4.p4(y1)),
    Equal.trans(Nat, Nat.add(Nat.mul(f4.p4(1n+gp), f4.p4(x1)), Nat.mul(f4.p4(1n+gp), f4.p4(y1))), Nat.add(f4.p4(1n+xv), Nat.mul(f4.p4(1n+gp), f4.p4(y1))), Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)),
      Equal.cong(Nat, Nat, z => Nat.add(z, Nat.mul(f4.p4(1n+gp), f4.p4(y1))), Nat.mul(f4.p4(1n+gp), f4.p4(x1)), f4.p4(1n+xv),
        Equal.trans(Nat, Nat.mul(f4.p4(1n+gp), f4.p4(x1)), f4.p4(Nat.mul(1n+gp, x1)), f4.p4(1n+xv),
          Equal.sym(Nat, f4.p4(Nat.mul(1n+gp, x1)), Nat.mul(f4.p4(1n+gp), f4.p4(x1)), f4.p4_mul(1n+gp, x1)),
          Equal.cong(Nat, Nat, z => f4.p4(z), Nat.mul(1n+gp, x1), 1n+xv, ex))),
      Equal.cong(Nat, Nat, z => Nat.add(f4.p4(1n+xv), z), Nat.mul(f4.p4(1n+gp), f4.p4(y1)), f4.p4(1n+yv),
        Equal.trans(Nat, Nat.mul(f4.p4(1n+gp), f4.p4(y1)), f4.p4(Nat.mul(1n+gp, y1)), f4.p4(1n+yv),
          Equal.sym(Nat, f4.p4(Nat.mul(1n+gp, y1)), Nat.mul(f4.p4(1n+gp), f4.p4(y1)), f4.p4_mul(1n+gp, y1)),
          Equal.cong(Nat, Nat, z => f4.p4(z), Nat.mul(1n+gp, y1), 1n+yv, ey)))))
# G^4 divides z^2, so (G^2)^2 | z^2 and therefore G^2 | z
def f4.gdvd(+xv: Nat, +yv: Nat, +zv: Nat, +gp: Nat, +x1: Nat, +y1: Nat, ex: {Nat.mul(1n+gp, x1) == 1n+xv : Nat}, ey: {Nat.mul(1n+gp, y1) == 1n+yv : Nat}, e: {Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)) == Nat.mul(1n+zv, 1n+zv) : Nat}) -> G.Dvd<Nat.mul(1n+gp, 1n+gp), 1n+zv>:
  G.gc.sq_dvd_pos(Nat.mul(1n+gp, 1n+gp), Unit{}, 1n+zv, Unit{},
    G.Dvd{Nat.add(f4.p4(x1), f4.p4(y1)),
      Equal.trans(Nat, Nat.mul(f4.p4(1n+gp), Nat.add(f4.p4(x1), f4.p4(y1))), Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)), Nat.mul(1n+zv, 1n+zv),
        f4.gsum(xv, yv, gp, x1, y1, ex, ey), e)})

def f4.no_coprime_pos(x: Nat, px: N.Nat.Lt(0n, x), y: Nat, py: N.Nat.Lt(0n, y), z: Nat, pz: N.Nat.Lt(0n, z), hcop: G.Cop<x, y>, e: {Nat.add(f4.p4(x), f4.p4(y)) == Nat.mul(z, z) : Nat}) -> Empty:
  match x:
    case 0n:
      match px:
    case 1n+xv:
      match y:
        case 0n:
          match py:
        case 1n+yv:
          match z:
            case 0n:
              match pz:
            case 1n+zv:
              f4.no_coprime(xv, yv, zv, hcop, e)

def f4.main_sq.fin(+xv: Nat, +yv: Nat, +zv: Nat, +gp: Nat, +x1: Nat, +y1: Nat, +z1: Nat, +ex: {Nat.mul(1n+gp, x1) == 1n+xv : Nat}, +ey: {Nat.mul(1n+gp, y1) == 1n+yv : Nat}, +ez: {Nat.mul(Nat.mul(1n+gp, 1n+gp), z1) == 1n+zv : Nat}, hcop: G.Cop<x1, y1>, +e: {Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)) == Nat.mul(1n+zv, 1n+zv) : Nat}) -> Empty:
  f4.no_coprime_pos(x1, G.gc.pos_of_mul_r(1n+gp, x1, xv, ex),
    y1, G.gc.pos_of_mul_r(1n+gp, y1, yv, ey),
    z1, G.gc.pos_of_mul_r(Nat.mul(1n+gp, 1n+gp), z1, zv, ez),
    hcop,
    G.gc.ncancel_pos(Nat.add(f4.p4(x1), f4.p4(y1)), Nat.mul(z1, z1), f4.p4(1n+gp), Unit{},
      Equal.trans(Nat, Nat.mul(f4.p4(1n+gp), Nat.add(f4.p4(x1), f4.p4(y1))), Nat.mul(1n+zv, 1n+zv), Nat.mul(f4.p4(1n+gp), Nat.mul(z1, z1)),
        Equal.trans(Nat, Nat.mul(f4.p4(1n+gp), Nat.add(f4.p4(x1), f4.p4(y1))), Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)), Nat.mul(1n+zv, 1n+zv),
          f4.gsum(xv, yv, gp, x1, y1, ex, ey), e),
        Equal.trans(Nat, Nat.mul(1n+zv, 1n+zv), Nat.mul(Nat.mul(Nat.mul(1n+gp, 1n+gp), z1), Nat.mul(Nat.mul(1n+gp, 1n+gp), z1)), Nat.mul(f4.p4(1n+gp), Nat.mul(z1, z1)),
          Equal.cong(Nat, Nat, w => Nat.mul(w, w), 1n+zv, Nat.mul(Nat.mul(1n+gp, 1n+gp), z1),
            Equal.sym(Nat, Nat.mul(Nat.mul(1n+gp, 1n+gp), z1), 1n+zv, ez)),
          G.gc.nsq_prod(Nat.mul(1n+gp, 1n+gp), z1)))))

def f4.main_sq.dz(+xv: Nat, +yv: Nat, +zv: Nat, +gp: Nat, +x1: Nat, +y1: Nat, +ex: {Nat.mul(1n+gp, x1) == 1n+xv : Nat}, +ey: {Nat.mul(1n+gp, y1) == 1n+yv : Nat}, hcop: G.Cop<x1, y1>, +e: {Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)) == Nat.mul(1n+zv, 1n+zv) : Nat}, dz: G.Dvd<Nat.mul(1n+gp, 1n+gp), 1n+zv>) -> Empty:
  match dz:
    case G.Dvd{z1, ez}:
      f4.main_sq.fin(xv, yv, zv, gp, x1, y1, z1, ex, ey, ez, hcop, e)

def f4.main_sq.op(+xv: Nat, +yv: Nat, +zv: Nat, +e: {Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)) == Nat.mul(1n+zv, 1n+zv) : Nat}, sp: G.Split<1n+xv, 1n+yv>) -> Empty:
  match sp:
    case G.Split{gp, x1, y1, ex, ey, hcop}:
      +ex2 = ex
      +ey2 = ey
      +x1c = x1
      +y1c = y1
      +gpc = gp
      f4.main_sq.dz(xv, yv, zv, gpc, x1c, y1c, ex2, ey2, hcop, e,
        f4.gdvd(xv, yv, zv, gpc, x1c, y1c, ex2, ey2, e))

# x^4 + y^4 is never a square, for positive x and y
def f4.main_sq(+xv: Nat, +yv: Nat, +zv: Nat, +e: {Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)) == Nat.mul(1n+zv, 1n+zv) : Nat}) -> Empty:
  f4.main_sq.op(xv, yv, zv, e, G.gc.split(xv, yv))

def f4.main_sq_pos(x: Nat, px: N.Nat.Lt(0n, x), y: Nat, py: N.Nat.Lt(0n, y), z: Nat, pz: N.Nat.Lt(0n, z), e: {Nat.add(f4.p4(x), f4.p4(y)) == Nat.mul(z, z) : Nat}) -> Empty:
  match x:
    case 0n:
      match px:
    case 1n+xv:
      match y:
        case 0n:
          match py:
        case 1n+yv:
          match z:
            case 0n:
              match pz:
            case 1n+zv:
              f4.main_sq(xv, yv, zv, e)

# and in particular x^4 + y^4 = z^4 is impossible: z^4 is the square of z^2
def f4.main(+xv: Nat, +yv: Nat, +zv: Nat, e: {Nat.add(f4.p4(1n+xv), f4.p4(1n+yv)) == f4.p4(1n+zv) : Nat}) -> Empty:
  f4.main_sq_pos(1n+xv, Unit{}, 1n+yv, Unit{}, Nat.mul(1n+zv, 1n+zv), Unit{}, e)

# ---------------------------------------------------------------------
# THE THEOREM, in mathlib's vocabulary
# ---------------------------------------------------------------------

# the equation carried from mathlib's operators to Base's
def f4.carry(+x2: Nat, +y2: Nat) -> {Nat.add(f4.p4(1n+x2), f4.p4(1n+y2)) == M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)) : Nat}:
  Equal.sym(Nat, M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), Nat.add(f4.p4(1n+x2), f4.p4(1n+y2)),
    Equal.trans(Nat, M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), Nat.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), Nat.add(f4.p4(1n+x2), f4.p4(1n+y2)),
      f4.add_bridge(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)),
      Equal.trans(Nat, Nat.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), Nat.add(f4.p4(1n+x2), M.pow(1n+y2, 4n)), Nat.add(f4.p4(1n+x2), f4.p4(1n+y2)),
        Equal.cong(Nat, Nat, w => Nat.add(w, M.pow(1n+y2, 4n)), M.pow(1n+x2, 4n), f4.p4(1n+x2), f4.pow4_bridge(1n+x2)),
        Equal.cong(Nat, Nat, w => Nat.add(f4.p4(1n+x2), w), M.pow(1n+y2, 4n), f4.p4(1n+y2), f4.pow4_bridge(1n+y2)))))

# x^4 + y^4 = z^4 has no solution in positive naturals (Fermat's own case,
# by infinite descent). The 1n+ makes each variable positive.
def fermat_n4(x: Nat, y: Nat, z: Nat) -> {M.nat_eq(M.add(M.pow(1n+x, 4n), M.pow(1n+y, 4n)), M.pow(1n+z, 4n)) == False{} : Bool}:
  +x2 = x
  +y2 = y
  +z2 = z
  f4.nat_eq_false(M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), M.pow(1n+z2, 4n),
    eq => f4.main(x2, y2, z2,
      Equal.trans(Nat, Nat.add(f4.p4(1n+x2), f4.p4(1n+y2)), M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), f4.p4(1n+z2),
        f4.carry(x2, y2),
        Equal.trans(Nat, M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), M.pow(1n+z2, 4n), f4.p4(1n+z2),
          eq, f4.pow4_bridge(1n+z2)))))

# the stronger statement the descent actually proves: x^4 + y^4 is never a
# perfect square, for positive x and y
def fermat_n4_square(x: Nat, y: Nat, z: Nat) -> {M.nat_eq(M.add(M.pow(1n+x, 4n), M.pow(1n+y, 4n)), M.mul(1n+z, 1n+z)) == False{} : Bool}:
  +x2 = x
  +y2 = y
  +z2 = z
  f4.nat_eq_false(M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), M.mul(1n+z2, 1n+z2),
    eq => f4.main_sq(x2, y2, z2,
      Equal.trans(Nat, Nat.add(f4.p4(1n+x2), f4.p4(1n+y2)), M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), Nat.mul(1n+z2, 1n+z2),
        f4.carry(x2, y2),
        Equal.trans(Nat, M.add(M.pow(1n+x2, 4n), M.pow(1n+y2, 4n)), M.mul(1n+z2, 1n+z2), Nat.mul(1n+z2, 1n+z2),
          eq, f4.mul_bridge(1n+z2, 1n+z2)))))