Theory Buffons_Needle

theory Buffons_Needle
imports Probability
(*
  File:    Buffons_Needle.thy
  Author:  Manuel Eberl <eberlm@in.tum.de>

  A formal solution of Buffon's needle problem.
*)
section ‹Buffon's Needle Problem›
theory Buffons_Needle
  imports "HOL-Probability.Probability"
begin

subsection ‹Auxiliary material›

lemma sin_le_zero': "sin x ≤ 0" if "x ≥ -pi" "x ≤ 0" for x
  by (metis minus_le_iff neg_0_le_iff_le sin_ge_zero sin_minus that(1) that(2))


subsection ‹Problem definition›

text ‹
  Consider a needle of length $l$ whose centre has the $x$-coordinate $x$. The following then
  defines the set of all $x$-coordinates that the needle covers 
  (i.e. the projection of the needle onto the $x$-axis.)
›
definition needle :: "real ⇒ real ⇒ real ⇒ real set" where
  "needle l x φ = closed_segment (x - l / 2 * sin φ) (x + l / 2 * sin φ)"

text ‹
  Buffon's Needle problem is then this: Assuming the needle's $x$ position is chosen uniformly
  at random in a strip of width $d$ centred at the origin, what is the probability that the 
  needle crosses at least one of the left/right boundaries of that strip (located at 
  $x = \pm\frac{1}{2}d$)?
›
definition buffon :: "real ⇒ real ⇒ bool measure" where
  "buffon l d = 
     do {
       (x, φ) ← uniform_measure lborel ({-d/2..d/2} × {-pi..pi});
       return (count_space UNIV) (needle l x φ ∩ {-d/2, d/2} ≠ {})
     }"


subsection ‹Derivation of the solution›

text ‹
  The following form is a bit easier to handle.
›
lemma buffon_altdef:
  "buffon l d =
     do {
       (x, φ) ← uniform_measure lborel ({-d/2..d/2} × {-pi..pi});
       return (count_space UNIV) 
         (let a = x - l / 2 * sin φ; b = x + l / 2 * sin φ
          in  min a b + d/2 ≤ 0 ∧ max a b + d/2 ≥ 0 ∨ min a b - d/2 ≤ 0 ∧ max a b - d/2 ≥ 0)
     }"
proof -
  note buffon_def[of l d]
  also {
    have "(λ(x,φ). needle l x φ ∩ {-d/2, d/2} ≠ {}) =
        (λ(x,φ). let a = x - l / 2 * sin φ; b = x + l / 2 * sin φ
                 in  -d/2 ≥ min a b ∧ -d/2 ≤ max a b ∨ min a b ≤ d/2 ∧ max a b ≥ d/2)"
      by (auto simp: needle_def Let_def closed_segment_eq_real_ivl min_def max_def)
    also have "… = 
      (λ(x,φ). let a = x - l / 2 * sin φ; b = x + l / 2 * sin φ
               in  min a b + d/2 ≤ 0 ∧ max a b + d/2 ≥ 0 ∨ min a b - d/2 ≤ 0 ∧ max a b - d/2 ≥ 0)"
      by (auto simp add: algebra_simps Let_def)
    finally have "(λ(x, φ). return (count_space UNIV) (needle l x φ ∩ {- d/2, d/2} ≠ {})) =
                  (λ(x,φ). return (count_space UNIV) 
                    (let a = x - l / 2 * sin φ; b = x + l / 2 * sin φ
                     in  min a b + d/2 ≤ 0 ∧ max a b + d/2 ≥ 0 ∨ min a b - d/2 ≤ 0 ∧ max a b - d/2 ≥ 0))"
      by (simp add: case_prod_unfold fun_eq_iff)
  }
  finally show ?thesis .
qed
    
text ‹
  It is obvious that the problem boils down to determining the measure of the following set:
›
definition buffon_set :: "real ⇒ real ⇒ (real × real) set" where
  "buffon_set l d = {(x,φ) ∈ {-d/2..d/2} × {-pi..pi}. abs x ≥ d / 2 - abs (sin φ) * l / 2}"

text ‹
  By using the symmetry inherent in the problem, we can reduce the problem to the following 
  set, which corresponds to one quadrant of the original set:
›
definition buffon_set' :: "real ⇒ real ⇒ (real × real) set" where
  "buffon_set' l d = {(x,φ) ∈ {0..d/2} × {0..pi}. x ≥ d / 2 - sin φ * l / 2}"

lemma closed_buffon_set [simp, intro, measurable]: "closed (buffon_set l d)"
proof -
  have "buffon_set l d = ({-d/2..d/2} × {-pi..pi}) ∩ 
          (λz. abs (fst z) + abs (sin (snd z)) * l / 2 - d / 2) -` {0..}" 
    (is "_ = ?A") unfolding buffon_set_def by auto
  also have "closed …"
    by (intro closed_Int closed_vimage closed_Times) (auto intro!: continuous_intros)
  finally show ?thesis by simp
qed

lemma closed_buffon_set' [simp, intro, measurable]: "closed (buffon_set' l d)"
proof -
  have "buffon_set' l d = ({0..d/2} × {0..pi}) ∩ 
          (λz. fst z + sin (snd z) * l / 2 - d / 2) -` {0..}" 
    (is "_ = ?A") unfolding buffon_set'_def by auto
  also have "closed …"
    by (intro closed_Int closed_vimage closed_Times) (auto intro!: continuous_intros)
  finally show ?thesis by simp
qed

lemma measurable_buffon_set [measurable]: "buffon_set l d ∈ sets borel" 
  by measurable

lemma measurable_buffon_set' [measurable]: "buffon_set' l d ∈ sets borel" 
  by measurable


context
  fixes d l :: real
  assumes d: "d > 0" and l: "l > 0"
begin

lemma buffon_altdef':
  "buffon l d = distr (uniform_measure lborel ({-d/2..d/2} × {-pi..pi}))
                  (count_space UNIV) (λz. z ∈ buffon_set l d)"
proof -
  let ?P = "λ(x,φ). let a = x - l / 2 * sin φ; b = x + l / 2 * sin φ
                    in  min a b + d/2 ≤ 0 ∧ max a b + d/2 ≥ 0 ∨ min a b - d/2 ≤ 0 ∧ max a b - d/2 ≥ 0"
  have "buffon l d = 
          uniform_measure lborel ({- d / 2..d / 2} × {-pi..pi}) ⤜
          (λz. return (count_space UNIV) (?P z))"
    unfolding buffon_altdef case_prod_unfold by simp
  also have "… = uniform_measure lborel ({- d / 2..d / 2} × {-pi..pi}) ⤜
          (λz. return (count_space UNIV) (z ∈ buffon_set l d))"
  proof (intro bind_cong_AE AE_uniform_measureI AE_I2 impI refl return_measurable, goal_cases)
    show "(λz. return (count_space UNIV) (?P z))
             ∈ uniform_measure lborel ({- d / 2..d / 2} × {- pi..pi}) →M
                 subprob_algebra (count_space UNIV)"
      unfolding Let_def case_prod_unfold lborel_prod [symmetric] by measurable
    show "(λz. return (count_space UNIV) (z ∈ buffon_set l d))
            ∈ uniform_measure lborel ({- d / 2..d / 2} × {- pi..pi}) →M
                subprob_algebra (count_space UNIV)" by simp
    
    case (4 z)
    hence "?P z ⟷ z ∈ buffon_set l d"
    proof (cases "snd z ≥ 0")
      case True
      with 4 have "fst z - l / 2 * sin (snd z) ≤ fst z + l / 2 * sin (snd z)" using l
        by (auto simp: sin_ge_zero)
      moreover from True and 4 have "sin (snd z) ≥ 0" by (auto simp: sin_ge_zero)
      ultimately show ?thesis using 4 True unfolding buffon_set_def
        by (force simp: field_simps Let_def min_def max_def case_prod_unfold abs_if)
    next
      case False
      with 4 have "fst z - l / 2 * sin (snd z) ≥ fst z + l / 2 * sin (snd z)" using l
        by (auto simp: sin_le_zero' mult_nonneg_nonpos)
      moreover from False and 4 have "sin (snd z) ≤ 0" by (auto simp: sin_le_zero')
      ultimately show ?thesis using 4 and False
        unfolding buffon_set_def using l d
        by (force simp: field_simps Let_def min_def max_def case_prod_unfold abs_if)
    qed
    thus ?case by (simp only: )
  qed (simp_all add: borel_prod [symmetric])
  also have "… = distr (uniform_measure lborel ({-d/2..d/2} × {-pi..pi})) 
                    (count_space UNIV) (λz. z ∈ buffon_set l d)"
    by (rule bind_return_distr') simp_all
  finally show ?thesis .
qed

lemma buffon_prob_aux:
  "emeasure (buffon l d) {True} = emeasure lborel (buffon_set l d) / ennreal (2 * d * pi)"
proof -
  have [measurable]: "A × B ∈ sets borel" if "A ∈ sets borel" "B ∈ sets borel" 
    for A B :: "real set" using that unfolding borel_prod [symmetric] by simp
    
  have "emeasure (buffon l d) {True} = 
          emeasure (uniform_measure lborel ({- (d / 2)..d / 2} × {-pi..pi}))
          ((λz. z ∈ buffon_set l d) -` {True})" (is "_ = emeasure ?M _")
    by (simp add: buffon_altdef' emeasure_distr)
  also have "(λz. z ∈ buffon_set l d) -` {True} = buffon_set l d" by auto
  also have "buffon_set l d ⊆ {-d/2..d/2} × {-pi..pi}"
    using l d by (auto simp: buffon_set_def)
  hence "emeasure ?M (buffon_set l d) = 
           emeasure lborel (buffon_set l d) / emeasure lborel ({- (d / 2)..d / 2} × {-pi..pi})"
    by (subst emeasure_uniform_measure) (simp_all add: Int_absorb1)
  also have "emeasure lborel ({- (d / 2)..d / 2} × {-pi..pi}) = ennreal (2 * pi * d)"
    using d by (simp add: lborel_prod [symmetric] lborel.emeasure_pair_measure_Times
                          ennreal_mult algebra_simps)
  finally show ?thesis by (simp add: mult_ac)
qed

lemma emeasure_buffon_set_conv_buffon_set':
  "emeasure lborel (buffon_set l d) = 4 * emeasure lborel (buffon_set' l d)"
proof -
  have distr_lborel [simp]: "distr M lborel f = distr M borel f" for M and f :: "real ⇒ real"
    by (rule distr_cong) simp_all
    
  define A where "A = buffon_set' l d"
  define B C D where "B = (λx. (-fst x, snd x)) -` A" and "C = (λx. (fst x, -snd x)) -` A" and
      "D = (λx. (-fst x, -snd x)) -` A"
  have meas [measurable]:
     "(λx::real × real. (-fst x, snd x)) ∈ borel_measurable borel"
     "(λx::real × real. (fst x, -snd x)) ∈ borel_measurable borel"
     "(λx::real × real. (-fst x, -snd x)) ∈ borel_measurable borel"
    unfolding borel_prod [symmetric] by measurable
  have meas' [measurable]: "A ∈ sets borel" "B ∈ sets borel" "C ∈ sets borel" "D ∈ sets borel"
    unfolding A_def B_def C_def D_def by (rule measurable_buffon_set' measurable_sets_borel meas)+
  
  have *: "buffon_set l d = A ∪ B ∪ C ∪ D"
  proof (intro equalityI subsetI, goal_cases)
    case (1 z)
    show ?case
    proof (cases "fst z ≥ 0"; cases "snd z ≥ 0")
      assume "fst z ≥ 0" "snd z ≥ 0"
      with 1 have "z ∈ A"
        by (auto split: prod.splits simp: buffon_set_def buffon_set'_def sin_ge_zero A_def)
      thus ?thesis by blast   
    next
      assume "¬(fst z ≥ 0)" "snd z ≥ 0"
      with 1 have "z ∈ B"
        by (auto split: prod.splits simp: buffon_set_def buffon_set'_def sin_ge_zero A_def B_def)
      thus ?thesis by blast
    next    
      assume "fst z ≥ 0" "¬(snd z ≥ 0)"
      with 1 have "z ∈ C"
        by (auto split: prod.splits simp: buffon_set_def buffon_set'_def sin_le_zero' A_def C_def)
      thus ?thesis by blast   
    next
      assume "¬(fst z ≥ 0)" "¬(snd z ≥ 0)"
      with 1 have "z ∈ D"
        by (auto split: prod.splits simp: buffon_set_def buffon_set'_def sin_le_zero' A_def D_def)
      thus ?thesis by blast
    qed
  qed (auto simp: buffon_set_def buffon_set'_def sin_ge_zero sin_le_zero'  A_def B_def C_def D_def)
  
  have "A ∩ B = {0} × ({0..pi} ∩ {φ. sin φ * l - d ≥ 0})"
    using d l by (auto simp: buffon_set'_def  A_def B_def C_def D_def)
  moreover have "emeasure lborel … = 0"
    unfolding lborel_prod [symmetric] by (subst lborel.emeasure_pair_measure_Times) simp_all
  ultimately have AB: "(A ∩ B) ∈ null_sets lborel"
    unfolding lborel_prod [symmetric] by (simp add: null_sets_def)
  
  have "C ∩ D = {0} × ({-pi..0} ∩ {φ. -sin φ * l - d ≥ 0})"
    using d l by (auto simp: buffon_set'_def  A_def B_def C_def D_def)
  moreover have "emeasure lborel … = 0"
    unfolding lborel_prod [symmetric] by (subst lborel.emeasure_pair_measure_Times) simp_all
  ultimately have CD: "(C ∩ D) ∈ null_sets lborel"
    unfolding lborel_prod [symmetric] by (simp add: null_sets_def)

  have "A ∩ D = {}" "B ∩ C = {}" using d l 
    by (auto simp: buffon_set'_def A_def D_def B_def C_def)
  moreover have "A ∩ C = {(d/2, 0)}" "B ∩ D = {(-d/2, 0)}"
    using d l by (auto simp: case_prod_unfold buffon_set'_def A_def B_def C_def D_def)
  ultimately have AD: "A ∩ D ∈ null_sets lborel" and BC: "B ∩ C ∈ null_sets lborel" and
    AC: "A ∩ C ∈ null_sets lborel" and BD: "B ∩ D ∈ null_sets lborel" by auto
  
  note *
  also have "emeasure lborel (A ∪ B ∪ C ∪ D) = emeasure lborel (A ∪ B ∪ C) + emeasure lborel D"
    using AB AC AD BC BD CD by (intro emeasure_Un') (auto simp: Int_Un_distrib2)
  also have "emeasure lborel (A ∪ B ∪ C) = emeasure lborel (A ∪ B) + emeasure lborel C"
    using AB AC BC using AB AC AD BC BD CD by (intro emeasure_Un') (auto simp: Int_Un_distrib2)
  also have "emeasure lborel (A ∪ B) = emeasure lborel A + emeasure lborel B"
    using AB using AB AC AD BC BD CD by (intro emeasure_Un') (auto simp: Int_Un_distrib2)
  also have "emeasure lborel B = emeasure (distr lborel lborel (λ(x,y). (-x, y))) A"
    (is "_ = emeasure ?M _") unfolding B_def 
    by (subst emeasure_distr) (simp_all add: case_prod_unfold)
  also have "?M = lborel" unfolding lborel_prod [symmetric]
    by (subst pair_measure_distr [symmetric]) (simp_all add: sigma_finite_lborel lborel_distr_uminus)
  also have "emeasure lborel C = emeasure (distr lborel lborel (λ(x,y). (x, -y))) A"
    (is "_ = emeasure ?M _") unfolding C_def 
    by (subst emeasure_distr) (simp_all add: case_prod_unfold)
  also have "?M = lborel" unfolding lborel_prod [symmetric]
    by (subst pair_measure_distr [symmetric]) (simp_all add: sigma_finite_lborel lborel_distr_uminus)
  also have "emeasure lborel D = emeasure (distr lborel lborel (λ(x,y). (-x, -y))) A"
    (is "_ = emeasure ?M _") unfolding D_def 
    by (subst emeasure_distr) (simp_all add: case_prod_unfold)
  also have "?M = lborel" unfolding lborel_prod [symmetric]
    by (subst pair_measure_distr [symmetric]) (simp_all add: sigma_finite_lborel lborel_distr_uminus)
  finally have "emeasure lborel (buffon_set l d) = 
                  of_nat (Suc (Suc (Suc (Suc 0)))) * emeasure lborel A"
    unfolding of_nat_Suc ring_distribs by simp
  also have "of_nat (Suc (Suc (Suc (Suc 0)))) = (4 :: ennreal)" by simp
  finally show ?thesis unfolding A_def .
qed 

text ‹
  It only remains now to compute the measure of @{const buffon_set'}. We first reduce this
  problem to a relatively simple integral:
›
lemma emeasure_buffon_set':
  "emeasure lborel (buffon_set' l d) = 
     ennreal (integral {0..pi} (λx. min (d / 2) (sin x * l / 2)))"
  (is "emeasure lborel ?A = _")
proof -  
  have "emeasure lborel ?A = nn_integral lborel (λx. indicator ?A x)"
    by (intro nn_integral_indicator [symmetric]) simp_all
  also have "(lborel :: (real × real) measure) = lborel ⨂M lborel" 
    by (simp only: lborel_prod)
  also have "nn_integral … (indicator ?A) = (∫+φ. ∫+x. indicator ?A (x, φ) ∂lborel ∂lborel)"
    by (subst lborel_pair.nn_integral_snd [symmetric]) (simp_all add: lborel_prod borel_prod)
  also have "… = (∫+φ. ∫+x. indicator {0..pi} φ * indicator {max 0 (d/2 - sin φ * l / 2) .. d/2} x ∂lborel ∂lborel)"
    using d l by (intro nn_integral_cong) (auto simp: indicator_def field_simps buffon_set'_def)
  also have "… = ∫+ φ. indicator {0..pi} φ * emeasure lborel {max 0 (d / 2 - sin φ * l / 2)..d / 2} ∂lborel"
    by (subst nn_integral_cmult) simp_all
  also have "… = ∫+ φ. ennreal (indicator {0..pi} φ * min (d / 2) (sin φ * l / 2)) ∂lborel"
    (is "_ = ?I") using d l by (intro nn_integral_cong) (auto simp: indicator_def sin_ge_zero max_def min_def)
  also have "integrable lborel (λφ. (d / 2) * indicator {0..pi} φ)" by simp
  hence int: "integrable lborel (λφ. indicator {0..pi} φ * min (d / 2) (sin φ * l / 2))"
    by (rule Bochner_Integration.integrable_bound)
       (insert l d, auto intro!: AE_I2 simp: indicator_def min_def sin_ge_zero)
  hence "?I = set_lebesgue_integral lborel {0..pi} (λφ. min (d / 2) (sin φ * l / 2))"
    by (subst nn_integral_eq_integral, assumption)
       (insert d l, auto intro!: AE_I2 simp: sin_ge_zero min_def indicator_def set_lebesgue_integral_def)
  also have "… = ennreal (integral {0..pi} (λx. min (d / 2) (sin x * l / 2)))"
    (is "_ = ennreal ?I") using int by (subst set_borel_integral_eq_integral) (simp_all add: set_integrable_def)
  finally show ?thesis by (simp add: lborel_prod)
qed

  
text ‹
  We now have to distinguish two cases: The first and easier one is that where the length 
  of the needle, $l$, is less than or equal to the strip width, $d$:
›
context
  assumes l_le_d: "l ≤ d"
begin

lemma emeasure_buffon_set'_short: "emeasure lborel (buffon_set' l d) = ennreal l"
proof -
  have "emeasure lborel (buffon_set' l d) =
          ennreal (integral {0..pi} (λx. min (d / 2) (sin x * l / 2)))" (is "_ = ennreal ?I")
    by (rule emeasure_buffon_set')
  also have *: "sin φ * l ≤ d" if "φ ≥ 0" "φ ≤ pi" for φ
    using mult_mono[OF l_le_d sin_le_one _ sin_ge_zero] that d by (simp add: algebra_simps)
  have "?I = integral {0..pi} (λx. (l / 2) * sin x)"
    using l d l_le_d  
    by (intro integral_cong) (auto dest: * simp: min_def sin_ge_zero)
  also have "… = l / 2 * integral {0..pi} sin" by simp
  also have "(sin has_integral (-cos pi - (- cos 0))) {0..pi}"
    by (intro fundamental_theorem_of_calculus)
       (auto intro!: derivative_eq_intros simp: has_field_derivative_iff_has_vector_derivative [symmetric])
  hence "integral {0..pi} sin = -cos pi - (-cos 0)"
    by (simp add: has_integral_iff)
  finally show ?thesis by (simp add: lborel_prod)
qed

lemma emeasure_buffon_set_short: "emeasure lborel (buffon_set l d) = 4 * ennreal l"
  by (simp add: emeasure_buffon_set_conv_buffon_set' emeasure_buffon_set'_short l_le_d)

theorem buffon_short: "emeasure (buffon l d) {True} = ennreal (2 * l / (d * pi))"
proof -
  have "emeasure (buffon l d) {True} = ennreal (4 * l) / ennreal (2 * d * pi)"
    using d l by (subst buffon_prob_aux) (simp add: emeasure_buffon_set_short ennreal_mult)
  also have "… = ennreal (4 * l / (2 * d * pi))"
    using d l by (subst divide_ennreal) simp_all
  also have "4 * l / (2 * d * pi) = 2 * l / (d * pi)" by simp
  finally show ?thesis .
qed

end


text ‹
  The other case where the needle is at least as long as the strip width is more complicated:
›
context
  assumes l_ge_d: "l ≥ d"
begin

lemma emeasure_buffon_set'_long: 
  "emeasure lborel (buffon_set' l d) =
     ennreal (l * (1 - sqrt (1 - (d / l)2)) + arccos (d / l) * d)"
proof -
  define φ' where "φ' = arcsin (d / l)"
  have φ'_nonneg: "φ' ≥ 0" unfolding φ'_def using d l l_ge_d arcsin_le_mono[of 0 "d/l"] 
    by (simp add: φ'_def)
  have φ'_le: "φ' ≤ pi / 2" unfolding φ'_def using arcsin_bounded[of "d/l"] d l l_ge_d
    by (simp add: field_simps)
  have ge_phi': "sin φ ≥ d / l" if "φ ≥ φ'" "φ ≤ pi / 2" for φ
    using arcsin_le_iff[of "d / l" "φ"] d l_ge_d that φ'_nonneg by (auto simp: φ'_def field_simps)
  have le_phi': "sin φ ≤ d / l" if "φ ≤ φ'" "φ ≥ 0" for φ
    using le_arcsin_iff[of "d / l" "φ"] d l_ge_d that φ'_le by (auto simp: φ'_def field_simps)
    
  let ?f = "(λx. min (d / 2) (sin x * l / 2))"
  have "emeasure lborel (buffon_set' l d) = ennreal (integral {0..pi} ?f)" (is "_ = ennreal ?I")
    by (rule emeasure_buffon_set')
  also have "?I = integral {0..pi/2} ?f + integral {pi/2..pi} ?f"
    by (rule Henstock_Kurzweil_Integration.integral_combine [symmetric]) (auto intro!: integrable_continuous_real continuous_intros)
  also have "integral {pi/2..pi} ?f = integral {-pi/2..0} (?f ∘ (λφ. φ + pi))"
    by (subst integral_shift) (auto intro!: continuous_intros)
  also have "… = integral {-(pi/2)..-0} (λx. min (d / 2) (sin (-x) * l / 2))" by (simp add: o_def)
  also have "… = integral {0..pi/2} ?f" (is "_ = ?I") by (subst Henstock_Kurzweil_Integration.integral_reflect_real) simp_all
  also have "… + … = 2 * …" by simp
  also have "?I = integral {0..φ'} ?f + integral {φ'..pi/2} ?f"
    using l d l_ge_d φ'_nonneg φ'_le
    by (intro Henstock_Kurzweil_Integration.integral_combine [symmetric]) (auto intro!: integrable_continuous_real continuous_intros)
  also have "integral {0..φ'} ?f = integral {0..φ'} (λx. l / 2 * sin x)"
    using l by (intro integral_cong) (auto simp: min_def field_simps dest: le_phi')
  also have "((λx. l / 2 * sin x) has_integral (- (l / 2 * cos φ') - (- (l / 2 * cos 0)))) {0..φ'}"
    using φ'_nonneg
    by (intro fundamental_theorem_of_calculus)
       (auto simp: has_field_derivative_iff_has_vector_derivative [symmetric] intro!: derivative_eq_intros)
  hence "integral {0..φ'} (λx. l / 2 * sin x) = (1 - cos φ') * l / 2"
    by (simp add: has_integral_iff algebra_simps)
  also have "integral {φ'..pi/2} ?f = integral {φ'..pi/2} (λ_. d / 2)"
    using l by (intro integral_cong) (auto simp: min_def field_simps dest: ge_phi')
  also have "… = arccos (d / l) * d / 2" using φ'_le d l l_ge_d 
    by (subst arccos_arcsin_eq) (auto simp: field_simps φ'_def)
  also have "cos φ' = sqrt (1 - (d / l)^2)"
    unfolding φ'_def by (rule cos_arcsin) (insert d l l_ge_d, auto simp: field_simps)
  also have "2 * ((1 - sqrt (1 - (d / l)2)) * l / 2 + arccos (d / l) * d / 2) = 
               l * (1 - sqrt (1 - (d / l)2)) + arccos (d / l) * d"
    using d l by (simp add: field_simps)
  finally show ?thesis .
qed

lemma emeasure_buffon_set_long: "emeasure lborel (buffon_set l d) = 
        4 * ennreal (l * (1 - sqrt (1 - (d / l)2)) + arccos (d / l) * d)"
  by (simp add: emeasure_buffon_set_conv_buffon_set' emeasure_buffon_set'_long l_ge_d)

theorem buffon_long: 
  "emeasure (buffon l d) {True} = 
     ennreal (2 / pi * ((l / d) - sqrt ((l / d)2 - 1) + arccos (d / l)))"
proof -
  have *: "l * sqrt ((l2 - d2) / l2) + 0 ≤ l + d * arccos (d / l)"
    using d l_ge_d by (intro add_mono mult_nonneg_nonneg arccos_lbound) (auto simp: field_simps)
  have "emeasure (buffon l d) {True} = 
          ennreal (4 * (l - l * sqrt (1 - (d / l)2) + arccos (d / l) * d)) / ennreal (2 * d * pi)"
    using d l l_ge_d * unfolding buffon_prob_aux emeasure_buffon_set_long ennreal_numeral [symmetric]
    by (subst ennreal_mult [symmetric])
       (auto intro!: add_nonneg_nonneg mult_nonneg_nonneg simp: field_simps)
  also have "… = ennreal ((4 * (l - l * sqrt (1 - (d / l)2) + arccos (d / l) * d)) / (2 * d * pi))"
    using d l * by (subst divide_ennreal) (auto simp: field_simps)
  also have "(4 * (l - l * sqrt (1 - (d / l)2) + arccos (d / l) * d)) / (2 * d * pi) =
               2 / pi * (l / d - l / d * sqrt ((d / l)^2 * ((l / d)^2 - 1)) + arccos (d / l))"
    using d l by (simp add: field_simps)
  also have "l / d * sqrt ((d / l)^2 * ((l / d)^2 - 1)) = sqrt ((l / d) ^ 2 - 1)"
    using d l l_ge_d unfolding real_sqrt_mult real_sqrt_abs by simp
  finally show ?thesis .
qed

end
end

end