LiDigamma.lean

Declarations

Source

/-
NguyenRH1/LiDigamma.lean — tier ii′ sub-stage ii′-3a, piece (A) core
(D-066 §4): the digamma partial-fraction series and the polygamma values
at 1/2.

Mathlib has no digamma series; we build it from BohrMollerup's
`logGammaSeq x n = x·log n + log n! − Σ_{m=0}^n log(x+m)` with
`tendsto_log_gamma`. The derivative sequence
  d/dx[logGammaSeq x n] = log n − Σ_{m=0}^n 1/(x+m)
  = (log n − H_n) + [−1/x + Σ_{m=1}^n (1/m − 1/(x+m))]
converges locally uniformly on (0,∞) (the harmonic piece to −γ by
`Real.tendsto_harmonic_sub_log`; the series piece by the Weierstrass
M-test with Σ 1/m²), so the deriv-limit interchange
(`hasDerivAt_of_tendstoLocallyUniformlyOn`) yields the real identity
  (log Γ)'(x) = −γ − 1/x + Σ_{m≥1} (1/m − 1/(x+m))   (x > 0).
The complex series converges locally uniformly on the right half-plane
{0 < re s}, hence is analytic there and agrees with Complex.digamma by
the identity principle (`AnalyticOnNhd.eqOn_of_preconnected_of_frequently_eq`
through the real segment); termwise iterated derivation
(`iteratedDerivWithin_tsum`) then gives the polygamma values
  digamma^{(j)}(1/2) = (−1)^{j+1}·j!·Σ_{m≥0} (m + 1/2)^{−(j+1)}   (j ≥ 1),
and the odd-m reindexing identifies the sum as 2^{j+1}·lagariasA (j+1).

Archive references: D-066 (route β4, piece A), journal-020.
-/
import Mathlib
import NguyenRH1.LiSwaps

open Zeta23 Complex Filter Topology Real
open scoped Nat

noncomputable section

namespace NguyenRH1

/-! ## 1. The real series: summability -/

/-- The digamma series summand: 1/(m+1) − 1/(x+m+1) (m : ℕ, x : ℝ). -/
def dpsiTermR (x : ℝ) (m : ℕ) : ℝ := 1 / (m + 1 : ℝ) - 1 / (x + m + 1)

/-- The real digamma series: −γ − 1/x + Σ_{m≥0} dpsiTermR x m. -/
def DinfR (x : ℝ) : ℝ :=
  -Real.eulerMascheroniConstant - 1 / x + ∑' m : ℕ, dpsiTermR x m

/-- For 0 < x the summand is x/((m+1)(x+m+1)) ≥ 0, ≤ x/(m+1)². -/
lemma dpsiTermR_eq {x : ℝ} (hx : 0 < x) (m : ℕ) :
    dpsiTermR x m = x / ((m + 1 : ℝ) * (x + m + 1)) := by
  have hm0 : (0 : ℝ) < (m : ℝ) + 1 := by positivity
  have hxm : (0 : ℝ) < x + m + 1 := by positivity
  rw [dpsiTermR]
  field_simp
  ring

lemma dpsiTermR_nonneg {x : ℝ} (hx : 0 < x) (m : ℕ) : 0 ≤ dpsiTermR x m := by
  rw [dpsiTermR_eq hx]
  positivity

lemma dpsiTermR_le {x : ℝ} (hx : 0 < x) (m : ℕ) :
    dpsiTermR x m ≤ x / ((m : ℝ) + 1) ^ 2 := by
  have hm0 : (0 : ℝ) < (m : ℝ) + 1 := by positivity
  have hxm : (0 : ℝ) < x + m + 1 := by positivity
  rw [dpsiTermR_eq hx]
  have hCB : ((m : ℝ) + 1) ^ 2 ≤ ((m : ℝ) + 1) * (x + m + 1) := by
    rw [pow_two]
    apply mul_le_mul_of_nonneg_left _ hm0.le
    linarith [Nat.cast_nonneg (α := ℝ) m]
  have h4 : 1 / ((↑m + 1) * (x + ↑m + 1)) ≤ 1 / ((↑m + 1) ^ 2) :=
    one_div_le_one_div_of_le (by positivity) hCB
  rw [one_div, one_div] at h4
  rw [div_eq_mul_inv, div_eq_mul_inv]
  exact mul_le_mul_of_nonneg_left h4 hx.le

/-- The dominator family is summable (p-series at p = 2, shifted). -/
lemma summable_div_add_one_sq (c : ℝ) :
    Summable (fun m : ℕ => c / ((m : ℝ) + 1) ^ 2) := by
  have h := Real.summable_one_div_nat_pow.mpr (show 1 < 2 from by norm_num)
  have h2 : Summable (fun m : ℕ => (1 : ℝ) / ((m : ℝ) + 1) ^ 2) := by
    have h3 : Summable (fun m : ℕ => (1 : ℝ) / ((m + 1 : ℕ) : ℝ) ^ 2) :=
      (summable_nat_add_iff 1).mpr h
    apply h3.congr
    intro m
    push_cast
    ring_nf
  exact (h2.mul_left c).congr (fun m => by ring)

/-- The real digamma series is summable for 0 < x. -/
lemma dpsiTermR_summable {x : ℝ} (hx : 0 < x) : Summable (dpsiTermR x) :=
  Summable.of_nonneg_of_le (dpsiTermR_nonneg hx) (dpsiTermR_le hx) (summable_div_add_one_sq x)

/-! ## 2. The deriv-limit interchange: (log Γ)' = DinfR on (0,∞) -/

/-- The derivative of logGammaSeq x n in x: log n − Σ_{m=0}^n 1/(x+m). -/
def dpsiSeqR (n : ℕ) (x : ℝ) : ℝ :=
  Real.log n - ∑ m ∈ Finset.range (n + 1), 1 / (x + m)

/-- The derivative sequence decomposes as (log n − H_n) + (−1/x + partial series). -/
lemma dpsiSeqR_eq (n : ℕ) (x : ℝ) :
    dpsiSeqR n x = (Real.log n - (harmonic n : ℝ)) +
      (-1 / x + ∑ m ∈ Finset.range n, dpsiTermR x m) := by
  rw [dpsiSeqR]
  have hh : (harmonic n : ℝ) = ∑ i ∈ Finset.range n, 1 / ((i : ℝ) + 1) := by
    rw [harmonic]
    push_cast
    apply Finset.sum_congr rfl
    intro i hi
    rw [one_div]
  have h2 : ∑ m ∈ Finset.range (n + 1), 1 / (x + (m : ℝ)) =
      1 / x + ∑ i ∈ Finset.range n, 1 / (x + i + 1) := by
    rw [Finset.sum_range_succ']
    have e : ∀ i : ℕ, 1 / (x + ((i + 1 : ℕ) : ℝ)) = 1 / (x + i + 1 : ℝ) := fun i => by
      push_cast
      ring_nf
    rw [Finset.sum_congr rfl (fun i _ => e i)]
    have e0 : 1 / (x + ((0 : ℕ) : ℝ)) = 1 / x := by simp
    rw [e0]
    ring
  have h3 : ∑ i ∈ Finset.range n, dpsiTermR x i =
      ∑ i ∈ Finset.range n, 1 / ((i : ℝ) + 1) - ∑ i ∈ Finset.range n, 1 / (x + i + 1) := by
    rw [← Finset.sum_sub_distrib]
    apply Finset.sum_congr rfl
    intro i hi
    rw [dpsiTermR]
  rw [hh, h2, h3]
  ring

/-- The series part converges uniformly on compacts inside (0,∞). -/
lemma dpsiTermR_summableLocallyUniformlyOn :
    SummableLocallyUniformlyOn (fun m : ℕ => fun x : ℝ => dpsiTermR x m) (Set.Ioi 0) := by
  apply SummableLocallyUniformlyOn_of_locally_bounded isOpen_Ioi
  intro K hK hKc
  by_cases hKe : K = ∅
  · subst hKe
    exact ⟨0, summable_zero, fun m x hx => by simp at hx⟩
  · obtain ⟨x0, hx0⟩ := Set.nonempty_iff_ne_empty.mpr hKe
    obtain ⟨C, hCK, hCle⟩ := hKc.exists_isMaxOn ⟨x0, hx0⟩ (continuous_id.continuousOn)
    refine ⟨fun m => C / ((m : ℝ) + 1) ^ 2, summable_div_add_one_sq C, fun m x hx => ?_⟩
    have hx0 : 0 < x := hK hx
    have hxC : x ≤ C := hCle hx
    rw [Real.norm_eq_abs, abs_of_nonneg (dpsiTermR_nonneg hx0 m)]
    calc dpsiTermR x m ≤ x / ((m : ℝ) + 1) ^ 2 := dpsiTermR_le hx0 m
      _ ≤ C / ((m : ℝ) + 1) ^ 2 :=
          (div_le_div_iff_of_pos_right (by positivity : (0 : ℝ) < ((m : ℝ) + 1) ^ 2)).mpr hxC

/-- The derivative sequence converges locally uniformly to DinfR on (0,∞). -/
lemma dpsiSeqR_tendstoLocallyUniformlyOn :
    TendstoLocallyUniformlyOn dpsiSeqR DinfR atTop (Set.Ioi 0) := by
  have hA : TendstoLocallyUniformlyOn (fun n : ℕ => fun _ : ℝ => Real.log n - (harmonic n : ℝ))
      (fun _ => -Real.eulerMascheroniConstant) atTop (Set.Ioi 0) := by
    have h := Real.tendsto_harmonic_sub_log
    have h2 : Tendsto (fun n : ℕ => Real.log n - (harmonic n : ℝ)) atTop
        (𝓝 (-Real.eulerMascheroniConstant)) := by
      have h3 : ∀ n : ℕ, -((harmonic n : ℝ) - Real.log n) =
          Real.log n - (harmonic n : ℝ) := fun n => by ring
      exact Real.tendsto_harmonic_sub_log.neg.congr'
        (Filter.Eventually.of_forall h3)
    rw [tendstoLocallyUniformlyOn_iff_forall_isCompact isOpen_Ioi]
    intro K hK hKc
    rw [Metric.tendstoUniformlyOn_iff]
    intro ε hε
    rw [Metric.tendsto_atTop] at h2
    obtain ⟨N, hN⟩ := h2 ε hε
    rw [Filter.eventually_atTop]
    exact ⟨N, fun n hn x hx => by rw [dist_comm]; exact hN n hn⟩
  have hB : TendstoLocallyUniformlyOn (fun _ : ℕ => fun x : ℝ => -1 / x) (fun x => -1 / x)
      atTop (Set.Ioi 0) := by
    rw [tendstoLocallyUniformlyOn_iff_forall_isCompact isOpen_Ioi]
    intro K hK hKc
    rw [Metric.tendstoUniformlyOn_iff]
    intro ε hε
    rw [Filter.eventually_atTop]
    exact ⟨0, fun n _ x _ => by simpa using hε⟩
  have hC : TendstoLocallyUniformlyOn (fun n : ℕ => fun x : ℝ => ∑ m ∈ Finset.range n, dpsiTermR x m)
      (fun x => ∑' m : ℕ, dpsiTermR x m) atTop (Set.Ioi 0) := by
    rw [tendstoLocallyUniformlyOn_iff_forall_isCompact isOpen_Ioi]
    intro K hK hKc
    by_cases hKe : K = ∅
    · subst hKe
      rw [Metric.tendstoUniformlyOn_iff]
      intro ε hε
      rw [Filter.eventually_atTop]
      exact ⟨0, fun n _ x hx => by simp at hx⟩
    · obtain ⟨x0, hx0⟩ := Set.nonempty_iff_ne_empty.mpr hKe
      obtain ⟨C, hCK, hCle⟩ := hKc.exists_isMaxOn ⟨x0, hx0⟩ (continuous_id.continuousOn)
      apply tendstoUniformlyOn_tsum_nat (summable_div_add_one_sq C)
      intro m x hx
      have hx0 : 0 < x := hK hx
      have hxC : x ≤ C := hCle hx
      rw [Real.norm_eq_abs, abs_of_nonneg (dpsiTermR_nonneg hx0 m)]
      calc dpsiTermR x m ≤ x / ((m : ℝ) + 1) ^ 2 := dpsiTermR_le hx0 m
        _ ≤ C / ((m : ℝ) + 1) ^ 2 :=
            (div_le_div_iff_of_pos_right (by positivity : (0 : ℝ) < ((m : ℝ) + 1) ^ 2)).mpr hxC
  have hdecomp : ∀ n : ℕ, (Set.Ioi 0).EqOn (dpsiSeqR n)
      (fun x => (Real.log n - (harmonic n : ℝ)) +
        (-1 / x + ∑ m ∈ Finset.range n, dpsiTermR x m)) :=
    fun n x _ => dpsiSeqR_eq n x
  have hlimit : (Set.Ioi 0).EqOn
      (fun x => -Real.eulerMascheroniConstant + (-1 / x + ∑' m : ℕ, dpsiTermR x m)) DinfR :=
    fun x _ => by rw [DinfR]; ring
  exact ((hA.add (hB.add hC)).congr (fun n x _ => (dpsiSeqR_eq n x).symm)).congr_right hlimit

/-- (log Γ)'(x) = DinfR(x) for x > 0 — the interchange. -/
lemma deriv_log_Gamma_real {x : ℝ} (hx : 0 < x) :
    deriv (Real.log ∘ Real.Gamma) x = DinfR x := by
  have hfg : ∀ y ∈ Set.Ioi 0, Tendsto (fun n => BohrMollerup.logGammaSeq y n) atTop
      (𝓝 ((Real.log ∘ Real.Gamma) y)) := by
    intro y hy
    exact BohrMollerup.tendsto_log_gamma hy
  have hf : ∀ᶠ n in atTop, ∀ y ∈ Set.Ioi 0,
      HasDerivAt (fun y => BohrMollerup.logGammaSeq y n) (dpsiSeqR n y) y := by
    apply Filter.Eventually.of_forall
    intro n y hy
    have h1 : HasDerivAt (fun y => y * Real.log n) (Real.log n) y := by
      simpa using (hasDerivAt_id y).mul_const (Real.log n)
    have h2 : HasDerivAt (fun _ : ℝ => Real.log (n !)) 0 y := hasDerivAt_const _ _
    have h3 : HasDerivAt (fun y => ∑ m ∈ Finset.range (n + 1), Real.log (y + (m : ℝ)))
        (∑ m ∈ Finset.range (n + 1), (1 / (y + (m : ℝ)))) y := by
      rw [show (fun y => ∑ m ∈ Finset.range (n + 1), Real.log (y + (m : ℝ))) =
        ∑ m ∈ Finset.range (n + 1), (fun y => Real.log (y + (m : ℝ))) from by
        funext y; simp [Finset.sum_apply]]
      apply HasDerivAt.sum
      intro m hm
      have hym : (0 : ℝ) < y + m := by
        have := hy
        rw [Set.mem_Ioi] at this
        positivity
      have h4 : HasDerivAt (fun y : ℝ => y + (m : ℝ)) 1 y := by
        simpa using (hasDerivAt_id y).add_const (m : ℝ)
      have h5 := (Real.hasDerivAt_log (ne_of_gt hym)).comp y h4
      rw [mul_one] at h5
      have heq : (Real.log ∘ fun y : ℝ => y + (m : ℝ)) = fun y => Real.log (y + (m : ℝ)) := rfl
      rw [heq] at h5
      have hinv : (1 : ℝ) / (y + (m : ℝ)) = (y + (m : ℝ))⁻¹ := by rw [one_div]
      rw [hinv]
      exact h5
    have hmain : HasDerivAt (fun y => y * Real.log n + Real.log (n !) -
        ∑ m ∈ Finset.range (n + 1), Real.log (y + (m : ℝ)))
        (Real.log n + 0 - ∑ m ∈ Finset.range (n + 1), 1 / (y + (m : ℝ))) y :=
      (h1.add h2).sub h3
    have hzero : Real.log n + 0 - ∑ m ∈ Finset.range (n + 1), 1 / (y + (m : ℝ)) =
        dpsiSeqR n y := by
      rw [dpsiSeqR]
      ring
    rw [hzero] at hmain
    have heq : (fun y => BohrMollerup.logGammaSeq y n) =
        fun y => y * Real.log n + Real.log (n !) -
          ∑ m ∈ Finset.range (n + 1), Real.log (y + (m : ℝ)) := by
      funext y
      rw [BohrMollerup.logGammaSeq]
    rw [← heq] at hmain
    exact hmain
  have h := hasDerivAt_of_tendstoLocallyUniformlyOn isOpen_Ioi dpsiSeqR_tendstoLocallyUniformlyOn
    hf hfg (Set.mem_Ioi.mpr hx)
  exact h.deriv

/-- ψ_real(x) := Γ'(x)/Γ(x) = DinfR(x) for x > 0. -/
lemma psiReal_eq_DinfR {x : ℝ} (hx : 0 < x) :
    deriv Real.Gamma x / Real.Gamma x = DinfR x := by
  have h := deriv_log_Gamma_real hx
  have hG : DifferentiableAt ℝ Real.Gamma x := by
    apply Real.differentiableAt_Gamma
    intro m hmeq
    have hxx : (0 : ℝ) < x := hx
    rw [hmeq] at hxx
    have hm : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
    linarith
  have hG0 : Real.Gamma x ≠ 0 := (Real.Gamma_pos_of_pos hx).ne'
  have hchain : deriv (Real.log ∘ Real.Gamma) x = deriv Real.Gamma x / Real.Gamma x := by
    have h1 : HasDerivAt Real.Gamma (deriv Real.Gamma x) x := hG.hasDerivAt
    have h2 : HasDerivAt Real.log (Real.Gamma x)⁻¹ (Real.Gamma x) :=
      Real.hasDerivAt_log hG0
    have h3 := h2.comp x h1
    have heq : (Real.log ∘ Real.Gamma) = fun y => Real.log (Real.Gamma y) := rfl
    rw [heq] at h3
    rw [heq]
    rw [h3.deriv]
    field_simp [hG0]
  rw [hchain] at h
  exact h

/-! ## 3. The complex series: summability, analyticity, cast bridge -/

/-- The complex digamma series summand. -/
def dpsiTermC (s : ℂ) (m : ℕ) : ℂ := 1 / ((m : ℂ) + 1) - 1 / (s + m + 1)

/-- The complex digamma series: −γ − 1/s + Σ_{m≥0} dpsiTermC s m. -/
def DinfC (s : ℂ) : ℂ :=
  -(Real.eulerMascheroniConstant : ℂ) - 1 / s + ∑' m : ℕ, dpsiTermC s m

/-- The domain: the ball of radius 1/2 around 1/2 (avoids 0, −1, −2, …;
convex, contains the real segment (0,1)). -/
def dpsiDomain : Set ℂ := Metric.ball ((1 : ℂ) / 2) (1 / 2)

lemma dpsiDomain_convex : Convex ℝ dpsiDomain := convex_ball _ _
lemma dpsiDomain_isOpen : IsOpen dpsiDomain := Metric.isOpen_ball
lemma dpsiDomain_half : (1 : ℂ) / 2 ∈ dpsiDomain :=
  Metric.mem_ball_self (by norm_num : (0 : ℝ) < 1 / 2)

/-- Points of dpsiDomain have positive real part. -/
lemma dpsiDomain_re_pos {s : ℂ} (hs : s ∈ dpsiDomain) : 0 < s.re := by
  have h1 : ‖s - (1 : ℂ) / 2‖ < 1 / 2 := by
    have h := hs
    rw [dpsiDomain, Metric.mem_ball, dist_eq_norm] at h
    exact h
  have h2 : |s.re - 1 / 2| ≤ ‖s - (1 : ℂ) / 2‖ := by
    have e : s.re - 1 / 2 = (s - (1 : ℂ) / 2).re := by
      simp [Complex.sub_re]
    rw [e]
    exact Complex.abs_re_le_norm _
  have h3 : |s.re - 1 / 2| < 1 / 2 := lt_of_le_of_lt h2 h1
  have h4 : -(1 / 2) < s.re - 1 / 2 := (abs_lt.mp h3).1
  linarith

lemma dpsiTermC_pole_free {s : ℂ} (hs : s ∈ dpsiDomain) (m : ℕ) : s + m + 1 ≠ 0 := by
  intro h
  have hre : 0 < s.re := dpsiDomain_re_pos hs
  have h2 : (s + (m : ℂ) + 1).re = 0 := by rw [h]; simp
  simp only [Complex.add_re, Complex.one_re] at h2
  have h3 : (↑m : ℂ).re = (m : ℝ) := by simp
  rw [h3] at h2
  have h4 : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
  linarith

/-- The summand in factored form: dpsiTermC s m = s/((m+1)(s+m+1)). -/
lemma dpsiTermC_eq {s : ℂ} (m : ℕ) (hsm : s + m + 1 ≠ 0) :
    dpsiTermC s m = s / ((↑m + 1) * (s + ↑m + 1)) := by
  have hm0 : ((m : ℂ) + 1) ≠ 0 := by
    intro h
    have hre : ((m : ℂ) + 1).re = 0 := by rw [h]; simp
    simp [Complex.add_re, Complex.one_re] at hre
    have : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
    linarith
  rw [dpsiTermC]
  field_simp [hm0, hsm]
  ring

/-- ‖(m : ℂ) + 1‖ = m + 1. -/
lemma norm_natCast_add_one (m : ℕ) : ‖((m : ℂ) + 1)‖ = (m : ℝ) + 1 := by
  have e : ((m : ℂ) + 1) = ↑((m : ℝ) + 1) := by push_cast; ring
  rw [e, Complex.norm_real, Real.norm_eq_abs]
  exact abs_of_nonneg (by positivity)

/-- The complex series is summable on dpsiDomain. -/
lemma dpsiTermC_summable {s : ℂ} (hs : s ∈ dpsiDomain) : Summable (dpsiTermC s) := by
  have hnorm : ∀ m : ℕ, ‖dpsiTermC s m‖ ≤ ‖s‖ / ((m : ℝ) + 1) ^ 2 := by
    intro m
    have hsm : s + m + 1 ≠ 0 := dpsiTermC_pole_free hs m
    rw [dpsiTermC_eq m hsm]
    rw [norm_div]
    have hden : ((m : ℝ) + 1) * ((m : ℝ) + 1) ≤ ‖((m : ℂ) + 1) * (s + ↑m + 1)‖ := by
      rw [norm_mul, norm_natCast_add_one]
      apply mul_le_mul_of_nonneg_left _ (by positivity : (0 : ℝ) ≤ (m : ℝ) + 1)
      have hre : 0 < s.re := dpsiDomain_re_pos hs
      have h5 : s.re + m + 1 ≤ ‖s + ↑m + 1‖ := by
        have h6 : (s + ↑m + 1).re = s.re + m + 1 := by
          have h3 : (↑m : ℂ).re = (m : ℝ) := by simp
          simp [Complex.add_re, Complex.one_re, h3]
        rw [← h6]
        exact Complex.re_le_norm _
      linarith [h5]
    have hpos : (0 : ℝ) < ((m : ℝ) + 1) * ((m : ℝ) + 1) := by positivity
    calc ‖s‖ / ‖((m : ℂ) + 1) * (s + ↑m + 1)‖ ≤ ‖s‖ / (((m : ℝ) + 1) * ((m : ℝ) + 1)) :=
          div_le_div_of_nonneg_left (norm_nonneg _) hpos hden
      _ = ‖s‖ / ((m : ℝ) + 1) ^ 2 := by ring
  exact Summable.of_norm_bounded (summable_div_add_one_sq ‖s‖) hnorm

/-! ## 4. Analyticity and the identity with Complex.digamma -/

/-- The series converges uniformly on dpsiDomain (M-test, Σ 1/(m+1)²). -/
lemma dpsiTermC_hasSumUniformlyOn :
    HasSumUniformlyOn (fun m => fun s : ℂ => dpsiTermC s m)
      (fun s => ∑' m : ℕ, dpsiTermC s m) dpsiDomain := by
  apply HasSumUniformlyOn.of_norm_le_summable (summable_div_add_one_sq 1)
  intro m s hs
  have hsm : s + m + 1 ≠ 0 := dpsiTermC_pole_free hs m
  have hre : 0 < s.re := dpsiDomain_re_pos hs
  rw [dpsiTermC_eq m hsm, norm_div]
  have hden : ((m : ℝ) + 1) * ((m : ℝ) + 1) ≤ ‖((m : ℂ) + 1) * (s + ↑m + 1)‖ := by
    rw [norm_mul, norm_natCast_add_one]
    apply mul_le_mul_of_nonneg_left _ (by positivity : (0 : ℝ) ≤ (m : ℝ) + 1)
    have h5 : s.re + m + 1 ≤ ‖s + ↑m + 1‖ := by
      have h6 : (s + ↑m + 1).re = s.re + m + 1 := by
        have h3 : (↑m : ℂ).re = (m : ℝ) := by simp
        simp [Complex.add_re, Complex.one_re, h3]
      rw [← h6]
      exact Complex.re_le_norm _
    linarith [h5]
  have hs1 : ‖s‖ ≤ 1 := by
    have h := hs
    rw [dpsiDomain, Metric.mem_ball, dist_eq_norm] at h
    have h12 : ‖((1 : ℂ) / 2)‖ = 1 / 2 := by
      rw [norm_div, norm_one]
      norm_num
    calc ‖s‖ ≤ ‖s - (1 : ℂ) / 2‖ + ‖((1 : ℂ) / 2)‖ := by
          have e : s = (s - (1 : ℂ) / 2) + (1 : ℂ) / 2 := by ring
          nth_rewrite 1 [e]
          exact norm_add_le _ _
      _ ≤ 1 / 2 + 1 / 2 := by linarith [h, h12]
      _ = 1 := by ring
  have hpos : (0 : ℝ) < ((m : ℝ) + 1) * ((m : ℝ) + 1) := by positivity
  calc ‖s‖ / ‖((m : ℂ) + 1) * (s + ↑m + 1)‖ ≤ ‖s‖ / (((m : ℝ) + 1) * ((m : ℝ) + 1)) :=
        div_le_div_of_nonneg_left (norm_nonneg _) hpos hden
    _ = ‖s‖ / ((m : ℝ) + 1) ^ 2 := by ring
    _ ≤ 1 / ((m : ℝ) + 1) ^ 2 :=
        (div_le_div_iff_of_pos_right (by positivity : (0 : ℝ) < ((m : ℝ) + 1) ^ 2)).mpr hs1

/-- The series function is differentiable on dpsiDomain (Weierstrass convergence). -/
lemma dpsiTermC_tsum_differentiableOn :
    DifferentiableOn ℂ (fun s => ∑' m : ℕ, dpsiTermC s m) dpsiDomain := by
  have htend : TendstoLocallyUniformlyOn
      (fun t : Finset ℕ => fun s => ∑ m ∈ t, dpsiTermC s m)
      (fun s => ∑' m : ℕ, dpsiTermC s m) atTop dpsiDomain :=
    TendstoUniformlyOn.tendstoLocallyUniformlyOn
      (hasSumUniformlyOn_iff_tendstoUniformlyOn.mp dpsiTermC_hasSumUniformlyOn)
  apply htend.differentiableOn _ dpsiDomain_isOpen
  apply Filter.Eventually.of_forall
  intro t s hs
  rw [show (fun s => ∑ m ∈ t, dpsiTermC s m) = ∑ m ∈ t, (fun s => dpsiTermC s m) from by
    funext s; simp [Finset.sum_apply]]
  apply DifferentiableAt.differentiableWithinAt
  apply DifferentiableAt.sum
  intro m hm
  have hsm : s + m + 1 ≠ 0 := dpsiTermC_pole_free hs m
  have h1 : DifferentiableAt ℂ (fun s : ℂ => 1 / ((m : ℂ) + 1)) s := differentiableAt_const _
  have h2 : DifferentiableAt ℂ (fun s : ℂ => 1 / (s + ↑m + 1)) s := by
    have h3 : DifferentiableAt ℂ (fun s : ℂ => s + ↑m + 1) s :=
      (differentiableAt_id.add (differentiableAt_const _)).add (differentiableAt_const _)
    exact (differentiableAt_const (1 : ℂ)).div h3 hsm
  exact h1.sub h2

/-- DinfC is analytic on dpsiDomain (Weierstrass convergence theorem). -/
lemma DinfC_analyticOn : AnalyticOnNhd ℂ DinfC dpsiDomain := by
  have hconst : DifferentiableOn ℂ (fun s : ℂ => -(Real.eulerMascheroniConstant : ℂ) - 1 / s)
      dpsiDomain := by
    intro s hs
    have hs0 : s ≠ 0 := by
      have h := dpsiDomain_re_pos hs
      intro hsz
      rw [hsz] at h
      simp at h
    have h1 : DifferentiableAt ℂ (fun s : ℂ => 1 / s) s :=
      (differentiableAt_const (1 : ℂ)).div differentiableAt_id hs0
    exact (differentiableAt_const _).sub h1 |>.differentiableWithinAt
  have hfull : DifferentiableOn ℂ DinfC dpsiDomain := by
    have e : DinfC = (fun s : ℂ => -(Real.eulerMascheroniConstant : ℂ) - 1 / s) +
        (fun s => ∑' m : ℕ, dpsiTermC s m) := by
      funext s
      rw [DinfC]
      simp [Pi.add_apply]
    rw [e]
    exact hconst.add dpsiTermC_tsum_differentiableOn
  exact hfull.analyticOnNhd dpsiDomain_isOpen

/-- Complex.digamma is analytic on dpsiDomain (Gamma nonzero + differentiable there). -/
lemma digamma_analyticOn : AnalyticOnNhd ℂ Complex.digamma dpsiDomain := by
  have hpole : ∀ z ∈ dpsiDomain, ∀ m : ℕ, z ≠ -(m : ℂ) := by
    intro z hz m hmeq
    have hre : 0 < z.re := dpsiDomain_re_pos hz
    rw [hmeq] at hre
    simp at hre
    have hm : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
    linarith
  have hGammaAn : AnalyticOnNhd ℂ Complex.Gamma dpsiDomain := by
    apply DifferentiableOn.analyticOnNhd _ dpsiDomain_isOpen
    intro z hz
    exact (Complex.differentiableAt_Gamma z (hpole z hz)).differentiableWithinAt
  have hdig : Complex.digamma = fun s => deriv Complex.Gamma s / Complex.Gamma s := by
    funext s
    rw [Complex.digamma_def, logDeriv_apply]
  rw [hdig]
  intro s hs
  exact (hGammaAn.deriv s hs).div (hGammaAn s hs)
    (Complex.Gamma_ne_zero (hpole s hs))

/-- The cast bridge: DinfC(↑x) = ↑(DinfR x) for 0 < x. -/
lemma DinfC_ofReal {x : ℝ} (_hx : 0 < x) : DinfC (↑x) = ↑(DinfR x) := by
  rw [DinfC, DinfR]
  push_cast
  congr 1
  apply tsum_congr
  intro m
  rw [dpsiTermR, dpsiTermC]
  push_cast
  ring

/-- Complex.digamma agrees with DinfC on the real segment inside the domain. -/
lemma digamma_eq_DinfC_ofReal {x : ℝ} (hx : x ∈ Set.Ioo 0 1) :
    Complex.digamma (↑x) = DinfC (↑x) := by
  have hx0 : 0 < x := hx.1
  rw [DinfC_ofReal hx0]
  -- digamma ↑x = deriv Gamma ↑x / Gamma ↑x = ↑(Γ'(x)/Γ(x)) = ↑(DinfR x)
  have hGx : Complex.Gamma (↑x : ℂ) = ↑(Real.Gamma x) := Complex.Gamma_ofReal x
  have hderiv : deriv Complex.Gamma (↑x : ℂ) = ↑(deriv Real.Gamma x) := by
    have h2 : HasDerivAt Real.Gamma (deriv Real.Gamma x) x := by
      have hd : DifferentiableAt ℝ Real.Gamma x := by
        apply Real.differentiableAt_Gamma
        intro m hmeq
        rw [hmeq] at hx0
        simp at hx0
        have hm : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
        linarith
      exact hd.hasDerivAt
    have h1 : HasDerivAt (fun y : ℝ => (↑(Real.Gamma y) : ℂ)) (↑(deriv Real.Gamma x) : ℂ) x := by
      have h3 := Complex.ofRealCLM.hasFDerivAt.comp_hasDerivAt x h2
      rw [Complex.ofRealCLM_apply] at h3
      have heq : (Complex.ofRealCLM ∘ Real.Gamma) = fun y : ℝ => ↑(Real.Gamma y) := by
        funext y
        rw [Function.comp_apply, Complex.ofRealCLM_apply]
      rw [heq] at h3
      exact h3
    -- deriv (Gamma ∘ ofReal) x = deriv Gamma ↑x · 1
    have h6 : deriv (fun y : ℝ => Complex.Gamma (↑y)) x = deriv Complex.Gamma (↑x : ℂ) := by
      have h7 : HasDerivAt (fun y : ℝ => (↑y : ℂ)) (1 : ℂ) x := by
        have h8 := Complex.ofRealCLM.hasDerivAt (x := x)
        rw [Complex.ofRealCLM_apply] at h8
        exact h8
      have h9 : HasDerivAt Complex.Gamma (deriv Complex.Gamma (↑x : ℂ)) (↑x : ℂ) := by
        apply (Complex.differentiableAt_Gamma _ _).hasDerivAt
        intro m hmeq
        have h2 : (x : ℝ) = -(m : ℝ) := by
          have h3 : (↑x : ℂ).re = (-(m : ℂ)).re := by rw [hmeq]
          simp at h3
          exact h3
        have hm : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
        linarith
      have h10 := h9.comp x h7
      -- comp gives deriv Gamma ↑x * 1
      rw [mul_one] at h10
      have heq : (Complex.Gamma ∘ fun y : ℝ => (↑y : ℂ)) = fun y : ℝ => Complex.Gamma (↑y) := rfl
      rw [heq] at h10
      exact h10.deriv
    have h5 : (fun y : ℝ => Complex.Gamma (↑y)) = fun y : ℝ => (↑(Real.Gamma y) : ℂ) :=
      funext fun y => Complex.Gamma_ofReal y
    have h1d : deriv (fun y : ℝ => (↑(Real.Gamma y) : ℂ)) x = ↑(deriv Real.Gamma x) := h1.deriv
    rw [← h5] at h1d
    rw [h6] at h1d
    exact h1d
  rw [Complex.digamma_def, logDeriv_apply, hderiv, hGx]
  rw [← Complex.ofReal_div]
  congr 1
  exact psiReal_eq_DinfR hx0

/-! ## 5. Congruence helpers and the inverse-power iterated derivative -/

/-- If f = g on a nbd of x, then deriv f = deriv g on a (smaller) nbd of x. -/
lemma eventuallyEq_deriv {f g : ℂ → ℂ} {x : ℂ} (h : f =ᶠ[𝓝 x] g) :
    deriv f =ᶠ[𝓝 x] deriv g := by
  have hS : {y | f y = g y} ∈ 𝓝 x := h
  have hxint : x ∈ interior {y | f y = g y} := by
    rw [mem_interior_iff_mem_nhds]
    exact hS
  filter_upwards [IsOpen.mem_nhds isOpen_interior hxint] with y hy
  have h3 : f =ᶠ[𝓝 y] g := by
    apply Filter.eventually_of_mem (IsOpen.mem_nhds isOpen_interior hy)
    intro z hz
    have h4 : z ∈ {y | f y = g y} := interior_subset hz
    exact h4
  exact h3.deriv_eq

lemma eventuallyEq_iteratedDeriv {f g : ℂ → ℂ} {x : ℂ} (h : f =ᶠ[𝓝 x] g) (n : ℕ) :
    iteratedDeriv n f =ᶠ[𝓝 x] iteratedDeriv n g := by
  induction n with
  | zero => simpa [iteratedDeriv_zero] using h
  | succ n ih =>
    rw [iteratedDeriv_succ, iteratedDeriv_succ]
    exact eventuallyEq_deriv ih

/-- Pointwise congruence of iteratedDeriv under eventual equality. -/
lemma iteratedDeriv_congr_of_eventuallyEq {f g : ℂ → ℂ} {x : ℂ} (h : f =ᶠ[𝓝 x] g) (n : ℕ) :
    iteratedDeriv n f x = iteratedDeriv n g x :=
  (eventuallyEq_iteratedDeriv h n).self_of_nhds

/-- The iterated derivative of s ↦ −(s+c)⁻¹:
(d/ds)^k [−(s+c)⁻¹] = −(−1)^k·k!·(s+c)^{−(k+1)}. -/
lemma iteratedDeriv_neg_zpow_neg_one_add (k : ℕ) (c : ℂ) {s : ℂ} (hs : s + c ≠ 0) :
    iteratedDeriv k (fun s : ℂ => -((s + c) ^ (-1 : ℤ))) s =
      -((-1 : ℂ)) ^ k * (k ! : ℂ) * (s + c) ^ (-(k : ℤ) - 1) := by
  induction k generalizing s with
  | zero =>
    simp
  | succ k ih =>
    rw [iteratedDeriv_succ]
    -- replace iteratedDeriv k F by its closed form on a nbd of s
    have hopen : IsOpen {t : ℂ | t + c ≠ 0} := by
      have hcont : Continuous (fun t : ℂ => t + c) := continuous_id.add continuous_const
      exact isOpen_compl_singleton.preimage hcont
    have hev : (fun t => iteratedDeriv k (fun u : ℂ => -((u + c) ^ (-1 : ℤ))) t) =ᶠ[𝓝 s]
        (fun t => -((-1 : ℂ)) ^ k * (k ! : ℂ) * (t + c) ^ (-(k : ℤ) - 1)) := by
      filter_upwards [hopen.mem_nhds hs] with t ht
      exact ih ht
    rw [hev.deriv_eq]
    -- compute deriv of the closed form at s
    have h1 : HasDerivAt (fun t : ℂ => (t + c) ^ (-(k : ℤ) - 1))
        (((-(k : ℤ) - 1 : ℤ) : ℂ) * (s + c) ^ (-(k : ℤ) - 1 - 1)) s := by
      have h2 : HasDerivAt (fun t : ℂ => t + c) 1 s := by
        simpa using (hasDerivAt_id s).add_const c
      have h3 := (hasDerivAt_zpow (-(k : ℤ) - 1) (s + c) (Or.inl hs)).comp s h2
      rw [mul_one] at h3
      have heq : ((fun x : ℂ => x ^ (-(k : ℤ) - 1)) ∘ fun t : ℂ => t + c) =
          fun t => (t + c) ^ (-(k : ℤ) - 1) := rfl
      rw [heq] at h3
      exact h3
    have h4 : deriv (fun t : ℂ => -((-1 : ℂ)) ^ k * (k ! : ℂ) * (t + c) ^ (-(k : ℤ) - 1)) s =
        -((-1 : ℂ)) ^ k * (k ! : ℂ) *
          (((-(k : ℤ) - 1 : ℤ) : ℂ) * (s + c) ^ (-(k : ℤ) - 1 - 1)) :=
      (h1.const_mul _).deriv
    rw [h4]
    -- normalize the exponents to −(k:ℤ)−2, then the constants
    have hexp1 : (-(k : ℤ) - 1 - 1 : ℤ) = (-(k : ℤ) - 2 : ℤ) := by ring
    have hexp2 : (-((k + 1 : ℕ) : ℤ) - 1 : ℤ) = (-(k : ℤ) - 2 : ℤ) := by push_cast; ring
    rw [hexp1, hexp2]
    have hconst : ((-(k : ℤ) - 1 : ℤ) : ℂ) = -((k : ℂ) + 1) := by push_cast; ring
    have hfact : (((k + 1) ! : ℕ) : ℂ) = ((k : ℂ) + 1) * (k ! : ℂ) := by
      rw [Nat.factorial_succ]
      push_cast
      ring
    rw [hconst, hfact, pow_succ]
    ring

/-! ## 6. The identity digamma = DinfC on the domain (identity principle) -/

/-- digamma and DinfC agree on dpsiDomain. -/
lemma digamma_eq_DinfC_on : Set.EqOn Complex.digamma DinfC dpsiDomain := by
  apply AnalyticOnNhd.eqOn_of_preconnected_of_frequently_eq digamma_analyticOn DinfC_analyticOn
    dpsiDomain_convex.isPreconnected dpsiDomain_half
  rw [frequently_iff]
  intro U hU
  rw [mem_nhdsWithin_iff_exists_mem_nhds_inter] at hU
  obtain ⟨V, hV, hVU⟩ := hU
  rw [Metric.mem_nhds_iff] at hV
  obtain ⟨δ, hδ, hδV⟩ := hV
  have hpos : (0 : ℝ) < min (δ / 2) (1 / 4) := lt_min (by positivity) (by norm_num)
  have hxI : (1 / 2 : ℝ) + min (δ / 2) (1 / 4) ∈ Set.Ioo 0 1 := by
    constructor
    · linarith
    · have h1 : min (δ / 2) (1 / 4) ≤ 1 / 4 := min_le_right _ _
      linarith
  refine ⟨↑((1 / 2 : ℝ) + min (δ / 2) (1 / 4)), hVU ?_, ?_⟩
  · constructor
    · apply hδV
      rw [Metric.mem_ball, dist_eq_norm]
      have e : (↑((1 / 2 : ℝ) + min (δ / 2) (1 / 4)) : ℂ) - (1 : ℂ) / 2 =
          ↑(min (δ / 2) (1 / 4)) := by push_cast; ring
      rw [e, Complex.norm_real, Real.norm_eq_abs, abs_of_nonneg hpos.le]
      calc min (δ / 2) (1 / 4) ≤ δ / 2 := min_le_left _ _
        _ < δ := by linarith [hδ]
    · intro hcontra
      have h2 : (1 / 2 : ℝ) + min (δ / 2) (1 / 4) = 1 / 2 := by
        have h3 : (↑((1 / 2 : ℝ) + min (δ / 2) (1 / 4)) : ℂ) = (1 : ℂ) / 2 :=
          Set.mem_singleton_iff.mp hcontra
        have h4 : ((1 : ℂ) / 2) = ↑((1 : ℝ) / 2) := by push_cast; ring
        rw [h4] at h3
        exact Complex.ofReal_inj.mp h3
      linarith
  · exact digamma_eq_DinfC_ofReal hxI

/-! ## 7. The iterated derivatives at 1/2 -/

/-- The shifted p-series is summable for any exponent ≥ 2. -/
lemma summable_div_add_one_pow (c : ℝ) (p : ℕ) (hp : 2 ≤ p) :
    Summable (fun m : ℕ => c / ((m : ℝ) + 1) ^ p) := by
  have h := Real.summable_one_div_nat_pow.mpr (show 1 < p by omega)
  have h2 : Summable (fun m : ℕ => (1 : ℝ) / ((m : ℝ) + 1) ^ p) := by
    have h3 : Summable (fun m : ℕ => (1 : ℝ) / ((m + 1 : ℕ) : ℝ) ^ p) :=
      (summable_nat_add_iff 1).mpr h
    apply h3.congr
    intro m
    push_cast
    ring_nf
  exact (h2.mul_left c).congr (fun m => by ring)

/-- The one-div form of the zpow iterated derivative. -/
lemma iteratedDeriv_neg_one_div_add_const (k : ℕ) (c : ℂ) {s : ℂ} (hs : s + c ≠ 0) :
    iteratedDeriv k (fun s : ℂ => -(1 / (s + c))) s =
      -((-1 : ℂ)) ^ k * (k ! : ℂ) * (s + c) ^ (-(k : ℤ) - 1) := by
  have e : (fun s : ℂ => -(1 / (s + c))) = fun s => -((s + c) ^ (-1 : ℤ)) := by
    funext s
    rw [zpow_neg_one, one_div]
  rw [e]
  exact iteratedDeriv_neg_zpow_neg_one_add k c hs

/-- The k-th iterated derivative (k ≥ 1) of a summand on the domain. -/
lemma iteratedDerivWithin_dpsiTermC (k : ℕ) (hk : 1 ≤ k) (m : ℕ) {s : ℂ}
    (hs : s ∈ dpsiDomain) :
    iteratedDerivWithin k (fun s : ℂ => dpsiTermC s m) dpsiDomain s =
      -((-1 : ℂ)) ^ k * (k ! : ℂ) * (s + ↑m + 1) ^ (-(k : ℤ) - 1) := by
  have hsm : s + m + 1 ≠ 0 := dpsiTermC_pole_free hs m
  rw [iteratedDerivWithin_of_isOpen dpsiDomain_isOpen hs]
  have heq : (fun s : ℂ => dpsiTermC s m) =
      (fun s : ℂ => 1 / ((m : ℂ) + 1)) + fun s : ℂ => -(1 / (s + ↑m + 1)) := by
    funext s
    simp only [Pi.add_apply]
    rw [dpsiTermC]
    ring
  rw [heq]
  have hcf : ContDiffAt ℂ k (fun s : ℂ => 1 / ((m : ℂ) + 1)) s := contDiff_const.contDiffAt
  have hcg : ContDiffAt ℂ k (fun s : ℂ => -(1 / (s + ↑m + 1))) s := by
    have h2 : AnalyticAt ℂ (fun s : ℂ => s + ↑m + 1) s :=
      (analyticAt_id.add analyticAt_const).add analyticAt_const
    have h3 : AnalyticAt ℂ (fun s : ℂ => 1 / (s + ↑m + 1)) s :=
      analyticAt_const.div h2 hsm
    exact h3.neg.contDiffAt
  rw [iteratedDeriv_add hcf hcg]
  rw [iteratedDeriv_const, if_neg (by omega : k ≠ 0)]
  have h4 : iteratedDeriv k (fun s : ℂ => -(1 / (s + ↑m + 1))) s =
      -((-1 : ℂ)) ^ k * (k ! : ℂ) * (s + ↑m + 1) ^ (-(k : ℤ) - 1) := by
    have h5 : (fun s : ℂ => -(1 / (s + ↑m + 1))) =
        fun s : ℂ => -(1 / (s + ((m : ℂ) + 1))) := by
      funext s
      rw [add_assoc]
    rw [h5]
    rw [iteratedDeriv_neg_one_div_add_const k ((m : ℂ) + 1) (by rw [← add_assoc]; exact hsm)]
    rw [← add_assoc]
  rw [h4]
  rw [zero_add]

/-- The k-th derivative family (k ≥ 1) is summable locally uniformly on the domain. -/
lemma dpsiTermC_deriv_summableLocallyUniformlyOn (k : ℕ) (hk : 1 ≤ k) :
    SummableLocallyUniformlyOn
      (fun m => (iteratedDerivWithin k (fun s : ℂ => dpsiTermC s m) dpsiDomain)) dpsiDomain := by
  apply SummableLocallyUniformlyOn_of_locally_bounded dpsiDomain_isOpen
  intro K hK hKc
  refine ⟨fun m => (k ! : ℝ) / ((m : ℝ) + 1) ^ (k + 1),
    summable_div_add_one_pow _ (k + 1) (by omega), fun m s hs => ?_⟩
  have hsd : s ∈ dpsiDomain := hK hs
  rw [iteratedDerivWithin_dpsiTermC k hk m hsd]
  have hre : 0 < s.re := dpsiDomain_re_pos hsd
  have hnormS : (m : ℝ) + 1 ≤ ‖s + ↑m + 1‖ := by
    have h6 : (s + ↑m + 1).re = s.re + m + 1 := by
      have h3 : (↑m : ℂ).re = (m : ℝ) := by simp
      simp [Complex.add_re, Complex.one_re, h3]
    have h5 : s.re + m + 1 ≤ ‖s + ↑m + 1‖ := by
      rw [← h6]
      exact Complex.re_le_norm _
    linarith [h5, hre]
  rw [norm_mul, norm_mul]
  have hsign : ‖-((-1 : ℂ)) ^ k‖ = 1 := by
    have h7 : ‖((-1 : ℂ))‖ = 1 := by rw [norm_neg, norm_one]
    rw [norm_neg, norm_pow, h7, one_pow]
  rw [hsign, one_mul, RCLike.norm_natCast, norm_zpow]
  rw [show (-(k : ℤ) - 1 : ℤ) = -(((k + 1 : ℕ) : ℤ)) by push_cast; ring, zpow_neg, zpow_natCast]
  show (↑k ! : ℝ) * (‖s + ↑m + 1‖ ^ (k + 1))⁻¹ ≤ (k ! : ℝ) / ((m : ℝ) + 1) ^ (k + 1)
  rw [div_eq_mul_inv]
  apply mul_le_mul_of_nonneg_left _ (Nat.cast_nonneg _)
  rw [← one_div, ← one_div]
  apply one_div_le_one_div_of_le (by positivity)
  exact pow_le_pow_left₀ (by positivity : (0 : ℝ) ≤ (m : ℝ) + 1) hnormS (k + 1)

/-- The iterated derivative of the series function at 1/2 (termwise). -/
lemma iteratedDerivWithin_tsum_dpsiTermC (k : ℕ) (hk : 1 ≤ k) :
    iteratedDerivWithin k (fun s : ℂ => ∑' m : ℕ, dpsiTermC s m) dpsiDomain ((1 : ℂ) / 2) =
      ∑' m : ℕ, -((-1 : ℂ)) ^ k * (k ! : ℂ) * ((1 / 2 : ℂ) + ↑m + 1) ^ (-(k : ℤ) - 1) := by
  rw [iteratedDerivWithin_tsum k dpsiDomain_isOpen dpsiDomain_half
    (fun s hs => dpsiTermC_summable hs)
    (fun j hj1 _ => dpsiTermC_deriv_summableLocallyUniformlyOn j hj1)
    ?_]
  · apply tsum_congr
    intro m
    exact iteratedDerivWithin_dpsiTermC k hk m dpsiDomain_half
  · intro n j r hjm hr
    rcases Nat.eq_zero_or_pos j with hj0 | hj0
    · subst hj0
      rw [iteratedDerivWithin_zero]
      have hsm : r + n + 1 ≠ 0 := dpsiTermC_pole_free hr n
      have heq : (fun s : ℂ => dpsiTermC s n) = fun s => 1 / ((n : ℂ) + 1) - 1 / (s + ↑n + 1) := by
        funext s; rw [dpsiTermC]
      rw [heq]
      apply DifferentiableAt.sub (differentiableAt_const _)
      have h3 : DifferentiableAt ℂ (fun s : ℂ => s + ↑n + 1) r :=
        (differentiableAt_id.add (differentiableAt_const _)).add (differentiableAt_const _)
      exact (differentiableAt_const (1 : ℂ)).div h3 hsm
    · have h5 : (fun r => iteratedDerivWithin j (fun z : ℂ => dpsiTermC z n) dpsiDomain r) =ᶠ[𝓝 r]
          fun r => -((-1 : ℂ)) ^ j * (j ! : ℂ) * (r + ↑n + 1) ^ (-(j : ℤ) - 1) := by
        filter_upwards [dpsiDomain_isOpen.mem_nhds hr] with t ht
        exact iteratedDerivWithin_dpsiTermC j (by omega) n ht
      apply h5.differentiableAt_iff.mpr
      have hsm : r + n + 1 ≠ 0 := dpsiTermC_pole_free hr n
      have h4 : DifferentiableAt ℂ (fun s : ℂ => s + ↑n + 1) r :=
        (differentiableAt_id.add (differentiableAt_const _)).add (differentiableAt_const _)
      have h6 : DifferentiableAt ℂ (fun s : ℂ => (s + ↑n + 1) ^ (-(j : ℤ) - 1)) r := by
        have h7 := ((hasDerivAt_zpow (-(j : ℤ) - 1) (r + ↑n + 1) (Or.inl hsm)).differentiableAt).comp r h4
        have heq : ((fun x : ℂ => x ^ (-(j : ℤ) - 1)) ∘ fun s : ℂ => s + ↑n + 1) =
            fun s : ℂ => (s + ↑n + 1) ^ (-(j : ℤ) - 1) := rfl
        rw [heq] at h7
        exact h7
      exact (differentiableAt_const _).mul h6


/-! ## 9. The values at 1/2 and the transfer to digamma -/

/-- The (1/2 + m)-power series is summable (norm ≤ 2^{k+1}/(m+1)^{k+1}). -/
lemma summable_zpow_half_add_nat (k : ℕ) (hk : 1 ≤ k) :
    Summable (fun m : ℕ => ((1 / 2 : ℂ) + ↑m) ^ (-(k : ℤ) - 1)) := by
  apply Summable.of_norm_bounded (summable_div_add_one_pow (2 ^ (k + 1)) (k + 1) (by omega))
  intro m
  have hpos : (0 : ℝ) < (m : ℝ) + 1 / 2 := by positivity
  have h1 : ‖((1 / 2 : ℂ) + ↑m)‖ = (m : ℝ) + 1 / 2 := by
    have e : ((1 / 2 : ℂ) + ↑m) = ↑((m : ℝ) + 1 / 2) := by push_cast; ring
    rw [e, Complex.norm_real, Real.norm_eq_abs, abs_of_nonneg hpos.le]
  rw [norm_zpow, h1]
  rw [show (-(k : ℤ) - 1 : ℤ) = -(((k + 1 : ℕ) : ℤ)) by push_cast; ring, zpow_neg, zpow_natCast]
  rw [← one_div]
  show (1 : ℝ) / ((m : ℝ) + 1 / 2) ^ (k + 1) ≤ 2 ^ (k + 1) / ((m : ℝ) + 1) ^ (k + 1)
  rw [div_le_div_iff₀ (by positivity : (0 : ℝ) < ((m : ℝ) + 1 / 2) ^ (k + 1))
    (by positivity : (0 : ℝ) < ((m : ℝ) + 1) ^ (k + 1))]
  rw [one_mul]
  rw [show (2 : ℝ) ^ (k + 1) * ((m : ℝ) + 1 / 2) ^ (k + 1) =
    (2 * ((m : ℝ) + 1 / 2)) ^ (k + 1) by rw [mul_pow]]
  apply pow_le_pow_left₀ (by positivity : (0 : ℝ) ≤ (m : ℝ) + 1)
  ring_nf
  linarith [Nat.cast_nonneg (α := ℝ) m]

/-- The polygamma values of DinfC at 1/2 (j ≥ 1). -/
lemma iteratedDeriv_DinfC_half (k : ℕ) (hk : 1 ≤ k) :
    iteratedDeriv k DinfC ((1 : ℂ) / 2) =
      ((-1 : ℂ)) ^ (k + 1) * (k ! : ℂ) * ∑' m : ℕ, ((1 / 2 : ℂ) + ↑m) ^ (-(k : ℤ) - 1) := by
  have hsign : ∀ k : ℕ, -((-1 : ℂ)) ^ k = (-1 : ℂ) ^ (k + 1) := fun k => by
    rw [pow_succ]; ring
  have hsplit : DinfC = (fun s : ℂ => -(Real.eulerMascheroniConstant : ℂ) - 1 / s) +
      (fun s => ∑' m : ℕ, dpsiTermC s m) := by
    funext s
    rw [DinfC]
    simp [Pi.add_apply]
  rw [hsplit]
  have hcfA : ContDiffAt ℂ k (fun s : ℂ => -(Real.eulerMascheroniConstant : ℂ) - 1 / s)
      ((1 : ℂ) / 2) := by
    have h1 : AnalyticAt ℂ (fun s : ℂ => 1 / s) ((1 : ℂ) / 2) :=
      analyticAt_const.div analyticAt_id (by norm_num)
    exact (analyticAt_const.sub h1).contDiffAt
  have hcfB : ContDiffAt ℂ k (fun s => ∑' m : ℕ, dpsiTermC s m) ((1 : ℂ) / 2) := by
    have h2 : AnalyticOnNhd ℂ (fun s => ∑' m : ℕ, dpsiTermC s m) dpsiDomain :=
      dpsiTermC_tsum_differentiableOn.analyticOnNhd dpsiDomain_isOpen
    exact (h2 _ dpsiDomain_half).contDiffAt
  rw [iteratedDeriv_add hcfA hcfB]
  have hA : iteratedDeriv k (fun s : ℂ => -(Real.eulerMascheroniConstant : ℂ) - 1 / s)
      ((1 : ℂ) / 2) = ((-1 : ℂ)) ^ (k + 1) * (k ! : ℂ) * ((1 : ℂ) / 2) ^ (-(k : ℤ) - 1) := by
    have heq : (fun s : ℂ => -(Real.eulerMascheroniConstant : ℂ) - 1 / s) =
        (fun _ : ℂ => -(Real.eulerMascheroniConstant : ℂ)) + fun s => -(1 / (s + 0)) := by
      funext s
      simp [Pi.add_apply]
      ring
    rw [heq]
    have hcf1 : ContDiffAt ℂ k (fun _ : ℂ => -(Real.eulerMascheroniConstant : ℂ)) ((1 : ℂ) / 2) :=
      contDiff_const.contDiffAt
    have hcf2 : ContDiffAt ℂ k (fun s : ℂ => -(1 / (s + 0))) ((1 : ℂ) / 2) := by
      have h2 : AnalyticAt ℂ (fun s : ℂ => s + 0) ((1 : ℂ) / 2) :=
        analyticAt_id.add analyticAt_const
      have h3 : AnalyticAt ℂ (fun s : ℂ => 1 / (s + 0)) ((1 : ℂ) / 2) :=
        analyticAt_const.div h2 (by simp)
      exact h3.neg.contDiffAt
    rw [iteratedDeriv_add hcf1 hcf2]
    rw [iteratedDeriv_const, if_neg (by omega : k ≠ 0)]
    rw [iteratedDeriv_neg_one_div_add_const k 0 (by simp)]
    rw [add_zero, zero_add, hsign]
  have hB : iteratedDeriv k (fun s => ∑' m : ℕ, dpsiTermC s m) ((1 : ℂ) / 2) =
      ((-1 : ℂ)) ^ (k + 1) * (k ! : ℂ) * ∑' m : ℕ, ((1 / 2 : ℂ) + ↑m + 1) ^ (-(k : ℤ) - 1) := by
    rw [← iteratedDerivWithin_of_isOpen dpsiDomain_isOpen dpsiDomain_half]
    rw [iteratedDerivWithin_tsum_dpsiTermC k hk]
    rw [tsum_mul_left]
    rw [hsign]
  rw [hA, hB]
  rw [← mul_add]
  congr 2
  have hsum := summable_zpow_half_add_nat k hk
  conv_rhs => rw [hsum.tsum_eq_zero_add]
  have e0 : ((1 / 2 : ℂ) + ((0 : ℕ) : ℂ)) ^ (-(k : ℤ) - 1) = ((1 : ℂ) / 2) ^ (-(k : ℤ) - 1) := by
    simp
  rw [e0]
  apply add_left_cancel_iff.mpr
  apply tsum_congr
  intro m
  push_cast
  ring

/-- The polygamma values of Complex.digamma at 1/2 (j ≥ 1). -/
lemma iteratedDeriv_digamma_half (j : ℕ) (hj : 1 ≤ j) :
    iteratedDeriv j Complex.digamma ((1 : ℂ) / 2) =
      ((-1 : ℂ)) ^ (j + 1) * (j ! : ℂ) * ∑' m : ℕ, ((1 / 2 : ℂ) + ↑m) ^ (-(j : ℤ) - 1) := by
  have hev : Complex.digamma =ᶠ[𝓝 ((1 : ℂ) / 2)] DinfC := by
    apply Filter.eventually_of_mem (dpsiDomain_isOpen.mem_nhds dpsiDomain_half)
    intro z hz
    exact digamma_eq_DinfC_on hz
  rw [iteratedDeriv_congr_of_eventuallyEq hev]
  exact iteratedDeriv_DinfC_half j hj

/-- The odd-shift equivalence: m ↦ 2m+3 from ℕ to {m odd, m ≥ 3}. -/
def oddGeThreeEquiv : ℕ ≃ {m : ℕ // Odd m ∧ 3 ≤ m} where
  toFun m := ⟨2 * m + 3, by
    constructor
    · exact Nat.odd_iff.mpr (by omega)
    · omega⟩
  invFun := fun ⟨m, _⟩ => (m - 3) / 2
  left_inv m := by simp
  right_inv := fun ⟨m, hm⟩ => by
    obtain ⟨hodd, h3⟩ := hm
    rw [Nat.odd_iff] at hodd
    apply Subtype.ext
    simp only []
    omega

/-- The series over odd m ≥ 3 of m^{-(j+1)} is summable (real, the
lagariasA tail's summand family). -/
lemma summable_odd_tail_inv_pow (j : ℕ) (hj : 1 ≤ j) :
    Summable (fun m : {m : ℕ // Odd m ∧ 3 ≤ m} => (1 / ((m : ℝ) : ℝ)) ^ (j + 1)) := by
  have h := oddSlice_summable (j + 1) (by omega)
  apply h.congr
  intro m
  rw [div_pow, one_pow]

/-- The reindexing: Σ_{m≥0} (2m+1)^{-(j+1)} = lagariasA (j+1) (as ℂ). -/
lemma tsum_two_mul_add_one_inv_pow (j : ℕ) (hj : 1 ≤ j) :
    (∑' m : ℕ, ((2 : ℂ) * ↑m + 1) ^ (-(j : ℤ) - 1)) = (lagariasA (j + 1) : ℝ) := by
  -- peel m = 0 (value 1), shift, then the odd≥3 equiv
  have hsumm : Summable (fun m : ℕ => ((2 : ℂ) * ↑m + 1) ^ (-(j : ℤ) - 1)) := by
    apply Summable.of_norm_bounded (summable_div_add_one_pow (4 ^ (j + 1)) (j + 1) (by omega))
    intro m
    have hpos : (0 : ℝ) < 2 * (m : ℝ) + 1 := by positivity
    have h1 : ‖((2 : ℂ) * ↑m + 1)‖ = 2 * (m : ℝ) + 1 := by
      have e : ((2 : ℂ) * ↑m + 1) = ↑(2 * (m : ℝ) + 1) := by push_cast; ring
      rw [e, Complex.norm_real, Real.norm_eq_abs, abs_of_nonneg hpos.le]
    rw [norm_zpow, h1]
    rw [show (-(j : ℤ) - 1 : ℤ) = -(((j + 1 : ℕ) : ℤ)) by push_cast; ring, zpow_neg, zpow_natCast]
    rw [← one_div]
    rw [div_le_div_iff₀ (by positivity : (0 : ℝ) < (2 * (m : ℝ) + 1) ^ (j + 1))
      (by positivity : (0 : ℝ) < ((m : ℝ) + 1) ^ (j + 1))]
    rw [one_mul]
    rw [show (4 : ℝ) ^ (j + 1) * (2 * (m : ℝ) + 1) ^ (j + 1) =
      (4 * (2 * (m : ℝ) + 1)) ^ (j + 1) by rw [mul_pow]]
    apply pow_le_pow_left₀ (by positivity : (0 : ℝ) ≤ (m : ℝ) + 1)
    linarith [Nat.cast_nonneg (α := ℝ) m]
  rw [hsumm.tsum_eq_zero_add]
  have hf0 : ((2 : ℂ) * ((0 : ℕ) : ℂ) + 1) ^ (-(j : ℤ) - 1) = 1 := by simp
  rw [hf0]
  -- the tail: Σ_{m≥0} (2↑m+3)^{-(j+1)} = Σ_{m' odd ≥ 3} (↑m')^{-(j+1)}
  have heq2 : ∀ m : ℕ, ((2 : ℂ) * ↑(m + 1) + 1) ^ (-(j : ℤ) - 1) =
      ((2 : ℂ) * ↑m + 3) ^ (-(j : ℤ) - 1) := fun m => by push_cast; ring_nf
  rw [tsum_congr heq2]
  have heq3 : (fun m : ℕ => ((2 : ℂ) * ↑m + 3) ^ (-(j : ℤ) - 1)) =
      (fun m' : {m : ℕ // Odd m ∧ 3 ≤ m} => (((m' : ℕ) : ℂ)) ^ (-(j : ℤ) - 1)) ∘ oddGeThreeEquiv := by
    funext m
    rw [Function.comp_apply]
    show ((2 : ℂ) * ↑m + 3) ^ (-(j : ℤ) - 1) = (((oddGeThreeEquiv m : ℕ) : ℂ)) ^ (-(j : ℤ) - 1)
    congr 1
    show (2 : ℂ) * ↑m + 3 = ((2 * m + 3 : ℕ) : ℂ)
    push_cast
    ring
  rw [heq3]
  rw [show (∑' m : ℕ, (((fun m' : {m : ℕ // Odd m ∧ 3 ≤ m} => (((m' : ℕ) : ℂ)) ^ (-(j : ℤ) - 1)) ∘
      oddGeThreeEquiv) m)) =
    ∑' m : ℕ, (fun m' : {m : ℕ // Odd m ∧ 3 ≤ m} => (((m' : ℕ) : ℂ)) ^ (-(j : ℤ) - 1))
      (oddGeThreeEquiv m) from rfl]
  rw [Equiv.tsum_eq oddGeThreeEquiv
    (fun m' : {m : ℕ // Odd m ∧ 3 ≤ m} => (((m' : ℕ) : ℂ)) ^ (-(j : ℤ) - 1))]
  -- now: 1 + Σ_{m' odd ≥ 3} (↑m')^{-(j+1)} = ↑(lagariasA (j+1))
  have hj1 : j + 1 ≠ 1 := by omega
  rw [lagariasA, if_neg hj1]
  push_cast
  congr 1
  apply tsum_congr
  intro m
  have hexp : (-(j : ℤ) - 1 : ℤ) = -(((j + 1 : ℕ) : ℤ)) := by push_cast; ring
  rw [hexp]
  rw [show (((m : ℕ) : ℂ)) ^ (-(((j + 1 : ℕ) : ℤ))) = ((((m : ℕ) : ℂ) ^ (j + 1))⁻¹) from by
    rw [zpow_neg, zpow_natCast]]
  rw [one_div_pow, inv_eq_one_div]

/-- AnalyticAt is preserved by iteratedDeriv. -/
lemma analyticAt_iteratedDeriv {F : ℂ → ℂ} {x : ℂ} (hF : AnalyticAt ℂ F x) (n : ℕ) :
    AnalyticAt ℂ (iteratedDeriv n F) x := by
  induction n with
  | zero => simpa [iteratedDeriv_zero] using hF
  | succ n ih =>
    rw [iteratedDeriv_succ]
    exact ih.deriv

/-- The affine comp: iterated derivatives of F((1+s)/2) at s. -/
lemma iteratedDeriv_one_add_div_two (F : ℂ → ℂ) (j : ℕ) {s : ℂ}
    (hF : AnalyticAt ℂ F ((1 + s) / 2))
    (hFn : ∀ᶠ t in 𝓝 s, AnalyticAt ℂ F ((1 + t) / 2)) :
    iteratedDeriv j (fun t : ℂ => F ((1 + t) / 2)) s =
      (1 / 2 : ℂ) ^ j * iteratedDeriv j F ((1 + s) / 2) := by
  induction j generalizing s with
  | zero => simp [iteratedDeriv_zero]
  | succ j ih =>
    rw [iteratedDeriv_succ, iteratedDeriv_succ]
    have hev : (fun t => iteratedDeriv j (fun u : ℂ => F ((1 + u) / 2)) t) =ᶠ[𝓝 s]
        (fun t => (1 / 2 : ℂ) ^ j * iteratedDeriv j F ((1 + t) / 2)) := by
      filter_upwards [hFn] with t ht
      have hcont : Tendsto (fun u : ℂ => (1 + u) / 2) (𝓝 t) (𝓝 ((1 + t) / 2)) := by
        have hc : ContinuousAt (fun u : ℂ => (1 + u) / 2) t :=
          (continuousAt_const.add continuousAt_id).div continuousAt_const (by norm_num)
        exact hc.tendsto
      exact ih ht (hcont.eventually ht.eventually_analyticAt)
    rw [hev.deriv_eq]
    have hcomp : HasDerivAt (fun t : ℂ => iteratedDeriv j F ((1 + t) / 2))
        (iteratedDeriv (j + 1) F ((1 + s) / 2) * (1 / 2)) s := by
      have hc : HasDerivAt (fun t : ℂ => (1 + t) / 2) (1 / 2) s := by
        have h := ((hasDerivAt_const s (1 : ℂ)).add (hasDerivAt_id s)).div_const (2 : ℂ)
        simpa using h
      have hd : HasDerivAt (iteratedDeriv j F) (iteratedDeriv (j + 1) F ((1 + s) / 2))
          ((1 + s) / 2) := by
        rw [iteratedDeriv_succ]
        exact (analyticAt_iteratedDeriv hF j).differentiableAt.hasDerivAt
      have hc2 := hd.comp s hc
      have heq : (iteratedDeriv j F ∘ fun t : ℂ => (1 + t) / 2) =
          fun t => iteratedDeriv j F ((1 + t) / 2) := rfl
      rw [heq] at hc2
      exact hc2
    rw [deriv_const_mul _ hcomp.differentiableAt, hcomp.deriv]
    rw [pow_succ, ← iteratedDeriv_succ]
    ring

/-- The scaling: (1/2 + m)^{-(j+1)} = 2^{j+1}·(2m+1)^{-(j+1)}. -/
lemma zpow_half_add_nat_scale (j : ℕ) (m : ℕ) :
    ((1 / 2 : ℂ) + ↑m) ^ (-(j : ℤ) - 1) =
      (2 : ℂ) ^ (j + 1) * ((2 : ℂ) * ↑m + 1) ^ (-(j : ℤ) - 1) := by
  have e : ((1 / 2 : ℂ) + ↑m) = ((2 : ℂ) * ↑m + 1) / 2 := by ring
  rw [e, div_zpow]
  rw [show ((2 : ℂ)) ^ (-(j : ℤ) - 1) = ((2 : ℂ) ^ (j + 1))⁻¹ from by
    rw [show (-(j : ℤ) - 1 : ℤ) = -(((j + 1 : ℕ) : ℤ)) by push_cast; ring, zpow_neg, zpow_natCast]]
  rw [div_eq_mul_inv, inv_inv, mul_comm]

/-- The (1/2+m)-series equals 2^{j+1}·lagariasA(j+1) (as ℂ). -/
lemma tsum_half_add_nat_zpow (j : ℕ) (hj : 1 ≤ j) :
    (∑' m : ℕ, ((1 / 2 : ℂ) + ↑m) ^ (-(j : ℤ) - 1)) = (2 : ℂ) ^ (j + 1) * (lagariasA (j + 1) : ℝ) := by
  rw [tsum_congr (fun m => zpow_half_add_nat_scale j m)]
  rw [tsum_mul_left, tsum_two_mul_add_one_inv_pow j hj]

/-- **The archimedean Taylor coefficients of digamma at 1/2** (j ≥ 1):
iteratedDeriv j digamma (1/2) = (−1)^{j+1}·j!·2^{j+1}·lagariasA(j+1). -/
lemma iteratedDeriv_digamma_half_value (j : ℕ) (hj : 1 ≤ j) :
    iteratedDeriv j Complex.digamma ((1 : ℂ) / 2) =
      ((-1 : ℂ)) ^ (j + 1) * (j ! : ℂ) * ((2 : ℂ) ^ (j + 1) * (lagariasA (j + 1) : ℝ)) := by
  rw [iteratedDeriv_digamma_half j hj, tsum_half_add_nat_zpow j hj]

end NguyenRH1