/-
Copyright 2026 The Monstrous Moonshine formalization authors.
Released under the Apache License 2.0.

# The analytic Jacobi product *is* the formal one (task M1-L1f)

`Definitions/MMLevelOneCusp.lean`'s footer names the one missing step of the level-one story: no
lemma in Mathlib (nor in this repository) identified the Taylor coefficients at the cusp of the
**analytic** Jacobi product `MM.levelOneDiscTprod = ∏' n, (1 - eta_q n ·)^24` with the coefficients
of the **formal** product over `ℤ`,

  `MM.formalDiscTprod = ∏' n, ((1 : PowerSeries ℤ) - X ^ (n + 1)) ^ 24`.

This file proves that comparison:

  `MM.qExpansion_levelOneDiscTprod_eq_formal :
     qExpansion 1 levelOneDiscTprod = PowerSeries.map (Int.castRingHom ℂ) formalDiscTprod`.

The content is: **the Fourier coefficients of the analytic product are the coefficients of the
formal product**, i.e. `Δ = q - 24q² + 252q³ - ⋯` is integral because the formal Jacobi product is
integral by construction (`MM.formalDiscTprod_coeff_isInt`).  With
`MM.qExpansion_levelOneDiscTprod` (`qExpansion 1 P = MM.levelOneDiscUnit`) this discharges the
integrality of the unit part `P = Δ/q`, and hence the Δ-side input of the level-one `hq`.

## The route actually taken, and what the brief got wrong

The brief proposed (i) a locally-uniform-limit route iterating `TendstoLocallyUniformlyOn.deriv`
by hand, or (ii) a coefficient-stabilisation route.  Both are viable and the proof below is a
hybrid, but the *engine* of the formal-to-analytic comparison is neither: it is

  **`PowerSeries.eval₂` / `PowerSeries.aeval` (`Mathlib/RingTheory/PowerSeries/Evaluation.lean`),
  the evaluation of a formal power series at a topologically nilpotent point, which is a
  `Continuous` ring homomorphism.**

That is the lemma the footer said did not exist ("no lemma identifies the Taylor coefficients at
`0` of an infinite product of analytic functions with the coefficients of the corresponding formal
product"): although no such comparison is stated for analytic products, `aeval` *is* a continuous
ring hom on `ℤ⟦X⟧` for the `X`-adic (`WithPiTopology`) topology, so it carries the `tprod`
`formalDiscTprod` to the corresponding `tprod` in `ℂ` — the analytic product — with no derivative
bookkeeping at all.  The route below therefore is:

1. `MM.discTprodFun` is the product as a function of `q`, and `MM.levelOneDiscTprod z` is
   `discTprodFun (𝕢 1 z)`;
2. `MM.cuspFunction_levelOneDiscTprod_eventuallyEq`: the cusp function of `P` agrees near `0` with
   `discTprodFun`, so the two have the same `iteratedDeriv` at `0`;
3. `MM.tendstoLocallyUniformlyOn_discTprodPartial`: the partial products converge locally uniformly
   to `discTprodFun` on `ball 0 (1/2)` — this is where the local-uniform input lives;
4. `MM.tendstoLocallyUniformlyOn_iteratedDeriv_discTprodPartial`: `TendstoLocallyUniformlyOn.deriv`
   iterated by hand (Mathlib has no iterated version), which is the step the file footer calls "not
   small" — it is about 25 lines;
5. `MM.iteratedDeriv_polynomial_eval`: for a polynomial `p`, `iteratedDeriv N (p.eval ·) 0 =
   N! * p.coeff N` — the stabilisation, done through `Polynomial.eval_eq_sum_range` and
   `iteratedDeriv_fun_pow_zero`;
6. `MM.coeff_prod_range_eq_formalDiscTprod`: the `N`-th coefficient of the formal product is
   already the `N`-th coefficient of its first `N + 1` factors (`MM.coeff_formalDiscTprod_eq_partial`
   of `Definitions/MMLevelOneIntegrality.lean`, used through `Finset.prod_range_add`).

**A fact that is missing from Mathlib and had to be worked around.**  `PowerSeries.eval₂Hom`
requires the *target* to carry a linear topology (`IsLinearTopology S S`).  There is no
`IsLinearTopology ℂ ℂ` (the submodules of the field `ℂ` are `⊥` and `⊤`, so a linear topology on
`ℂ` is discrete), so `PowerSeries.aeval` **cannot** be applied with target `ℂ`.  It applies with
target `ℤ⟦X⟧`-like rings but not to the analytic product.  That is why the route above uses the
locally uniform limit after all; the two ingredients that make it short are
`Summable.hasProdLocallyUniformlyOn_nat_one_add` (for step 3) and `Polynomial.eval_eq_sum_range`
(for step 5). -/

import Definitions.MMLevelOneIntegrality
import Mathlib.Analysis.Complex.LocallyUniformLimit
import Mathlib.Analysis.Normed.Module.MultipliableUniformlyOn

noncomputable section

set_option autoImplicit false
set_option maxHeartbeats 1000000

open scoped Classical Topology PowerSeries.WithPiTopology
open Finset Filter UpperHalfPlane

namespace MM

local notation "𝕢" => Periodic.qParam
local notation "Δ" => ModularForm.discriminant

/-! ## The analytic product as a function of the local parameter `q` -/

/-- **The analytic Jacobi product as a function of `q`**: `∏' n, (1 - q^(n+1))^24`, the value of
`MM.levelOneDiscTprod` at the point of `ℍ` whose local parameter is `q`. -/
noncomputable def discTprodFun (q : ℂ) : ℂ :=
  ∏' n : ℕ, (1 - q ^ (n + 1)) ^ 24

/-- **The `K`-th partial product** of the analytic Jacobi product, as a function of `q`. -/
noncomputable def discTprodPartial (K : ℕ) (q : ℂ) : ℂ :=
  ∏ n ∈ Finset.range K, (1 - q ^ (n + 1)) ^ 24

/-- The partial product, bundled as a polynomial, so that its Taylor coefficients at `0` can be
read off from `Polynomial.eval_eq_sum_range`. -/
noncomputable def discPoly (K : ℕ) : Polynomial ℂ :=
  ∏ n ∈ Finset.range K, ((1 : Polynomial ℂ) - Polynomial.X ^ (n + 1)) ^ 24

theorem levelOneDiscTprod_eq_discTprodFun (z : UpperHalfPlane) :
    levelOneDiscTprod z = discTprodFun (𝕢 1 (z : ℂ)) := by
  simp only [levelOneDiscTprod, discTprodFun, ModularForm.eta_q]

theorem discTprodPartial_eq_discPoly_eval (K : ℕ) :
    discTprodPartial K = fun q : ℂ => (discPoly K).eval q := by
  funext q
  simp only [discTprodPartial, discPoly]
  rw [Polynomial.eval_prod]
  exact Finset.prod_congr rfl fun n _ => by simp

/-! ## The stabilisation at the finite level

The `N`-th Taylor coefficient at `0` of a polynomial function is `N!` times the `N`-th coefficient
of the polynomial.  This is the analytic counterpart of `MM.coeff_formalDiscTprod_eq_partial` of
`Definitions/MMLevelOneIntegrality.lean`; it is what makes the analytic partial products comparable
with the formal product without enumerating coefficients. -/

/-- **The `N`-th Taylor coefficient of a polynomial function.**  For `p : ℂ[X]`,
`iteratedDeriv N (fun q => p.eval q) 0 = N! * p.coeff N`. -/
theorem iteratedDeriv_polynomial_eval (p : Polynomial ℂ) (N : ℕ) :
    iteratedDeriv N (fun q : ℂ => p.eval q) 0
      = (N.factorial : ℂ) * (p : PowerSeries ℂ).coeff N := by
  have hfun : (fun q : ℂ => p.eval q)
      = fun q : ℂ => ∑ i ∈ Finset.range (p.natDegree + 1), p.coeff i * q ^ i := by
    funext q
    exact Polynomial.eval_eq_sum_range q
  have hderiv : iteratedDeriv N
        (fun q : ℂ => ∑ i ∈ Finset.range (p.natDegree + 1), p.coeff i * q ^ i) 0
      = ∑ i ∈ Finset.range (p.natDegree + 1),
          iteratedDeriv N (fun q : ℂ => p.coeff i * q ^ i) 0 :=
    iteratedDeriv_fun_sum fun _ _ => by fun_prop
  rw [hfun, hderiv]
  simp only [iteratedDeriv_const_mul_field, iteratedDeriv_fun_pow_zero, Polynomial.coeff_coe]
  by_cases hN : N ≤ p.natDegree
  · refine (Finset.sum_eq_single N ?_ ?_).trans ?_
    · intro i _ hi
      rw [if_neg (fun h => hi h.symm), mul_zero]
    · intro h
      exact absurd (Finset.mem_range.mpr (Nat.lt_succ_of_le hN)) h
    · rw [if_pos rfl]; ring
  · push_neg at hN
    have hzero : ∀ i ∈ Finset.range (p.natDegree + 1),
        p.coeff i * (if N = i then (i.factorial : ℂ) else 0) = 0 := by
      intro i hi
      have hi' : i ≤ p.natDegree := by have := Finset.mem_range.mp hi; omega
      rw [if_neg (by omega), mul_zero]
    rw [Finset.sum_eq_zero hzero, Polynomial.coeff_eq_zero_of_natDegree_lt hN, mul_zero]

/-- **The formal stabilisation for an arbitrary cutoff.**  For `K ≥ N + 1` the `N`-th coefficient of
the `K`-th formal partial product is already the `N`-th coefficient of the full product
`MM.formalDiscTprod` (`MM.coeff_formalDiscTprod_eq_partial` is the case `K = N + 1`). -/
theorem coeff_prod_range_eq_formalDiscTprod {K N : ℕ} (hK : N + 1 ≤ K) :
    (∏ n ∈ Finset.range K, ((1 : PowerSeries ℤ) - PowerSeries.X ^ (n + 1)) ^ 24).coeff N
      = formalDiscTprod.coeff N := by
  obtain ⟨j, rfl⟩ : ∃ j, K = N + 1 + j := ⟨K - (N + 1), by omega⟩
  rw [Finset.prod_range_add,
    coeff_mul_eq_of_coeff_eq _ _ N (fun i hi => coeff_prod_range_tail N j i hi),
    coeff_formalDiscTprod_eq_partial]

/-- **The stabilisation, in the form the comparison uses.**  For `K ≥ N + 1`, the `N`-th Taylor
coefficient at `0` of the `K`-th analytic partial product is `N!` times the `N`-th coefficient of
the formal product. -/
theorem iteratedDeriv_discTprodPartial_eq (K N : ℕ) (hK : N + 1 ≤ K) :
    iteratedDeriv N (discTprodPartial K) 0
      = (N.factorial : ℂ) * (PowerSeries.map (Int.castRingHom ℂ) formalDiscTprod).coeff N := by
  rw [discTprodPartial_eq_discPoly_eval, iteratedDeriv_polynomial_eval]
  congr 1
  rw [← PowerSeries.coeff_map]
  congr 1
  rw [← coeff_prod_range_eq_formalDiscTprod hK]
  push_cast
  rw [Polynomial.coeToPowerSeries.ringHom_apply, map_prod]
  refine Finset.prod_congr rfl fun n _ => ?_
  rw [map_pow, map_sub, map_one, map_mul, Polynomial.coeToPowerSeries.ringHom_apply]
  simp

/-! ## Local uniform convergence of the partial products

The engine of the comparison: the partial products converge locally uniformly to the analytic
product on the disc of radius `1/2`.  The radius is chosen so that the geometric domination
`‖(1 - q^(n+1))^24 - 1‖ ≤ 24·2^24·‖q‖^(n+1)` is uniform, which is exactly the hypothesis of
`Summable.hasProdLocallyUniformlyOn_nat_one_add`. -/

theorem norm_one_sub_pow_pow_sub_one_le (m : ℕ) {q : ℂ} (hq : ‖q‖ < 1) :
    ‖(1 - q ^ m) ^ 24 - 1‖ ≤ 24 * 2 ^ 24 * ‖q‖ ^ m := by
  have hz : ‖q ^ m‖ < 1 := by
    rw [norm_pow]
    exact lt_of_le_of_lt (pow_le_one₀ (norm_nonneg q) hq.le) (by norm_num)
  have hfac : (1 - q ^ m) ^ 24 - 1
      = (-(q ^ m)) * ∑ i ∈ Finset.range 24, (1 - q ^ m) ^ i := by
    rw [← geom_sum_mul (1 - q ^ m) 24,
      show (1 - q ^ m) - 1 = -(q ^ m) by ring]
  have h2 : ‖1 - q ^ m‖ ≤ 2 := by
    have h := norm_sub_le (1 : ℂ) (q ^ m)
    rw [norm_one] at h
    linarith [norm_nonneg (q ^ m), hz.le]
  have hnorm : ‖∑ i ∈ Finset.range 24, (1 - q ^ m) ^ i‖ ≤ 24 * 2 ^ 24 := by
    refine (norm_sum_le _ _).trans ?_
    calc ∑ i ∈ Finset.range 24, ‖(1 - q ^ m) ^ i‖
        = ∑ i ∈ Finset.range 24, ‖1 - q ^ m‖ ^ i := by simp only [norm_pow]
      _ ≤ ∑ _i ∈ Finset.range 24, (2 : ℝ) ^ 24 := by
          refine Finset.sum_le_sum fun i hi => ?_
          calc ‖1 - q ^ m‖ ^ i ≤ (2 : ℝ) ^ i := pow_le_pow_left₀ (norm_nonneg _) h2 i
            _ ≤ (2 : ℝ) ^ 24 := pow_le_pow_right₀ (by norm_num) (by have := Finset.mem_range.mp hi; omega)
      _ = 24 * 2 ^ 24 := by simp
  rw [hfac, norm_mul, norm_neg, norm_pow]
  calc ‖q‖ ^ m * ‖∑ i ∈ Finset.range 24, (1 - q ^ m) ^ i‖
      ≤ ‖q‖ ^ m * (24 * 2 ^ 24) := mul_le_mul_of_nonneg_left hnorm (by positivity)
    _ ≤ (1 : ℝ) * (24 * 2 ^ 24) := by
        exact mul_le_mul_of_nonneg_right (pow_le_one₀ (norm_nonneg q) hq.le) (by positivity)
    _ = 24 * 2 ^ 24 * ‖q‖ ^ m := by rw [one_mul]; ring

theorem tendstoLocallyUniformlyOn_discTprodPartial :
    TendstoLocallyUniformlyOn (fun K : ℕ => discTprodPartial K) discTprodFun atTop
      (Metric.ball (0 : ℂ) (1 / 2)) := by
  have hu : Summable (fun n : ℕ => (24 * 2 ^ 24 : ℝ) * (1 / 2 : ℝ) ^ (n + 1)) :=
    ((summable_nat_add_iff 1).mpr
      (summable_geometric_of_lt_one (by norm_num) (by norm_num))).mul_left _
  have hbase : HasProdLocallyUniformlyOn
      (fun n : ℕ => fun x : ℂ => (1 - x ^ (n + 1)) ^ 24)
      discTprodFun (Metric.ball (0 : ℂ) (1 / 2)) := by
    have h := hu.hasProdLocallyUniformlyOn_nat_one_add (K := Metric.ball (0 : ℂ) (1 / 2))
      (f := fun n : ℕ => fun x : ℂ => (1 - x ^ (n + 1)) ^ 24 - 1) Metric.isOpen_ball ?_ ?_
    · simpa only [add_sub_cancel, discTprodFun] using h
    · filter_upwards with n
      intro x hx
      have hx1 : ‖x‖ < 1 / 2 := by simpa [Metric.mem_ball, dist_zero_right] using hx
      have hx2 : ‖x ^ (n + 1)‖ < 1 := by
        rw [norm_pow]
        exact lt_of_le_of_lt (pow_le_one₀ (norm_nonneg x) hx1.le) (by norm_num)
      refine (norm_one_sub_pow_pow_sub_one_le (n + 1) hx2).trans ?_
      rw [norm_pow]
      calc ‖x‖ ^ (n + 1) * (24 * 2 ^ 24) ≤ (1 / 2 : ℝ) ^ (n + 1) * (24 * 2 ^ 24) :=
            mul_le_mul_of_nonneg_right (pow_le_pow_left₀ (norm_nonneg x) hx1.le _) (by positivity)
        _ = (24 * 2 ^ 24 : ℝ) * (1 / 2 : ℝ) ^ (n + 1) := by ring
    · intro n
      fun_prop
  intro u hu' x hx
  obtain ⟨t, ht, hFt⟩ := hbase u hu' x hx
  exact ⟨t, ht, tendsto_finset_range.eventually hFt⟩

/-! ## Iterated derivatives of the limit

`Complex.TendstoLocallyUniformlyOn.deriv` says that locally uniform convergence implies locally
uniform convergence of the derivatives.  Mathlib has no iterated version, so it is iterated here by
hand; note that no differentiability of the *limit* is needed, only of the approximants, which are
polynomials. -/

theorem tendstoLocallyUniformlyOn_iteratedDeriv_discTprodPartial (N : ℕ) :
    TendstoLocallyUniformlyOn (fun K : ℕ => iteratedDeriv N (discTprodPartial K))
      (iteratedDeriv N discTprodFun) atTop (Metric.ball (0 : ℂ) (1 / 2)) := by
  induction N with
  | zero => simpa only [iteratedDeriv_zero] using tendstoLocallyUniformlyOn_discTprodPartial
  | succ N IH =>
    have hdiff : ∀ᶠ K : ℕ in atTop,
        DifferentiableOn ℂ (iteratedDeriv N (discTprodPartial K)) (Metric.ball (0 : ℂ) (1 / 2)) :=
      Eventually.of_forall fun K => by
        rw [discTprodPartial_eq_discPoly_eval]
        exact ((contDiff_aeval (discPoly K) (N + 1 : WithTop ℕ∞)).differentiable_iteratedDeriv'
          N).differentiableOn
    simpa only [Function.comp_apply, ← iteratedDeriv_succ] using
      IH.deriv hdiff Metric.isOpen_ball

/-- **The limit of the Taylor coefficients is the Taylor coefficient of the limit**, in the form
needed: the `N`-th derivative at `0` of the analytic product is `N!` times the `N`-th coefficient of
the formal product. -/
theorem iteratedDeriv_discTprodFun (N : ℕ) :
    iteratedDeriv N discTprodFun 0
      = (N.factorial : ℂ) * (PowerSeries.map (Int.castRingHom ℂ) formalDiscTprod).coeff N := by
  refine tendsto_nhds_unique
    ((tendstoLocallyUniformlyOn_iteratedDeriv_discTprodPartial N).tendsto_at
      (Metric.mem_ball_self (by norm_num)) |>.congr' ?_) tendsto_const_nhds
  filter_upwards [eventually_atTop.mpr ⟨N + 1, fun K hK => hK⟩] with K hK
  exact iteratedDeriv_discTprodPartial_eq K N hK

/-! ## The cusp function of `P` is the analytic product -/

/-- **The cusp function of `P` is the analytic product.**  For `q` in the unit disc,
`cuspFunction 1 levelOneDiscTprod q = ∏' n, (1 - q^(n+1))^24`, because `q` is the local parameter:
`q = 𝕢 1 τ` for `τ = invQParam 1 q ∈ ℍ`.  Only the germ at `0` is needed. -/
theorem cuspFunction_levelOneDiscTprod_eventuallyEq :
    cuspFunction 1 levelOneDiscTprod =ᶠ[𝓝 (0 : ℂ)] discTprodFun := by
  filter_upwards [Metric.ball_mem_nhds (0 : ℂ) one_pos] with q hq
  by_cases hq0 : q = 0
  · subst hq0
    rw [cuspFunction_apply_zero one_pos analyticAt_cuspFunction_levelOneDiscTprod
      periodic_levelOneDiscTprod, Tendsto.limUnder_eq tendsto_levelOneDiscTprod]
    simp [discTprodFun]
  · have him := Periodic.im_invQParam_pos_of_norm_lt_one one_pos (by simpa using hq) hq0
    have hper := periodic_levelOneDiscTprod
    simp only [cuspFunction, Periodic.cuspFunction_eq_of_nonzero 1 _ hq0, Function.comp_apply,
      ofComplex_apply_of_im_pos him, levelOneDiscTprod_eq_discTprodFun,
      Periodic.qParam_right_inv one_ne_zero hq0]

/-! ## The comparison -/

/-- **The analytic `q`-expansion of the Jacobi product is the formal product over `ℤ`.**  This is
the missing step named at the end of `Definitions/MMLevelOneCusp.lean`: the Taylor coefficients at
the cusp of `P = Δ/q = ∏' n, (1 - eta_q n ·)^24` are the coefficients of
`∏' n, ((1 : ℤ⟦X⟧) - X^(n+1))^24`, so they are integers by construction. -/
theorem qExpansion_levelOneDiscTprod_eq_formal :
    qExpansion 1 levelOneDiscTprod
      = PowerSeries.map (Int.castRingHom ℂ) formalDiscTprod := by
  ext N
  rw [UpperHalfPlane.qExpansion_coeff,
    Filter.EventuallyEq.iteratedDeriv_eq N cuspFunction_levelOneDiscTprod_eventuallyEq,
    iteratedDeriv_discTprodFun]
  exact inv_mul_cancel_left₀ (Nat.cast_ne_zero.mpr (Nat.factorial_ne_zero N)) _

/-! ## The `Δ`-side consequences -/

/-- `Δ = 𝕢 · P`, as `q`-expansions: `qExpansion 1 Δ = X * qExpansion 1 P`. -/
theorem qExpansion_discriminant_eq_X_mul :
    qExpansion 1 Δ = PowerSeries.X * qExpansion 1 levelOneDiscTprod := by
  have hfun : Δ = (fun τ : UpperHalfPlane => 𝕢 1 (τ : ℂ)) * levelOneDiscTprod := by
    funext z
    rw [Pi.mul_apply]
    exact ModularForm.discriminant_eq_q_prod z
  rw [hfun, UpperHalfPlane.qExpansion_mul analyticAt_cuspFunction_qParam
    analyticAt_cuspFunction_levelOneDiscTprod, qExpansion_qParam_one]

/-- **The `q`-expansion of `Δ` is `X` times the formal Jacobi product.**  This is the Δ-side of the
comparison: `Δ = q·∏' n, (1-q^(n+1))^24` with the product's coefficients integral. -/
theorem qExpansion_discriminant_eq_X_mul_formal :
    qExpansion 1 Δ = PowerSeries.X * PowerSeries.map (Int.castRingHom ℂ) formalDiscTprod := by
  rw [qExpansion_discriminant_eq_X_mul, qExpansion_levelOneDiscTprod_eq_formal]

/-- **The unit part of `Δ` has integer coefficients.**  `MM.levelOneDiscUnit.coeff n` is the
coefficient of `q^(n+1)` in `qExpansion 1 Δ`, hence the coefficient of `X^n` in the formal Jacobi
product over `ℤ`. -/
theorem levelOneDiscUnit_coeff_isInt (n : ℕ) :
    ∃ z : ℤ, levelOneDiscUnit.coeff n = (z : ℂ) := by
  obtain ⟨z, hz⟩ := formalDiscTprod_coeff_isInt n
  exact ⟨z, by rw [← qExpansion_levelOneDiscTprod, qExpansion_levelOneDiscTprod_eq_formal]; exact hz⟩

end MM
