-rw-r--r-- 16038 leangoppa-20230726/Goppadecoding/approximant.lean raw
import Mathlib.RingTheory.Polynomial.Basic
import Mathlib.LinearAlgebra.FiniteDimensional
import Goppadecoding.nat
import Goppadecoding.dvd
import Goppadecoding.sum
import Goppadecoding.lineardependence
import Goppadecoding.zcoeff

open Polynomial BigOperators Finset

namespace Goppadecoding

section approximants

variable {k: Type _} [Field k]

def approximant (A B: k[X]) (t: ℕ) (a b: k[X]): Prop :=
  IsCoprime a b ∧ a.degree ≤ t ∧ t + (a*B-b*A).degree < A.degree

def small_approximant (A B: k[X]) (t: ℕ) (a b: k[X]): Prop :=
  approximant A B t a b ∧ b.degree < t

-- this proof follows the "explicitly" line
-- alternative is to avoid coeffs for k-module {(a,b): deg a ≤ t, deg b < t}:
--   "simply" set up k-linear map to the interesting coeffs of (a*B-b*A)
theorem approximant_exists_lemma {A B: k[X]} (degAB: A.degree > B.degree) (t: ℕ)
: ∃ (a b:k[X]), (a ≠ 0 ∨ b ≠ 0) ∧ a.degree ≤ t ∧ b.degree < t ∧ t + (a*B-b*A).degree < A.degree :=
by
  have degAnat: degree A ≠ ⊥ := LT.lt.ne_bot degAB
  let n := (degree A).unbot degAnat
  have degAn: degree A = n := Eq.symm (WithBot.coe_unbot (degree A) degAnat)
  let V := Fin (2*t) → k
  let M: (Fin (2*t+1) → V) := fun j ↦ fun i ↦ if (j:ℕ) < t+1 then (zcoeff B (n-t+(i:ℤ)-(j:ℤ))) else (-zcoeff A (n+1+(i:ℤ)-(j:ℤ)))
  have explicitdependency: ∃ c : Fin (2*t+1) → k, ∑ j, c j • M j = 0 ∧ ∃ j, c j ≠ 0 := by
    have rankVbelowMnum: FiniteDimensional.finrank k V < Fintype.card (Fin (2*t+1)) := by
      calc
        FiniteDimensional.finrank k V = 2*t := FiniteDimensional.finrank_fin_fun k 
        _ < 2*t+1 := by simp only [lt_add_iff_pos_right]
        _ = Fintype.card (Fin (2*t+1)) := by rw[eq_comm,Fintype.card_fin (2*t+1)]
    have dependency: ¬LinearIndependent k M := by
      apply dependency_if_finrank_space_lt_card_family
      apply rankVbelowMnum
    exact Fintype.not_linearIndependent_iff.mp dependency
  cases explicitdependency with | intro c cdep =>
    have cdepleft: (∑ j: Fin (2*t+1), c j • M j = 0) := by simp only [cdep.left]
    have cdepleft2: (∀ i: Fin (2*t), ((∑ j: Fin (2*t+1), c j * M j i) = 0)) := by
      intro i
      calc
        ∑ j: Fin (2*t+1), c j * M j i
          = ∑ j: Fin (2*t+1), (c j • M j) i := rfl
        _ = (∑ j: Fin (2*t+1), (c j • M j)) i := by simp only [Pi.smul_apply, smul_eq_mul, mul_ite, mul_neg, Finset.sum_apply]
        _ = 0 := by rw[cdepleft]; simp only [Pi.zero_apply]
    let a := ∑ j in range (t+1), (monomial: ℕ → k →ₗ[k] k[X]) j (c j) 
    let b := ∑ j in range t, (monomial: ℕ → k →ₗ[k] k[X]) j (c (j+t+1))
    have nonzero: a ≠ 0 ∨ b ≠ 0 := by
      cases cdep.right with | intro x cxnonzero =>
        let xn: ℕ := x
        by_cases xnt: xn < t+1 
        . left 
          have axn: coeff a xn ≠ 0 := by 
            calc
              coeff a xn
                = coeff (∑ j in range (t+1), (monomial: ℕ → k →ₗ[k] k[X]) j (c j)) xn := by rfl
              _ = c xn := by rw[coeff_sum_range_monomial]; simp only [xnt, ite_true]
              _ = c x := by simp only; rw[Fin.cast_val_eq_self]
              _ ≠ 0 := by apply cxnonzero
          apply nonzero_if_coeff_nonzero axn
        right
        let i: ℕ := xn-(t+1)
        have xn_ge_t1: xn ≥ t+1 := Nat.not_lt.mp xnt
        have ixn: xn = i+(t+1) := (Nat.sub_eq_iff_eq_add xn_ge_t1).mp rfl
        have xn_lt_2t1: xn < 2*t+1 := Fin.prop x
        have i_lt_t: i < t := by linarith
        have bxn: coeff b i ≠ 0 := by
          calc
            coeff b i
              = coeff (∑ j in range t, (monomial: ℕ → k →ₗ[k] k[X]) j (c (j+t+1))) i := by rfl
            _ = c (i+t+1) := by rw[coeff_sum_range_monomial]; simp only [ge_iff_le, i_lt_t, ite_true]
            _ = c (i+(t+1)) := by group
            _ = c xn := by rw[ixn]; simp only [ge_iff_le, Nat.cast_add, Nat.cast_one]
            _ = c x := by simp only; rw[Fin.cast_val_eq_self]
            _ ≠ 0 := by apply cxnonzero
        apply nonzero_if_coeff_nonzero bxn
    use a 
    use b
    constructor
    . apply nonzero
    have degale: degree a ≤ t := by
      have degalt: degree a < t+1 := by apply degree_sum_range_lt 
      apply Order.lt_succ_iff.mp
      apply degalt
    constructor
    . apply degale
    have degblt: degree b < t := by apply degree_sum_range_lt 
    constructor
    . apply degblt
    have aBbAcoeff: ∀ d, coeff (a*B-b*A) d = (∑ j in range (t+1), (c j * zcoeff B (d-j))) + (∑ j in range t, (c (j+t+1) * -zcoeff A (d-j))) := by
      intro d
      calc
        coeff (a*B-b*A) d
          = coeff (a*B) d - coeff (b*A) d := by rw[coeff_sub]
        _ = (∑ j in range (t+1), (c j * zcoeff B (d-j))) - (∑ j in range t, (c (j+t+1) * zcoeff A (d-j))) := by repeat rw[sum_monomials_mul_coeff]
        _ = (∑ j in range (t+1), (c j * zcoeff B (d-j))) + -(∑ j in range t, (c (j+t+1) * zcoeff A (d-j))) := by apply sub_eq_add_neg
        _ = (∑ j in range (t+1), (c j * zcoeff B (d-j))) + (∑ j in range t, c (j+t+1) * -zcoeff A (d-j)) := by simp only [mul_neg, sum_neg_distrib]
    have abbacoeff0: ∀ d:ℕ, (n-t:ℤ) ≤ (d:ℤ) → coeff (a*B-b*A) d = 0 := by
      intro d dlower
      by_cases dupper: t+n ≤ d
      . have aBbAdeg: degree (a*B-b*A) < t+n := by
          have tdegB: t + degree B < t + n := by
            by_cases bot: degree B = ⊥
            . simp only [bot, WithBot.add_bot]
              apply WithBot.bot_lt_coe
            let m := (degree B).unbot bot
            have degBm: degree B = m := Eq.symm (WithBot.coe_unbot (degree B) bot)
            rw[degBm]
            apply WithBot.add_lt_add_left
            exact WithBot.nat_ne_bot t
            rw[← degBm,← degAn]
            exact degAB
          have aBdeg: degree (a*B) < t+n := by calc
            degree (a*B) = degree a + degree B := degree_mul
            _ ≤ t + degree B := add_le_add_right degale (degree B)
            _ < t + n := tdegB
          have bAdeg: degree (b*A) < t+n := by calc
            degree (b*A) = degree b + degree A := degree_mul
            _ < t + degree A := WithBot.add_lt_add_right degAnat degblt
            _ ≤ t+n := Eq.le (congrArg (HAdd.hAdd ↑t) degAn)
          calc
            degree (a*B-b*A) ≤ max (degree (a*B)) (degree (b*A)) := degree_sub_le (a*B) (b*A)
            _ < t+n := max_lt aBdeg bAdeg
        have degd: degree (a*B-b*A) < d := by calc
          degree (a*B-b*A) < t+n := aBbAdeg
          _ ≤ d := by norm_cast
        exact coeff_eq_zero_of_degree_lt degd
      let i:ℤ := (d:ℤ)-(n-t:ℤ)
      have di: d = n-t+i := by group
      let inat:ℕ := i.toNat
      have inatupper: inat < 2*t := by
        have iupper: i < 2*t := by
          have dntupper: (d:ℤ)-(n-t:ℤ) ≤ 2*t-1 := by linarith
          exact Int.le_sub_one_iff.mp dntupper
        have dntlower: 0 ≤ (d:ℤ)-(n-t:ℤ) := by linarith
        exact (Int.toNat_lt dntlower).mpr iupper
      let ifin: Fin (2*t) := ⟨inat,inatupper⟩
      have iifin: i = (ifin:ℤ) := by
        have iinat: inat = i := Int.toNat_sub_of_le dlower
        exact id (Eq.symm iinat)
      calc
        coeff (a*B-b*A) d
          = (∑ j in range (t+1), (c j * zcoeff B (d-j))) + (∑ j in range t, c (j+t+1) * -zcoeff A (d-j))
            := aBbAcoeff d
        _ = (∑ j in range (t+1), (c j * zcoeff B (n-t+i-j))) + (∑ j in range t, c (j+t+1) * -zcoeff A (n-t+i-j))
            := by rw[di]
        _ = (∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c (j:ℕ) * zcoeff B (n-t+i-(j:ℕ))) else 0) + (∑ j in range t, c (j+t+1) * -zcoeff A (n-t+i-j))
            := by
                have t12t1: t+1 ≤ 2*t+1 := by linarith
                rw[sum_range_eq_sum_fin_ite t12t1]
        _ = (∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c (j:ℕ) * zcoeff B (n-t+i-(j:ℕ))) else 0) + (∑ j: Fin (2*t+1), if (j:ℕ) < t then c ((j:ℕ)+t+1) * -zcoeff A (n-t+i-(j:ℕ)) else 0)
            := by
                have t2t1: t ≤ 2*t+1 := by linarith
                rw[sum_range_eq_sum_fin_ite t2t1]
        _ = (∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c (j:ℕ) * zcoeff B (n-t+i-(j:ℕ))) else 0) + (∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then 0 else c (((j:ℕ)-(t+1):ℕ)+t+1) * -zcoeff A (n-t+i-((j:ℕ)-(t+1):ℕ)))
            := by
                have tt1: 2*t+1 = t+(t+1) := by group
                rw[sum_fin_ite_eq_sum_fin_ite_sub tt1 (fun j ↦ c ((j:ℕ)+t+1) * -zcoeff A (n-t+i-(j:ℕ)))]
        _ = ∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c (j:ℕ) * zcoeff B (n-t+i-(j:ℕ))) else c (((j:ℕ)-(t+1):ℕ)+t+1) * -zcoeff A (n-t+i-((j:ℕ)-(t+1):ℕ))
            := by apply sum_ite_add_sum_ite
        _ = ∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c (j:ℕ) * zcoeff B (n-t+i-(j:ℕ))) else (c j * -zcoeff A (n-t+i-((j:ℕ)-(t+1):ℕ)))
            := by
                apply Finset.sum_congr
                simp only
                intro x _
                simp only [add_sub_cancel'_right, ge_iff_le, mul_neg]
                by_cases xt1: (x:ℕ) < t+1
                . simp only [xt1, ge_iff_le, ite_true]
                . simp only [xt1, ge_iff_le, ite_false, neg_inj, mul_eq_mul_right_iff]
                  left; apply congr_arg
                  rw[fin_nat_sub_eq_sub (Nat.not_lt.mp xt1)]
                  rw[Nat.cast_add]
                  group
        _ = ∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c j * zcoeff B (n-t+i-(j:ℕ))) else (c j * -zcoeff A (n-t+i-((j:ℕ)-(t+1):ℕ)))
            := by
                apply Finset.sum_congr
                simp only
                intro x _
                simp only [add_sub_cancel'_right, ge_iff_le, mul_neg]
                by_cases xt1: (x:ℕ) < t+1
                . simp only [xt1, ge_iff_le, ite_true, mul_eq_mul_right_iff]
                  left; apply congr_arg; exact Fin.cast_val_eq_self x
                . simp only [xt1, ge_iff_le, ite_false]
        _ = ∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c j * zcoeff B (n-t+i-j)) else (c j * -zcoeff A (n-t+i-((j:ℕ)-(t+1):ℕ)))
            := by simp only [add_sub_cancel'_right, ge_iff_le, mul_neg]
        _ = ∑ j: Fin (2*t+1), if (j:ℕ) < t+1 then (c j * zcoeff B (n-t+i-j)) else (c j * -zcoeff A (n+1+i-j))
            := by
                apply Finset.sum_congr
                simp only
                intro x _
                simp only [add_sub_cancel'_right, ge_iff_le, mul_neg]
                by_cases xt1: (x:ℕ) < t+1
                . simp only [xt1, ge_iff_le, ite_true]
                . simp only [xt1, ge_iff_le, ite_false, neg_inj, mul_eq_mul_left_iff]
                  left; apply congr_arg
                  have xsum: ((x:ℕ):ℤ) = ((((x:ℕ)-(t+1):ℕ) + (t+1):ℕ):ℤ) := by
                    rw[Nat.sub_add_cancel (Nat.not_lt.mp xt1)]
                  rw[xsum]
                  group
        _ = ∑ j: Fin (2*t+1), c j * (if (j:ℕ) < t+1 then zcoeff B (n-t+i-j) else -zcoeff A (n+1+i-j))
            := by simp only [add_sub_cancel'_right, mul_neg, mul_ite]
        _ = ∑ j: Fin (2*t+1), c j * (if (j:ℕ) < t+1 then zcoeff B (n-t+ifin-j) else -zcoeff A (n+1+ifin-j))
            := by rw[iifin]
        _ = ∑ j: Fin (2*t+1), c j * M j ifin
            := rfl
        _ = 0 := cdepleft2 ifin
    rw[degAn]
    exact degree_add_lt_if_coeff_zero (a*B-b*A) t n abbacoeff0
    
theorem small_approximant_exists {A B: k[X]} (degAB: A.degree > B.degree) {t: ℕ}
: ∃ a b, small_approximant A B t a b :=
by classical
  unfold small_approximant
  unfold approximant
  have partway: ∃ a b, (a ≠ 0 ∨ b ≠ 0) ∧ degree a ≤ t ∧ degree b < t ∧ t + degree (a*B-b*A) < degree A := approximant_exists_lemma degAB t
  cases partway with | intro a acond =>
    cases acond with | intro b abcond =>
      let g := EuclideanDomain.gcd a b
      have existsag: ∃ ag, a = g*ag := EuclideanDomain.gcd_dvd_left a b
      have existsbg: ∃ bg, b = g*bg := EuclideanDomain.gcd_dvd_right a b
      cases existsag with | intro ag agformula =>
        cases existsbg with | intro bg bgformula =>
          use ag
          use bg
          have g0: g ≠ 0 := by
            apply gcd_ne_zero_iff.mpr
            simp only [ne_eq, abcond]
          have degg: 0 ≤ degree g := zero_le_degree_iff.mpr g0
          constructor
          . constructor
            . let x := EuclideanDomain.gcdA a b
              let y := EuclideanDomain.gcdB a b
              have axby: g = a*x+b*y := EuclideanDomain.gcd_eq_gcd_ab a b
              have gagxbgy: g*(ag*x+bg*y) = g := by rw[mul_add,← mul_assoc,← mul_assoc,← agformula,← bgformula,axby]
              have agxbgy: ag*x+bg*y = 1 := eq_one_if_mul_eq_self g0 gagxbgy
              unfold IsCoprime
              use x
              use y
              rw[← agxbgy]; ring
            constructor
            . have degagag: degree a = degree g + degree ag := by rw[agformula,degree_mul]
              have degag: degree ag ≤ degree a := by rw[degagag]; apply le_add_of_nonneg_left degg
              apply le_trans degag abcond.2.1
            calc
              t + degree (ag*B-bg*A)
                ≤ degree g + (t + degree (ag*B-bg*A)) := le_add_self_if_le degg
              _ = t + degree (a*B-b*A) := by
                    have h: a*B-b*A = g*(ag*B-bg*A) := by rw[mul_sub,← mul_assoc,← mul_assoc,agformula,bgformula]
                    have h2: degree (a*B-b*A) = degree g + degree (ag*B-bg*A) := by rw[h,degree_mul]
                    rw[h2]
                    group
              _ < degree A := by simp only [abcond]
          have degbgbg: degree b = degree g + degree bg := by rw[bgformula,degree_mul]
          have degbg: degree bg ≤ degree b := by rw[degbgbg]; apply le_add_of_nonneg_left degg
          apply lt_of_le_of_lt degbg abcond.2.2.1

theorem approximant_best {A B a b c d: k[X]} {t: ℕ}
  (appr: approximant A B t a b)
  (cdeg: c.degree ≤ t)
  (cBdAdeg: t + (c*B-d*A).degree < A.degree)
: ∃ (L: k[X]), (c = a*L ∧ d = b*L) :=
by
  unfold approximant at appr
  let coprime := appr.1
  let dega := appr.2.1
  let degaBbA := appr.2.2
  have D: a*d = b*c := by
    have D1: degree (a*d-c*b) + degree A < degree A := by
      have degc2: degree (c*(a*B-b*A)) < degree A := by
        calc
          _ = degree c + degree (a*B-b*A) := degree_mul
          _ ≤ t + degree (a*B-b*A) := add_le_add_right cdeg (degree (a*B-b*A))
          _ < degree A := degaBbA
      have dega2: degree (a*(c*B-d*A)) < degree A := by
        calc
          _ = degree a + degree (c*B-d*A) := degree_mul
          _ ≤ t + degree (c*B-d*A) := add_le_add_right dega (degree (c*B-d*A))
          _ < degree A := cBdAdeg
      calc
        _ = degree ((a*d-c*b)*A) := by rw[← degree_mul]
        _ = degree (c*(a*B-b*A)-a*(c*B-d*A)) := by ring_nf
        _ ≤ max (degree (c*(a*B-b*A))) (degree (a*(c*B-d*A))) := degree_sub_le (c*(a*B-b*A)) (a*(c*B-d*A))
        _ < degree A := max_lt degc2 dega2
    have D2: degree (a*d-c*b) < 0 := neg_of_add_lt_left D1
    have D3: degree (a*d-c*b) = ⊥ := (Nat.WithBot.lt_zero_iff (degree (a*d-c*b))).mp D2
    have D4: a*d-c*b = 0 := degree_eq_bot.mp D3
    calc
      _ = a*d - 0 := by group
      _ = a*d - (a*d-c*b) := by rw[← D4]
      _ = b*c := by ring
  have cdiv: a ∣ b*c := by
    apply Dvd.intro d
    rw[D]
  have cdiv2: a ∣ c := by apply dvd_if_dvd_mul_left coprime cdiv
  have ddiv: b ∣ a*d := by
    apply Dvd.intro c
    rw[← D]
  have ddiv2: b ∣ d := by apply dvd_if_dvd_mul_left (IsCoprime.symm coprime) ddiv
  cases cdiv2 with | intro L caL =>
    by_cases a0: a = 0
    . cases ddiv2 with | intro M dbM =>
        use M
        constructor
        . calc
            c = a*L := caL
            _ = 0*L := by rw[a0]
            _ = 0*M := by ring
            _ = a*M := by rw[← a0]
        . apply dbM
    . use L
      constructor
      . apply caL
      . have h: a*d = a*(b*L) := by calc
          a*d = b*c := D
          _ = b*(a*L) := by rw[caL]
          _ = a*(b*L) := by ring
        apply mul_left_cancel₀ a0 h

end approximants