-- Executable Euclidean GCD with a proved sufficient iteration bound
-- Part of Cumulus, the Grace17 mathematics library.
-- Copyright (c) 2026 Grant Passmore.
-- Released under the Apache License 2.0. See LICENSE in the repository root.
-- Authors: Grant Passmore

module Grace.Cumulus.NumberTheory.Euclidean
import GCD
import Nat.Division

-- The second input strictly decreases on each nonterminal Euclidean step.
-- The public function uses b steps; correctness below proves this sufficient.
def gcdSteps (fuel a b : Nat) → Nat
  | zero, a, b => a
  | succ(fuel), a, b => chooseNat(isZero(b), a, gcdSteps(fuel, b, a % b))

def gcd (a b : Nat) : Nat := gcdSteps(b, a, b)

theorem gcdSteps.zero_right (a fuel : Nat)
  shows gcdSteps(fuel, a, 0) = a
proof
  induction fuel
  | zero => settle with gcdSteps
  | succ(k) using ih => settle with gcdSteps, isZero, chooseNat

theorem gcdSteps.successor (fuel a b : Nat)
  shows gcdSteps(succ(fuel), a, succ(b)) = gcdSteps(fuel, succ(b), a % succ(b))
proof
  simp [gcdSteps, isZero, chooseNat]

theorem gcdSteps.natural (fuel a b : Nat)
  shows Nat(gcdSteps(fuel, a, b))
proof
  induction fuel generalizing a b
  | zero =>
    intros a b aNat bNat
    simp [gcdSteps, isZero, chooseNat]
  | succ(k) using ih =>
    intros a b aNat bNat
    have remainderNat := mod.natural at b, a with *
    have recursive := ih with b, a % b, bNat, remainderNat
    have boolean := isZero.boolean at b with *
    have natural := chooseNat.natural at isZero(b), a, gcdSteps(k, b, a % b) with *
    simpa [gcdSteps] using natural

theorem gcdSteps.correct (fuel a b : Nat)
  shows b ≤ fuel → IsGCD(a, b, gcdSteps(fuel, a, b))
proof
  induction fuel generalizing a b
  | zero =>
    intros a b aNat bNat bound
    have bZero : b = 0
    . have lower := Grace.Cumulus.Nat.Order.LessEqual.zero at b with *
      exact LessEqual.antisymmetric at b, 0 with *
    rw bZero
    have terminal := IsGCD.zero_right at a with *
    simpa [gcdSteps, isZero, chooseNat] using terminal
  | succ(k) using ih =>
    intros a b aNat bNat bound
    have shape := natural.zero_or_successor at b with *
    cases shape with zero | successor
    . rw zero
      have terminal := IsGCD.zero_right at a with *
      have value := gcdSteps.zero_right at a, succ(k) with *
      rwa value
    . obtain ⟨j, jNat, bShape⟩ := successor
      rw bShape at bound
      rw bShape
      have nonzero : succ(j) ≠ 0
      . intro equal
        have distinct := Naturals.disjoint.zero.succ at j
        symm at equal
        exact distinct with equal
      have remainderNat := mod.natural at succ(j), a with *
      have remainderBound := mod.less at a, succ(j) with *
      have enough : a % succ(j) ≤ k
      . arith
      have recursive := ih with succ(j), a % succ(j), assumption, remainderNat, enough
      obtain ⟨q, qNat, equation⟩ := mod.equation at a, succ(j) with *
      have preserved := IsGCD.euclidean_backward at a, succ(j), q, a % succ(j), gcdSteps(k, succ(j), a % succ(j)) with *
      have value := gcdSteps.successor at k, a, j with *
      rwa value

theorem gcd.natural (a b : Nat)
  shows Nat(gcd(a, b))
proof
  unfold gcd
  exact gcdSteps.natural with *

theorem gcd.correct (a b : Nat)
  shows IsGCD(a, b, gcd(a, b))
proof
  have bound := LessEqual.reflexive at b with *
  have correct := gcdSteps.correct at b, a, b with *, bound
  simpa [gcd] using correct

theorem gcd.unique (a b g : Nat)
  requires certificate : IsGCD(a, b, g)
  shows gcd(a, b) = g
proof
  have correct := gcd.correct at a, b with *
  exact IsGCD.unique with *

theorem IsGCD.exists (a b : Nat)
  shows ∃ g, Nat(g) ∧ IsGCD(a, b, g)
proof
  have natural := gcd.natural at a, b with *
  have correct := gcd.correct at a, b with *
  refine ⟨gcd(a, b), natural, correct⟩

theorem gcd.zero_right (a : Nat)
  shows gcd(a, 0) = a
proof
  have certificate := IsGCD.zero_right with *
  exact gcd.unique with *

theorem gcd.zero_left (a : Nat)
  shows gcd(0, a) = a
proof
  have certificate := IsGCD.zero_left with *
  exact gcd.unique with *

theorem gcd.commutative (a b : Nat)
  shows gcd(a, b) = gcd(b, a)
proof
  have reverse := gcd.correct at b, a with *
  have certificate := IsGCD.symmetric with *
  exact gcd.unique with *

theorem gcd.divides_left (a b : Nat)
  shows gcd(a, b) ∣ a
proof
  have correct := gcd.correct at a, b with *
  exact IsGCD.divides_left with *

theorem gcd.divides_right (a b : Nat)
  shows gcd(a, b) ∣ b
proof
  have correct := gcd.correct at a, b with *
  exact IsGCD.divides_right with *

theorem gcd.greatest (a b d : Nat)
  requires left : d ∣ a
  requires right : d ∣ b
  shows d ∣ gcd(a, b)
proof
  have correct := gcd.correct at a, b with *
  exact IsGCD.greatest with *

theorem gcd.euclidean (a b : Nat)
  requires nonzero : b ≠ 0
  shows gcd(a, b) = gcd(b, a % b)
proof
  have recursive := gcd.correct at b, a % b with *
  obtain ⟨q, qNat, equation⟩ := mod.equation at a, b with *
  have certificate := IsGCD.euclidean_backward with *
  exact gcd.unique with *

theorem gcd.of_divides (a b : Nat)
  requires divides : a ∣ b
  shows gcd(a, b) = a
proof
  have certificate := IsGCD.of_divides with *
  exact gcd.unique with *

theorem gcd.zero_iff (a b : Nat)
  shows gcd(a, b) = 0 ↔ (a = 0 ∧ b = 0)
proof
  constructor
  . intro zero
    have left := gcd.divides_left at a, b with *
    have right := gcd.divides_right at a, b with *
    rw zero at left
    rw zero at right
    have aZero := zero.divides_only_zero at a with *
    have bZero := zero.divides_only_zero at b with *
    tauto
  . intro zeros
    rcases zeros with ⟨aZero, bZero⟩
    rw [aZero, bZero]
    exact gcd.zero_right with *