LiArchimedean.lean

Declarations

Source

/-
NguyenRH1/LiArchimedean.lean — tier ii′ sub-stage ii′-3a, piece (A)
(D-066 §4): the archimedean coefficients of the (4.10) split.

Target (D-066): the Taylor coefficients of s ↦ logDeriv Gammaℝ (1+s) at 0:

  (1/j!)·iteratedDeriv j (fun s => logDeriv Gammaℝ (1+s)) 0 = (−1)^{j+1} a_{j+1},

i.e. `iteratedDeriv j ... 0 = (−1)^{j+1}·j!·lagariasA (j+1)`.

Part 1 (this file, complete): the Gammaℝ logDeriv split
  logDeriv Gammaℝ (1+s) = −(log π)/2 + (1/2)·digamma((1+s)/2)
(`logDeriv_Gammaℝ_one_add`) and the j = 0 evaluation
  logDeriv Gammaℝ 1 = −(log 4π + γ)/2 = −lagariasA1
(`logDeriv_Gammaℝ_one`) — via Mathlib's `Complex.digamma_one_half`.

Part 2 (next): the digamma series from BohrMollerup's logGammaSeq and the
j ≥ 1 polygamma values at 1/2, then the odd-m reindexing to lagariasA.

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

open Zeta23 Complex Filter Topology
open scoped Nat

noncomputable section

namespace NguyenRH1

/-! ## 1. The Gammaℝ logDeriv split -/

/-- The cpow factor as an exponential: π^(−w/2) = exp(−(log π)/2·w). -/
lemma pi_cpow_neg_half_eq_exp (w : ℂ) :
    (↑Real.pi : ℂ) ^ (-w / 2) = Complex.exp (-(Real.log Real.pi : ℂ) / 2 * w) := by
  rw [Complex.cpow_def_of_ne_zero (by exact_mod_cast Real.pi_ne_zero)]
  congr 1
  rw [← Complex.ofReal_log Real.pi_nonneg]
  ring_nf

/-- The Gammaℝ logDeriv split: for (1+s)/2 off the poles of Γ,
  logDeriv Gammaℝ (1+s) = −(log π)/2 + (1/2)·ψ((1+s)/2). -/
lemma logDeriv_Gammaℝ_one_add {s : ℂ}
    (hs : Complex.Gamma ((1 + s) / 2) ≠ 0)
    (hs2 : ∀ m : ℕ, ((1 + s) / 2) ≠ -m) :
    logDeriv Complex.Gammaℝ (1 + s) =
      -(Real.log Real.pi : ℂ) / 2 + (1 / 2) * Complex.digamma ((1 + s) / 2) := by
  have hfun : Complex.Gammaℝ = fun w : ℂ =>
      Complex.exp (-(Real.log Real.pi : ℂ) / 2 * w) * Complex.Gamma (w / 2) := by
    funext w
    show Complex.Gammaℝ w = _
    rw [Complex.Gammaℝ_def]
    congr 1
    exact pi_cpow_neg_half_eq_exp w
  rw [hfun]
  have hdif1 : DifferentiableAt ℂ (fun w : ℂ => Complex.exp (-(Real.log Real.pi : ℂ) / 2 * w))
      (1 + s) := by
    have h1 : DifferentiableAt ℂ (fun w : ℂ => -(Real.log Real.pi : ℂ) / 2 * w) (1 + s) :=
      differentiableAt_id.const_mul _
    exact Complex.differentiableAt_exp.comp (1 + s) h1
  have hdif2 : DifferentiableAt ℂ (fun w : ℂ => Complex.Gamma (w / 2)) (1 + s) := by
    have h2 : DifferentiableAt ℂ (fun w : ℂ => w / 2) (1 + s) :=
      differentiableAt_id.div_const 2
    exact (Complex.differentiableAt_Gamma _ hs2).comp (1 + s) h2
  rw [logDeriv_mul (f := fun w : ℂ => Complex.exp (-(Real.log Real.pi : ℂ) / 2 * w))
    (g := fun w => Complex.Gamma (w / 2)) (1 + s) (Complex.exp_ne_zero _) hs hdif1 hdif2]
  have h1 : logDeriv (fun w : ℂ => Complex.exp (-(Real.log Real.pi : ℂ) / 2 * w)) (1 + s) =
      -(Real.log Real.pi : ℂ) / 2 := by
    rw [logDeriv_apply]
    have hd : deriv (fun w : ℂ => Complex.exp (-(Real.log Real.pi : ℂ) / 2 * w)) (1 + s) =
        Complex.exp (-(Real.log Real.pi : ℂ) / 2 * (1 + s)) * (-(Real.log Real.pi : ℂ) / 2) := by
      have hh1 : HasDerivAt (fun w : ℂ => -(Real.log Real.pi : ℂ) / 2 * w)
          (-(Real.log Real.pi : ℂ) / 2) (1 + s) := by
        simpa using (hasDerivAt_id (1 + s)).const_mul (-(Real.log Real.pi : ℂ) / 2)
      have hh2 := (Complex.hasDerivAt_exp (-(Real.log Real.pi : ℂ) / 2 * (1 + s))).comp (1 + s) hh1
      have heq : (Complex.exp ∘ fun w : ℂ => -(Real.log Real.pi : ℂ) / 2 * w) =
          fun w : ℂ => Complex.exp (-(Real.log Real.pi : ℂ) / 2 * w) := rfl
      rw [heq] at hh2
      exact hh2.deriv
    rw [hd]
    field_simp [Complex.exp_ne_zero (-(Real.log Real.pi : ℂ) / 2 * (1 + s))]
  have h2 : logDeriv (fun w : ℂ => Complex.Gamma (w / 2)) (1 + s) =
      (1 / 2) * Complex.digamma ((1 + s) / 2) := by
    rw [logDeriv_apply]
    have hd : deriv (fun w : ℂ => Complex.Gamma (w / 2)) (1 + s) =
        deriv Complex.Gamma ((1 + s) / 2) * (1 / 2) := by
      have hh1 : HasDerivAt (fun w : ℂ => w / 2) (1 / 2) (1 + s) := by
        simpa using (hasDerivAt_id (1 + s)).div_const (2 : ℂ)
      have hh2 := ((Complex.differentiableAt_Gamma _ hs2).hasDerivAt).comp (1 + s) hh1
      have heq : (Complex.Gamma ∘ fun w : ℂ => w / 2) = fun w : ℂ => Complex.Gamma (w / 2) := rfl
      rw [heq] at hh2
      exact hh2.deriv
    rw [hd, Complex.digamma_def, logDeriv_apply]
    field_simp [hs]
  rw [h1, h2]

/-! ## 2. The j = 0 evaluation: τ₀ = −a₁ -/

/-- log(4π) = 2·log 2 + log π. -/
lemma log_four_pi : Real.log (4 * Real.pi) = 2 * Real.log 2 + Real.log Real.pi := by
  rw [Real.log_mul (by norm_num) (by positivity)]
  congr 1
  rw [show (4 : ℝ) = 2 ^ 2 from by norm_num, Real.log_pow]
  ring

/-- **τ₀** (D-066 §4, j = 0): logDeriv Gammaℝ 1 = −(log 4π + γ)/2 = −a₁. -/
theorem logDeriv_Gammaℝ_one :
    logDeriv Complex.Gammaℝ 1 = -(lagariasA1 : ℂ) := by
  have hs : Complex.Gamma ((1 + (0 : ℂ)) / 2) ≠ 0 := by
    have e : (1 + (0 : ℂ)) / 2 = 1 / 2 := by ring
    rw [e, Complex.Gamma_one_half_eq]
    have hpi : (↑Real.pi : ℂ) ≠ 0 := by exact_mod_cast Real.pi_ne_zero
    rw [Complex.cpow_def_of_ne_zero hpi]
    exact Complex.exp_ne_zero _
  have hs2 : ∀ m : ℕ, ((1 + (0 : ℂ)) / 2) ≠ -m := by
    intro m hm
    have e : (1 + (0 : ℂ)) / 2 = 1 / 2 := by ring
    rw [e] at hm
    have hre : ((1 : ℂ) / 2).re = (-(m : ℂ)).re := by rw [hm]
    simp at hre
    have hm0 : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
    linarith
  have h := logDeriv_Gammaℝ_one_add (s := 0) hs hs2
  rw [add_zero] at h
  rw [h]
  rw [Complex.digamma_one_half]
  rw [show Complex.log (2 : ℂ) = (Real.log 2 : ℂ) by
    rw [show (2 : ℂ) = ((2 : ℝ) : ℂ) by simp]
    exact (Complex.ofReal_log (by norm_num : (0 : ℝ) ≤ 2)).symm]
  rw [lagariasA1, log_four_pi]
  push_cast
  ring

/-! ## 3. The full archimedean coefficient evaluation (all j) -/

/-- The split's hypotheses hold on a neighbourhood of 0. -/
lemma eventually_Gamma_ne_zero_and_pole_free :
    (∀ᶠ s in 𝓝 (0 : ℂ), Complex.Gamma ((1 + s) / 2) ≠ 0) ∧
      (∀ᶠ s in 𝓝 (0 : ℂ), ∀ m : ℕ, (1 + s) / 2 ≠ -(m : ℂ)) := by
  have hcont1 : ContinuousAt (fun s : ℂ => Complex.Gamma ((1 + s) / 2)) 0 := by
    have haffine : ContinuousAt (fun s : ℂ => (1 + s) / 2) 0 := by
      have h := (continuousAt_const.add continuousAt_id : ContinuousAt (fun s : ℂ => 1 + s) 0)
      simpa using h.div_const (2 : ℂ)
    have hG : ContinuousAt Complex.Gamma ((1 : ℂ) / 2) := by
      apply Complex.continuousAt_Gamma
      intro m hcontra
      have h2 : ((1 : ℂ) / 2).re = (-(m : ℂ)).re := by rw [hcontra]
      simp at h2
      have hm : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
      linarith
    have hcomp := hG.comp_of_eq haffine (by show (1 + (0 : ℂ)) / 2 = 1 / 2; ring)
    have heq : (Complex.Gamma ∘ fun s : ℂ => (1 + s) / 2) =
        fun s => Complex.Gamma ((1 + s) / 2) := rfl
    rw [heq] at hcomp
    exact hcomp
  have h1 : ∀ᶠ s in 𝓝 (0 : ℂ), Complex.Gamma ((1 + s) / 2) ≠ 0 := by
    have hG1 : Complex.Gamma ((1 + (0 : ℂ)) / 2) ≠ 0 := by
      have e : (1 + (0 : ℂ)) / 2 = 1 / 2 := by ring
      rw [e, Complex.Gamma_one_half_eq]
      have hpi : (↑Real.pi : ℂ) ≠ 0 := by exact_mod_cast Real.pi_ne_zero
      rw [Complex.cpow_def_of_ne_zero hpi]
      exact Complex.exp_ne_zero _
    exact hcont1.tendsto.eventually_ne hG1
  have h2 : ∀ᶠ s in 𝓝 (0 : ℂ), ∀ m : ℕ, (1 + s) / 2 ≠ -(m : ℂ) := by
    have hopen : IsOpen {s : ℂ | 0 < (((1 + s) / 2).re)} := by
      have hcont : Continuous (fun s : ℂ => (((1 + s) / 2).re)) := by
        have h1 : Continuous (fun s : ℂ => (1 + s) / 2) :=
          (continuous_const.add continuous_id).div_const (2 : ℂ)
        exact Complex.continuous_re.comp h1
      exact isOpen_Ioi.preimage hcont
    have hmem : (0 : ℂ) ∈ {s : ℂ | 0 < (((1 + s) / 2).re)} := by
      show 0 < (((1 + (0 : ℂ)) / 2).re)
      have e : (1 + (0 : ℂ)) / 2 = 1 / 2 := by ring
      rw [e]
      simp
    filter_upwards [hopen.mem_nhds hmem] with s hs
    intro m hmeq
    rw [hmeq] at hs
    simp at hs
    have hm : (0 : ℝ) ≤ (m : ℝ) := Nat.cast_nonneg m
    linarith
  exact ⟨h1, h2⟩

/-- **Piece (A), the archimedean Taylor coefficients (D-066 §4)**: for all j,
  iteratedDeriv j (fun s => logDeriv Gammaℝ (1+s)) 0 = (−1)^{j+1}·j!·↑lagariasA(j+1). -/
theorem iteratedDeriv_logDeriv_Gammaℝ (j : ℕ) :
    iteratedDeriv j (fun s : ℂ => logDeriv Complex.Gammaℝ (1 + s)) 0 =
      ((-1 : ℂ)) ^ (j + 1) * (j ! : ℂ) * (lagariasA (j + 1) : ℝ) := by
  rcases Nat.eq_zero_or_pos j with hj0 | hj0
  · subst hj0
    rw [iteratedDeriv_zero]
    rw [lagariasA, if_pos rfl]
    simp [logDeriv_Gammaℝ_one]
  · have hj : 1 ≤ j := hj0
    obtain ⟨hG, hpole⟩ := eventually_Gamma_ne_zero_and_pole_free
    have hev : (fun s : ℂ => logDeriv Complex.Gammaℝ (1 + s)) =ᶠ[𝓝 (0 : ℂ)]
        (fun s => -(Real.log Real.pi : ℂ) / 2 + (1 / 2) * Complex.digamma ((1 + s) / 2)) := by
      filter_upwards [hG, hpole] with s h1 h2
      exact logDeriv_Gammaℝ_one_add h1 h2
    rw [iteratedDeriv_congr_of_eventuallyEq hev]
    -- the constant dies for j ≥ 1
    have hsplit2 : (fun s : ℂ => -(Real.log Real.pi : ℂ) / 2 +
        (1 / 2) * Complex.digamma ((1 + s) / 2)) =
        (fun _ : ℂ => -(Real.log Real.pi : ℂ) / 2) +
          (fun s => (1 / 2) * Complex.digamma ((1 + s) / 2)) := rfl
    rw [hsplit2]
    have hcd0 : ContDiffAt ℂ j (fun s => Complex.digamma ((1 + s) / 2)) 0 := by
      have h1 : AnalyticAt ℂ Complex.digamma ((1 : ℂ) / 2) :=
        digamma_analyticOn _ dpsiDomain_half
      have h2 : AnalyticAt ℂ (fun s : ℂ => (1 + s) / 2) 0 :=
        (analyticAt_const.add analyticAt_id).div analyticAt_const (by norm_num)
      have h3 := AnalyticAt.comp_of_eq h1 h2 (by show (1 + (0 : ℂ)) / 2 = 1 / 2; ring)
      have heq : (Complex.digamma ∘ (fun s : ℂ => (1 + s) / 2)) =
          fun s => Complex.digamma ((1 + s) / 2) := rfl
      rw [heq] at h3
      exact h3.contDiffAt
    have hcd : ContDiffAt ℂ j (fun s => (1 / 2) * Complex.digamma ((1 + s) / 2)) 0 :=
      (contDiff_const (𝕜 := ℂ) (c := (1 / 2 : ℂ))).contDiffAt.mul hcd0
    rw [iteratedDeriv_add contDiff_const.contDiffAt hcd]
    rw [iteratedDeriv_const, if_neg (by omega : j ≠ 0), zero_add]
    rw [iteratedDeriv_const_mul (1 / 2 : ℂ) hcd0]
    rw [iteratedDeriv_one_add_div_two Complex.digamma j (by
      have e : (1 + (0 : ℂ)) / 2 = 1 / 2 := by ring
      rw [e]
      exact digamma_analyticOn _ dpsiDomain_half) (by
      have hcont : Tendsto (fun s : ℂ => (1 + s) / 2) (𝓝 0) (𝓝 ((1 + (0 : ℂ)) / 2)) := by
        have hc : ContinuousAt (fun s : ℂ => (1 + s) / 2) 0 := by
          have h := (continuousAt_const.add continuousAt_id : ContinuousAt (fun s : ℂ => 1 + s) 0)
          simpa using h.div_const (2 : ℂ)
        exact hc.tendsto
      have h1 : ((1 + (0 : ℂ)) / 2) = 1 / 2 := by ring
      rw [h1] at hcont
      exact hcont.eventually ((digamma_analyticOn _ dpsiDomain_half).eventually_analyticAt))]
    rw [show ((1 : ℂ) + 0) / 2 = 1 / 2 from by ring]
    rw [iteratedDeriv_digamma_half_value j hj]
    have h2pow : ((1 : ℂ) / 2) * ((1 : ℂ) / 2) ^ j * (2 : ℂ) ^ (j + 1) = 1 := by
      rw [show ((1 : ℂ) / 2) * ((1 : ℂ) / 2) ^ j = ((1 : ℂ) / 2) ^ (j + 1) from by
        rw [pow_succ]; ring]
      rw [← mul_pow]
      simp
    have hkey : (1 / 2 : ℂ) * ((1 / 2 : ℂ) ^ j *
        (((-1 : ℂ)) ^ (j + 1) * (j ! : ℂ) * ((2 : ℂ) ^ (j + 1) * (lagariasA (j + 1) : ℝ)))) =
        ((-1 : ℂ)) ^ (j + 1) * (j ! : ℂ) * (lagariasA (j + 1) : ℝ) := by
      rw [show (1 / 2 : ℂ) * ((1 / 2 : ℂ) ^ j * (((-1 : ℂ)) ^ (j + 1) * (j ! : ℂ) *
          ((2 : ℂ) ^ (j + 1) * (lagariasA (j + 1) : ℝ)))) =
        (((1 : ℂ) / 2) * ((1 : ℂ) / 2) ^ j * (2 : ℂ) ^ (j + 1)) *
          (((-1 : ℂ)) ^ (j + 1) * (j ! : ℂ) * (lagariasA (j + 1) : ℝ)) by ring]
      rw [h2pow]
      ring
    exact hkey