Explicit second-order Euler approximants and saddle models

import Mathlib.Analysis.SpecialFunctions.Trigonometric.Basic
import Mathlib.NumberTheory.Harmonic.EulerMascheroni
import Mathlib.Data.Nat.Choose.Basic

open scoped BigOperators

namespace EulerMascheroni.P2

/-- The published p=2 binomial summand, over the rationals. -/
def coefficient (n k : ℕ) : ℚ :=
  (n.choose k : ℚ)^2 * ((n+k).choose k : ℚ)^2 / (k.factorial : ℚ)

def Q (n : ℕ) : ℚ := ∑ k ∈ Finset.range (n+1), coefficient n k

def P (n : ℕ) : ℚ := - ∑ k ∈ Finset.range (n+1),
  coefficient n k * (2 * harmonic (n+k) + 2 * harmonic (n-k) - 5 * harmonic k)

noncomputable def F (n : ℕ) : ℝ :=
  (Q n : ℝ) * Real.eulerMascheroniConstant - (P n : ℝ)

noncomputable def scale (n : ℕ) : ℝ := (n : ℝ) ^ (4 / 5 : ℝ)
noncomputable def subscale (n : ℕ) : ℝ := (n : ℝ) ^ (2 / 5 : ℝ)
noncomputable def theta : ℝ := 2 * Real.pi / 5
noncomputable def rate : ℝ := 5 * (1 - Real.cos theta)

noncomputable def phase (n : ℕ) : ℝ :=
  5 * scale n * Real.sin theta - (2/3 : ℝ) * subscale n * Real.sin (3*theta)
    - 2*theta

noncomputable def qModel (n : ℕ) : ℝ :=
  Real.exp (5 * scale n - (2/3 : ℝ) * subscale n + 1/5) /
    (4 * Real.pi^2 * Real.sqrt 5 * scale n^2)

noncomputable def fModel (n : ℕ) : ℝ :=
  Real.exp (5 * scale n * Real.cos theta -
    (2/3 : ℝ) * subscale n * Real.cos (3*theta) + 1/5) /
    (Real.pi * Real.sqrt 5 * scale n^2)

/-- The two analytic assertions required by the relative-error reduction.
They are a proposition, not axioms or claimed theorems. Indices start at 1. -/
def SaddleLimits : Prop :=
  Filter.Tendsto (fun n : ℕ => (Q (n+1) : ℝ) / qModel (n+1) - 1)
    Filter.atTop (nhds 0) ∧
  Filter.Tendsto (fun n : ℕ => F (n+1) / fModel (n+1) - Real.sin (phase (n+1)))
    Filter.atTop (nhds 0)

/-- A zero-safe formulation of the predicted sharp exponential rate.
Unlike a real-valued log convention, this does not assign log(0)=0. -/
def SharpRate : Prop := ∀ ε : ℝ, 0 < ε →
  (∀ᶠ n : ℕ in Filter.atTop,
    |Real.eulerMascheroniConstant - (P (n+1) : ℝ) / (Q (n+1) : ℝ)| ≤
      Real.exp ((-rate + ε) * scale (n+1))) ∧
  (∃ᶠ n : ℕ in Filter.atTop,
    Real.exp ((-rate - ε) * scale (n+1)) ≤
      |Real.eulerMascheroniConstant - (P (n+1) : ℝ) / (Q (n+1) : ℝ)|)

end EulerMascheroni.P2