-- Generated by the TypeDensity exporter from TypeDensity/Reciprocity/Defs.lean (namespace TypeDensity); do not edit.
import Definitions.Def_TypeDensity_Basic
import Mathlib.FieldTheory.RatFunc.Basic
import Mathlib.MeasureTheory.Measure.Haar.Unique
import Mathlib.NumberTheory.ArithmeticFunction.Moebius
import Mathlib.NumberTheory.Padics.PadicIntegers
import Mathlib.RingTheory.IntegralClosure.IntegralRestrict
import Mathlib.RingTheory.MvPowerSeries.Basic
import Mathlib.RingTheory.RamificationInertia.Basic
import Mathlib.Topology.Algebra.Nonarchimedean.AdicTopology
import Mathlib.FieldTheory.RatFunc.AsPolynomial
import Mathlib.RingTheory.PowerSeries.Expand
import Mathlib.RingTheory.PowerSeries.Log
set_option autoImplicit false

noncomputable section

namespace TypeDensity

namespace Expr

/-- Rational well-formedness: every inverted subexpression has a nonzero numerator polynomial, so
that `Expr.rational` is a ring homomorphism on the expression (its unreduced denominators never
vanish). -/
def wf : Expr → Prop
  | .const _ => True
  | .var => True
  | .add x y => wf x ∧ wf y
  | .mul x y => wf x ∧ wf y
  | .inv x => wf x ∧ x.fraction.1 ≠ 0

end Expr

/-- Polynomials evaluated at X⁻¹ in ℚ(X). -/
def recipPoly : Polynomial ℚ →+* RatFunc ℚ :=
  Polynomial.eval₂RingHom (RatFunc.C : ℚ →+* RatFunc ℚ) (RatFunc.X)⁻¹

/-- Polynomials evaluated at X^δ in ℚ(X). -/
def adamsPoly (δ : ℕ) : Polynomial ℚ →+* RatFunc ℚ :=
  Polynomial.eval₂RingHom (RatFunc.C : ℚ →+* RatFunc ℚ) (RatFunc.X ^ δ)

/-- Evaluating a polynomial at `X` in `ℚ(X)` is the canonical embedding. -/
theorem eval₂_C_X_eq_algebraMap (g : Polynomial ℚ) :
    Polynomial.eval₂ (RatFunc.C : ℚ →+* RatFunc ℚ) (RatFunc.X : RatFunc ℚ) g =
      algebraMap (Polynomial ℚ) (RatFunc ℚ) g := by
  induction g using Polynomial.induction_on' with
  | add p q hp hq => simp [hp, hq]
  | monomial n a => simp [RatFunc.algebraMap_monomial]

/-- A nonzero polynomial stays nonzero after the substitution `X ↦ X⁻¹`. -/
theorem recipPoly_ne_zero (p : Polynomial ℚ) (hp : p ≠ 0) : recipPoly p ≠ 0 := by
  let _ : Invertible (RatFunc.X : RatFunc ℚ) := invertibleOfNonzero RatFunc.X_ne_zero
  have key := Polynomial.eval₂_reflect_eq_zero_iff (RatFunc.C : ℚ →+* RatFunc ℚ)
    (RatFunc.X : RatFunc ℚ) p.natDegree (p.reflect p.natDegree)
    (by simpa using Polynomial.natDegree_reflect_le (N := p.natDegree) (p := p))
  rw [Polynomial.reflect_reflect, invOf_eq_inv] at key
  intro h
  have h3 := key.mp h
  rw [eval₂_C_X_eq_algebraMap] at h3
  exact hp (Polynomial.reflect_eq_zero_iff.mp
    (not_not.1 fun hne ↦ RatFunc.algebraMap_ne_zero hne h3))

/-- `recipPoly` is injective. -/
theorem recipPoly_injective : Function.Injective recipPoly :=
  (injective_iff_map_eq_zero _).2 fun p hp ↦ not_not.1 fun hne ↦ recipPoly_ne_zero p hne hp

/-- A nonzero polynomial stays nonzero after the substitution `X ↦ X^δ`, `δ > 0`. -/
theorem adamsPoly_ne_zero (δ : ℕ) (hδ : 0 < δ) (p : Polynomial ℚ) (hp : p ≠ 0) :
    adamsPoly δ p ≠ 0 := by
  intro h
  have h2 : Polynomial.eval₂ (RatFunc.C : ℚ →+* RatFunc ℚ) (RatFunc.X : RatFunc ℚ)
      ((Polynomial.X : Polynomial ℚ) ^ δ) = (RatFunc.X : RatFunc ℚ) ^ δ := by simp
  change Polynomial.eval₂ _ ((RatFunc.X : RatFunc ℚ) ^ δ) p = 0 at h
  rw [← h2, ← Polynomial.eval₂_comp, eval₂_C_X_eq_algebraMap] at h
  rcases Polynomial.comp_eq_zero_iff.mp
      (not_not.1 fun hne ↦ RatFunc.algebraMap_ne_zero hne h) with h4 | ⟨_, h5⟩
  · exact hp h4
  · have h6 := congrArg Polynomial.natDegree h5
    rw [Polynomial.natDegree_X_pow, Polynomial.natDegree_C] at h6
    omega

/-- `adamsPoly δ` is injective for `δ > 0`. -/
theorem adamsPoly_injective (δ : ℕ) (hδ : 0 < δ) : Function.Injective (adamsPoly δ) :=
  (injective_iff_map_eq_zero _).2 fun p hp ↦ not_not.1 fun hne ↦ adamsPoly_ne_zero δ hδ p hne hp

/-- The automorphism q ↦ q⁻¹ of ℚ(X) (reciprocal substitution, the proof's `v`). -/
def recipRat : RatFunc ℚ →+* RatFunc ℚ :=
  RatFunc.liftRingHom recipPoly
    (nonZeroDivisors_le_comap_nonZeroDivisors_of_injective recipPoly recipPoly_injective)

/-- The Adams substitution q ↦ q^δ on ℚ(X) (δ > 0; the identity for δ = 0, never used). -/
def adamsRat (δ : ℕ) : RatFunc ℚ →+* RatFunc ℚ :=
  if hδ : 0 < δ then
    RatFunc.liftRingHom (adamsPoly δ)
      (nonZeroDivisors_le_comap_nonZeroDivisors_of_injective (adamsPoly δ)
        (adamsPoly_injective δ hδ))
  else RingHom.id (RatFunc ℚ)

/-- The coefficient algebra ℚ(q)[Z_{e,f}]: finitely supported functions on types, multiplied by
adding the multisets (the product of Z-monomials). A law becomes an element of this commutative
ring. -/
abbrev Coef := AddMonoidAlgebra (RatFunc ℚ) FactorType

/-- The scalar embedding ℚ(q) → Coef (the empty type is the monomial 1). -/
def scalar (f : RatFunc ℚ) : Coef := AddMonoidAlgebra.single 0 f

/-- q^k for an integer k, as a scalar of Coef. -/
def qpow (k : ℤ) : Coef := scalar (RatFunc.X ^ k)

/-- A law as an element of Coef: the sum over its entries of (rational value) · Z^{type}. -/
def lawCoef (L : Law) : Coef := (L.map fun z ↦ AddMonoidAlgebra.single z.1 z.2.rational).sum

/-- scaleType δ as an additive monoid homomorphism of types (multisets add). -/
def scaleTypeHom (δ : ℕ) : FactorType →+ FactorType where
  toFun := scaleType δ
  map_zero' := by simp [scaleType]
  map_add' := fun σ τ ↦ by simp [scaleType, Multiset.map_add]

/-- The Adams operation ψ_δ on Coef: q ↦ q^δ on every coefficient and Z_{e,f} ↦ Z_{e,δf} on
types. -/
def adamsCoef (δ : ℕ) : Coef →+* Coef :=
  (AddMonoidAlgebra.mapDomainRingHom (RatFunc ℚ) (scaleTypeHom δ)).comp
    (AddMonoidAlgebra.mapRingHom FactorType (adamsRat δ))

/-- Reciprocal substitution v on Coef: q ↦ q⁻¹ on every coefficient, types unchanged. -/
def recipCoef : Coef →+* Coef := AddMonoidAlgebra.mapRingHom FactorType recipRat

/-- Every entry of the law is a well-formed expression. -/
def LawWf (L : Law) : Prop := ∀ z ∈ L, z.2.wf

/-- Series in the cluster-degree variable x with coefficients in Coef. -/
abbrev Ser := PowerSeries Coef

/-- Agreement of two series in every degree ≤ n (congruence modulo x^{n+1}). -/
def SCongr (n : ℕ) (f g : Ser) : Prop :=
  ∀ m : ℕ, m ≤ n → PowerSeries.coeff m f = PowerSeries.coeff m g

/-- The Adams operation ψ_k on series: ψ_k on every coefficient and x ↦ x^k (zero for k = 0,
never used). -/
def adamsSer (k : ℕ) (f : Ser) : Ser :=
  if hk : k = 0 then 0 else PowerSeries.expand k hk (PowerSeries.map (adamsCoef k) f)

/-- Σ_{k ≥ 1} ψ_k(f)/k, coefficientwise finite: the coefficient of xⁿ is
Σ_{k ∣ n} (1/k) ψ_k(f_{n/k}). -/
def plethSum (f : Ser) : Ser :=
  PowerSeries.mk fun n ↦ ∑ k ∈ n.divisors, ((k : ℚ)⁻¹) • adamsCoef k (PowerSeries.coeff (n / k) f)

/-- Σ_{k ≥ 1} μ(k) ψ_k(g)/k, coefficientwise finite. -/
def moebiusSum (g : Ser) : Ser :=
  PowerSeries.mk fun n ↦ ∑ k ∈ n.divisors,
    ((((ArithmeticFunction.moebius k : ℤ) : ℚ)) / (k : ℚ)) •
      adamsCoef k (PowerSeries.coeff (n / k) g)

/-- Σ_{m ≥ 1} C_m ψ_m(L) with C_m = availablePoints ε m: the coefficient of xⁿ is
Σ_{m ∣ n} C_m ψ_m(L_{n/m}). -/
def countSum (ε : ℤ) (L : Ser) : Ser :=
  PowerSeries.mk fun n ↦ ∑ m ∈ n.divisors,
    (availablePoints ε m).rational • adamsCoef m (PowerSeries.coeff (n / m) L)

/-- Ordinary exponential of a series with zero constant term (composition with Σ xⁿ/n!). -/
def expSer (f : Ser) : Ser := (PowerSeries.exp Coef).subst f

/-- Ordinary logarithm of a series with constant term one. -/
def logSer (g : Ser) : Ser := PowerSeries.logOf g

/-- The formal power g^c = exp(c · log g) for a scalar exponent c ∈ ℚ(q). -/
def formalPow (c : RatFunc ℚ) (g : Ser) : Ser := expSer (PowerSeries.C (scalar c) * logSer g)

/-- Plethystic exponential Exp f = exp(Σ_k ψ_k(f)/k). -/
def pExp (f : Ser) : Ser := expSer (plethSum f)

/-- Plethystic logarithm Log g = Σ_k μ(k) ψ_k(log g)/k. -/
def pLog (g : Ser) : Ser := moebiusSum (logSer g)

/-- The plethystic power P_w(g) = Exp(w · Log g) for a scalar w ∈ ℚ(q). -/
def pPow (w : RatFunc ℚ) (g : Ser) : Ser := pExp (PowerSeries.C (scalar w) * pLog g)

/-- Reciprocal substitution on every coefficient of a series (the proof's f ↦ f^v). -/
def recipSer (f : Ser) : Ser := PowerSeries.map recipCoef f

/-- Dilation S: x ↦ qx. -/
def dilate (f : Ser) : Ser := PowerSeries.rescale (scalar RatFunc.X) f

/-- Inverse dilation S⁻¹: x ↦ x/q. -/
def codilate (f : Ser) : Ser := PowerSeries.rescale (scalar RatFunc.X⁻¹) f

/-- The cluster series B_a = 1 + Σ_{n ≥ 1} F_{S(a;n)} xⁿ of vertical multiplicity a. -/
def clusterSer (a : ℕ) : Ser :=
  PowerSeries.mk fun n ↦ if n = 0 then 1 else lawCoef (stateLaw (State.single a n))

/-- The Euler transform P_{q+ε}(B_a) with q + ε rational points: G_a (ε = 0), P_a (ε = −1),
Q_a (ε = 1). -/
def eulerSer (ε : ℤ) (a : ℕ) : Ser := pPow (RatFunc.X + (ε : RatFunc ℚ)) (clusterSer a)

/-- The off-boundary pattern sum of weight w in ambient dimension N for vertical multiplicity a and
offset ε:
Σ over patterns m of weight w of divisorCount ε m · Π_{(i,j)} ψ_{i+1}(F_{S(a;j+1)})^{m(i,j)}. -/
def patternSum (ε : ℤ) (a N w : ℕ) : Coef :=
  ∑ m ∈ (Finset.univ : Finset (Pattern N)).filter (fun m ↦ patternWeight m = w),
    (divisorCount ε m).rational • ∏ ij : Fin N × Fin N,
      (adamsCoef (ij.1.val + 1) (lawCoef (stateLaw (State.single a (ij.2.val + 1))))) ^ (m ij).val

/-- Contact series: coefficient functions on contact pairs (c,h), i.e. Σ U(c,h) E_{c,h} in the
normal basis. -/
abbrev QSer := ℕ → ℕ → Coef

/-- The twisted product E_{c,h} E_{d,k} = q^{-ck} E_{c+d,h+k}, as a finite convolution. -/
def qmul (U V : QSer) : QSer := fun c h ↦
  ∑ p ∈ Finset.HasAntidiagonal.antidiagonal c, ∑ r ∈ Finset.HasAntidiagonal.antidiagonal h,
    qpow (-((p.1 : ℤ) * (r.2 : ℤ))) * (U p.1 r.1 * V p.2 r.2)

/-- The unit E_{0,0}. -/
def qone : QSer := fun c h ↦ if c = 0 ∧ h = 0 then 1 else 0

/-- Reciprocal substitution on every coefficient (the proof's V^v). -/
def qrecip (U : QSer) : QSer := fun c h ↦ recipCoef (U c h)

/-- The anti-involution †: (f E_{c,h})† = v(f) q^{c+h−ch} E_{c,h}. -/
def qdagger (U : QSer) : QSer := fun c h ↦
  qpow ((c : ℤ) + (h : ℤ) - (c : ℤ) * (h : ℤ)) * recipCoef (U c h)

/-- ι₊: E_{c,h} ↦ q^{-C(c,2)} E_{c,c+h}. -/
def qplus (U : QSer) : QSer := fun c k ↦
  if c ≤ k then qpow (-(chooseTwo c : ℤ)) * U c (k - c) else 0

/-- ι₋: E_{c,h} ↦ q^{-C(h,2)} E_{c+h,h}. -/
def qminus (U : QSer) : QSer := fun k h ↦
  if h ≤ k then qpow (-(chooseTwo h : ℤ)) * U (k - h) h else 0

/-- Diagonal substitution x ↦ J = E_{1,1}: Σ aₘ xᵐ ↦ Σ aₘ Jᵐ with Jᵐ = q^{-C(m,2)} E_{m,m}. -/
def qdiag (f : Ser) : QSer := fun c h ↦
  if c = h then qpow (-(chooseTwo c : ℤ)) * PowerSeries.coeff c f else 0

/-- U_{a,b} = 1 + (q−1) Σ_{c,h ≥ 1} F_{T(a,b;c,h)} E_{c,h}; zero on the two axes. -/
def crossingQ (a b : ℕ) : QSer := fun c h ↦
  if c = 0 ∧ h = 0 then 1
  else if 0 < c ∧ 0 < h then scalar (RatFunc.X - 1) * lawCoef (stateLaw (State.crossing a b c h))
  else 0

/-- Agreement through total contact degree j. -/
def QCongr (j : ℕ) (U V : QSer) : Prop := ∀ c h : ℕ, c + h ≤ j → U c h = V c h

/-- Vanishing on the two axes away from the constant term (support in the positive quadrant plus
`(0,0)`). -/
def QPositive (U : QSer) : Prop := ∀ c h : ℕ, (c = 0 ∨ h = 0) → ¬ (c = 0 ∧ h = 0) → U c h = 0

/-- The representation ρ: ρ(E_{c,h}) f = q^{-tri c} x^{c+h} f(x/q^c); coefficient n sums over
(c+h) + m = n. -/
def qact (U : QSer) (f : Ser) : Ser := PowerSeries.mk fun n ↦
  ∑ p ∈ Finset.HasAntidiagonal.antidiagonal n, ∑ r ∈ Finset.HasAntidiagonal.antidiagonal p.1,
    qpow (-(tri r.1 : ℤ) - (r.1 : ℤ) * (p.2 : ℤ)) * (U r.1 r.2 * PowerSeries.coeff p.2 f)

/-- The same representation with q replaced by q⁻¹ (all Gaussian exponents negated). -/
def qactRecip (U : QSer) (f : Ser) : Ser := PowerSeries.mk fun n ↦
  ∑ p ∈ Finset.HasAntidiagonal.antidiagonal n, ∑ r ∈ Finset.HasAntidiagonal.antidiagonal p.1,
    qpow ((tri r.1 : ℤ) + (r.1 : ℤ) * (p.2 : ℤ)) * (U r.1 r.2 * PowerSeries.coeff p.2 f)

/-- Φ: xᵐ ↦ q^{-tri m} xᵐ. -/
def gaussPhi (f : Ser) : Ser := PowerSeries.mk fun n ↦ qpow (-(tri n : ℤ)) * PowerSeries.coeff n f

/-- Φ⁻¹: xᵐ ↦ q^{tri m} xᵐ. -/
def gaussPhiInv (f : Ser) : Ser := PowerSeries.mk fun n ↦ qpow (tri n : ℤ) * PowerSeries.coeff n f

/-- Transpose of a contact series. -/
def qtranspose (U : QSer) : QSer := fun c h ↦ U h c

/-- The summand of the quadruple sum computing `(ι₊(A) · qdiag(D) · ι₋(C))_{c,h}`: the outer
product is split at `(a, b)` and the inner one at `(a', b')`. -/
def qtripleTerm (A Cq : QSer) (D : Ser) (c h a b a' b' : ℕ) : Coef :=
  if a' ≤ a ∧ b' ≤ b then
    qpow (-((a : ℤ) * ((h - b : ℕ) : ℤ))) *
      ((qpow (-((a' : ℤ) * ((b - b' : ℕ) : ℤ))) * (qplus A a' b' * qdiag D (a - a') (b - b'))) *
        qminus Cq (c - a) (h - b))
  else 0

end TypeDensity