2607.19337v1 / ChampernowneA.lean

all files

import Mathlib

/-!
# The Champernowne-type element of the finite adeles

This file formalizes the main theorem of Seki, *On the A-transcendence of a
Champernowne-type constant*.  The only external inputs are isolated in the two
propositions `PrimeNumberTheorem` and `MaynardTaoBoundedGaps` below.
-/

open Filter Finset
open scoped BigOperators Topology Polynomial Polynomial.Bivariate
open Asymptotics

namespace ChampernowneA

/-- The zero-indexed sequence `2, 3, 5, ...` of prime numbers. -/
noncomputable def nthPrime (n : ℕ) : ℕ := Nat.nth Nat.Prime n

/-- The prime number theorem, in precisely the form used in the paper. -/
def PrimeNumberTheorem : Prop :=
  (fun n : ℕ ↦ (nthPrime n : ℝ)) ~[atTop]
    (fun n : ℕ ↦ (n : ℝ) * Real.log n)

/-- The bounded-gap consequence of the Maynard--Tao theorem used in the paper.

For each `k`, infinitely many blocks of `k + 1` consecutive primes have bounded
diameter.  This is equivalent to the finiteness of the relevant liminf. -/
def MaynardTaoBoundedGaps : Prop :=
  ∀ k : ℕ, ∃ B : ℕ,
    {n : ℕ | nthPrime (n + k) - nthPrime n ≤ B}.Infinite

lemma nthPrime_strictMono : StrictMono nthPrime :=
  Nat.nth_strictMono Nat.infinite_setOfPred_prime

lemma nthPrime_mono : Monotone nthPrime := nthPrime_strictMono.monotone

lemma nthPrime_prime (n : ℕ) : Nat.Prime (nthPrime n) :=
  Nat.prime_nth_prime n

lemma nthPrime_ge (n : ℕ) : n + 2 ≤ nthPrime n :=
  Nat.add_two_le_nth_prime n

/-- A fixed pattern of `k + 1` consecutive prime offsets which occurs
infinitely often.  This is the paper's pigeonhole step following Maynard--Tao. -/
lemma fixed_prime_cluster (hMT : MaynardTaoBoundedGaps) (k : ℕ) :
    ∃ (h : Fin (k + 1) → ℕ) (M : Set ℕ),
      M.Infinite ∧ h 0 = 0 ∧ StrictMono h ∧
        ∀ n ∈ M, ∀ j : Fin (k + 1),
          nthPrime (n + j) = nthPrime n + h j := by
  obtain ⟨B, hB⟩ := hMT k
  let S : Set ℕ := {n : ℕ | nthPrime (n + k) - nthPrime n ≤ B}
  have hS : S.Infinite := hB
  let offset : S → (Fin (k + 1) → Fin (B + 1)) := fun n j ↦
    ⟨nthPrime (n.1 + j) - nthPrime n.1, by
      have hjk : (j : ℕ) ≤ k := Nat.le_of_lt_succ j.isLt
      have hmono : nthPrime (n.1 + j) ≤ nthPrime (n.1 + k) :=
        nthPrime_mono (Nat.add_le_add_left hjk n.1)
      have hbase : nthPrime n.1 ≤ nthPrime (n.1 + j) :=
        nthPrime_mono (Nat.le_add_right n.1 j)
      have hdiam : nthPrime (n.1 + k) - nthPrime n.1 ≤ B := n.2
      omega⟩
  letI : Infinite S := hS.to_subtype
  obtain ⟨a, ha⟩ := Finite.exists_infinite_fiber offset
  let M : Set ℕ := {n | ∃ hn : n ∈ S, offset ⟨n, hn⟩ = a}
  have hM : M.Infinite := by
    letI : Infinite {n : S | offset n = a} := ha
    let e : {n : S | offset n = a} ↪ M :=
      ⟨fun n ↦ ⟨n.1.1, n.1.2, n.2⟩, fun x y hxy ↦ by
        apply Subtype.ext
        apply Subtype.ext
        exact congrArg (fun z : M ↦ z.1) hxy⟩
    exact Set.infinite_coe_iff.mp (Infinite.of_injective e e.injective)
  let h : Fin (k + 1) → ℕ := fun j ↦ a j
  refine ⟨h, M, hM, ?_, ?_, ?_⟩
  · obtain ⟨n, hn⟩ := hM.nonempty
    obtain ⟨hnS, hna⟩ := hn
    have hz := congrFun hna (0 : Fin (k + 1))
    simpa [offset, h] using hz.symm
  · intro i j hij
    obtain ⟨n, hn⟩ := hM.nonempty
    obtain ⟨hnS, hna⟩ := hn
    have hi := congrFun hna i
    have hj := congrFun hna j
    simp only [offset] at hi hj
    dsimp [h]
    rw [← hi, ← hj]
    exact Nat.sub_lt_sub_right
      (nthPrime_mono (Nat.le_add_right n i))
      (nthPrime_strictMono (Nat.add_lt_add_left hij n))
  · intro n hn j
    obtain ⟨hnS, hna⟩ := hn
    have hj := congrFun hna j
    simp only [offset] at hj
    dsimp [h]
    rw [← hj]
    exact (Nat.add_sub_of_le
      (nthPrime_mono (Nat.le_add_right n j))).symm

section PrimeAsymptotics

/-- The slowly growing factor `pₙ / n`.  Its value at zero is immaterial. -/
noncomputable def primeRatio (n : ℕ) : ℝ := (nthPrime n : ℝ) / n

lemma pnt_primeRatio (hPNT : PrimeNumberTheorem) :
    primeRatio ~[atTop] (fun n : ℕ ↦ Real.log n) := by
  have hn : (fun n : ℕ ↦ (n : ℝ)) ~[atTop] (fun n : ℕ ↦ (n : ℝ)) :=
    IsEquivalent.refl
  have hdiv := hPNT.div hn
  refine (hdiv.congr_left (Eventually.of_forall fun n ↦ rfl)).congr_right ?_
  filter_upwards [eventually_ne_atTop (0 : ℕ)] with n hn0
  simp [hn0]

lemma tendsto_primeRatio_atTop (hPNT : PrimeNumberTheorem) :
    Tendsto primeRatio atTop atTop := by
  apply (pnt_primeRatio hPNT).tendsto_atTop_iff.mpr
  exact Real.tendsto_log_atTop.comp tendsto_natCast_atTop_atTop

/-- Every fixed power of `pₙ/n` is `o(n)`.  This is the only upper-growth
consequence of PNT needed for polynomial nonvanishing. -/
lemma primeRatio_pow_isLittleO (hPNT : PrimeNumberTheorem) (m : ℕ) :
    (fun n : ℕ ↦ primeRatio n ^ m) =o[atTop] (fun n : ℕ ↦ (n : ℝ)) := by
  have hlog :
      (fun n : ℕ ↦ Real.log (n : ℝ) ^ m) =o[atTop] (fun n : ℕ ↦ (n : ℝ)) :=
    Real.isLittleO_pow_log_id_atTop.comp_tendsto tendsto_natCast_atTop_atTop
  exact (pnt_primeRatio hPNT).pow m |>.trans_isLittleO hlog

private lemma natPow_isBigO_natPow {a b : ℕ} (hab : a ≤ b) :
    (fun n : ℕ ↦ (n : ℝ) ^ a) =O[atTop] (fun n : ℕ ↦ (n : ℝ) ^ b) := by
  apply IsBigO.of_bound' (eventually_ge_atTop (1 : ℕ) |>.mono fun n hn ↦ ?_)
  rw [Real.norm_eq_abs, Real.norm_eq_abs, abs_of_nonneg (by positivity),
    abs_of_nonneg (by positivity)]
  exact pow_le_pow_right₀ (by exact_mod_cast hn) hab

private lemma natPow_isBigO_natPow_mul_logPow (D J : ℕ) :
    (fun n : ℕ ↦ (n : ℝ) ^ D) =O[atTop]
      (fun n : ℕ ↦ (n : ℝ) ^ D * Real.log n ^ J) := by
  apply IsBigO.of_bound' (eventually_ge_atTop (3 : ℕ) |>.mono fun n hn ↦ ?_)
  rw [Real.norm_eq_abs, Real.norm_eq_abs, abs_of_nonneg (by positivity),
    abs_of_nonneg (mul_nonneg (by positivity)
      (pow_nonneg (Real.log_nonneg (by norm_cast; omega)) J))]
  nth_rw 1 [← mul_one ((n : ℝ) ^ D)]
  gcongr
  apply one_le_pow₀
  apply (Real.le_log_iff_exp_le (by positivity)).mpr
  exact Real.exp_one_lt_three.le.trans (by exact_mod_cast hn)

/-- Lexicographic dominance of the monomials `n^D (log n)^J`, where total
degree is compared first and log-degree second. -/
lemma logMonomial_isLittleO {t j D J : ℕ}
    (h : t < D ∨ t = D ∧ j < J) :
    (fun n : ℕ ↦ (n : ℝ) ^ t * Real.log n ^ j) =o[atTop]
      (fun n : ℕ ↦ (n : ℝ) ^ D * Real.log n ^ J) := by
  rcases h with ht | ⟨rfl, hj⟩
  · have hlog :
        (fun n : ℕ ↦ Real.log (n : ℝ) ^ j) =o[atTop] (fun n : ℕ ↦ (n : ℝ)) :=
      Real.isLittleO_pow_log_id_atTop.comp_tendsto tendsto_natCast_atTop_atTop
    have h₁ := (isBigO_refl (fun n : ℕ ↦ (n : ℝ) ^ t) atTop).mul_isLittleO hlog
    have h₂ :
        (fun n : ℕ ↦ (n : ℝ) ^ t * (n : ℝ)) =O[atTop]
          (fun n : ℕ ↦ (n : ℝ) ^ D) := by
      have := natPow_isBigO_natPow (Nat.succ_le_iff.mpr ht)
      exact this.congr' (Eventually.of_forall fun n ↦ by simp [pow_succ]) (Eventually.of_forall fun _ ↦ rfl)
    exact h₁.trans_isBigO (h₂.trans (natPow_isBigO_natPow_mul_logPow D J))
  · have hp :
        (fun x : ℝ ↦ x ^ j) =o[atTop] (fun x : ℝ ↦ x ^ J) :=
      isLittleO_pow_pow_atTop_of_lt hj
    have hlog := hp.comp_tendsto
      (Real.tendsto_log_atTop.comp tendsto_natCast_atTop_atTop)
    exact (isBigO_refl (fun n : ℕ ↦ (n : ℝ) ^ t) atTop).mul_isLittleO hlog

/-- PNT converts a monomial in `(n,pₙ)` to its corresponding monomial in
`(n,log n)`. -/
lemma primeMonomial_isEquivalent (hPNT : PrimeNumberTheorem) (i j : ℕ) :
    (fun n : ℕ ↦ (n : ℝ) ^ i * (nthPrime n : ℝ) ^ j) ~[atTop]
      (fun n : ℕ ↦ (n : ℝ) ^ (i + j) * Real.log n ^ j) := by
  have hi :
      (fun n : ℕ ↦ (n : ℝ) ^ i) ~[atTop] (fun n : ℕ ↦ (n : ℝ) ^ i) :=
    IsEquivalent.refl
  have h := hi.mul (hPNT.pow j)
  refine h.congr_right (Eventually.of_forall fun n ↦ ?_)
  simp only [Pi.mul_apply, Pi.pow_apply]
  rw [mul_pow, ← mul_assoc, ← pow_add]

lemma primeMonomial_isLittleO (hPNT : PrimeNumberTheorem)
    {i j I J : ℕ} (h : i + j < I + J ∨ i + j = I + J ∧ j < J) :
    (fun n : ℕ ↦ (n : ℝ) ^ i * (nthPrime n : ℝ) ^ j) =o[atTop]
      (fun n : ℕ ↦ (n : ℝ) ^ I * (nthPrime n : ℝ) ^ J) := by
  have hsmall := logMonomial_isLittleO (t := i + j) (j := j)
    (D := I + J) (J := J) h
  have hactual :
      (fun n : ℕ ↦ (n : ℝ) ^ i * (nthPrime n : ℝ) ^ j) =o[atTop]
        (fun n : ℕ ↦ (n : ℝ) ^ (I + J) * Real.log n ^ J) :=
    (primeMonomial_isEquivalent hPNT i j).trans_isLittleO hsmall
  exact hactual.trans_isBigO
    (primeMonomial_isEquivalent hPNT I J).isBigO_symm

/-- The ordering on bivariate monomials used in the paper: total degree first,
then the exponent of the prime variable. -/
def monomialKey (q : ℕ × ℕ) : Lex (ℕ × ℕ) :=
  toLex (q.1 + q.2, q.2)

lemma monomialKey_injective : Function.Injective monomialKey := by
  rintro ⟨i, j⟩ ⟨I, J⟩ h
  change toLex (i + j, j) = toLex (I + J, J) at h
  have h' : (i + j, j) = (I + J, J) := congrArg ofLex h
  have hsum : i + j = I + J := congrArg Prod.fst h'
  have hj : j = J := congrArg Prod.snd h'
  apply Prod.ext
  · simp only
    omega
  · simpa only using hj

/-- A nonempty finite linear combination of distinct monomials in `(n,pₙ)` is
eventually nonzero.  This packages equation (2.3) of the paper. -/
lemma finite_bivariate_sum_eventually_ne_zero (hPNT : PrimeNumberTheorem)
    (s : Finset (ℕ × ℕ)) (a : ℕ × ℕ → ℝ)
    (hs : s.Nonempty) (ha : ∀ q ∈ s, a q ≠ 0) :
    ∀ᶠ n : ℕ in atTop,
      ∑ q ∈ s, a q * ((n : ℝ) ^ q.1 * (nthPrime n : ℝ) ^ q.2) ≠ 0 := by
  let keys : Finset (Lex (ℕ × ℕ)) := s.image monomialKey
  have hkeys : keys.Nonempty := hs.image _
  let K := keys.max' hkeys
  have hKmem : K ∈ keys := keys.max'_mem hkeys
  obtain ⟨Q, hQs, hQK⟩ := Finset.mem_image.mp hKmem
  have hdom : ∀ q ∈ s, q ≠ Q →
      q.1 + q.2 < Q.1 + Q.2 ∨
        q.1 + q.2 = Q.1 + Q.2 ∧ q.2 < Q.2 := by
    intro q hqs hqQ
    have hkq : monomialKey q ∈ keys := Finset.mem_image.mpr ⟨q, hqs, rfl⟩
    have hle : monomialKey q ≤ K := Finset.le_max' keys _ hkq
    rw [← hQK] at hle
    have hne : monomialKey q ≠ monomialKey Q :=
      fun h ↦ hqQ (monomialKey_injective h)
    have hlt : monomialKey q < monomialKey Q := lt_of_le_of_ne hle hne
    change toLex (q.1 + q.2, q.2) < toLex (Q.1 + Q.2, Q.2) at hlt
    simpa only [Prod.Lex.lt_iff, ofLex_toLex] using hlt
  let lead : ℕ → ℝ := fun n ↦
    a Q * ((n : ℝ) ^ Q.1 * (nthPrime n : ℝ) ^ Q.2)
  have hrest :
      (fun n : ℕ ↦
        ∑ q ∈ s.erase Q,
          a q * ((n : ℝ) ^ q.1 * (nthPrime n : ℝ) ^ q.2)) =o[atTop] lead := by
    apply IsLittleO.sum
    intro q hq
    simp only [mem_erase] at hq
    have hmono := primeMonomial_isLittleO hPNT (hdom q hq.2 hq.1)
    exact (hmono.const_mul_left (a q)).const_mul_right (ha Q hQs)
  have hequiv :
      (fun n : ℕ ↦
        ∑ q ∈ s, a q * ((n : ℝ) ^ q.1 * (nthPrime n : ℝ) ^ q.2)) ~[atTop] lead := by
    have h₁ := hrest.add_isEquivalent (IsEquivalent.refl : lead ~[atTop] lead)
    refine h₁.congr_left (Eventually.of_forall fun n ↦ ?_)
    simpa only [Pi.add_apply, lead] using (Finset.sum_erase_add s _ hQs)
  obtain ⟨φ, hφ, heq⟩ := hequiv.exists_eq_mul
  have hφpos : ∀ᶠ n in atTop, (0 : ℝ) < φ n :=
    (tendsto_order.mp hφ).1 0 zero_lt_one
  filter_upwards [heq, hφpos, eventually_ne_atTop (0 : ℕ)] with n heq hφpos hn
  rw [heq]
  apply mul_ne_zero hφpos.ne'
  dsimp [lead]
  apply mul_ne_zero (ha Q hQs)
  apply mul_ne_zero
  · exact pow_ne_zero _ (Nat.cast_ne_zero.mpr hn)
  · exact pow_ne_zero _ (Nat.cast_ne_zero.mpr (nthPrime_prime n).ne_zero)

section BivariatePolynomial

variable {R : Type*} [CommSemiring R]

/-- The pairs `(i,j)` for which the coefficient of `X^i Y^j` is nonzero. -/
noncomputable def bivariateSupport (P : R[X][Y]) : Finset (ℕ × ℕ) :=
  P.support.biUnion fun j ↦
    (P.coeff j).support.image fun i ↦ (i, j)

lemma mem_bivariateSupport {P : R[X][Y]} {q : ℕ × ℕ} :
    q ∈ bivariateSupport P ↔ (P.coeff q.2).coeff q.1 ≠ 0 := by
  constructor
  · intro hq
    simp only [bivariateSupport, mem_biUnion, Polynomial.mem_support_iff, mem_image] at hq
    obtain ⟨j, hj, i, hi, hij⟩ := hq
    subst q
    simpa only [Prod.fst, Prod.snd] using hi
  · intro hq
    have hj : q.2 ∈ P.support := by
      rw [Polynomial.mem_support_iff]
      intro hz
      rw [hz] at hq
      simp at hq
    rw [bivariateSupport, mem_biUnion]
    refine ⟨q.2, hj, ?_⟩
    rw [mem_image]
    exact ⟨q.1, Polynomial.mem_support_iff.mpr hq, by simp⟩

lemma bivariateSupport_nonempty {P : R[X][Y]} :
    (bivariateSupport P).Nonempty ↔ P ≠ 0 := by
  constructor
  · rintro ⟨q, hq⟩ rfl
    simpa using (mem_bivariateSupport.mp hq)
  · intro hP
    obtain ⟨j, hj⟩ := Polynomial.support_nonempty.mpr hP
    obtain ⟨i, hi⟩ := Polynomial.support_nonempty.mpr (Polynomial.mem_support_iff.mp hj)
    exact ⟨(i, j), mem_bivariateSupport.mpr (Polynomial.mem_support_iff.mp hi)⟩

lemma evalEval_eq_bivariate_sum (P : R[X][Y]) (x y : R) :
    P.evalEval x y =
      ∑ q ∈ bivariateSupport P,
        (P.coeff q.2).coeff q.1 * (x ^ q.1 * y ^ q.2) := by
  rw [← Polynomial.eval₂_evalRingHom x, Polynomial.eval₂_eq_sum]
  simp only [Polynomial.sum]
  rw [bivariateSupport, sum_biUnion]
  · apply sum_congr rfl
    intro j hj
    rw [sum_image]
    · change Polynomial.eval x (P.coeff j) * y ^ j = _
      rw [Polynomial.eval_eq_sum]
      simp only [Polynomial.sum]
      rw [Finset.sum_mul]
      apply sum_congr rfl
      intro i hi
      change (P.coeff j).coeff i * x ^ i * y ^ j =
        (P.coeff j).coeff i * (x ^ i * y ^ j)
      ac_rfl
    · intro i₁ hi₁ i₂ hi₂ h
      simpa using congrArg Prod.fst h
  · intro j₁ hj₁ j₂ hj₂ hj
    change Disjoint
      ((P.coeff j₁).support.image fun i ↦ (i, j₁))
      ((P.coeff j₂).support.image fun i ↦ (i, j₂))
    rw [Finset.disjoint_left]
    intro q hq₁ hq₂
    simp only [mem_image] at hq₁ hq₂
    obtain ⟨i₁, hi₁, rfl⟩ := hq₁
    obtain ⟨i₂, hi₂, hpair⟩ := hq₂
    exact hj (congrArg Prod.snd hpair).symm

end BivariatePolynomial

/-- Equation (2.3), stated directly for a nonzero bivariate real polynomial. -/
lemma bivariate_eventually_ne_zero (hPNT : PrimeNumberTheorem)
    (P : ℝ[X][Y]) (hP : P ≠ 0) :
    ∀ᶠ n : ℕ in atTop, P.evalEval (n : ℝ) (nthPrime n : ℝ) ≠ 0 := by
  have hsum := finite_bivariate_sum_eventually_ne_zero hPNT
    (bivariateSupport P)
    (fun q ↦ (P.coeff q.2).coeff q.1)
    (bivariateSupport_nonempty.mpr hP)
    (fun q hq ↦ mem_bivariateSupport.mp hq)
  filter_upwards [hsum] with n hn
  rw [evalEval_eq_bivariate_sum]
  exact hn

end PrimeAsymptotics

section LinearAlgebra

/-- Clear a finite family of rational denominators without losing
nontriviality. -/
lemma clear_denominators {ι : Type*} [Fintype ι] [DecidableEq ι]
    (c : ι → ℚ) (hc : c ≠ 0) :
    ∃ z : ι → ℤ, z ≠ 0 ∧
      ∃ D : ℕ, 0 < D ∧ ∀ i, (z i : ℚ) = (D : ℚ) * c i := by
  let D : ℕ := ∏ i : ι, (c i).den
  have hDpos : 0 < D := Finset.prod_pos fun i _ ↦ (c i).den_pos
  have hden (i : ι) : (c i).den ∣ D :=
    Finset.dvd_prod_of_mem (fun i ↦ (c i).den) (Finset.mem_univ i)
  choose t ht using hden
  let z : ι → ℤ := fun i ↦ (c i).num * (t i : ℤ)
  have hzcast (i : ι) : (z i : ℚ) = (D : ℚ) * c i := by
    have hnum : ((c i).den : ℚ) * c i = ((c i).num : ℚ) := by
      calc
        ((c i).den : ℚ) * c i =
            ((c i).den : ℚ) * (((c i).num : ℚ) / ((c i).den : ℚ)) := by
              congr 1
              exact (Rat.num_div_den (c i)).symm
        _ = ((c i).num : ℚ) := by
          field_simp
    have htQ : (D : ℚ) = ((c i).den : ℚ) * (t i : ℚ) := by
      exact_mod_cast ht i
    simp only [z, Int.cast_mul, Int.cast_natCast]
    rw [htQ]
    calc
      ((c i).num : ℚ) * (t i : ℚ) = (t i : ℚ) * (c i).num := mul_comm _ _
      _ = (t i : ℚ) * (((c i).den : ℚ) * c i) := by rw [hnum]
      _ = ((c i).den : ℚ) * (t i : ℚ) * c i := by ring
  refine ⟨z, ?_, D, hDpos, hzcast⟩
  intro hz
  apply hc
  funext i
  have := hzcast i
  rw [hz, Pi.zero_apply, Int.cast_zero] at this
  exact (mul_eq_zero.mp this.symm).resolve_left (Nat.cast_ne_zero.mpr hDpos.ne')

/-- The integer solution of the paper's `d(d+1)`-by-`d(d+1)+1`
homogeneous system.  We use the equivalent coefficient system for the shifted
polynomials `f(X+j)`, which makes the later cancellation direct. -/
lemma integer_cancellation_coefficients (f : ℤ[X]) (d : ℕ)
    (_hdeg : f.natDegree = d) (h : Fin (d * (d + 1) + 1) → ℕ) :
    ∃ c : Fin (d * (d + 1) + 1) → ℤ, c ≠ 0 ∧
      ∀ (v : Fin d) (u : Fin (d + 1)),
        ∑ j, c j * (-(h j : ℤ)) ^ (v : ℕ) *
          (f.comp (Polynomial.X + Polynomial.C (j : ℤ))).coeff u = 0 := by
  let A : Matrix (Fin d × Fin (d + 1)) (Fin (d * (d + 1) + 1)) ℚ :=
    fun vu j ↦
      (((-(h j : ℤ)) ^ (vu.1 : ℕ) *
        (f.comp (Polynomial.X + Polynomial.C (j : ℤ))).coeff vu.2 : ℤ) : ℚ)
  let L := Matrix.mulVecLin A
  have hdim :
      Module.finrank ℚ (Fin d × Fin (d + 1) → ℚ) <
        Module.finrank ℚ (Fin (d * (d + 1) + 1) → ℚ) := by
    simp only [Module.finrank_pi, Fintype.card_prod, Fintype.card_fin]
    omega
  have hker : LinearMap.ker L ≠ ⊥ :=
    LinearMap.ker_ne_bot_of_finrank_lt hdim
  obtain ⟨cQ, hcQker, hcQne⟩ := (Submodule.ne_bot_iff _).mp hker
  have hmul : A.mulVec cQ = 0 := by
    simpa only [L, Matrix.mulVecLin_apply] using LinearMap.mem_ker.mp hcQker
  obtain ⟨c, hcne, D, hDpos, hcast⟩ := clear_denominators cQ hcQne
  refine ⟨c, hcne, ?_⟩
  intro v u
  apply (Int.cast_eq_zero (α := ℚ)).mp
  simp only [Int.cast_sum, Int.cast_mul, Int.cast_pow, Int.cast_neg,
    Int.cast_natCast]
  calc
    ∑ j, (c j : ℚ) * (-(h j : ℚ)) ^ (v : ℕ) *
        ((f.comp (Polynomial.X + Polynomial.C (j : ℤ))).coeff u : ℚ) =
      (D : ℚ) * ∑ j,
        cQ j * (-(h j : ℚ)) ^ (v : ℕ) *
          ((f.comp (Polynomial.X + Polynomial.C (j : ℤ))).coeff u : ℚ) := by
            simp_rw [hcast]
            rw [Finset.mul_sum]
            apply sum_congr rfl
            intro j hj
            ring
    _ = (D : ℚ) * (A.mulVec cQ (v, u)) := by
      congr 1
      simp only [Matrix.mulVec, dotProduct, A, Int.cast_mul, Int.cast_pow,
        Int.cast_neg, Int.cast_natCast]
      apply sum_congr rfl
      intro j hj
      ring
    _ = 0 := by rw [hmul]; simp

/-- The coefficient conditions imply the finite cancellation identity used
after expanding each reciprocal as a geometric series. -/
lemma cancellation_eval (f : ℤ[X]) (d : ℕ) (hdeg : f.natDegree = d)
    (h : Fin (d * (d + 1) + 1) → ℕ)
    (c : Fin (d * (d + 1) + 1) → ℤ)
    (hc : ∀ (v : Fin d) (u : Fin (d + 1)),
      ∑ j, c j * (-(h j : ℤ)) ^ (v : ℕ) *
        (f.comp (Polynomial.X + Polynomial.C (j : ℤ))).coeff u = 0)
    (v : Fin d) (n : ℤ) :
    ∑ j, c j * f.eval (n + (j : ℤ)) * (-(h j : ℤ)) ^ (v : ℕ) = 0 := by
  let Q : ℤ[X] := ∑ j,
    (c j * (-(h j : ℤ)) ^ (v : ℕ)) •
      f.comp (Polynomial.X + Polynomial.C (j : ℤ))
  have hQdeg : Q.natDegree ≤ d := by
    apply Polynomial.natDegree_sum_le_of_forall_le
    intro j hj
    apply (Polynomial.natDegree_smul_le _ _).trans
    rw [Polynomial.natDegree_comp]
    rw [hdeg, Polynomial.natDegree_X_add_C]
    simp
  have hQ : Q = 0 := by
    apply Polynomial.ext
    intro u
    change Q.coeff u = 0
    by_cases hu : u < d + 1
    · let uf : Fin (d + 1) := ⟨u, hu⟩
      dsimp [Q]
      change Polynomial.lcoeff ℤ u
        (∑ j, (c j * (-(h j : ℤ)) ^ (v : ℕ)) •
          f.comp (Polynomial.X + Polynomial.C (j : ℤ))) = 0
      rw [map_sum]
      simpa only [Polynomial.lcoeff_apply, Polynomial.coeff_smul, smul_eq_mul]
        using hc v uf
    · have hdu : d < u := by omega
      exact (Polynomial.coeff_eq_zero_of_natDegree_lt
        (lt_of_le_of_lt hQdeg hdu))
  have heval : Q.eval n = 0 := by rw [hQ, Polynomial.eval_zero]
  calc
    ∑ j, c j * f.eval (n + (j : ℤ)) * (-(h j : ℤ)) ^ (v : ℕ) =
        Q.eval n := by
      change _ = Polynomial.eval n
        (∑ j, (c j * (-(h j : ℤ)) ^ (v : ℕ)) •
          f.comp (Polynomial.X + Polynomial.C (j : ℤ)))
      rw [Polynomial.eval_finsetSum]
      apply sum_congr rfl
      intro j hj
      simp only [Polynomial.eval_smul, smul_eq_mul, Polynomial.eval_comp,
        Polynomial.eval_add, Polynomial.eval_X, Polynomial.eval_C]
      ring
    _ = 0 := heval

end LinearAlgebra

section AuxiliaryPolynomial

/-- The paper's auxiliary polynomial
`Φ(x,y)=Σ c_j f(x+j) ∏_{i≠j}(y+h_i)`. -/
noncomputable def auxiliaryPhi {m : ℕ} (f : ℤ[X])
    (h : Fin (m + 1) → ℕ) (c : Fin (m + 1) → ℤ) : ℤ[X][Y] :=
  ∑ j,
    Polynomial.C
        (Polynomial.C (c j) *
          f.comp (Polynomial.X + Polynomial.C (j : ℤ))) *
      ∏ i ∈ Finset.univ.erase j,
        (Polynomial.X + Polynomial.C (Polynomial.C (h i : ℤ)))

lemma auxiliaryPhi_eval {m : ℕ} (f : ℤ[X])
    (h : Fin (m + 1) → ℕ) (c : Fin (m + 1) → ℤ) (x y : ℤ) :
    (auxiliaryPhi f h c).evalEval x y =
      ∑ j, c j * f.eval (x + (j : ℤ)) *
        ∏ i ∈ Finset.univ.erase j, (y + (h i : ℤ)) := by
  rw [auxiliaryPhi, Polynomial.evalEval_finsetSum]
  apply sum_congr rfl
  intro j hj
  simp only [Polynomial.evalEval_mul, Polynomial.evalEval_C,
    Polynomial.eval_mul, Polynomial.eval_C, Polynomial.eval_comp,
    Polynomial.eval_add, Polynomial.eval_X, Polynomial.evalEval_prod,
    Polynomial.evalEval_add, Polynomial.evalEval_X]

lemma auxiliaryPhi_ne_zero {m : ℕ} (f : ℤ[X]) (hf : f ≠ 0)
    (h : Fin (m + 1) → ℕ) (hh : StrictMono h)
    (c : Fin (m + 1) → ℤ) (hc : c ≠ 0) :
    auxiliaryPhi f h c ≠ 0 := by
  obtain ⟨j, hcj⟩ : ∃ j, c j ≠ 0 := by
    simpa only [ne_eq, funext_iff, not_forall, not_not, Pi.zero_apply] using hc
  intro hzero
  have heval :
      Polynomial.eval (Polynomial.C (-(h j : ℤ))) (auxiliaryPhi f h c) = 0 := by
    rw [hzero, Polynomial.eval_zero]
  rw [auxiliaryPhi, Polynomial.eval_finsetSum] at heval
  have hsingle :
      ∑ l,
          Polynomial.eval (Polynomial.C (-(h j : ℤ)))
            (Polynomial.C
                (Polynomial.C (c l) *
                  f.comp (Polynomial.X + Polynomial.C (l : ℤ))) *
              ∏ i ∈ Finset.univ.erase l,
                (Polynomial.X + Polynomial.C (Polynomial.C (h i : ℤ)))) =
        Polynomial.eval (Polynomial.C (-(h j : ℤ)))
          (Polynomial.C
              (Polynomial.C (c j) *
                f.comp (Polynomial.X + Polynomial.C (j : ℤ))) *
            ∏ i ∈ Finset.univ.erase j,
              (Polynomial.X + Polynomial.C (Polynomial.C (h i : ℤ)))) := by
    apply Finset.sum_eq_single j
    · intro l hl hlj
      rw [Polynomial.eval_mul, Polynomial.eval_prod]
      apply mul_eq_zero_of_right
      apply Finset.prod_eq_zero (i := j)
      · exact Finset.mem_erase.mpr ⟨Ne.symm hlj, Finset.mem_univ _⟩
      · simp
    · simp
  rw [hsingle] at heval
  have hshift :
      f.comp (Polynomial.X + Polynomial.C (j : ℤ)) ≠ 0 :=
    Polynomial.comp_X_add_C_ne_zero_iff.mpr hf
  have hprod :
      Polynomial.eval (Polynomial.C (-(h j : ℤ)))
          (∏ i ∈ Finset.univ.erase j,
            (Polynomial.X + Polynomial.C (Polynomial.C (h i : ℤ)))) ≠ 0 := by
    rw [Polynomial.eval_prod]
    apply Finset.prod_ne_zero_iff.mpr
    intro i hi
    simp only [Finset.mem_erase] at hi
    simp only [Polynomial.eval_add, Polynomial.eval_X, Polynomial.eval_C,
      ne_eq]
    have hne : h i ≠ h j := fun hij ↦ hi.1 (hh.injective hij)
    intro hz
    have hz0 := congrArg (fun p : ℤ[X] ↦ p.coeff 0) hz
    simp only [Polynomial.coeff_add, Polynomial.coeff_C_zero,
      Polynomial.coeff_zero] at hz0
    apply hne
    omega
  have hleft :
      Polynomial.eval (Polynomial.C (-(h j : ℤ)))
        (Polynomial.C
          (Polynomial.C (c j) *
            f.comp (Polynomial.X + Polynomial.C (j : ℤ)))) ≠ 0 := by
    simp only [Polynomial.eval_C]
    exact mul_ne_zero (Polynomial.C_ne_zero.mpr hcj) hshift
  rw [Polynomial.eval_mul] at heval
  exact (mul_ne_zero hleft hprod) heval

end AuxiliaryPolynomial

section SmallInteger

lemma finite_geometric_reciprocal (d : ℕ) {p h : ℚ}
    (hp : p ≠ 0) (hph : p + h ≠ 0) :
    1 / (p + h) =
      ∑ v ∈ Finset.range d, (-h) ^ v / p ^ (v + 1) +
        (-h) ^ d / (p ^ d * (p + h)) := by
  have hsum :
      (∑ v ∈ Finset.range d, (-h) ^ v / p ^ (v + 1)) =
        (1 / p) * ∑ v ∈ Finset.range d, (-h / p) ^ v := by
    rw [Finset.mul_sum]
    apply Finset.sum_congr rfl
    intro v hv
    rw [div_pow, pow_succ]
    field_simp
  rw [hsum]
  have hrem :
      (-h) ^ d / (p ^ d * (p + h)) = (-h / p) ^ d / (p + h) := by
    rw [div_pow]
    field_simp
  rw [hrem]
  have hg := geom_sum_mul (-h / p) d
  have hg' :
      (∑ v ∈ Finset.range d, (-h / p) ^ v) * (1 - (-h / p)) +
          (-h / p) ^ d = 1 := by
    linarith
  have hr : 1 - (-h / p) = (p + h) / p := by
    field_simp
    ring
  rw [hr] at hg'
  field_simp [hp] at hg'
  field_simp [hp, hph]
  nlinarith

/-- After the `d` cancellation identities, only the remainder term in the
finite geometric expansion survives. -/
lemma reciprocal_sum_eq_remainder (f : ℤ[X]) (d : ℕ)
    (h : Fin (d * (d + 1) + 1) → ℕ)
    (c : Fin (d * (d + 1) + 1) → ℤ)
    (hcancel : ∀ (v : Fin d) (n : ℤ),
      ∑ j, c j * f.eval (n + (j : ℤ)) *
        (-(h j : ℤ)) ^ (v : ℕ) = 0)
    (n : ℕ) :
    ∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) /
          ((nthPrime n : ℚ) + h j) =
      ∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) *
        (-(h j : ℚ)) ^ d /
          ((nthPrime n : ℚ) ^ d * ((nthPrime n : ℚ) + h j)) := by
  have hp : (nthPrime n : ℚ) ≠ 0 :=
    Nat.cast_ne_zero.mpr (nthPrime_prime n).ne_zero
  have hph (j : Fin (d * (d + 1) + 1)) :
      (nthPrime n : ℚ) + h j ≠ 0 := by positivity
  calc
    ∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) /
          ((nthPrime n : ℚ) + h j) =
        ∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) *
          (∑ v ∈ Finset.range d,
              (-(h j : ℚ)) ^ v / (nthPrime n : ℚ) ^ (v + 1) +
            (-(h j : ℚ)) ^ d /
              ((nthPrime n : ℚ) ^ d * ((nthPrime n : ℚ) + h j))) := by
      apply sum_congr rfl
      intro j hj
      rw [div_eq_mul_inv]
      rw [show ((nthPrime n : ℚ) + h j)⁻¹ =
          ∑ v ∈ Finset.range d,
              (-(h j : ℚ)) ^ v / (nthPrime n : ℚ) ^ (v + 1) +
            (-(h j : ℚ)) ^ d /
              ((nthPrime n : ℚ) ^ d * ((nthPrime n : ℚ) + h j)) by
        simpa only [one_div] using
          finite_geometric_reciprocal d hp (hph j)]
    _ =
        (∑ v ∈ Finset.range d,
          (∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) *
            (-(h j : ℚ)) ^ v) /
              (nthPrime n : ℚ) ^ (v + 1)) +
        ∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) *
          (-(h j : ℚ)) ^ d /
            ((nthPrime n : ℚ) ^ d * ((nthPrime n : ℚ) + h j)) := by
      simp_rw [mul_add]
      rw [Finset.sum_add_distrib]
      apply congrArg₂ (· + ·)
      · calc
        ∑ x, (c x : ℚ) * ((f.eval ((n : ℤ) + (x : ℤ)) : ℤ) : ℚ) *
            ∑ v ∈ Finset.range d,
              (-(h x : ℚ)) ^ v / (nthPrime n : ℚ) ^ (v + 1) =
          ∑ x, ∑ v ∈ Finset.range d,
            ((c x : ℚ) * ((f.eval ((n : ℤ) + (x : ℤ)) : ℤ) : ℚ)) *
              ((-(h x : ℚ)) ^ v / (nthPrime n : ℚ) ^ (v + 1)) := by
                apply sum_congr rfl
                intro x hx
                rw [Finset.mul_sum]
        _ = ∑ v ∈ Finset.range d, ∑ x,
            ((c x : ℚ) * ((f.eval ((n : ℤ) + (x : ℤ)) : ℤ) : ℚ)) *
              ((-(h x : ℚ)) ^ v / (nthPrime n : ℚ) ^ (v + 1)) := by
                rw [Finset.sum_comm]
          _ = _ := by
            apply sum_congr rfl
            intro v hv
            rw [Finset.sum_div]
            apply sum_congr rfl
            intro l hl
            ring
      · apply sum_congr rfl
        intro j hj
        ring
    _ = _ := by
      have hzero :
          ∑ v ∈ Finset.range d,
            (∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) *
              (-(h j : ℚ)) ^ v) /
                (nthPrime n : ℚ) ^ (v + 1) = 0 := by
        apply Finset.sum_eq_zero
        intro v hv
        have hvd : v < d := Finset.mem_range.mp hv
        let vf : Fin d := ⟨v, hvd⟩
        have hz := hcancel vf (n : ℤ)
        have hzQ :
            ∑ j, (c j : ℚ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℚ) *
              (-(h j : ℚ)) ^ v = 0 := by
          exact_mod_cast hz
        rw [hzQ, zero_div]
      rw [hzero, zero_add]

lemma tendsto_nthPrime_real :
    Tendsto (fun n : ℕ ↦ (nthPrime n : ℝ)) atTop atTop :=
  tendsto_natCast_atTop_iff.mpr nthPrime_strictMono.tendsto_atTop

lemma shifted_polynomial_eval_isBigO (f : ℤ[X]) (d : ℕ)
    (hdeg : f.natDegree = d) (j : ℕ) :
    (fun n : ℕ ↦ ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ)) =O[atTop]
      (fun n : ℕ ↦ (n : ℝ) ^ d) := by
  let P : ℝ[X] :=
    (f.map (Int.castRingHom ℝ)).comp
      (Polynomial.X + Polynomial.C (j : ℝ))
  have hPdeg : P.natDegree = d := by
    dsimp [P]
    rw [Polynomial.natDegree_comp,
      Polynomial.natDegree_map_eq_of_injective Int.cast_injective,
      hdeg, Polynomial.natDegree_X_add_C]
    simp
  have hOreal :
      (fun x : ℝ ↦ P.eval x) =O[atTop] (fun x : ℝ ↦ x ^ P.natDegree) :=
    P.isEquivalent_atTop_lead.isBigO.trans
      (isBigO_const_mul_self P.leadingCoeff (fun x : ℝ ↦ x ^ P.natDegree) atTop)
  have hOnat := hOreal.comp_tendsto (tendsto_natCast_atTop_atTop (R := ℝ))
  refine hOnat.congr' (Eventually.of_forall fun n ↦ ?_)
    (Eventually.of_forall fun n ↦ by simp [Function.comp_apply, hPdeg])
  simp only [Function.comp_apply, P, Polynomial.eval_comp, Polynomial.eval_add,
    Polynomial.eval_X, Polynomial.eval_C, Polynomial.eval_map]
  rw [show (n : ℝ) + (j : ℝ) = (((n : ℤ) + (j : ℤ) : ℤ) : ℝ) by norm_cast]
  exact Polynomial.eval₂_at_apply _ _

lemma prime_add_isEquivalent (h : ℕ) :
    (fun n : ℕ ↦ (nthPrime n : ℝ) + h) ~[atTop]
      (fun n : ℕ ↦ (nthPrime n : ℝ)) := by
  have hconst :
      (fun _ : ℕ ↦ (h : ℝ)) =o[atTop] (fun n : ℕ ↦ (nthPrime n : ℝ)) :=
    (isLittleO_const_id_atTop (h : ℝ)).comp_tendsto tendsto_nthPrime_real
  exact (IsEquivalent.refl :
    (fun n : ℕ ↦ (nthPrime n : ℝ)) ~[atTop]
      (fun n : ℕ ↦ (nthPrime n : ℝ))).add_isLittleO hconst

lemma residual_term_tendsto_zero (hPNT : PrimeNumberTheorem)
    (f : ℤ[X]) (d : ℕ) (hdeg : f.natDegree = d)
    (c : ℤ) (h : ℕ) :
    Tendsto
      (fun n : ℕ ↦
        (c : ℝ) * ((f.eval ((n : ℤ)) : ℤ) : ℝ) * (-(h : ℝ)) ^ d /
          ((nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h)))
      atTop (𝓝 0) := by
  have hnumO :
      (fun n : ℕ ↦
        (c : ℝ) * ((f.eval (n : ℤ) : ℤ) : ℝ) * (-(h : ℝ)) ^ d) =O[atTop]
        (fun n : ℕ ↦ (n : ℝ) ^ d) := by
    have hfO := shifted_polynomial_eval_isBigO f d hdeg 0
    have hfO' :
        (fun n : ℕ ↦ ((f.eval (n : ℤ) : ℤ) : ℝ)) =O[atTop]
          (fun n : ℕ ↦ (n : ℝ) ^ d) := by
      refine hfO.congr' (Eventually.of_forall fun n ↦ ?_)
        (Eventually.of_forall fun _ ↦ rfl)
      norm_num
    have := hfO'.const_mul_left ((c : ℝ) * (-(h : ℝ)) ^ d)
    exact this.congr' (Eventually.of_forall fun n ↦ by ring)
      (Eventually.of_forall fun _ ↦ rfl)
  have hpowSmall :
      (fun n : ℕ ↦ (n : ℝ) ^ d) =o[atTop]
        (fun n : ℕ ↦ (nthPrime n : ℝ) ^ (d + 1)) := by
    have hm := primeMonomial_isLittleO hPNT
      (i := d) (j := 0) (I := 0) (J := d + 1) (by omega)
    exact hm.congr'
      (Eventually.of_forall fun n ↦ by simp)
      (Eventually.of_forall fun n ↦ by simp)
  have hden :
      (fun n : ℕ ↦
        (nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h)) ~[atTop]
        (fun n : ℕ ↦ (nthPrime n : ℝ) ^ (d + 1)) := by
    have hpRefl :
        (fun n : ℕ ↦ (nthPrime n : ℝ) ^ d) ~[atTop]
          (fun n : ℕ ↦ (nthPrime n : ℝ) ^ d) := IsEquivalent.refl
    have hm := hpRefl.mul (prime_add_isEquivalent h)
    refine hm.congr_right (Eventually.of_forall fun n ↦ ?_)
    simp only [Pi.mul_apply]
    rw [pow_succ]
  have hsmall :
      (fun n : ℕ ↦
        (c : ℝ) * ((f.eval (n : ℤ) : ℤ) : ℝ) * (-(h : ℝ)) ^ d) =o[atTop]
        (fun n : ℕ ↦
          (nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h)) :=
    (hnumO.trans_isLittleO hpowSmall).trans_isBigO hden.isBigO_symm
  apply (isLittleO_iff_tendsto' ?_).mp hsmall
  filter_upwards with n
  intro hz
  exfalso
  have hp : (0 : ℝ) < nthPrime n := by
    exact_mod_cast (nthPrime_prime n).pos
  have hph : (0 : ℝ) < nthPrime n + h :=
    add_pos_of_pos_of_nonneg hp (Nat.cast_nonneg _)
  exact (mul_pos (pow_pos hp _) hph).ne' hz

/-- The complete remainder sum in (2.4) tends to zero. -/
lemma remainder_sum_tendsto_zero (hPNT : PrimeNumberTheorem)
    (f : ℤ[X]) (d : ℕ) (hdeg : f.natDegree = d)
    (h : Fin (d * (d + 1) + 1) → ℕ)
    (c : Fin (d * (d + 1) + 1) → ℤ) :
    Tendsto
      (fun n : ℕ ↦
        ∑ j, (c j : ℝ) *
          ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) *
          (-(h j : ℝ)) ^ d /
            ((nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h j)))
      atTop (𝓝 0) := by
  have hj (j : Fin (d * (d + 1) + 1)) :
      Tendsto
        (fun n : ℕ ↦
          (c j : ℝ) * ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) *
            (-(h j : ℝ)) ^ d /
              ((nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h j)))
        atTop (𝓝 0) := by
    let fj := f.comp (Polynomial.X + Polynomial.C (j : ℤ))
    have hfjdeg : fj.natDegree = d := by
      dsimp [fj]
      rw [Polynomial.natDegree_comp, hdeg,
        Polynomial.natDegree_X_add_C]
      simp
    have ht := residual_term_tendsto_zero hPNT fj d hfjdeg (c j) (h j)
    refine ht.congr' (Eventually.of_forall fun n ↦ ?_)
    simp only [fj, Polynomial.eval_comp, Polynomial.eval_add,
      Polynomial.eval_X, Polynomial.eval_C]
  simpa using tendsto_finsetSum Finset.univ
    (fun j _ ↦ hj j)

end SmallInteger

section MainTheorem

/-- The zero-indexed form of the paper's central theorem. -/
theorem infinitely_many_nthPrime_not_dvd (hPNT : PrimeNumberTheorem)
    (hMT : MaynardTaoBoundedGaps) (f : ℤ[X]) (hf : f ≠ 0) :
    {n : ℕ | ¬(nthPrime n : ℤ) ∣ f.eval (n : ℤ)}.Infinite := by
  by_contra hfinite
  have hbad : {n : ℕ | ¬(nthPrime n : ℤ) ∣ f.eval (n : ℤ)}.Finite :=
    Set.not_infinite.mp hfinite
  have hdiv : ∀ᶠ n : ℕ in atTop, (nthPrime n : ℤ) ∣ f.eval (n : ℤ) := by
    apply atTop_le_cofinite
    have heq :
        {n : ℕ | ¬(nthPrime n : ℤ) ∣ f.eval (n : ℤ)}ᶜ =
          {n : ℕ | (nthPrime n : ℤ) ∣ f.eval (n : ℤ)} := by
      ext n
      simp
    rw [← heq]
    exact hbad.compl_mem_cofinite
  let d := f.natDegree
  let k := d * (d + 1)
  obtain ⟨h, M, hM, hh0, hhmono, hcluster⟩ :=
    fixed_prime_cluster hMT k
  have hk : k = d * (d + 1) := rfl
  subst k
  obtain ⟨c, hcne, hcoeff⟩ :=
    integer_cancellation_coefficients f d rfl h
  have hcancel : ∀ (v : Fin d) (n : ℤ),
      ∑ j, c j * f.eval (n + (j : ℤ)) *
        (-(h j : ℤ)) ^ (v : ℕ) = 0 :=
    fun v n ↦ cancellation_eval f d rfl h c hcoeff v n
  let ΦZ : ℤ[X][Y] := auxiliaryPhi f h c
  have hΦZ : ΦZ ≠ 0 := auxiliaryPhi_ne_zero f hf h hhmono c hcne
  let castZH : ℤ →+* ℝ := Int.castRingHom ℝ
  let ΦR : ℝ[X][Y] := ΦZ.map (Polynomial.mapRingHom castZH)
  have hinner :
      Function.Injective (Polynomial.mapRingHom castZH) :=
    Polynomial.map_injective castZH Int.cast_injective
  have hΦR : ΦR ≠ 0 :=
    (Polynomial.map_ne_zero_iff hinner).mpr hΦZ
  have hΦevent :
      ∀ᶠ n : ℕ in atTop,
        ΦR.evalEval (n : ℝ) (nthPrime n : ℝ) ≠ 0 :=
    bivariate_eventually_ne_zero hPNT ΦR hΦR
  have hrem :=
    remainder_sum_tendsto_zero hPNT f d rfl h c
  have hsmall :
      ∀ᶠ n : ℕ in atTop,
        ‖∑ j, (c j : ℝ) *
          ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) *
          (-(h j : ℝ)) ^ d /
            ((nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h j))‖ < 1 := by
    have hm := (Metric.tendsto_nhds.mp hrem) 1 zero_lt_one
    simpa only [dist_zero_right] using hm
  obtain ⟨Ndiv, hNdiv⟩ := eventually_atTop.mp hdiv
  have hdivAll :
      ∀ᶠ n : ℕ in atTop, ∀ j : Fin (d * (d + 1) + 1),
        (nthPrime (n + j) : ℤ) ∣ f.eval ((n : ℤ) + (j : ℤ)) := by
    apply eventually_atTop.mpr
    refine ⟨Ndiv, fun n hn j ↦ ?_⟩
    simpa only [Nat.cast_add] using
      hNdiv (n + j) (hn.trans (Nat.le_add_right n j))
  have hall :
      ∀ᶠ n : ℕ in atTop,
        ΦR.evalEval (n : ℝ) (nthPrime n : ℝ) ≠ 0 ∧
        ‖∑ j, (c j : ℝ) *
          ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) *
          (-(h j : ℝ)) ^ d /
            ((nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h j))‖ < 1 ∧
        ∀ j : Fin (d * (d + 1) + 1),
          (nthPrime (n + j) : ℤ) ∣ f.eval ((n : ℤ) + (j : ℤ)) := by
    filter_upwards [hΦevent, hsmall, hdivAll] with n hnΦ hnsmall hndiv
    exact ⟨hnΦ, hnsmall, hndiv⟩
  obtain ⟨N, hN⟩ := eventually_atTop.mp hall
  obtain ⟨n, hnM, hnN⟩ := hM.exists_gt N
  obtain ⟨hnΦ, hnsmall, hndiv⟩ := hN n hnN.le
  have hncluster (j : Fin (d * (d + 1) + 1)) :
      nthPrime (n + j) = nthPrime n + h j :=
    hcluster n hnM j
  have hdvd (j : Fin (d * (d + 1) + 1)) :
      ((nthPrime n + h j : ℕ) : ℤ) ∣
        f.eval ((n : ℤ) + (j : ℤ)) := by
    rw [← hncluster j]
    exact hndiv j
  choose q hq using hdvd
  let qsum : ℤ := ∑ j, c j * q j
  have hrecip := reciprocal_sum_eq_remainder f d h c hcancel n
  have hrecipR :
      ∑ j, (c j : ℝ) *
          ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) /
            ((nthPrime n : ℝ) + h j) =
        ∑ j, (c j : ℝ) *
          ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) *
          (-(h j : ℝ)) ^ d /
            ((nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h j)) := by
    exact_mod_cast hrecip
  have hqsumR :
      (qsum : ℝ) =
        ∑ j, (c j : ℝ) *
          ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) *
          (-(h j : ℝ)) ^ d /
            ((nthPrime n : ℝ) ^ d * ((nthPrime n : ℝ) + h j)) := by
    calc
      (qsum : ℝ) = ∑ j, (c j : ℝ) * (q j : ℝ) := by
        simp only [qsum, Int.cast_sum, Int.cast_mul]
      _ = ∑ j, (c j : ℝ) *
          ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) /
            ((nthPrime n : ℝ) + h j) := by
        apply sum_congr rfl
        intro j hj
        have hqR :
            ((f.eval ((n : ℤ) + (j : ℤ)) : ℤ) : ℝ) =
              ((nthPrime n : ℝ) + h j) * (q j : ℝ) := by
          exact_mod_cast hq j
        rw [hqR]
        have hpR : (0 : ℝ) < nthPrime n := by
          exact_mod_cast (nthPrime_prime n).pos
        have hden :
            (nthPrime n : ℝ) + h j ≠ 0 :=
          (add_pos_of_pos_of_nonneg hpR (Nat.cast_nonneg _)).ne'
        field_simp [hden]
      _ = _ := hrecipR
  have hqsumSmall : ‖(qsum : ℝ)‖ < 1 := by
    rw [hqsumR]
    exact hnsmall
  have hqsum0 : qsum = 0 := by
    apply Int.abs_lt_one_iff.mp
    rw [← Int.cast_lt (R := ℝ), Int.cast_one, Int.cast_abs]
    simpa only [Real.norm_eq_abs] using hqsumSmall
  have hΦzero : ΦZ.evalEval (n : ℤ) (nthPrime n : ℤ) = 0 := by
    rw [auxiliaryPhi_eval]
    calc
      ∑ j, c j * f.eval ((n : ℤ) + (j : ℤ)) *
          ∏ i ∈ Finset.univ.erase j, ((nthPrime n : ℤ) + h i) =
        ∑ j, c j * q j *
          ∏ i, ((nthPrime n : ℤ) + h i) := by
            apply sum_congr rfl
            intro j hj
            have hq' :
                f.eval ((n : ℤ) + (j : ℤ)) =
                  ((nthPrime n : ℤ) + (h j : ℤ)) * q j := by
              exact_mod_cast hq j
            rw [← Finset.mul_prod_erase Finset.univ
              (fun i ↦ (nthPrime n : ℤ) + h i) (Finset.mem_univ j)]
            rw [hq']
            ring
      _ = (∑ j, c j * q j) *
          ∏ i, ((nthPrime n : ℤ) + h i) := by
            rw [Finset.sum_mul]
      _ = 0 := by
        change qsum * ∏ i, ((nthPrime n : ℤ) + h i) = 0
        rw [hqsum0, zero_mul]
  have hmapEval :=
    Polynomial.map_mapRingHom_evalEval castZH ΦZ (n : ℤ) (nthPrime n : ℤ)
  have hΦcast :
      ΦR.evalEval (n : ℝ) (nthPrime n : ℝ) =
        ((ΦZ.evalEval (n : ℤ) (nthPrime n : ℤ) : ℤ) : ℝ) := by
    exact hmapEval
  rw [hΦcast, hΦzero, Int.cast_zero] at hnΦ
  exact hnΦ rfl

/-- The paper uses one-indexed primes: `paperPrime 1 = 2`.  The value at zero
is irrelevant because the theorem explicitly restricts to positive indices. -/
noncomputable def paperPrime (n : ℕ) : ℕ := nthPrime (n - 1)

@[simp] lemma paperPrime_succ (n : ℕ) : paperPrime (n + 1) = nthPrime n := by
  simp [paperPrime]

/-- **Theorem 1.1 (Seki).**  For every nonzero integer polynomial `f`, there
are infinitely many positive integers `n` for which the `n`-th prime does not
divide `f(n)`. -/
theorem theorem_1_1 (hPNT : PrimeNumberTheorem)
    (hMT : MaynardTaoBoundedGaps) (f : ℤ[X]) (hf : f ≠ 0) :
    {n : ℕ | 0 < n ∧ ¬(paperPrime n : ℤ) ∣ f.eval (n : ℤ)}.Infinite := by
  let g := f.comp (Polynomial.X + Polynomial.C (1 : ℤ))
  have hg : g ≠ 0 := Polynomial.comp_X_add_C_ne_zero_iff.mpr hf
  have hzero := infinitely_many_nthPrime_not_dvd hPNT hMT g hg
  let S : Set ℕ :=
    {n : ℕ | 0 < n ∧ ¬(paperPrime n : ℤ) ∣ f.eval (n : ℤ)}
  let e :
      {m : ℕ | ¬(nthPrime m : ℤ) ∣ g.eval (m : ℤ)} ↪ S :=
    ⟨fun m ↦ ⟨m.1 + 1, by
        constructor
        · omega
        · have hm := m.2
          change ¬(nthPrime m.1 : ℤ) ∣ g.eval (m.1 : ℤ) at hm
          simpa only [paperPrime_succ, g, Polynomial.eval_comp,
            Polynomial.eval_add, Polynomial.eval_X, Polynomial.eval_C,
            Nat.cast_add, Nat.cast_one] using hm⟩,
      fun a b hab ↦ by
        apply Subtype.ext
        have hv := congrArg (fun z : S ↦ z.1) hab
        change a.1 + 1 = b.1 + 1 at hv
        omega⟩
  haveI : Infinite {m : ℕ | ¬(nthPrime m : ℤ) ∣ g.eval (m : ℤ)} :=
    Set.infinite_coe_iff.mpr hzero
  have hSI : Infinite S := Infinite.of_injective e e.injective
  exact Set.infinite_coe_iff.mp hSI

end MainTheorem

end ChampernowneA