feat (special_functions/gamma): doubling formula (#19136)

This PR adds the proof of the Legendre doubling formula, relating `Gamma s * Gamma (s + 1 / 2)` to `Gamma (2 * s)`.

src/analysis/analytic/isolated_zeros.lean

@@ -195,4 +195,15 @@ theorem eq_on_of_preconnected_of_mem_closure (hf : analytic_on 𝕜 f U) (hg : a
   eq_on f g U :=
 hf.eq_on_of_preconnected_of_frequently_eq hg hU h₀ (mem_closure_ne_iff_frequently_within.mp hfg)
 
+/-- The *identity principle* for analytic functions, global version: if two functions on a normed
+field `𝕜` are analytic everywhere and coincide at points which accumulate to a point `z₀`, then
+they coincide globally.
+For higher-dimensional versions requiring that the functions coincide in a neighborhood of `z₀`,
+see `eq_of_eventually_eq`. -/
+theorem eq_of_frequently_eq [connected_space 𝕜]
+  (hf : analytic_on 𝕜 f univ) (hg : analytic_on 𝕜 g univ)
+  (hfg : ∃ᶠ z in 𝓝[≠] z₀, f z = g z) : f = g :=
+funext (λ x, eq_on_of_preconnected_of_frequently_eq hf hg is_preconnected_univ
+    (mem_univ z₀) hfg (mem_univ x))
+
 end analytic_on

src/analysis/analytic/uniqueness.lean

@@ -89,7 +89,7 @@ begin
   exact uniform_space.completion.coe_injective F this,
 end
 
-/-- The *identity principle* for analytic functions: If two analytic function coincide in a whole
+/-- The *identity principle* for analytic functions: If two analytic functions coincide in a whole
 neighborhood of a point `z₀`, then they coincide globally along a connected set.
 For a one-dimensional version assuming only that the functions coincide at some points
 arbitrarily close to `z₀`, see `eq_on_of_preconnected_of_frequently_eq`. -/
@@ -103,4 +103,14 @@ begin
     λ z hz, (hf.sub hg).eq_on_zero_of_preconnected_of_eventually_eq_zero hU h₀ hfg' hz,
 end
 
+/-- The *identity principle* for analytic functions: If two analytic functions on a normed space
+coincide in a neighborhood of a point `z₀`, then they coincide everywhere.
+For a one-dimensional version assuming only that the functions coincide at some points
+arbitrarily close to `z₀`, see `eq_of_frequently_eq`. -/
+theorem eq_of_eventually_eq {f g : E → F} [preconnected_space E]
+  (hf : analytic_on 𝕜 f univ) (hg : analytic_on 𝕜 g univ) {z₀ : E} (hfg : f =ᶠ[𝓝 z₀] g) :
+  f = g :=
+funext (λ x, eq_on_of_preconnected_of_eventually_eq hf hg is_preconnected_univ
+    (mem_univ z₀) hfg (mem_univ x))
+
 end analytic_on

src/analysis/special_functions/gamma/beta.lean

@@ -5,7 +5,8 @@ Authors: David Loeffler
 -/
 import analysis.convolution
 import analysis.special_functions.trigonometric.euler_sine_prod
-import analysis.special_functions.gamma.basic
+import analysis.special_functions.gamma.bohr_mollerup
+import analysis.analytic.isolated_zeros
 
 /-!
 # The Beta function, and further properties of the Gamma function
@@ -28,8 +29,10 @@ refined properties of the Gamma function using these relations.
 * `complex.Gamma_mul_Gamma_one_sub`: Euler's reflection formula
   `Gamma s * Gamma (1 - s) = π / sin π s`.
 * `complex.differentiable_one_div_Gamma`: the function `1 / Γ(s)` is differentiable everywhere.
+* `complex.Gamma_mul_Gamma_add_half`: Legendre's duplication formula
+  `Gamma s * Gamma (s + 1 / 2) = Gamma (2 * s) * 2 ^ (1 - 2 * s) * sqrt π`.
 * `real.Gamma_ne_zero`, `real.Gamma_seq_tendsto_Gamma`,
-  `real.Gamma_mul_Gamma_one_sub`: real versions of the above results.
+  `real.Gamma_mul_Gamma_one_sub`, `real.Gamma_mul_Gamma_add_half`: real versions of the above.
 -/
 
 noncomputable theory
@@ -560,3 +563,64 @@ end
 end complex
 
 end inv_gamma
+
+section doubling
+/-!
+## The doubling formula for Gamma
+
+We prove the doubling formula for arbitrary real or complex arguments, by analytic continuation from
+the positive real case. (Knowing that `Γ⁻¹` is analytic everywhere makes this much simpler, since we
+do not have to do any special-case handling for the poles of `Γ`.)
+-/
+
+namespace complex
+
+theorem Gamma_mul_Gamma_add_half (s : ℂ) :
+  Gamma s * Gamma (s + 1 / 2) = Gamma (2 * s) * 2 ^ (1 - 2 * s) * ↑(real.sqrt π) :=
+begin
+  suffices : (λ z, (Gamma z)⁻¹ * (Gamma (z + 1 / 2))⁻¹) =
+    (λ z, (Gamma (2 * z))⁻¹  * 2 ^ (2 * z - 1) / ↑(real.sqrt π)),
+  { convert congr_arg has_inv.inv (congr_fun this s) using 1,
+    { rw [mul_inv, inv_inv, inv_inv] },
+    { rw [div_eq_mul_inv, mul_inv, mul_inv, inv_inv, inv_inv, ←cpow_neg, neg_sub] } },
+  have h1 : analytic_on ℂ (λ z : ℂ, (Gamma z)⁻¹ * (Gamma (z + 1 / 2))⁻¹) univ,
+  { refine differentiable_on.analytic_on _ is_open_univ,
+    refine (differentiable_one_div_Gamma.mul _).differentiable_on,
+    exact differentiable_one_div_Gamma.comp (differentiable_id.add (differentiable_const _)) },
+  have h2 : analytic_on ℂ (λ z, (Gamma (2 * z))⁻¹  * 2 ^ (2 * z - 1) / ↑(real.sqrt π)) univ,
+  { refine differentiable_on.analytic_on _ is_open_univ,
+    refine (differentiable.mul _ (differentiable_const _)).differentiable_on,
+    apply differentiable.mul,
+    { exact differentiable_one_div_Gamma.comp (differentiable_id'.const_mul _) },
+    { refine λ t, differentiable_at.const_cpow _ (or.inl two_ne_zero),
+      refine differentiable_at.sub_const (differentiable_at_id.const_mul _) _ } },
+  have h3 : tendsto (coe : ℝ → ℂ) (𝓝[≠] 1) (𝓝[≠] 1),
+  { rw tendsto_nhds_within_iff, split,
+    { exact tendsto_nhds_within_of_tendsto_nhds continuous_of_real.continuous_at },
+    { exact eventually_nhds_within_iff.mpr (eventually_of_forall $ λ t ht, of_real_ne_one.mpr ht)}},
+  refine analytic_on.eq_of_frequently_eq h1 h2 (h3.frequently _),
+  refine ((eventually.filter_mono nhds_within_le_nhds) _).frequently,
+  refine (eventually_gt_nhds zero_lt_one).mp (eventually_of_forall $ λ t ht, _),
+  rw [←mul_inv, Gamma_of_real, (by push_cast : (t : ℂ) + 1 / 2 = ↑(t + 1 / 2)), Gamma_of_real,
+    ←of_real_mul, Gamma_mul_Gamma_add_half_of_pos ht, of_real_mul, of_real_mul, ←Gamma_of_real,
+    mul_inv, mul_inv, (by push_cast : 2 * (t : ℂ) = ↑(2 * t)), Gamma_of_real,
+    of_real_cpow zero_le_two, of_real_bit0, of_real_one, ←cpow_neg, of_real_sub, of_real_one,
+    neg_sub, ←div_eq_mul_inv]
+end
+
+end complex
+
+namespace real
+open complex
+
+lemma Gamma_mul_Gamma_add_half (s : ℝ) :
+  Gamma s * Gamma (s + 1 / 2) = Gamma (2 * s) * 2 ^ (1 - 2 * s) * sqrt π :=
+begin
+  rw [←of_real_inj],
+  simpa only [←Gamma_of_real, of_real_cpow zero_le_two, of_real_mul, of_real_add, of_real_div,
+    of_real_bit0, of_real_one, of_real_sub] using complex.Gamma_mul_Gamma_add_half ↑s
+end
+
+end real
+
+end doubling

src/analysis/special_functions/gamma/bohr_mollerup.lean

@@ -4,6 +4,7 @@ Released under Apache 2.0 license as described in the file LICENSE.
 Authors: David Loeffler
 -/
 import analysis.special_functions.gamma.basic
+import analysis.special_functions.gaussian
 
 /-! # Convexity properties of the Gamma function
 
@@ -25,11 +26,65 @@ context of this proof, we place them in a separate namespace `real.bohr_mollerup
 general form of the Euler limit formula valid for any real or complex `x`; see
 `real.Gamma_seq_tendsto_Gamma` and `complex.Gamma_seq_tendsto_Gamma` in the file
 `analysis.special_functions.gamma.beta`.)
+
+As an application of the Bohr-Mollerup theorem we prove the Legendre doubling formula for the
+Gamma function for real positive `s` (which will be upgraded to a proof for all complex `s` in a
+later file).
+
+TODO: This argument can be extended to prove the general `k`-multiplication formula (at least up
+to a constant, and it should be possible to deduce the value of this constant using Stirling's
+formula).
 -/
 
 noncomputable theory
 open filter set measure_theory
-open_locale nat ennreal topology big_operators
+open_locale nat ennreal topology big_operators real
+
+section convexity
+
+-- Porting note: move the following lemmas to `Analysis.Convex.Function`
+variables {𝕜 E β : Type*} {s : set E} {f g : E → β}
+  [ordered_semiring 𝕜] [has_smul 𝕜 E] [add_comm_monoid E] [ordered_add_comm_monoid β]
+
+lemma convex_on.congr [has_smul 𝕜 β] (hf : convex_on 𝕜 s f) (hfg : eq_on f g s) :
+  convex_on 𝕜 s g :=
+⟨hf.1, λ x hx y hy a b ha hb hab,
+  by simpa only [←hfg hx, ←hfg hy, ←hfg (hf.1 hx hy ha hb hab)] using hf.2 hx hy ha hb hab⟩
+
+lemma concave_on.congr [has_smul 𝕜 β](hf : concave_on 𝕜 s f) (hfg : eq_on f g s) :
+  concave_on 𝕜 s g :=
+⟨hf.1, λ x hx y hy a b ha hb hab,
+  by simpa only [←hfg hx, ←hfg hy, ←hfg (hf.1 hx hy ha hb hab)] using hf.2 hx hy ha hb hab⟩
+
+lemma strict_convex_on.congr [has_smul 𝕜 β] (hf : strict_convex_on 𝕜 s f) (hfg : eq_on f g s) :
+  strict_convex_on 𝕜 s g :=
+⟨hf.1, λ x hx y hy hxy a b ha hb hab, by simpa only
+  [←hfg hx, ←hfg hy, ←hfg (hf.1 hx hy ha.le hb.le hab)] using hf.2 hx hy hxy ha hb hab⟩
+
+lemma strict_concave_on.congr [has_smul 𝕜 β] (hf : strict_concave_on 𝕜 s f) (hfg : eq_on f g s) :
+  strict_concave_on 𝕜 s g :=
+⟨hf.1, λ x hx y hy hxy a b ha hb hab, by simpa only
+  [←hfg hx, ←hfg hy, ←hfg (hf.1 hx hy ha.le hb.le hab)] using hf.2 hx hy hxy ha hb hab⟩
+
+lemma convex_on.add_const [module 𝕜 β] (hf : convex_on 𝕜 s f) (b : β) :
+  convex_on 𝕜 s (f + (λ _, b)) :=
+hf.add (convex_on_const _ hf.1)
+
+lemma concave_on.add_const [module 𝕜 β] (hf : concave_on 𝕜 s f) (b : β) :
+  concave_on 𝕜 s (f + (λ _, b)) :=
+hf.add (concave_on_const _ hf.1)
+
+lemma strict_convex_on.add_const {γ : Type*} {f : E → γ}
+  [ordered_cancel_add_comm_monoid γ] [module 𝕜 γ] (hf : strict_convex_on 𝕜 s f) (b : γ) :
+  strict_convex_on 𝕜 s (f + (λ _, b)) :=
+hf.add_convex_on (convex_on_const _ hf.1)
+
+lemma strict_concave_on.add_const {γ : Type*} {f : E → γ}
+  [ordered_cancel_add_comm_monoid γ] [module 𝕜 γ] (hf : strict_concave_on 𝕜 s f) (b : γ) :
+  strict_concave_on 𝕜 s (f + (λ _, b)) :=
+hf.add_concave_on (concave_on_const _ hf.1)
+
+end convexity
 
 namespace real
 
@@ -112,17 +167,12 @@ end
 
 lemma convex_on_Gamma : convex_on ℝ (Ioi 0) Gamma :=
 begin
-  refine ⟨convex_Ioi 0, λ x hx y hy a b ha hb hab, _⟩,
-  have := convex_on.comp (convex_on_exp.subset (subset_univ _) _) convex_on_log_Gamma
-    (λ u hu v hv huv, exp_le_exp.mpr huv),
-  convert this.2 hx hy ha hb hab,
-  { rw [function.comp_app, exp_log (Gamma_pos_of_pos $ this.1 hx hy ha hb hab)] },
-  { rw [function.comp_app, exp_log (Gamma_pos_of_pos hx)] },
-  { rw [function.comp_app, exp_log (Gamma_pos_of_pos hy)] },
-  { rw convex_iff_is_preconnected,
-    refine is_preconnected_Ioi.image _ (λ x hx, continuous_at.continuous_within_at _),
-    refine (differentiable_at_Gamma (λ m, _)).continuous_at.log (Gamma_pos_of_pos hx).ne',
-    exact (neg_lt_iff_pos_add.mpr (add_pos_of_pos_of_nonneg hx (nat.cast_nonneg m))).ne' }
+  refine ((convex_on_exp.subset (subset_univ _) _).comp convex_on_log_Gamma
+    (exp_monotone.monotone_on _)).congr (λ x hx, exp_log (Gamma_pos_of_pos hx)),
+  rw convex_iff_is_preconnected,
+  refine is_preconnected_Ioi.image _ (λ x hx, continuous_at.continuous_within_at _),
+  refine (differentiable_at_Gamma (λ m, _)).continuous_at.log (Gamma_pos_of_pos hx).ne',
+  exact (neg_lt_iff_pos_add.mpr (add_pos_of_pos_of_nonneg hx (nat.cast_nonneg m))).ne',
 end
 
 end convexity
@@ -377,4 +427,89 @@ end
 
 end strict_mono
 
+section doubling
+
+/-!
+## The Gamma doubling formula
+
+As a fun application of the Bohr-Mollerup theorem, we prove the Gamma-function doubling formula
+(for positive real `s`). The idea is that `2 ^ s * Gamma (s / 2) * Gamma (s / 2 + 1 / 2)` is
+log-convex and satisfies the Gamma functional equation, so it must actually be a constant
+multiple of `Gamma`, and we can compute the constant by specialising at `s = 1`. -/
+
+/-- Auxiliary definition for the doubling formula (we'll show this is equal to `Gamma s`) -/
+def doubling_Gamma (s : ℝ) : ℝ := Gamma (s / 2) * Gamma (s / 2 + 1 / 2) * 2 ^ (s - 1) / sqrt π
+
+lemma doubling_Gamma_add_one (s : ℝ) (hs : s ≠ 0) :
+  doubling_Gamma (s + 1) = s * doubling_Gamma s :=
+begin
+  rw [doubling_Gamma, doubling_Gamma, (by abel : s + 1 - 1 = s - 1 + 1), add_div, add_assoc,
+    add_halves (1 : ℝ), Gamma_add_one (div_ne_zero hs two_ne_zero), rpow_add two_pos, rpow_one],
+  ring,
+end
+
+lemma doubling_Gamma_one : doubling_Gamma 1 = 1 :=
+by simp_rw [doubling_Gamma, Gamma_one_half_eq, add_halves (1 : ℝ), sub_self, Gamma_one, mul_one,
+  rpow_zero, mul_one, div_self (sqrt_ne_zero'.mpr pi_pos)]
+
+lemma log_doubling_Gamma_eq :
+  eq_on (log ∘ doubling_Gamma) (λ s, log (Gamma (s / 2)) + log (Gamma (s / 2 + 1 / 2))
+    + s * log 2 - log (2 * sqrt π)) (Ioi 0) :=
+begin
+  intros s hs,
+  have h1 : sqrt π ≠ 0, from sqrt_ne_zero'.mpr pi_pos,
+  have h2 : Gamma (s / 2) ≠ 0, from (Gamma_pos_of_pos $ div_pos hs two_pos).ne',
+  have h3 : Gamma (s / 2 + 1 / 2) ≠ 0,
+    from (Gamma_pos_of_pos $ add_pos (div_pos hs two_pos) one_half_pos).ne',
+  have h4 : (2 : ℝ) ^ (s - 1) ≠ 0, from (rpow_pos_of_pos two_pos _).ne',
+  rw [function.comp_app, doubling_Gamma, log_div (mul_ne_zero (mul_ne_zero h2 h3) h4) h1,
+    log_mul (mul_ne_zero h2 h3) h4, log_mul h2 h3, log_rpow two_pos, log_mul two_ne_zero h1],
+  ring_nf,
+end
+
+lemma doubling_Gamma_log_convex_Ioi : convex_on ℝ (Ioi (0:ℝ)) (log ∘ doubling_Gamma) :=
+begin
+  refine (((convex_on.add _ _).add _).add_const _).congr log_doubling_Gamma_eq.symm,
+  { convert convex_on_log_Gamma.comp_affine_map
+      (distrib_mul_action.to_linear_map ℝ ℝ (1 / 2 : ℝ)).to_affine_map,
+    { simpa only [zero_div] using (preimage_const_mul_Ioi (0 : ℝ) one_half_pos).symm, },
+    { ext1 x,
+      change log (Gamma (x / 2)) = log (Gamma ((1 / 2 : ℝ) • x)),
+      rw [smul_eq_mul, mul_comm, mul_one_div] } },
+  { refine convex_on.subset _ (Ioi_subset_Ioi $ neg_one_lt_zero.le) (convex_Ioi _),
+    convert convex_on_log_Gamma.comp_affine_map ((distrib_mul_action.to_linear_map ℝ ℝ
+      (1 / 2 : ℝ)).to_affine_map + affine_map.const _ _ (1 / 2 : ℝ)),
+    { change Ioi (-1 : ℝ) = ((λ x : ℝ, x + 1 / 2) ∘ (λ x : ℝ, (1 / 2 : ℝ) * x)) ⁻¹' (Ioi 0),
+      rw [preimage_comp, preimage_add_const_Ioi, zero_sub, preimage_const_mul_Ioi (_ : ℝ)
+        one_half_pos, neg_div, div_self (@one_half_pos ℝ _).ne'] },
+    { ext1 x,
+      change log (Gamma (x / 2 + 1 / 2)) = log (Gamma ((1 / 2 : ℝ) • x + 1 / 2)),
+      rw [smul_eq_mul, mul_comm, mul_one_div] } },
+  { simpa only [mul_comm _ (log _)]
+      using (convex_on_id (convex_Ioi (0 : ℝ))).smul (log_pos one_lt_two).le }
+end
+
+lemma doubling_Gamma_eq_Gamma {s : ℝ} (hs : 0 < s) : doubling_Gamma s = Gamma s :=
+begin
+  refine eq_Gamma_of_log_convex doubling_Gamma_log_convex_Ioi
+    (λ y hy, doubling_Gamma_add_one y hy.ne') (λ y hy, _) doubling_Gamma_one hs,
+  apply_rules [mul_pos, Gamma_pos_of_pos, add_pos, inv_pos_of_pos,
+    rpow_pos_of_pos, two_pos, one_pos, sqrt_pos_of_pos pi_pos]
+end
+
+/-- Legendre's doubling formula for the Gamma function, for positive real arguments. Note that
+we shall later prove this for all `s` as `real.Gamma_mul_Gamma_add_half` (superseding this result)
+but this result is needed as an intermediate step. -/
+lemma Gamma_mul_Gamma_add_half_of_pos {s : ℝ} (hs : 0 < s) :
+  Gamma s * Gamma (s + 1 / 2) = Gamma (2 * s) * 2 ^ (1 - 2 * s) * sqrt π :=
+begin
+  rw [←(doubling_Gamma_eq_Gamma (mul_pos two_pos hs)),
+    doubling_Gamma, mul_div_cancel_left _ (two_ne_zero' ℝ),
+    (by abel : 1 - 2 * s = -(2 * s - 1)), rpow_neg zero_le_two],
+  field_simp [(sqrt_pos_of_pos pi_pos).ne', (rpow_pos_of_pos two_pos (2 * s - 1)).ne'],
+  ring,
+end
+
+end doubling
+
 end real