Mathlib.Analysis.Special­Functions.Jacobi

theorem jacobiShifted_zero (k : ℕ) (α β : ℝ) : jacobiShifted k α β 0 = 1 := by unfold jacobiShifted calc _ = (-1 : ℝ) ^ 0 * (k.choose 0 : ℝ) * (rising ((k : ℝ) + α + β + 1) 0 / rising (α + 1) 0) * (0 : ℝ) ^ 0 := by apply Finset.sum_eq_single · intro b hb hne simp [hne] · simp _ = 1 := by simp [rising]

theorem jacobiShifted_coefficient_recurrence (k m : ℕ) (α β : ℝ) (hα : -1 < α) (hm : m < k) : ((m : ℝ) + 1) * (α + (m : ℝ) + 1) * ((-1 : ℝ) ^ (m + 1) * (k.choose (m + 1) : ℝ) * (rising ((k : ℝ) + α + β + 1) (m + 1) / rising (α + 1) (m + 1))) = -((k : ℝ) - (m : ℝ)) * ((k : ℝ) + α + β + (m : ℝ) + 1) * ((-1 : ℝ) ^ m * (k.choose m : ℝ) * (rising ((k : ℝ) + α + β + 1) m / rising (α + 1) m)) := by have hden : rising (α + 1) m ≠ 0 := ne_of_gt (rising_pos (α + 1) m (by linarith)) have hstep : α + (m : ℝ) + 1 ≠ 0 := by have hmnonneg : (0 : ℝ) ≤ m := Nat.cast_nonneg m linarith have hchoose : (k.choose (m + 1) : ℝ) * ((m : ℝ) + 1) = (k.choose m : ℝ) * ((k : ℝ) - (m : ℝ)) := by have hc := congrArg (fun n : ℕ => (n : ℝ)) (Nat.choose_succ_right_eq k m) push_cast at hc rw [Nat.cast_sub hm.le] at hc exact hc simp only [rising, ascPochhammer_succ_eval, pow_succ] change _ * _ * (_ * _ * (((ascPochhammer ℝ m).eval ((k : ℝ) + α + β + 1) * (((k : ℝ) + α + β + 1) + (m : ℝ))) / ((ascPochhammer ℝ m).eval (α + 1) * ((α + 1) + (m : ℝ))))) = _ rw [show (α + 1) + (m : ℝ) = α + (m : ℝ) + 1 by ring] change (ascPochhammer ℝ m).eval (α + 1) ≠ 0 at hden field_simp [hden, hstep] linear_combination -((ascPochhammer ℝ m).eval ((k : ℝ) + α + β + 1)) * ((k : ℝ) + α + β + (m : ℝ) + 1) * hchoose

theorem jacobiShifted_ode (k : ℕ) (α β x : ℝ) (hα : -1 < α) : x * (1 - x) * deriv (deriv (jacobiShifted k α β)) x + (α + 1 - (α + β + 2) * x) * deriv (jacobiShifted k α β) x + (k : ℝ) * ((k : ℝ) + α + β + 1) * jacobiShifted k α β x = 0 := by let c : ℕ → ℝ := fun m => (-1 : ℝ) ^ m * (k.choose m : ℝ) * (rising ((k : ℝ) + α + β + 1) m / rising (α + 1) m) let p : Polynomial ℝ := ∑ m ∈ range (k + 1), C (c m) * X ^ m have hp (m : ℕ) : p.coeff m = if m ≤ k then c m else 0 := by simp [p, Finset.mem_range] have hc (m : ℕ) (hm : m < k) : ((m : ℝ) + 1) * (α + (m : ℝ) + 1) * c (m + 1) = -((k : ℝ) - (m : ℝ)) * ((k : ℝ) + α + β + (m : ℝ) + 1) * c m := by exact jacobiShifted_coefficient_recurrence k m α β hα hm have hcoeff (m : ℕ) (K : ℝ) : (X * (1 - X) * derivative (derivative p) + (C (α + 1) - C (α + β + 2) * X) * derivative p + C K * p).coeff m = ((m : ℝ) + 1) * (α + (m : ℝ) + 1) * p.coeff (m + 1) + (K - (m : ℝ) * ((m : ℝ) + α + β + 1)) * p.coeff m := by have hform : X * (1 - X) * derivative (derivative p) + (C (α + 1) - C (α + β + 2) * X) * derivative p + C K * p = X * derivative (derivative p) - X * (X * derivative (derivative p)) + C (α + 1) * derivative p - X * (C (α + β + 2) * derivative p) + C K * p := by ring rw [hform] cases m with | zero => simp only [coeff_add, coeff_sub, coeff_X_mul_zero, coeff_C_mul, coeff_derivative, Nat.cast_zero, zero_add, sub_zero] ring | succ m => cases m with | zero => simp only [coeff_add, coeff_sub, coeff_X_mul_zero, coeff_X_mul, coeff_C_mul, coeff_derivative, Nat.cast_succ, Nat.cast_zero] ring | succ m => simp only [coeff_add, coeff_sub, coeff_X_mul, coeff_C_mul, coeff_derivative, Nat.cast_succ] ring have hpoly : X * (1 - X) * derivative (derivative p) + (C (α + 1) - C (α + β + 2) * X) * derivative p + C ((k : ℝ) * ((k : ℝ) + α + β + 1)) * p = 0 := by ext m rw [hcoeff, coeff_zero, hp, hp] by_cases hm : m < k · rw [if_pos hm.le, if_pos (Nat.succ_le_iff.mpr hm)] linear_combination hc m hm · by_cases heq : m = k · subst m simp only [if_pos le_rfl, if_neg (Nat.not_succ_le_self k), mul_zero, zero_add] ring · have hk : k < m := by omega simp only [if_neg (by omega : ¬ m + 1 ≤ k), if_neg (by omega : ¬ m ≤ k), mul_zero, add_zero] have heval (y : ℝ) : p.eval y = jacobiShifted k α β y := by simp [p, jacobiShifted, c, eval_finsetSum, eval_mul, eval_pow] have hfun : jacobiShifted k α β = fun y => p.eval y := funext fun y => (heval y).symm rw [hfun] have hderivfun : deriv (fun y => p.eval y) = fun y => (derivative p).eval y := funext fun y => Polynomial.deriv p rw [hderivfun] simp only [Polynomial.deriv] have := congrArg (eval x) hpoly simpa [eval_add, eval_mul, eval_sub] using this

theorem jacobiShifted_one (k : ℕ) (α β : ℝ) (hα : -1 < α) : jacobiShifted k α β 1 = (-1 : ℝ) ^ k * rising (β + 1) k / rising (α + 1) k := by have hneg (a : ℝ) (n : ℕ) : (descPochhammer ℝ n).eval (-a) = (-1 : ℝ) ^ n * rising a n := by induction n with | zero => simp [rising] | succ n ih => rw [descPochhammer_succ_eval] simp only [rising, ascPochhammer_succ_eval] at * rw [ih] ring have hcast (n : ℕ) (x : ℝ) : (descPochhammer ℤ n).smeval x = (descPochhammer ℝ n).eval x := by rw [descPochhammer_smeval_eq_ascPochhammer, ascPochhammer_smeval_eq_eval, descPochhammer_eval_eq_ascPochhammer] let A : ℝ := (k : ℝ) + α + β + 1 let C : ℝ := α + 1 have hC : 0 < C := by dsimp [C]; linarith have hden : rising C k ≠ 0 := ne_of_gt (rising_pos C k hC) have hmul (m : ℕ) (hm : m ≤ k) : rising C m * rising (C + m) (k - m) = rising C k := by have h := congrArg (eval C) (ascPochhammer_mul ℝ m (k - m)) simp only [eval_mul, eval_comp, eval_add, eval_X, eval_natCast] at h have heq : m + (k - m) = k := Nat.add_sub_of_le hm simpa only [rising, heq] using h have hv : (descPochhammer ℝ k).eval (-A + (C + (k : ℝ) - 1)) = ∑ m ∈ range (k + 1), (k.choose m : ℝ) * ((descPochhammer ℝ m).eval (-A) * (descPochhammer ℝ (k - m)).eval (C + (k : ℝ) - 1)) := by have h := Ring.descPochhammer_smeval_add k (Commute.all (-A) (C + (k : ℝ) - 1)) simp_rw [hcast] at h simpa only [Finset.Nat.sum_antidiagonal_eq_sum_range_succ_mk, Nat.succ_eq_add_one] using h have hv' : (-1 : ℝ) ^ k * rising (β + 1) k = ∑ m ∈ range (k + 1), (k.choose m : ℝ) * ((-1 : ℝ) ^ m * rising A m * rising (C + m) (k - m)) := by convert hv using 1 · have heq : -A + (C + (k : ℝ) - 1) = -(β + 1) := by dsimp [A, C]; ring rw [heq, hneg] · apply Finset.sum_congr rfl intro m hm have hm' : m ≤ k := by have := Finset.mem_range.mp hm omega rw [hneg] have heq : C + (k : ℝ) - 1 - (k - m : ℕ) + 1 = C + m := by rw [Nat.cast_sub hm'] ring rw [descPochhammer_eval_eq_ascPochhammer, heq] dsimp [rising] apply (eq_div_iff hden).2 unfold jacobiShifted simp only [one_pow, mul_one, Finset.sum_mul] rw [hv'] apply Finset.sum_congr rfl intro m hm have hm' : m ≤ k := by have := Finset.mem_range.mp hm omega rw [← hmul m hm'] have hmden : rising C m ≠ 0 := ne_of_gt (rising_pos C m hC) have hA : (k : ℝ) + α + β + 1 = A := rfl have hC' : α + 1 = C := rfl rw [hA, hC'] field_simp [hmden]

theorem jacobiShifted_abs_le_one (k : ℕ) (α β x : ℝ) (hβα : β ≤ α) (hβ : -1 < β) (hα : -(1 / 2 : ℝ) ≤ α) (hx : x ∈ Icc (0 : ℝ) 1) : |jacobiShifted k α β x| ≤ 1 := by have hα' : -1 < α := by linarith have hmono (n : ℕ) (a b : ℝ) (ha : 0 < a) (hab : a ≤ b) : rising a n ≤ rising b n := by induction n with | zero => simp [rising] | succ n ih => rw [rising, ascPochhammer_succ_eval, rising, ascPochhammer_succ_eval] have hn : (0 : ℝ) ≤ n := Nat.cast_nonneg n have hp : 0 ≤ rising a n := (rising_pos a n ha).le calc rising a n * (a + n) ≤ rising b n * (a + n) := mul_le_mul_of_nonneg_right ih (by linarith) _ ≤ rising b n * (b + n) := mul_le_mul_of_nonneg_left (by linarith) (rising_pos b n (by linarith)).le have hright : |jacobiShifted k α β 1| ≤ 1 := by rw [jacobiShifted_one k α β hα', abs_div, abs_mul, abs_pow, abs_neg, abs_one, one_pow, one_mul] have hpβ := rising_pos (β + 1) k (by linarith) have hpα := rising_pos (α + 1) k (by linarith) rw [abs_of_pos hpβ, abs_of_pos hpα] exact (div_le_iff₀ hpα).2 (by simpa using hmono k (β + 1) (α + 1) (by linarith) (by linarith)) have hleft : |jacobiShifted k α β 0| ≤ 1 := by rw [jacobiShifted_zero] norm_num by_cases hk : k = 0 · subst k simp [jacobiShifted, rising] have hkpos : 0 < k := Nat.pos_of_ne_zero hk let f : ℝ → ℝ := jacobiShifted k α β let L : ℝ := (k : ℝ) * ((k : ℝ) + α + β + 1) have hL : 0 < L := by dsimp [L] have hk' : (0 : ℝ) < k := Nat.cast_pos.mpr hkpos have : 0 < (k : ℝ) + α + β + 1 := by have : (1 : ℝ) ≤ k := by exact_mod_cast hkpos linarith positivity let E : ℝ → ℝ := fun y => f y ^ 2 + y * (1 - y) * (deriv f y) ^ 2 / L let p : Polynomial ℝ := ∑ m ∈ range (k + 1), C ((-1 : ℝ) ^ m * (k.choose m : ℝ) * (rising ((k : ℝ) + α + β + 1) m / rising (α + 1) m)) * X ^ m have hp : f = fun y => p.eval y := by funext y simp [f, p, jacobiShifted, eval_finsetSum, eval_mul, eval_pow] have hf (y : ℝ) : DifferentiableAt ℝ f y := by rw [hp] exact p.differentiableAt have hdp : deriv f = fun y => p.derivative.eval y := by funext y rw [hp] exact Polynomial.deriv p have hdf (y : ℝ) : DifferentiableAt ℝ (deriv f) y := by rw [hdp] exact p.derivative.differentiableAt have hfc : Continuous f := continuous_iff_continuousAt.mpr fun y => (hf y).continuousAt have hdfc : Continuous (deriv f) := continuous_iff_continuousAt.mpr fun y => (hdf y).continuousAt have hE : ContinuousOn E (Icc (0 : ℝ) 1) := by apply Continuous.continuousOn dsimp [E] exact (hfc.pow 2).add (((continuous_id.mul (continuous_const.sub continuous_id)).mul (hdfc.pow 2)).div_const L) have hd (y : ℝ) (hy : y ∈ Ioo (0 : ℝ) 1) : DifferentiableAt ℝ E y := by dsimp [E] fun_prop have hder (y : ℝ) (hy : y ∈ Ioo (0 : ℝ) 1) : deriv E y = (-(2 * α + 1) + 2 * (α + β + 1) * y) * ((deriv f y) ^ 2 / L) := by have ho := jacobiShifted_ode k α β y hα' change y * (1 - y) * deriv (deriv f) y + (α + 1 - (α + β + 2) * y) * deriv f y + L * f y = 0 at ho have hcalc : deriv E y = 2 * f y * deriv f y + ((1 - 2 * y) * (deriv f y) ^ 2 + 2 * y * (1 - y) * deriv f y * deriv (deriv f) y) / L := by have hpow : deriv (fun z => f z ^ 2) y = 2 * f y * deriv f y := by rw [show (fun z => f z ^ 2) = f ^ 2 from rfl, deriv_pow (hf y) 2] ring have hdpow : deriv (fun z => (deriv f z) ^ 2) y = 2 * deriv f y * deriv (deriv f) y := by rw [show (fun z => (deriv f z) ^ 2) = (deriv f) ^ 2 from rfl, deriv_pow (hdf y) 2] ring have hlin : deriv (fun z : ℝ => z * (1 - z)) y = 1 - 2 * y := by rw [show (fun z : ℝ => z * (1 - z)) = (fun z : ℝ => z) * (fun z : ℝ => 1 - z) from rfl, deriv_mul (by fun_prop) (by fun_prop)] convert_to (1 : ℝ) * (1 - y) + y * (0 - 1) = 1 - 2 * y · rw [deriv_id'', deriv_const_sub, deriv_id''] norm_num · ring have hprod : deriv (fun z => z * (1 - z) * (deriv f z) ^ 2) y = (1 - 2 * y) * (deriv f y) ^ 2 + y * (1 - y) * (2 * deriv f y * deriv (deriv f) y) := by rw [show (fun z => z * (1 - z) * (deriv f z) ^ 2) = (fun z : ℝ => z * (1 - z)) * (fun z => (deriv f z) ^ 2) from rfl, deriv_mul (by fun_prop) (by fun_prop)] rw [hlin, hdpow] change deriv ((fun z => f z ^ 2) + (fun z => z * (1 - z) * (deriv f z) ^ 2 / L)) y = _ rw [deriv_add (by fun_prop) (by fun_prop)] rw [hpow, deriv_div_const, hprod] ring rw [hcalc] apply (div_left_inj' hL.ne').mp field_simp linear_combination 2 * deriv f y * ho have hbound := affine_deriv_energy_le_max E (2 * α + 1) (2 * (α + β + 1)) x (by linarith) hE hd (by intro y hy refine ⟨(deriv f y) ^ 2 / L, div_nonneg (sq_nonneg _) hL.le, ?_⟩ exact hder y hy) hx have hEx : (f x) ^ 2 ≤ E x := by dsimp [E] have : 0 ≤ x * (1 - x) := mul_nonneg hx.1 (by linarith [hx.2]) have : 0 ≤ x * (1 - x) * (deriv f x) ^ 2 / L := div_nonneg (mul_nonneg this (sq_nonneg _)) hL.le linarith have hE0 : E 0 = 1 := by simp [E, f, jacobiShifted_zero] have hE1 : E 1 ≤ 1 := by simp only [E, sub_self, mul_zero, zero_mul, zero_div, add_zero] have := (sq_le_sq₀ (abs_nonneg (f 1)) (by norm_num : (0 : ℝ) ≤ 1)).mpr hright simpa [f] using this rw [hE0] at hbound have hfx : f x ^ 2 ≤ 1 := le_trans hEx (le_trans hbound (max_le (le_refl 1) hE1)) have : |f x| ≤ 1 := by nlinarith [sq_nonneg (f x), sq_abs (f x)] exact this

theorem rodKernel_deriv_expansion (k i : ℕ) (α x : ℝ) (hx : 0 < x) : (deriv^[i] (rodKernel k α)) x = ∑ m ∈ range (k + 1), (-1 : ℝ) ^ m * (k.choose m : ℝ) * (descPochhammer ℝ i).eval (α + (k : ℝ) + (m : ℝ)) * x ^ (α + (k : ℝ) + (m : ℝ) - (i : ℝ)) := by let c (m : ℕ) : ℝ := (-1 : ℝ) ^ m * (k.choose m : ℝ) let p (m : ℕ) : ℝ := α + (k : ℝ) + (m : ℝ) have hbin (y : ℝ) : (1 - y) ^ k = ∑ m ∈ range (k + 1), c m * y ^ m := by have hh := add_pow (-y) (1 : ℝ) k simp only [one_pow, mul_one] at hh convert hh using 1 · ring · apply Finset.sum_congr rfl intro m hm dsimp [c] rw [neg_pow y] ring have hfun : EqOn (rodKernel k α) (fun y : ℝ => ∑ m ∈ range (k + 1), c m * y ^ p m) (Ioi 0) := by intro y hy have hy' : 0 < y := hy rw [rodKernel, hbin, Finset.mul_sum] apply Finset.sum_congr rfl intro m hm simp only [p] have hr : y ^ (α + (k : ℝ) + (m : ℝ)) = y ^ (α + (k : ℝ)) * y ^ m := by rw [show α + (k : ℝ) + (m : ℝ) = (α + (k : ℝ)) + (m : ℝ) by ring, Real.rpow_add hy', Real.rpow_natCast] rw [hr] ring have hder := (hfun.iteratedDeriv_of_isOpen isOpen_Ioi i) hx rw [iteratedDeriv_eq_iterate, iteratedDeriv_eq_iterate] at hder rw [hder, ← iteratedDeriv_eq_iterate] rw [iteratedDeriv_fun_sum] · apply Finset.sum_congr rfl intro m hm rw [iteratedDeriv_eq_iterate] simp only [p, c] rw [← iteratedDeriv_eq_iterate, iteratedDeriv_const_mul_field, iteratedDeriv_eq_iterate, Real.iter_deriv_rpow_const] ring · intro m hm exact contDiffAt_const.mul (Real.contDiffAt_rpow_const_of_ne (ne_of_gt hx))

theorem rodKernel_expansion_eq_h (k : ℕ) (α x : ℝ) (hα : 0 < α) (hx : 0 < x) : (∑ m ∈ range (k + 1), (-1 : ℝ) ^ m * (k.choose m : ℝ) * (descPochhammer ℝ k).eval (α + (k : ℝ) + (m : ℝ)) * x ^ (α + (k : ℝ) + (m : ℝ) - (k : ℝ))) = rising (α + 1) k * (x ^ α * h k α x) := by have hmul (n m : ℕ) : rising (α + 1) n * rising (α + 1 + (n : ℝ)) m = rising (α + 1) (n + m) := by have hh := congrArg (eval (α + 1)) (ascPochhammer_mul ℝ n m) simpa only [rising, Polynomial.eval_mul, Polynomial.eval_comp, Polynomial.eval_add, Polynomial.eval_X, Polynomial.eval_natCast] using hh have hcoef (m : ℕ) : (descPochhammer ℝ k).eval (α + (k : ℝ) + (m : ℝ)) = rising (α + 1) k * (rising ((k : ℝ) + α + 1) m / rising (α + 1) m) := by have hmpos : rising (α + 1) m ≠ 0 := ne_of_gt (rising_pos (α + 1) m (by linarith)) have hleft := hmul m k have hright := hmul k m have heq : rising (α + 1 + (m : ℝ)) k * rising (α + 1) m = rising (α + 1) k * rising ((k : ℝ) + α + 1) m := by calc _ = rising (α + 1) (m + k) := by rw [← hleft]; ring _ = rising (α + 1) (k + m) := by rw [Nat.add_comm] _ = _ := by rw [← hright]; congr 1; ring rw [descPochhammer_eval_eq_ascPochhammer] change rising (α + (k : ℝ) + (m : ℝ) - (k : ℝ) + 1) k = _ have hb : α + (k : ℝ) + (m : ℝ) - (k : ℝ) + 1 = α + 1 + (m : ℝ) := by ring rw [hb] rw [← mul_div_assoc] apply (eq_div_iff hmpos).2 exact heq unfold h jacobiShifted rw [Finset.mul_sum, Finset.mul_sum] apply Finset.sum_congr rfl intro m hm rw [hcoef m] have hpow : x ^ (α + (k : ℝ) + (m : ℝ) - (k : ℝ)) = x ^ α * x ^ m := by convert Real.rpow_add hx α (m : ℝ) using 1 · ring · simp rw [hpow] ring

theorem rodKernel_deriv_one (k i : ℕ) (α : ℝ) (hik : i < k) : (deriv^[i] (rodKernel k α)) 1 = 0 := by change (deriv^[i] (fun x : ℝ => x ^ (α + (k : ℝ)) * (1 - x) ^ k)) 1 = 0 apply iter_deriv_mul_one_sub_pow_at_one · exact Real.contDiffAt_rpow_const_of_ne (by norm_num) · exact hik

theorem rodKernel_deriv_tendsto_zero (k i : ℕ) (α : ℝ) (hik : i < k) (hα : 0 < α) : Tendsto (fun x : ℝ => (deriv^[i] (rodKernel k α)) x) (nhdsWithin 0 (Ioi 0)) (nhds 0) := by have hterm (m : ℕ) : Tendsto (fun x : ℝ => (-1 : ℝ) ^ m * (k.choose m : ℝ) * (descPochhammer ℝ i).eval (α + (k : ℝ) + (m : ℝ)) * x ^ (α + (k : ℝ) + (m : ℝ) - (i : ℝ))) (nhdsWithin 0 (Ioi 0)) (nhds 0) := by have hi : (i : ℝ) < (k : ℝ) := by exact_mod_cast hik have hm : 0 ≤ (m : ℝ) := Nat.cast_nonneg m have hp : 0 < α + (k : ℝ) + (m : ℝ) - (i : ℝ) := by linarith have hid : Tendsto (fun x : ℝ => x) (nhdsWithin 0 (Ioi 0)) (nhds 0) := (tendsto_id : Tendsto (fun x : ℝ => x) (nhds 0) (nhds 0)).mono_left inf_le_left have hpow := hid.rpow_const (inr hp.le) simpa [Real.zero_rpow (ne_of_gt hp)] using hpow.const_mul ((-1 : ℝ) ^ m * (k.choose m : ℝ) * (descPochhammer ℝ i).eval (α + (k : ℝ) + (m : ℝ))) have hsum : Tendsto (fun x : ℝ => ∑ m ∈ range (k + 1), (-1 : ℝ) ^ m * (k.choose m : ℝ) * (descPochhammer ℝ i).eval (α + (k : ℝ) + (m : ℝ)) * x ^ (α + (k : ℝ) + (m : ℝ) - (i : ℝ))) (nhdsWithin 0 (Ioi 0)) (nhds 0) := by simpa using (tendsto_finsetSum (range (k + 1)) (fun m _ => hterm m)) apply hsum.congr' filter_upwards [self_mem_nhdsWithin] with x hx exact (rodKernel_deriv_expansion k i α x hx).symm

theorem iter_deriv_mul_one_sub_pow_at_one (g : ℝ → ℝ) (k i : ℕ) (hg : ContDiffAt ℝ ⊤ g 1) (hi : i < k) : (deriv^[i] (fun x : ℝ => g x * (1 - x) ^ k)) 1 = 0 := by rw [← iteratedDeriv_eq_iterate] have hp : ContDiffAt ℝ i (fun x : ℝ => (1 - x) ^ k) 1 := by fun_prop change iteratedDeriv i (g * (fun x : ℝ => (1 - x) ^ k)) 1 = 0 rw [iteratedDeriv_mul (hg.of_le (by simp)) hp] apply Finset.sum_eq_zero intro j hj simp only [Finset.mem_range] at hj have hjle : j ≤ i := by omega have hlt : i - j < k := by omega have hpow : iteratedDeriv (i - j) (fun x : ℝ => (1 - x) ^ k) 1 = 0 := by rw [iteratedDeriv_comp_const_sub (i - j) (fun y : ℝ => y ^ k) 1] simp only [sub_self, iteratedDeriv_pow] have hpos : 0 < k - (i - j) := by omega rw [zero_pow (Nat.ne_of_gt hpos), mul_zero, smul_zero] rw [hpow] ring

theorem integral_mul_iterated_deriv_with_boundary (k : ℕ) (a b : ℝ) (p f : ℝ → ℝ) (hp : ∀ x ∈ Icc a b, ContDiffAt ℝ ⊤ p x) (hf : ∀ x ∈ Icc a b, ContDiffAt ℝ ⊤ f x) (hab : a < b) : (∫ x in a..b, p x * (deriv^[k] f) x) = (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) b * (deriv^[k - 1 - i] f) b - (deriv^[i] p) a * (deriv^[k - 1 - i] f) a)) + (-1 : ℝ) ^ k * (∫ x in a..b, (deriv^[k] p) x * f x) := by induction k generalizing p f with | zero => simp | succ k ih => have hfd : ∀ x ∈ Icc a b, ContDiffAt ℝ ⊤ (deriv f) x := by simpa only [Function.iterate_one] using contDiffAt_iterate_deriv_on a b f hf 1 have hpk := contDiffAt_iterate_deriv_on a b p hp k have hpk' := contDiffAt_iterate_deriv_on a b p hp (k + 1) have hintp : IntervalIntegrable (deriv^[k + 1] p) volume a b := by apply ContinuousOn.intervalIntegrable rw [Set.uIcc_of_le hab.le] intro x hx exact (hpk' x hx).continuousAt.continuousWithinAt have hintf : IntervalIntegrable (deriv f) volume a b := by apply ContinuousOn.intervalIntegrable rw [Set.uIcc_of_le hab.le] intro x hx exact (hfd x hx).continuousAt.continuousWithinAt have hip := intervalIntegral.integral_mul_deriv_eq_deriv_mul (u := deriv^[k] p) (v := f) (u' := deriv^[k + 1] p) (v' := deriv f) (by intro x hx have hx' : x ∈ Icc a b := by simpa only [Set.uIcc_of_le hab.le] using hx simpa only [Function.iterate_succ_apply'] using ((hpk x hx').differentiableAt (by simp)).hasDerivAt) (by intro x hx have hx' : x ∈ Icc a b := by simpa only [Set.uIcc_of_le hab.le] using hx exact ((hf x hx').differentiableAt (by simp)).hasDerivAt) hintp hintf have hih := ih p (deriv f) hp hfd simp_rw [← Function.iterate_succ_apply] at hih have hsum : (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) b * (deriv^[(k - 1 - i).succ] f) b - (deriv^[i] p) a * (deriv^[(k - 1 - i).succ] f) a)) = (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) b * (deriv^[k - i] f) b - (deriv^[i] p) a * (deriv^[k - i] f) a)) := by apply Finset.sum_congr rfl intro i hi have hlt : i < k := Finset.mem_range.mp hi have heq : (k - 1 - i).succ = k - i := by omega rw [heq] rw [hih, hip, hsum, Finset.sum_range_succ] simp only [Nat.add_sub_cancel_right, Nat.sub_self, Function.iterate_zero, id_eq] rw [pow_succ] ring

theorem tendsto_integral_one_div_nat (f : ℝ → ℝ) (hf : IntervalIntegrable f volume 0 1) : Tendsto (fun n : ℕ => ∫ x in (1 / ((n + 1 : ℕ) : ℝ))..1, f x) atTop (nhds (∫ x in (0 : ℝ)..1, f x)) := by have hcont : ContinuousOn (fun x : ℝ => ∫ t in x..(1 : ℝ), f t) (Icc (0 : ℝ) 1) := by have hi : IntegrableOn f (uIcc (0 : ℝ) 1) volume := by simpa only [Set.uIcc_of_le (by norm_num : (0 : ℝ) ≤ 1)] using ((intervalIntegrable_iff_integrableOn_Icc_of_le (by norm_num : (0 : ℝ) ≤ 1)).mp hf) simpa only [Set.uIcc_of_le (by norm_num : (0 : ℝ) ≤ 1)] using (intervalIntegral.continuousOn_primitive_interval_left hi) have hseq : Tendsto (fun n : ℕ => (1 / ((n + 1 : ℕ) : ℝ))) atTop (nhds (0 : ℝ)) := by simpa only [Nat.cast_add, Nat.cast_one] using (tendsto_one_div_add_atTop_nhds_zero_nat (𝕜 := ℝ)) have hmem : ∀ n : ℕ, (1 / ((n + 1 : ℕ) : ℝ)) ∈ Icc (0 : ℝ) 1 := by intro n have hn : (1 : ℝ) ≤ ((n + 1 : ℕ) : ℝ) := by exact_mod_cast Nat.succ_le_succ (zero_le n) constructor · positivity · apply (div_le_iff₀ (by positivity : (0 : ℝ) < ((n + 1 : ℕ) : ℝ))).2 simpa only [one_mul] using hn exact (hcont.continuousWithinAt (by norm_num : (0 : ℝ) ∈ Icc 0 1)).tendsto.comp (tendsto_nhdsWithin_iff.mpr ⟨hseq, Filter.Eventually.of_forall hmem⟩)

lemma jacobiShifted_zero_zero_eq_shiftedLegendre (n : ℕ) (x : ℝ) : jacobiShifted n 0 0 x = aeval x (shiftedLegendre n) := by unfold jacobiShifted rising shiftedLegendre simp only [map_sum, map_mul, map_pow, Polynomial.aeval_C, Polynomial.aeval_X] apply Finset.sum_congr rfl intro m hm rw [show (n : ℝ) + 0 + 0 + 1 = ((n + 1 : ℕ) : ℝ) by norm_num, show (0 : ℝ) + 1 = ((1 : ℕ) : ℝ) by norm_num, ascPochhammer_nat_eq_natCast_ascFactorial, ascPochhammer_nat_eq_natCast_ascFactorial, Nat.ascFactorial_eq_factorial_mul_choose, Nat.one_ascFactorial] norm_cast field_simp rw [Nat.choose_symm_add (a := n) (b := m)] push_cast ring

lemma iter_deriv_aeval_int (p : Polynomial ℤ) (n : ℕ) : deriv^[n] (fun x : ℝ => aeval x p) = fun x : ℝ => aeval x (derivative^[n] p) := by induction n with | zero => rfl | succ n ih => rw [Function.iterate_succ_apply', ih, Function.iterate_succ_apply'] funext x exact Polynomial.deriv_aeval _

lemma rod_boundary_zero (k i : ℕ) (hi : i < k) : (deriv^[i] (fun x : ℝ => x^k * (1-x)^k)) 0 = 0 := by rw [← iteratedDeriv_eq_iterate] have hp : ContDiffAt ℝ i (fun x : ℝ => x^k) 0 := by fun_prop have hq : ContDiffAt ℝ i (fun x : ℝ => (1-x)^k) 0 := by fun_prop change iteratedDeriv i ((fun x : ℝ => x^k) * (fun x : ℝ => (1-x)^k)) 0 = 0 rw [iteratedDeriv_mul hp hq] apply Finset.sum_eq_zero intro j hj have hjle : j ≤ i := by simpa using Finset.mem_range.mp hj have hjk : j < k := lt_of_le_of_lt hjle hi have hpow : iteratedDeriv j (fun x : ℝ => x^k) 0 = 0 := by simp [iteratedDeriv_pow, zero_pow (Nat.sub_ne_zero_of_lt hjk)] rw [hpow] ring

lemma shiftedLegendre_orthogonal_unit_of_lt {j k : ℕ} (hjk : j < k) : (∫ x in (0 : ℝ)..1, aeval x (shiftedLegendre j) * aeval x (shiftedLegendre k)) = 0 := by let p : ℝ → ℝ := fun x => aeval x (shiftedLegendre j) let f : ℝ → ℝ := fun x => x^k * (1-x)^k have hibp := integral_mul_iterated_deriv_eq k 0 1 p f (by intro x hx exact ((shiftedLegendre j).contDiff_aeval ⊤).contDiffAt) (by intro x hx; dsimp [f]; fun_prop) (by intro i hi refine ⟨rod_boundary_zero k i hi, ?_⟩ exact iter_deriv_mul_one_sub_pow_at_one (fun x : ℝ => x^k) k i (by fun_prop) hi) (by norm_num) have hpzero : deriv^[k] p = 0 := by change deriv^[k] (fun x : ℝ => aeval x (shiftedLegendre j)) = 0 rw [iter_deriv_aeval_int, Polynomial.iterate_derivative_eq_zero (by simpa using hjk)] funext x simp rw [hpzero] at hibp simp only [Pi.zero_apply, zero_mul, intervalIntegral.integral_zero, mul_zero] at hibp have hscale : (k.factorial : ℝ) * (∫ x in (0 : ℝ)..1, aeval x (shiftedLegendre j) * aeval x (shiftedLegendre k)) = ∫ x in (0 : ℝ)..1, p x * (deriv^[k] f) x := by calc _ = ∫ x in (0 : ℝ)..1, (k.factorial : ℝ) * (aeval x (shiftedLegendre j) * aeval x (shiftedLegendre k)) := by rw [intervalIntegral.integral_const_mul] _ = _ := by apply intervalIntegral.integral_congr intro x hx change (k.factorial : ℝ) * (aeval x (shiftedLegendre j) * aeval x (shiftedLegendre k)) = aeval x (shiftedLegendre j) * (deriv^[k] (fun y : ℝ => y^k * (1-y)^k)) x rw [← shiftedLegendre_rodrigues_eval k x] ring rw [hibp] at hscale exact (mul_eq_zero.mp hscale).resolve_left (by positivity)

lemma iterate_derivative_shiftedLegendre (k : ℕ) : derivative^[k] (shiftedLegendre k) = C ((k.factorial : ℤ) * ((-1 : ℤ)^k * ((2*k).choose k : ℤ))) := by calc _ = C ((derivative^[k] (shiftedLegendre k)).coeff 0) := Polynomial.eq_C_of_natDegree_le_zero ((Polynomial.natDegree_iterate_derivative _ _).trans (by simp)) _ = _ := by congr 1 rw [Polynomial.coeff_iterate_derivative] simp [Polynomial.coeff_shiftedLegendre, Nat.descFactorial_self] all_goals ring all_goals simp

theorem h_weight_intervalIntegrable (k j : ℕ) (α : ℝ) (hα : 0 < α) : IntervalIntegrable (fun x : ℝ => x ^ (α - 1 + (j : ℝ)) * h k α x) volume 0 1 := by have hp : -1 < α - 1 + (j : ℝ) := by have hj : (0 : ℝ) ≤ j := Nat.cast_nonneg j linarith have hh : Continuous (h k α) := by unfold h jacobiShifted fun_prop exact (intervalIntegrable_rpow' hp).mul_continuousOn hh.continuousOn

theorem h_weighted_moment_zero (k j : ℕ) (α : ℝ) (hα : 0 < α) (hj : 1 ≤ j) (hjk : j ≤ k) : (∫ x in (0 : ℝ)..1, x ^ (α - 1 + (j : ℝ)) * h k α x) = 0 := by let p : ℝ → ℝ := fun x => x ^ (j - 1) let f := rodKernel k α let c := rising (α + 1) k have hc : c ≠ 0 := ne_of_gt (rising_pos (α + 1) k (by linarith)) have hpow : ∀ x : ℝ, 0 < x → x ^ (α - 1 + (j : ℝ)) * h k α x = c⁻¹ * (p x * (deriv^[k] f) x) := by intro x hx have hexp : α - 1 + (j : ℝ) = ((j - 1 : ℕ) : ℝ) + α := by have hj' : ((j - 1 : ℕ) : ℝ) = (j : ℝ) - 1 := by simp [Nat.cast_sub hj] rw [hj']; ring rw [hexp, Real.rpow_add hx, Real.rpow_natCast] have hr := h_rodrigues k α x hα hx change x ^ α * h k α x = (deriv^[k] f) x / c at hr rw [mul_assoc, hr] dsimp [p] field_simp have hp : ∀ x : ℝ, ContDiffAt ℝ ⊤ p x := by intro x dsimp [p] fun_prop have hkzero : ∀ x : ℝ, (deriv^[k] p) x = 0 := by intro x rw [← iteratedDeriv_eq_iterate, show p = (fun x : ℝ => x ^ (j - 1)) from rfl, iteratedDeriv_pow] have hlt : j - 1 < k := by omega simp [Nat.descFactorial_eq_zero_iff_lt.mpr hlt] have hpleft (i : ℕ) : ContinuousAt (deriv^[i] p) 0 := by have heq : (deriv^[i] p) = fun x : ℝ => (j - 1).descFactorial i * x ^ (j - 1 - i) := by funext x rw [← iteratedDeriv_eq_iterate, show p = (fun x : ℝ => x ^ (j - 1)) from rfl, iteratedDeriv_pow] rw [heq] fun_prop have hleft (i : ℕ) (hik : i < k) : Tendsto (fun x : ℝ => (deriv^[i] f) x) (nhdsWithin 0 (Ioi 0)) (nhds 0) := rodKernel_deriv_tendsto_zero k i α hik hα have hboundary : Tendsto (fun a : ℝ => ∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) a * (deriv^[k - 1 - i] f) a)) (nhdsWithin 0 (Ioi 0)) (nhds 0) := by have ht := (tendsto_finsetSum (range k) (fun i hi => by have hlt : k - 1 - i < k := by have hii : i < k := Finset.mem_range.mp hi omega simpa using ((hpleft i).tendsto.mono_left inf_le_left |>.mul (hleft (k - 1 - i) hlt) |>.const_mul ((-1 : ℝ) ^ i)))) convert ht using 1 <;> first | rfl | simp let a : ℕ → ℝ := fun n => 1 / ((n + 1 : ℕ) : ℝ) have ha (n : ℕ) : 0 < a n := by dsimp [a]; positivity have ha1 (n : ℕ) : a n ≤ 1 := by dsimp [a] have hn : (1 : ℝ) ≤ ((n + 1 : ℕ) : ℝ) := by exact_mod_cast Nat.succ_le_succ (zero_le n) exact (div_le_iff₀ (by positivity)).2 (by simpa using hn) have hatend : Tendsto a atTop (nhdsWithin 0 (Ioi 0)) := by apply tendsto_nhdsWithin_iff.mpr constructor · simpa only [a, Nat.cast_add, Nat.cast_one] using (tendsto_one_div_add_atTop_nhds_zero_nat (𝕜 := ℝ)) · exact Filter.Eventually.of_forall (fun n => ha n) have htrunc (n : ℕ) (hn : 0 < n) : (∫ x in a n..1, x ^ (α - 1 + (j : ℝ)) * h k α x) = -c⁻¹ * (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n))) := by have hfp : ∀ x ∈ Icc (a n) 1, ContDiffAt ℝ ⊤ f x := by intro x hx have hxpos : 0 < x := lt_of_lt_of_le (ha n) hx.1 dsimp [f, rodKernel] exact (Real.contDiffAt_rpow_const_of_ne (ne_of_gt hxpos)).mul (by fun_prop) have han : a n < 1 := by dsimp [a] have hn' : (1 : ℝ) < ((n + 1 : ℕ) : ℝ) := by exact_mod_cast Nat.succ_lt_succ hn exact (div_lt_iff₀ (by positivity)).2 (by simpa using hn') have hibp := integral_mul_iterated_deriv_with_boundary k (a n) 1 p f (fun x hx => hp x) hfp han simp only [hkzero, zero_mul, intervalIntegral.integral_zero, mul_zero, add_zero] at hibp have hright : (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) 1 * (deriv^[k - 1 - i] f) 1 - (deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n))) = -(∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n))) := by rw [← Finset.sum_neg_distrib] apply Finset.sum_congr rfl intro i hi have hlt : k - 1 - i < k := by have hii : i < k := Finset.mem_range.mp hi omega rw [rodKernel_deriv_one k (k - 1 - i) α hlt] ring rw [hright] at hibp calc _ = ∫ x in a n..1, c⁻¹ * (p x * (deriv^[k] f) x) := by apply intervalIntegral.integral_congr intro x hx have hx' : x ∈ Icc (a n) 1 := by simpa only [Set.uIcc_of_le (ha1 n)] using hx exact hpow x (lt_of_lt_of_le (ha n) hx'.1) _ = c⁻¹ * (∫ x in a n..1, p x * (deriv^[k] f) x) := by rw [intervalIntegral.integral_const_mul] _ = _ := by rw [hibp]; ring have hlim1 := tendsto_integral_one_div_nat (fun x : ℝ => x ^ (α - 1 + (j : ℝ)) * h k α x) (h_weight_intervalIntegrable k j α hα) have hlim2 : Tendsto (fun n : ℕ => -c⁻¹ * (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n)))) atTop (nhds 0) := by simpa using (hboundary.comp hatend).const_mul (-c⁻¹) have heq : (∫ x in (0 : ℝ)..1, x ^ (α - 1 + (j : ℝ)) * h k α x) = 0 := tendsto_nhds_unique (hlim1.congr' (Filter.eventually_gt_atTop 0 |>.mono (fun n hn => htrunc n hn))) hlim2 exact heq

theorem h_zeroth_moment_rising (k : ℕ) (α : ℝ) (hα : 0 < α) : α * (∫ x in (0 : ℝ)..1, x ^ (α - 1) * h k α x) = ((Nat.factorial k : ℝ) / rising (α + 1) k) ^ 2 := by let p : ℝ → ℝ := fun x => x⁻¹ let f := rodKernel k α let c := rising (α + 1) k let g : ℝ → ℝ := fun x => x ^ (α - 1) * (1 - x) ^ k let q : ℝ → ℝ := fun x => x ^ (α - 1) * h k α x have hc : c ≠ 0 := ne_of_gt (rising_pos (α + 1) k (by linarith)) have hp (x : ℝ) (hx : 0 < x) : ContDiffAt ℝ ⊤ p x := by dsimp [p] exact contDiffAt_inv ℝ (ne_of_gt hx) have hfp (x : ℝ) (hx : 0 < x) : ContDiffAt ℝ ⊤ f x := by dsimp [f, rodKernel] exact (Real.contDiffAt_rpow_const_of_ne (ne_of_gt hx)).mul (by fun_prop) have hq (x : ℝ) (hx : 0 < x) : q x = c⁻¹ * (p x * (deriv^[k] f) x) := by have hr := h_rodrigues k α x hα hx change x ^ α * h k α x = (deriv^[k] f) x / c at hr dsimp [q, p] rw [show α - 1 = α + (-1 : ℝ) by ring, Real.rpow_add hx, Real.rpow_neg_one] calc x ^ α * x⁻¹ * h k α x = x⁻¹ * (x ^ α * h k α x) := by ring _ = _ := by rw [hr]; ring have hder (x : ℝ) (hx : 0 < x) : (-1 : ℝ) ^ k * ((deriv^[k] p) x * f x) = (k.factorial : ℝ) * g x := by dsimp [p, f, rodKernel, g] rw [iter_deriv_inv] have he : ((-1 - (k : ℤ) : ℤ) : ℝ) + (α + (k : ℝ)) = α - 1 := by push_cast; ring rw [← Real.rpow_intCast] have hs : (-1 : ℝ) ^ k * (-1 : ℝ) ^ k = 1 := by rw [← pow_add, ← two_mul, pow_mul] norm_num calc _ = ((-1 : ℝ) ^ k * (-1 : ℝ) ^ k) * (k.factorial : ℝ) * (x ^ ((-1 - (k : ℤ) : ℤ) : ℝ) * x ^ (α + (k : ℝ))) * (1 - x) ^ k := by ring _ = _ := by rw [hs, ← Real.rpow_add hx, he]; ring let a : ℕ → ℝ := fun n => 1 / ((n + 1 : ℕ) : ℝ) have ha (n : ℕ) : 0 < a n := by dsimp [a]; positivity have ha1 (n : ℕ) : a n ≤ 1 := by dsimp [a] have hn : (1 : ℝ) ≤ ((n + 1 : ℕ) : ℝ) := by exact_mod_cast Nat.succ_le_succ (zero_le n) exact (div_le_iff₀ (by positivity)).2 (by simpa using hn) have hatend : Tendsto a atTop (nhdsWithin 0 (Ioi 0)) := by apply tendsto_nhdsWithin_iff.mpr constructor · simpa only [a, Nat.cast_add, Nat.cast_one] using (tendsto_one_div_add_atTop_nhds_zero_nat (𝕜 := ℝ)) · exact Filter.Eventually.of_forall (fun n => ha n) have hboundary : Tendsto (fun x : ℝ => ∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) x * (deriv^[k - 1 - i] f) x)) (nhdsWithin 0 (Ioi 0)) (nhds 0) := by simpa [p, f] using (tendsto_finsetSum (range k) (fun i hi => (inv_mul_rodKernel_deriv_tendsto_zero k i α hα (Finset.mem_range.mp hi)).const_mul ((-1 : ℝ) ^ i))) have htrunc (n : ℕ) (hn : 0 < n) : (∫ x in a n..1, q x) = -(c⁻¹ * (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n)))) + c⁻¹ * (k.factorial : ℝ) * (∫ x in a n..1, g x) := by have han : a n < 1 := by dsimp [a] have hn' : (1 : ℝ) < ((n + 1 : ℕ) : ℝ) := by exact_mod_cast Nat.succ_lt_succ hn exact (div_lt_iff₀ (by positivity)).2 (by simpa using hn') have hibp := integral_mul_iterated_deriv_with_boundary k (a n) 1 p f (fun x hx => hp x (lt_of_lt_of_le (ha n) hx.1)) (fun x hx => hfp x (lt_of_lt_of_le (ha n) hx.1)) han have hright : (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) 1 * (deriv^[k - 1 - i] f) 1 - (deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n))) = -(∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n))) := by rw [← Finset.sum_neg_distrib] apply Finset.sum_congr rfl intro i hi have hlt : k - 1 - i < k := by have hii : i < k := Finset.mem_range.mp hi omega rw [show f = rodKernel k α from rfl, rodKernel_deriv_one k (k - 1 - i) α hlt] ring rw [hright] at hibp calc _ = c⁻¹ * (∫ x in a n..1, p x * (deriv^[k] f) x) := by rw [← intervalIntegral.integral_const_mul] apply intervalIntegral.integral_congr intro x hx have hx' : x ∈ Icc (a n) 1 := by simpa only [Set.uIcc_of_le (ha1 n)] using hx exact hq x (lt_of_lt_of_le (ha n) hx'.1) _ = _ := by rw [hibp, mul_add, ← intervalIntegral.integral_const_mul] have heq : (∫ x in a n..1, (-1 : ℝ) ^ k * ((deriv^[k] p) x * f x)) = (∫ x in a n..1, (k.factorial : ℝ) * g x) := by apply intervalIntegral.integral_congr intro x hx have hx' : x ∈ Icc (a n) 1 := by simpa only [Set.uIcc_of_le (ha1 n)] using hx exact hder x (lt_of_lt_of_le (ha n) hx'.1) rw [heq, intervalIntegral.integral_const_mul] ring have hlimq := tendsto_integral_one_div_nat q (by simpa only [q, Nat.cast_zero, add_zero] using (h_weight_intervalIntegrable k 0 α hα)) have hlimg := tendsto_integral_one_div_nat g (by dsimp [g] exact (intervalIntegrable_rpow' (by linarith : -1 < α - 1)).mul_continuousOn (by fun_prop)) have hlimb : Tendsto (fun n : ℕ => -(c⁻¹ * (∑ i ∈ range k, (-1 : ℝ) ^ i * ((deriv^[i] p) (a n) * (deriv^[k - 1 - i] f) (a n))))) atTop (nhds 0) := by simpa using ((hboundary.comp hatend).const_mul c⁻¹).neg have heq : (∫ x in (0 : ℝ)..1, q x) = c⁻¹ * (k.factorial : ℝ) * (∫ x in (0 : ℝ)..1, g x) := by apply tendsto_nhds_unique (hlimq.congr' (Filter.eventually_gt_atTop 0 |>.mono (fun n hn => htrunc n hn))) simpa only [a, Nat.cast_add, Nat.cast_one, zero_add, mul_assoc] using hlimb.add ((hlimg.const_mul (k.factorial : ℝ)).const_mul c⁻¹) rw [show (∫ x in (0 : ℝ)..1, x ^ (α - 1) * h k α x) = (∫ x in (0 : ℝ)..1, q x) from rfl, heq, show (∫ x in (0 : ℝ)..1, g x) = (k.factorial : ℝ) / (α * c) by simpa [g, c] using integral_rpow_mul_one_sub_pow k α hα] dsimp [c] have hr : rising (α + 1) k ≠ 0 := hc field_simp [hr]

theorem rising_eq_gamma_ratio (k : ℕ) (α : ℝ) (hα : 0 < α) : rising (α + 1) k = Gamma (α + (k : ℝ) + 1) / Gamma (α + 1) := by have hΓ : Gamma (α + 1) ≠ 0 := ne_of_gt (Real.Gamma_pos_of_pos (by linarith)) induction k with | zero => simp [rising, hΓ] | succ k ih => rw [rising, ascPochhammer_succ_eval] change rising (α + 1) k * (α + 1 + (k : ℝ)) = _ rw [ih] have hpos : α + (k : ℝ) + 1 ≠ 0 := by positivity have hrec := Real.Gamma_add_one hpos have heq : α + ((k + 1 : ℕ) : ℝ) + 1 = (α + (k : ℝ) + 1) + 1 := by push_cast; ring rw [heq, hrec] ring

theorem factorial_div_rising_power_bounds (α : ℝ) (hα : 0 < α) : ∃ c C : ℝ, 0 < c ∧ 0 < C ∧ ∀ k : ℕ, c * ((k + 1 : ℕ) : ℝ) ^ (-α) ≤ ((Nat.factorial k : ℝ) / rising (α + 1) k) ∧ ((Nat.factorial k : ℝ) / rising (α + 1) k) ≤ C * ((k + 1 : ℕ) : ℝ) ^ (-α) := by have hprod (n : ℕ) : (∏ j ∈ range (n + 1), (α + (j : ℝ))) = α * rising (α + 1) n := by induction n with | zero => simp [rising] | succ n ih => rw [Finset.prod_range_succ, ih] change α * rising (α + 1) n * (α + ↑(n + 1)) = α * (ascPochhammer ℝ (n + 1)).eval (α + 1) rw [ascPochhammer_succ_eval] change α * (ascPochhammer ℝ n).eval (α + 1) * (α + ↑(n + 1)) = α * ((ascPochhammer ℝ n).eval (α + 1) * (α + 1 + ↑n)) push_cast ring let f : ℕ → ℝ := fun k => Real.GammaSeq α (k + 1) have hfpos (k : ℕ) : 0 < f k := by simp only [f, Real.GammaSeq] apply div_pos · exact mul_pos (Real.rpow_pos_of_pos (by positivity) _) (by positivity) · exact Finset.prod_pos fun j hj => by have : 0 ≤ (j : ℝ) := Nat.cast_nonneg j positivity have hΓ : 0 < Gamma α := Real.Gamma_pos_of_pos hα have hlim : Tendsto f atTop (nhds (Gamma α)) := by exact (Real.GammaSeq_tendsto_Gamma α).comp (Filter.tendsto_add_atTop_nat 1) obtain ⟨B, hB⟩ := hlim.bddAbove_range obtain ⟨D, hD⟩ := (hlim.inv₀ hΓ.ne').bddAbove_range have hBpos : 0 < B := (hfpos 0).trans_le (hB ⟨0, rfl⟩) have hDpos : 0 < D := (inv_pos.mpr (hfpos 0)).trans_le (hD ⟨0, rfl⟩) have hflo (k : ℕ) : 1 / D ≤ f k := by apply (div_le_iff₀ hDpos).2 have hi : (f k)⁻¹ ≤ D := hD ⟨k, rfl⟩ have hu := mul_le_mul_of_nonneg_right hi (le_of_lt (hfpos k)) have he : (f k)⁻¹ * f k = 1 := inv_mul_cancel₀ (hfpos k).ne' nlinarith have hfhi (k : ℕ) : f k ≤ B := hB ⟨k, rfl⟩ have hscaled (k : ℕ) : ((k + 1 : ℕ) : ℝ) ^ α * ((Nat.factorial k : ℝ) / rising (α + 1) k) = f k * (α * (α + (k : ℝ) + 1) / ((k + 1 : ℕ) : ℝ)) := by have hr : 0 < rising (α + 1) k := rising_pos _ _ (by linarith) have hr' : 0 < rising (α + 1) (k + 1) := rising_pos _ _ (by linarith) have hk : (0 : ℝ) < (k + 1 : ℕ) := by positivity have ha : 0 < α + (k : ℝ) + 1 := by positivity simp only [f, Real.GammaSeq] rw [hprod (k + 1), Nat.factorial_succ, Nat.cast_mul, Nat.cast_add, Nat.cast_one] simp only [rising] rw [ascPochhammer_succ_eval] field_simp ring refine ⟨α / D, B * α * (α + 1), div_pos hα hDpos, mul_pos (mul_pos hBpos hα) (by linarith), ?_⟩ intro k have hk : (0 : ℝ) < ((k + 1 : ℕ) : ℝ) := by positivity have hratio : 0 < rising (α + 1) k := rising_pos _ _ (by linarith) have hfaclow : α ≤ α * (α + (k : ℝ) + 1) / ((k + 1 : ℕ) : ℝ) := by apply (le_div_iff₀ hk).2 push_cast nlinarith [(mul_nonneg (sq_nonneg α) (Nat.cast_nonneg k : (0 : ℝ) ≤ (k : ℝ)))] have hfachi : α * (α + (k : ℝ) + 1) / ((k + 1 : ℕ) : ℝ) ≤ α * (α + 1) := by apply (div_le_iff₀ hk).2 push_cast nlinarith [(mul_nonneg (sq_nonneg α) (Nat.cast_nonneg k : (0 : ℝ) ≤ (k : ℝ)))] have hfacpos : 0 ≤ α * (α + (k : ℝ) + 1) / ((k + 1 : ℕ) : ℝ) := le_trans hα.le hfaclow have hlo : α / D ≤ ((k + 1 : ℕ) : ℝ) ^ α * ((Nat.factorial k : ℝ) / rising (α + 1) k) := by rw [hscaled] calc α / D = (1 / D) * α := by ring _ ≤ f k * α := mul_le_mul_of_nonneg_right (hflo k) hα.le _ ≤ f k * (α * (α + (k : ℝ) + 1) / ((k + 1 : ℕ) : ℝ)) := mul_le_mul_of_nonneg_left hfaclow (hfpos k).le have hhi : ((k + 1 : ℕ) : ℝ) ^ α * ((Nat.factorial k : ℝ) / rising (α + 1) k) ≤ B * α * (α + 1) := by rw [hscaled] calc _ ≤ B * (α * (α + (k : ℝ) + 1) / ((k + 1 : ℕ) : ℝ)) := mul_le_mul_of_nonneg_right (hfhi k) hfacpos _ ≤ B * (α * (α + 1)) := mul_le_mul_of_nonneg_left hfachi hBpos.le _ = B * α * (α + 1) := by ring have hp : 0 < ((k + 1 : ℕ) : ℝ) ^ α := Real.rpow_pos_of_pos hk _ constructor · rw [Real.rpow_neg hk.le] change (α / D) / (((k + 1 : ℕ) : ℝ) ^ α) ≤ _ exact (div_le_iff₀ hp).2 (by simpa [mul_comm] using hlo) · rw [Real.rpow_neg hk.le] change _ ≤ (B * α * (α + 1)) / (((k + 1 : ℕ) : ℝ) ^ α) exact (le_div_iff₀ hp).2 (by simpa [mul_comm] using hhi)