Erdős problem 859

/-!
# Erdos 859: counterexample submission

Single-file proof body for conjectures.io's counterexample task
fc-8432eac9-erdos859-erdos-859-0087a1d40d-counterexample-v1.
The trusted submission wrapper supplies the imports and the Bounty namespace.
All helper declarations and the third-party license are included below.
-/

noncomputable section
open scoped Classical


/- BEGIN isolated base development. -/
section

open Filter Asymptotics Finset
open scoped Topology

namespace Erdos859Verification

/-- Every divisor used to represent `t` is at most `t`; hence membership for
positive integers depends only on their residue class modulo `t!`. -/
theorem divisorSumSet_iff_of_modEq_factorial
    {t n m : ℕ} (hn : n ≠ 0) (hm : m ≠ 0)
    (hmod : Nat.ModEq t.factorial n m) :
    n ∈ Erdos859.DivisorSumSet t ↔ m ∈ Erdos859.DivisorSumSet t := by
  have transfer {a b : ℕ} (hb : b ≠ 0)
      (hab : Nat.ModEq t.factorial a b)
      (ha : a ∈ Erdos859.DivisorSumSet t) :
      b ∈ Erdos859.DivisorSumSet t := by
    obtain ⟨s, hs, hsum⟩ := ha
    refine ⟨s, ?_, hsum⟩
    intro i hi
    have hid : i ∈ a.divisors := hs hi
    have hip : 0 < i := Nat.pos_of_mem_divisors hid
    have hit : i ≤ t := by
      rw [hsum]
      exact Finset.single_le_sum (fun j _ ↦ Nat.zero_le j) hi
    exact Nat.mem_divisors.mpr
      ⟨(hab.dvd_iff (Nat.dvd_factorial hip hit)).mp (Nat.mem_divisors.mp hid).1, hb⟩
  exact ⟨transfer hm hmod, transfer hn hmod.symm⟩

/-- Sum a periodic real sequence over a whole number of periods. -/
theorem sum_range_periodic_mul {f : ℕ → ℝ} {m : ℕ}
    (hp : Function.Periodic f m) (q : ℕ) :
    ∑ i ∈ range (q * m), f i = (q : ℝ) * ∑ i ∈ range m, f i := by
  induction q with
  | zero => simp
  | succ q ih =>
    rw [Nat.succ_mul, sum_range_add, ih]
    have heq : (∑ i ∈ range m, f (q * m + i)) = ∑ i ∈ range m, f i := by
      apply sum_congr rfl
      intro i _
      simpa [add_comm] using hp.nat_mul q i
    rw [heq, Nat.cast_add, Nat.cast_one]
    ring

/-- Quotient and remainder split the sum of a periodic sequence. -/
theorem sum_range_periodic_decompose {f : ℕ → ℝ} {m : ℕ}
    (hp : Function.Periodic f m) (n : ℕ) :
    ∑ i ∈ range n, f i =
      (n / m : ℕ) * (∑ i ∈ range m, f i) + ∑ i ∈ range (n % m), f i := by
  conv_lhs => rw [← Nat.div_add_mod n m, Nat.mul_comm m (n / m)]
  rw [sum_range_add, sum_range_periodic_mul hp]
  congr 1
  apply sum_congr rfl
  intro i _
  simpa [add_comm] using hp.nat_mul (n / m) i

/-- The Cesàro mean of a bounded periodic sequence is its mean over one period. -/
theorem tendsto_mean_periodic {f : ℕ → ℝ} {m : ℕ}
    (hm : 0 < m) (hp : Function.Periodic f m)
    (h0 : ∀ i, 0 ≤ f i) (h1 : ∀ i, f i ≤ 1) :
    Tendsto (fun n : ℕ ↦ (∑ i ∈ range n, f i) / n)
      atTop (𝓝 ((∑ i ∈ range m, f i) / m)) := by
  let s := ∑ i ∈ range m, f i
  have hmR : (0 : ℝ) < m := by exact_mod_cast hm
  have hs0 : 0 ≤ s := sum_nonneg (fun i _ ↦ h0 i)
  have hs1 : s ≤ m := by
    calc
      s ≤ ∑ _i ∈ range m, (1 : ℝ) := sum_le_sum (fun i _ ↦ h1 i)
      _ = m := by simp
  have hsdiv0 : 0 ≤ s / m := div_nonneg hs0 hmR.le
  have hsdiv1 : s / m ≤ 1 := (div_le_one hmR).mpr hs1
  have hbound : ∀ n : ℕ,
      |(∑ i ∈ range n, f i) - (n : ℝ) * (s / m)| ≤ m := by
    intro n
    have hr0 : 0 ≤ ∑ i ∈ range (n % m), f i := sum_nonneg (fun i _ ↦ h0 i)
    have hr1 : (∑ i ∈ range (n % m), f i) ≤ m := by
      calc
        _ ≤ ∑ _i ∈ range (n % m), (1 : ℝ) := sum_le_sum (fun i _ ↦ h1 i)
        _ = (n % m : ℕ) := by simp
        _ ≤ (m : ℝ) := by exact_mod_cast (Nat.mod_lt n hm).le
    have hrm0 : (0 : ℝ) ≤ (n % m : ℕ) * (s / m) :=
      mul_nonneg (Nat.cast_nonneg _) hsdiv0
    have hrm1 : (n % m : ℕ) * (s / m) ≤ (m : ℝ) := by
      calc
        _ ≤ (n % m : ℕ) * (1 : ℝ) :=
          mul_le_mul_of_nonneg_left hsdiv1 (Nat.cast_nonneg _)
        _ ≤ (m : ℝ) := by simpa using (show ((n % m : ℕ) : ℝ) ≤ m by
          exact_mod_cast (Nat.mod_lt n hm).le)
    have hn : (n : ℝ) = m * (n / m : ℕ) + (n % m : ℕ) := by
      exact_mod_cast (Nat.div_add_mod n m).symm
    rw [sum_range_periodic_decompose hp, hn]
    change |(n / m : ℕ) * s + (∑ i ∈ range (n % m), f i) -
      (m * (n / m : ℕ) + (n % m : ℕ)) * (s / m)| ≤ m
    have heq : (n / m : ℕ) * s + (∑ i ∈ range (n % m), f i) -
        (m * (n / m : ℕ) + (n % m : ℕ)) * (s / m) =
        (∑ i ∈ range (n % m), f i) - (n % m : ℕ) * (s / m) := by
      field_simp
      ring
    rw [heq, abs_le]
    constructor <;> linarith
  have hzero : Tendsto (fun n : ℕ ↦
      (∑ i ∈ range n, f i) / n - s / m) atTop (𝓝 0) := by
    apply squeeze_zero_norm' _ (tendsto_const_div_atTop_nhds_zero_nat (m : ℝ))
    filter_upwards [eventually_gt_atTop 0] with n hn
    have hnR : (0 : ℝ) < n := by exact_mod_cast hn
    have heq : (∑ i ∈ range n, f i) / n - s / m =
        ((∑ i ∈ range n, f i) - (n : ℝ) * (s / m)) / n := by
      field_simp
    rw [Real.norm_eq_abs, heq, abs_div, abs_of_pos hnR]
    exact div_le_div_of_nonneg_right (hbound n) hnR.le
  exact tendsto_sub_nhds_zero_iff.mp hzero

/-- Express the partial density as an average of a real-valued indicator. -/
theorem partialDensity_eq_indicator_mean (S : Set ℕ) (n : ℕ) :
    S.partialDensity Set.univ n =
      (∑ i ∈ range n, S.indicator (fun _ ↦ (1 : ℝ)) i) / n := by
  classical
  have hset : S ∩ Set.Iio n = (((range n).filter (· ∈ S) : Finset ℕ) : Set ℕ) := by
    ext i
    simp [and_comm]
  simp only [Set.partialDensity, Set.inter_univ, Set.univ_inter, hset,
    Set.ncard_coe_finset, Nat.ncard_Iio]
  congr 1
  simp [Set.indicator]

/-- An eventual equality of indicator functions preserves their natural density. -/
theorem hasDensity_of_indicator_mean {S : Set ℕ} {f : ℕ → ℝ} {a : ℝ}
    (hmean : Tendsto (fun n : ℕ ↦ (∑ i ∈ range n, f i) / n) atTop (𝓝 a))
    (heq : ∀ᶠ n in atTop, S.indicator (fun _ ↦ (1 : ℝ)) n = f n) :
    S.HasDensity a := by
  have hdiff : Tendsto (fun n ↦ S.indicator (fun _ ↦ (1 : ℝ)) n - f n)
      atTop (𝓝 0) := by
    apply tendsto_const_nhds.congr'
    filter_upwards [heq] with n hn
    simp [hn]
  have hcesaro := hdiff.cesaro
  have hcombined := hcesaro.add hmean
  change Tendsto (fun n ↦ S.partialDensity Set.univ n) atTop (𝓝 a)
  convert hcombined using 1
  · funext n
    rw [partialDensity_eq_indicator_mean]
    simp only [sum_sub_distrib, div_eq_mul_inv]
    ring
  · simp

/-- A set agreeing eventually with a bounded periodic indicator has a density. -/
theorem hasDensity_of_periodic_indicator {S : Set ℕ} {f : ℕ → ℝ} {m : ℕ}
    (hm : 0 < m) (hp : Function.Periodic f m)
    (h0 : ∀ i, 0 ≤ f i) (h1 : ∀ i, f i ≤ 1)
    (heq : ∀ᶠ n in atTop, S.indicator (fun _ ↦ (1 : ℝ)) n = f n) :
    S.HasDensity ((∑ i ∈ range m, f i) / m) :=
  hasDensity_of_indicator_mean (tendsto_mean_periodic hm hp h0 h1) heq

/-- A periodic indicator for the divisor-sum event, using positive representatives. -/
noncomputable def divisorIndicator (t n : ℕ) : ℝ :=
  (Erdos859.DivisorSumSet t).indicator (fun _ ↦ (1 : ℝ))
    (n % t.factorial + t.factorial)

theorem divisorIndicator_nonneg (t n : ℕ) : 0 ≤ divisorIndicator t n := by
  classical
  unfold divisorIndicator
  by_cases h : n % t.factorial + t.factorial ∈ Erdos859.DivisorSumSet t
  · simp [h]
  · simp [h]

theorem divisorIndicator_le_one (t n : ℕ) : divisorIndicator t n ≤ 1 := by
  classical
  unfold divisorIndicator
  by_cases h : n % t.factorial + t.factorial ∈ Erdos859.DivisorSumSet t
  · simp [h]
  · simp [h]

theorem divisorIndicator_periodic (t : ℕ) :
    Function.Periodic (divisorIndicator t) t.factorial := by
  intro n
  simp [divisorIndicator]

theorem divisorIndicator_eventually_eq (t : ℕ) :
    ∀ᶠ n in atTop, (Erdos859.DivisorSumSet t).indicator (fun _ ↦ (1 : ℝ)) n =
      divisorIndicator t n := by
  classical
  filter_upwards [eventually_gt_atTop 0] with n hn
  have hrep : n % t.factorial + t.factorial ≠ 0 :=
    ne_of_gt (lt_of_lt_of_le (Nat.factorial_pos t) (Nat.le_add_left _ _))
  have hmod : Nat.ModEq t.factorial n (n % t.factorial + t.factorial) := by
    simp [Nat.ModEq]
  have hiff := divisorSumSet_iff_of_modEq_factorial hn.ne' hrep hmod
  simp only [divisorIndicator, Set.indicator_apply, hiff]

/-- The exact density, as a finite average over residue classes modulo `t!`. -/
noncomputable def divisorDensity (t : ℕ) : ℝ :=
  (∑ i ∈ range t.factorial, divisorIndicator t i) / t.factorial

/-- Unconditional existence of the natural density of each divisor-sum set. -/
theorem divisorSumSet_hasDensity (t : ℕ) :
    (Erdos859.DivisorSumSet t).HasDensity (divisorDensity t) :=
  hasDensity_of_periodic_indicator (Nat.factorial_pos t) (divisorIndicator_periodic t)
    (divisorIndicator_nonneg t) (divisorIndicator_le_one t)
    (divisorIndicator_eventually_eq t)

theorem factorial_mem_divisorSumSet (t : ℕ) :
    t.factorial ∈ Erdos859.DivisorSumSet t := by
  by_cases ht : t = 0
  · subst t
    exact ⟨∅, by simp, by simp⟩
  · refine ⟨{t}, ?_, by simp⟩
    intro i hi
    have hi' : i = t := by simpa using hi
    subst i
    exact Nat.mem_divisors.mpr
      ⟨Nat.dvd_factorial (Nat.pos_of_ne_zero ht) le_rfl, Nat.factorial_ne_zero t⟩

theorem divisorIndicator_zero (t : ℕ) : divisorIndicator t 0 = 1 := by
  simp [divisorIndicator, factorial_mem_divisorSumSet]

/-- Each exact density is strictly positive. -/
theorem divisorDensity_pos (t : ℕ) : 0 < divisorDensity t := by
  have hs : 1 ≤ ∑ i ∈ range t.factorial, divisorIndicator t i := by
    rw [← divisorIndicator_zero t]
    exact single_le_sum (fun i _ ↦ divisorIndicator_nonneg t i)
      (mem_range.mpr (Nat.factorial_pos t))
  exact div_pos (lt_of_lt_of_le zero_lt_one hs) (by exact_mod_cast Nat.factorial_pos t)

theorem divisorDensity_le_one (t : ℕ) : divisorDensity t ≤ 1 := by
  exact (divisorSumSet_hasDensity t).mono (Set.subset_univ _) Set.HasDensity.univ

/-- The positive-density variant, proved without using its imported placeholder. -/
theorem divisorSumSet_hasPosDensity (t : ℕ) :
    (Erdos859.DivisorSumSet t).HasPosDensity :=
  ⟨divisorDensity t, divisorDensity_pos t, divisorSumSet_hasDensity t⟩

/-- A uniformly sampled residue is divisible by `d` with frequency `1/d`. -/
theorem tendsto_mean_dvd_indicator {d : ℕ} (hd : 0 < d) :
    Tendsto (fun N : ℕ ↦ (∑ n ∈ range N, if d ∣ n then (1 : ℝ) else 0) / N)
      atTop (𝓝 (1 / (d : ℝ))) := by
  have hp : Function.Periodic (fun n ↦ if d ∣ n then (1 : ℝ) else 0) d := by
    intro n
    simp only [← Nat.dvd_add_iff_left (dvd_refl d)]
  have h0 : ∀ n : ℕ, (0 : ℝ) ≤ if d ∣ n then 1 else 0 := by
    intro n
    split_ifs <;> norm_num
  have h1 : ∀ n : ℕ, (if d ∣ n then (1 : ℝ) else 0) ≤ 1 := by
    intro n
    split_ifs <;> norm_num
  have hs : (∑ n ∈ range d, if d ∣ n then (1 : ℝ) else 0) = 1 := by
    rw [sum_eq_single 0]
    · simp
    · intro n hn hn0
      have hnd : ¬ d ∣ n := by
        intro hdn
        have hle := Nat.le_of_dvd (Nat.pos_of_ne_zero hn0) hdn
        exact (not_le.mpr (mem_range.mp hn)) hle
      simp [hnd]
    · intro hnot
      exact False.elim (hnot (mem_range.mpr hd))
  simpa only [hs] using tendsto_mean_periodic hd hp h0 h1

/-- Sum the small positive divisors, with the zero sample harmlessly included. -/
noncomputable def smallDivisorSum (y n : ℕ) : ℝ :=
  ∑ d ∈ Icc 1 y, if d ∣ n then (d : ℝ) else 0

theorem smallDivisorSum_nonneg (y n : ℕ) : 0 ≤ smallDivisorSum y n := by
  apply sum_nonneg
  intro d _
  split_ifs <;> positivity

/-- Each small divisor contributes exactly one to the limiting mean. -/
theorem tendsto_mean_smallDivisorSum (y : ℕ) :
    Tendsto (fun N : ℕ ↦ (∑ n ∈ range N, smallDivisorSum y n) / N)
      atTop (𝓝 (y : ℝ)) := by
  have heach : ∀ d ∈ Icc 1 y,
      Tendsto (fun N : ℕ ↦ (∑ n ∈ range N, if d ∣ n then (d : ℝ) else 0) / N)
        atTop (𝓝 (1 : ℝ)) := by
    intro d hd
    have hd0 : 0 < d := (mem_Icc.mp hd).1
    have hdR : (d : ℝ) ≠ 0 := by exact_mod_cast hd0.ne'
    have h := (tendsto_mean_dvd_indicator hd0).const_mul (d : ℝ)
    convert h using 1
    · funext N
      rw [← mul_div_assoc, mul_sum]
      congr 1
      apply sum_congr rfl
      intro n _
      split_ifs <;> simp
    · field_simp
  have hsum := tendsto_finsetSum (Icc 1 y) heach
  convert hsum using 1
  · funext N
    simp only [smallDivisorSum, sum_div]
    rw [sum_comm]
  · simp

/-- Integers with a divisor in the half-open interval `(y,t]`. -/
def DivisorIntervalSet (y t : ℕ) : Set ℕ :=
  {n | ∃ d : ℕ, y < d ∧ d ≤ t ∧ d ∣ n}

theorem divisorInterval_periodic (y t : ℕ) :
    Function.Periodic ((DivisorIntervalSet y t).indicator (fun _ ↦ (1 : ℝ)))
      t.factorial := by
  classical
  intro n
  have heq : n + t.factorial ∈ DivisorIntervalSet y t ↔ n ∈ DivisorIntervalSet y t := by
    constructor
    · rintro ⟨d, hyd, hdt, hd⟩
      exact ⟨d, hyd, hdt,
        (Nat.dvd_add_iff_left (Nat.dvd_factorial (lt_of_le_of_lt (Nat.zero_le y) hyd)
          hdt)).mpr hd⟩
    · rintro ⟨d, hyd, hdt, hd⟩
      exact ⟨d, hyd, hdt,
        (Nat.dvd_add_iff_left (Nat.dvd_factorial (lt_of_le_of_lt (Nat.zero_le y) hyd)
          hdt)).mp hd⟩
  simp only [Set.indicator_apply, heq]

/-- The divisor-interval density as an exact residue-class average. -/
noncomputable def divisorIntervalDensity (y t : ℕ) : ℝ :=
  (∑ n ∈ range t.factorial, (DivisorIntervalSet y t).indicator (fun _ ↦ (1 : ℝ)) n) /
    t.factorial

theorem divisorInterval_hasDensity (y t : ℕ) :
    (DivisorIntervalSet y t).HasDensity (divisorIntervalDensity y t) := by
  classical
  apply hasDensity_of_periodic_indicator (Nat.factorial_pos t) (divisorInterval_periodic y t)
  · intro n
    by_cases hn : n ∈ DivisorIntervalSet y t <;> simp [hn]
  · intro n
    by_cases hn : n ∈ DivisorIntervalSet y t <;> simp [hn]
  · exact Filter.Eventually.of_forall (fun _ ↦ rfl)

/-- A representation avoiding large divisors must use at least `t` in small divisors. -/
theorem le_smallDivisorSum_of_no_large_divisor {t y n : ℕ}
    (hrep : n ∈ Erdos859.DivisorSumSet t) (hno : n ∉ DivisorIntervalSet y t) :
    (t : ℝ) ≤ smallDivisorSum y n := by
  classical
  obtain ⟨s, hs, hsum⟩ := hrep
  have hsub : s ⊆ Icc 1 y := by
    intro d hd
    have hdiv : d ∈ n.divisors := hs hd
    have hdt : d ≤ t := by
      rw [hsum]
      exact single_le_sum (fun i _ ↦ Nat.zero_le i) hd
    have hdy : d ≤ y := by
      by_contra h
      exact hno ⟨d, lt_of_not_ge h, hdt, (Nat.mem_divisors.mp hdiv).1⟩
    exact mem_Icc.mpr ⟨Nat.pos_of_mem_divisors hdiv, hdy⟩
  calc
    (t : ℝ) = ∑ d ∈ s, (d : ℝ) := by rw [hsum, Nat.cast_sum]
    _ = ∑ d ∈ s, if d ∣ n then (d : ℝ) else 0 := by
      apply sum_congr rfl
      intro d hd
      simp [(Nat.mem_divisors.mp (hs hd)).1]
    _ ≤ smallDivisorSum y n := by
      apply sum_le_sum_of_subset_of_nonneg hsub
      intro i _ _
      split_ifs <;> positivity

/-- The deterministic inequality behind the divisor-interval upper bound. -/
theorem divisor_sum_indicator_le (t y n : ℕ) :
    (t : ℝ) * (Erdos859.DivisorSumSet t).indicator (fun _ ↦ (1 : ℝ)) n ≤
      (t : ℝ) * (DivisorIntervalSet y t).indicator (fun _ ↦ (1 : ℝ)) n +
        smallDivisorSum y n := by
  classical
  by_cases hrep : n ∈ Erdos859.DivisorSumSet t
  · by_cases hlarge : n ∈ DivisorIntervalSet y t
    · simp only [Set.indicator_of_mem hrep, Set.indicator_of_mem hlarge, mul_one]
      linarith [smallDivisorSum_nonneg y n]
    · simpa only [Set.indicator_of_mem hrep, Set.indicator_of_notMem hlarge,
        mul_one, mul_zero, zero_add] using le_smallDivisorSum_of_no_large_divisor hrep hlarge
  · simp only [Set.indicator_of_notMem hrep, mul_zero]
    have hi : 0 ≤ (DivisorIntervalSet y t).indicator (fun _ ↦ (1 : ℝ)) n := by
      by_cases hn : n ∈ DivisorIntervalSet y t <;> simp [hn]
    exact add_nonneg (mul_nonneg (Nat.cast_nonneg t) hi) (smallDivisorSum_nonneg y n)

/-- The exact density bound used before invoking a divisor-interval theorem. -/
theorem divisorDensity_le_intervalDensity_add (t y : ℕ) (ht : 0 < t) :
    divisorDensity t ≤ divisorIntervalDensity y t + (y : ℝ) / t := by
  have hleft : Tendsto (fun N : ℕ ↦
      (t : ℝ) * ((Erdos859.DivisorSumSet t).partialDensity Set.univ N))
      atTop (𝓝 ((t : ℝ) * divisorDensity t)) :=
    (divisorSumSet_hasDensity t).const_mul (t : ℝ)
  have hright : Tendsto (fun N : ℕ ↦
      (t : ℝ) * ((DivisorIntervalSet y t).partialDensity Set.univ N) +
        (∑ n ∈ range N, smallDivisorSum y n) / N)
      atTop (𝓝 ((t : ℝ) * divisorIntervalDensity y t + y)) :=
    ((divisorInterval_hasDensity y t).const_mul (t : ℝ)).add
      (tendsto_mean_smallDivisorSum y)
  have hle : (t : ℝ) * divisorDensity t ≤ (t : ℝ) * divisorIntervalDensity y t + y := by
    apply le_of_tendsto_of_tendsto' hleft hright
    intro N
    rw [partialDensity_eq_indicator_mean, partialDensity_eq_indicator_mean,
      ← mul_div_assoc, ← mul_div_assoc, ← add_div]
    apply div_le_div_of_nonneg_right _ (Nat.cast_nonneg N)
    rw [mul_sum, mul_sum, ← sum_add_distrib]
    exact sum_le_sum (fun n _ ↦ divisor_sum_indicator_le t y n)
  have htR : (0 : ℝ) < t := by exact_mod_cast ht
  have hdiv := (div_le_div_iff_of_pos_right htR).mpr hle
  simpa [add_div, htR.ne'] using hdiv

/-- The coprime factors of a divisor product are uniquely determined. -/
theorem divisor_product_injective {b n d₁ d₂ e₁ e₂ : ℕ}
    (hbn : Nat.Coprime b n)
    (hd₁ : d₁ ∈ n.divisors) (hd₂ : d₂ ∈ n.divisors)
    (he₁ : e₁ ∈ b.divisors) (he₂ : e₂ ∈ b.divisors)
    (heq : d₁ * e₁ = d₂ * e₂) : d₁ = d₂ ∧ e₁ = e₂ := by
  have hc₁ : Nat.Coprime d₁ e₂ :=
    (hbn.symm.of_dvd_left (Nat.mem_divisors.mp hd₁).1).of_dvd_right
      (Nat.mem_divisors.mp he₂).1
  have hc₂ : Nat.Coprime d₂ e₁ :=
    (hbn.symm.of_dvd_left (Nat.mem_divisors.mp hd₂).1).of_dvd_right
      (Nat.mem_divisors.mp he₁).1
  have h₁ : d₁ ∣ d₂ := hc₁.dvd_of_dvd_mul_right (heq ▸ dvd_mul_right d₁ e₁)
  have h₂ : d₂ ∣ d₁ := hc₂.dvd_of_dvd_mul_right (heq.symm ▸ dvd_mul_right d₂ e₂)
  have hd : d₁ = d₂ := Nat.dvd_antisymm h₁ h₂
  refine ⟨hd, ?_⟩
  rw [hd] at heq
  exact Nat.eq_of_mul_eq_mul_left (Nat.pos_of_mem_divisors hd₂) heq

/-- A practical coprime factor replaces each bounded coefficient by a disjoint
sum of actual divisors. This verifies the lifting step of the manuscript. -/
theorem bounded_coefficients_mem_divisorSumSet
    {b n : ℕ} (hb : Nat.IsPractical b) (hb0 : b ≠ 0) (hn0 : n ≠ 0)
    (hbn : Nat.Coprime b n) (c : ℕ → ℕ)
    (hc : ∀ d ∈ n.divisors, c d ≤ b) :
    b * n ∈ Erdos859.DivisorSumSet (∑ d ∈ n.divisors, d * c d) := by
  classical
  have hchoice : ∀ d : ℕ, ∃ s : Finset ℕ, s ⊆ b.divisors ∧
      (d ∈ n.divisors → c d = ∑ e ∈ s, e) := by
    intro d
    by_cases hd : d ∈ n.divisors
    · obtain ⟨s, hs, heq⟩ := hb (c d) (hc d hd)
      exact ⟨s, hs, fun _ ↦ heq⟩
    · exact ⟨∅, by simp, fun h ↦ False.elim (hd h)⟩
  choose s hs hsum using hchoice
  let pairs : Finset (Σ _d : ℕ, ℕ) := n.divisors.sigma s
  let values : Finset ℕ := pairs.image (fun p ↦ p.1 * p.2)
  have hinj : ∀ p ∈ pairs, ∀ q ∈ pairs, p.1 * p.2 = q.1 * q.2 → p = q := by
    intro p hp q hq heq
    obtain ⟨hdp, hep⟩ := mem_sigma.mp hp
    obtain ⟨hdq, heq'⟩ := mem_sigma.mp hq
    obtain ⟨hfirst, hsecond⟩ := divisor_product_injective hbn hdp hdq
      (hs p.1 hep) (hs q.1 heq') heq
    exact Sigma.ext hfirst (heq_of_eq hsecond)
  refine ⟨values, ?_, ?_⟩
  · intro z hz
    obtain ⟨p, hp, rfl⟩ := mem_image.mp hz
    obtain ⟨hdp, hep⟩ := mem_sigma.mp hp
    apply Nat.mem_divisors.mpr
    constructor
    · rw [mul_comm b n]
      exact Nat.mul_dvd_mul (Nat.mem_divisors.mp hdp).1
        (Nat.mem_divisors.mp (hs p.1 hep)).1
    · exact Nat.mul_ne_zero hb0 hn0
  · change (∑ d ∈ n.divisors, d * c d) = ∑ z ∈ pairs.image (fun p ↦ p.1 * p.2), z
    rw [sum_image hinj]
    change (∑ d ∈ n.divisors, d * c d) = ∑ p ∈ n.divisors.sigma s, p.1 * p.2
    rw [sum_sigma]
    apply sum_congr rfl
    intro d hd
    rw [hsum d hd, mul_sum]

/-- Adjoining a new prime at most `b + 1` preserves practicality. -/
theorem practical_mul_prime {b p : ℕ} (hb : Nat.IsPractical b)
    (hp : Nat.Prime p) (hcop : Nat.Coprime b p) (hpb : p ≤ b + 1) :
    Nat.IsPractical (b * p) := by
  intro m hm
  have hb0 : b ≠ 0 := by
    have := hp.two_le
    omega
  have hq : m / p ≤ b := by
    calc
      m / p ≤ (b * p) / p := Nat.div_le_div_right hm
      _ = b := Nat.mul_div_cancel b hp.pos
  have hr : m % p ≤ b := by
    have := Nat.mod_lt m hp.pos
    omega
  let c : ℕ → ℕ := fun d ↦ if d = 1 then m % p else if d = p then m / p else 0
  have hc : ∀ d ∈ p.divisors, c d ≤ b := by
    intro d _
    dsimp [c]
    split_ifs <;> omega
  have hsum : (∑ d ∈ p.divisors, d * c d) = m := by
    rw [hp.sum_divisors]
    simpa [c, hp.ne_one, add_comm] using Nat.mod_add_div m p
  have hrep := bounded_coefficients_mem_divisorSumSet hb hb0 hp.ne_zero hcop c hc
  rw [hsum] at hrep
  exact hrep