Theory Buffons_Needle

(*
  File:    Buffons_Needle.thy
  Author:  Manuel Eberl <manuel@pruvisto.org>

  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_raw ‹
\begin{figure}
\begin{center}
\begin{tikzpicture}
\coordinate (lefttick) at (-3,0);
\coordinate (righttick) at (3,0);

\draw (lefttick) -- (righttick);
\draw [thick] (lefttick) ++ (0,0.4) -- ++(0,3);
\draw [thick] (righttick) ++ (0,0.4) -- ++(0,3);

\coordinate (needle) at (1,2);

\newcommand{\needleangle}{55}
\newcommand{\needlelength}{{1}}
\newcommand{\needlethickness}{0.6pt}

\draw ($(lefttick)+(0,4pt)$) -- ($(lefttick)-(0,4pt)$);
\draw ($(righttick)+(0,4pt)$) -- ($(righttick)-(0,4pt)$);
\draw (0,4pt) -- (0,-4pt);

\draw [densely dashed, thin] let \p1 = (needle) in (\x1, 0) -- (needle);
\draw [densely dashed, thin] let \p1 = (needle) in (needle) -- (3, \y1);
\draw (needle) ++ (15pt,0) arc(0:\needleangle:15pt);
\path (needle) -- ++(15pt,0) node [above, midway, yshift=-1.9pt, xshift=1.8pt] {$\scriptstyle\varphi$};

\node [below, xshift=-3.5pt] at ($(lefttick)-(0,4pt)$) {$-\nicefrac{d}{2}$};
\node [below] at ($(righttick)-(0,4pt)$) {$\nicefrac{d}{2}$};
\node [below,yshift=-1pt] at (0,-4pt) {$0$};
\node [below,yshift=-2pt] at (needle |- 0,-4pt) {$x$};

\draw[<->] (needle) ++({\needleangle+90}:5pt) ++(\needleangle:{-\needlelength}) -- ++(\needleangle:2) node [midway, above, rotate=\needleangle] {$\scriptstyle l$};

\draw [line width=0.7pt,fill=white] (needle) ++({\needleangle+90}:\needlethickness) -- ++(\needleangle:\needlelength) arc({\needleangle+90}:{\needleangle-90}:\needlethickness) 
  -- ++(\needleangle:-\needlelength) -- ++(\needleangle:-\needlelength) arc({\needleangle+270}:{\needleangle+90}:\needlethickness) -- ++(\needleangle:\needlelength);

\end{tikzpicture}
\end{center}
\caption{A sketch of the situation in Buffon's needle experiment. There is a needle of length $l$
with its centre at a certain $x$ coordinate, angled at an angle $\varphi$ off the horizontal axis.
The two vertical lines are a distance of $d$ apart, each being $\nicefrac{d}{2}$ away from the
origin.}
\label{fig:buffon}
\end{figure}

\definecolor{myred}{HTML}{cc2428}
\begin{figure}[h]
\begin{center}
\begin{tikzpicture}
  \begin{axis}[
          xmin=0, xmax=7, ymin=0, ymax=1,
          width=\textwidth, height=0.6\textwidth,
          xlabel={$l/d$}, ylabel={$\mathcal P$}, tick style={thin,black},
          ylabel style = {rotate=270,anchor=west},
  ] 
  \addplot [color=myred, line width=1pt, mark=none,domain=0:1,samples=200] ({x}, {2/pi*x}); 
  \addplot [color=myred, line width=1pt, mark=none,domain=1:7,samples=200] ({x}, {2/pi*(x-sqrt(x*x-1)+acos(1/x)/180*pi)}); 
  \end{axis}
\end{tikzpicture}
\caption{The probability $\mathcal P$ of the needle hitting one of the lines, as a function of the quotient $l/d$ (where $l$ is the length of the needle and $d$ the horizontal distance between the lines).}
\label{fig:buffonplot}
\end{center}
\end{figure}
›

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$)?

  We will show that, if we let $x := \nicefrac{l}{d}$, the probability of this is
  \[
  \mathcal P_{l,d} =
    \begin{cases}
      \nicefrac{2}{\pi} \cdot x & \text{if}\ l \leq d\\
      \nicefrac{2}{\pi}\cdot(x - \sqrt{x^2 - 1} + \arccos (\nicefrac{1}{x})) & \text{if}\ l \geq d
    \end{cases} 
  \]
  A plot of this function can be found in Figure~\ref{fig:buffonplot}.
›

locale Buffon =
  fixes d l :: real
  assumes d: "d > 0" and l: "l > 0"
begin

definition Buffon :: "(real × real) measure" where
  "Buffon = uniform_measure lborel ({-d/2..d/2} × {-pi..pi})"

lemma space_Buffon [simp]: "space Buffon = UNIV"
  by (simp add: Buffon_def)

definition Buffon_set :: "(real × real) set" where
  "Buffon_set = {(x,φ)  {-d/2..d/2} × {-pi..pi}. needle l x φ  {-d/2, d/2}  {}}"


subsection ‹Derivation of the solution›

text ‹
  The following form is a bit easier to handle.
›
lemma Buffon_set_altdef1:
  "Buffon_set =
     {(x,φ)  {-d/2..d/2} × {-pi..pi}.
         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 -
  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 show ?thesis unfolding Buffon_set_def case_prod_unfold
    by (intro Collect_cong conj_cong refl) meson
qed

lemma Buffon_set_altdef2:
  "Buffon_set = {(x,φ)  {-d/2..d/2} × {-pi..pi}. abs x  d / 2 - abs (sin φ) * l / 2}"
  unfolding Buffon_set_altdef1
proof (intro Collect_cong prod.case_cong refl conj_cong)
  fix x φ
  assume *: "(x, φ)  {-d/2..d/2} × {-pi..pi}"
  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"

  show "?P x φ  (d / 2 - ¦sin φ¦ * l / 2  ¦x¦)"
  proof (cases "φ  0")
    case True
    have "x - l / 2 * sin φ  x + l / 2 * sin φ" using l True *
      by (auto simp: sin_ge_zero)
    moreover from True and * have "sin φ  0" by (auto simp: sin_ge_zero)
    ultimately show ?thesis using * True
      by (force simp: field_simps Let_def min_def max_def case_prod_unfold abs_if)
  next
    case False
    with * have "x - l / 2 * sin φ  x + l / 2 * sin φ" using l
      by (auto simp: sin_le_zero' mult_nonneg_nonpos)
    moreover from False and * have "sin φ  0" by (auto simp: sin_le_zero')
    ultimately show ?thesis using * False l d
      by (force simp: field_simps Let_def min_def max_def case_prod_unfold abs_if)
  qed
qed
  
    
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) set" where
  "Buffon_set' = {(x,φ)  {0..d/2} × {0..pi}. x  d / 2 - sin φ * l / 2}"

lemma closed_buffon_set [simp, intro, measurable]: "closed Buffon_set"
  proof -
  have "Buffon_set = ({-d/2..d/2} × {-pi..pi})  
          (λz. abs (fst z) + abs (sin (snd z)) * l / 2 - d / 2) -` {0..}" 
    (is "_ = ?A") unfolding Buffon_set_altdef2 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'"
proof -
  have "Buffon_set' = ({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  sets borel" 
  by measurable

lemma measurable_buffon_set' [measurable]: "Buffon_set'  sets borel" 
  by measurable


sublocale prob_space Buffon
  unfolding Buffon_def
proof -
  have "emeasure lborel ({- d / 2..d / 2} × {- pi..pi}) = ennreal (2 * d * pi)"
    unfolding lborel_prod [symmetric] using d
    by (subst lborel.emeasure_pair_measure_Times)
       (auto simp: ennreal_mult mult_ac simp flip: ennreal_numeral)
  also have "  0    "
    using d by auto
  finally show "prob_space (uniform_measure lborel ({- d / 2..d / 2} × {- pi..pi}))"
    by (intro prob_space_uniform_measure) auto
qed

lemma buffon_prob_aux:
  "emeasure Buffon {(x,φ). needle l x φ  {-d/2, d/2}  {}} =
     emeasure lborel Buffon_set / 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 "{(x, φ). needle l x φ  {- d / 2, d / 2}  {}}  sets borel"
    by (intro pred_Collect_borel)
       (simp add: borel_prod [symmetric] needle_def closed_segment_eq_real_ivl case_prod_unfold)
  hence "emeasure Buffon {(x,φ). needle l x φ  {-d/2, d/2}  {}} =
           emeasure lborel (({-d/2..d/2} × {- pi..pi})  {(x,φ). needle l x φ  {-d/2, d/2}  {}}) /
           emeasure lborel ({-(d/2)..d/2} × {-pi..pi})"
    unfolding Buffon_def Buffon_set_def by (subst emeasure_uniform_measure) simp_all
  also have "({-d/2..d/2} × {- pi..pi})  {(x, φ). needle l x φ  {-d/2, d/2}  {}} = Buffon_set"
    unfolding Buffon_set_def by auto
  also have "emeasure lborel ({-(d/2)..d/2} × {-pi..pi}) = ennreal (2 * d * pi)"
    using d by (simp flip: lborel_prod ennreal_mult add: lborel.emeasure_pair_measure_Times)
  finally show ?thesis .
qed

lemma emeasure_buffon_set_conv_buffon_set':
  "emeasure lborel Buffon_set = 4 * emeasure lborel Buffon_set'"
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'"
  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 = 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_altdef2 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  B"
        by (auto split: prod.splits simp: Buffon_set_altdef2 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_altdef2 Buffon_set'_def sin_le_zero' A_def B_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_altdef2 Buffon_set'_def sin_le_zero' A_def B_def D_def)
      thus ?thesis by blast
    qed
  next
    case (2 z)
    thus ?case using d l
      by (auto simp: Buffon_set_altdef2 Buffon_set'_def sin_ge_zero sin_le_zero'  A_def B_def C_def D_def)
  qed
  
  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 = 
                  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' = 
     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' = ennreal l"
proof -
  have "emeasure lborel Buffon_set' =
          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_real_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 = 4 * ennreal l"
  by (simp add: emeasure_buffon_set_conv_buffon_set' emeasure_buffon_set'_short l_le_d)

lemma prob_short_aux:
  "Buffon {(x, φ). needle l x φ  {- d / 2, d / 2}  {}} = ennreal (2 * l / (d * pi))"
  unfolding buffon_prob_aux emeasure_buffon_set_short using d l
  by (simp flip: ennreal_mult ennreal_numeral add: divide_ennreal)

lemma prob_short: "𝒫((x,φ) in Buffon. needle l x φ  {-d/2, d/2}  {}) = 2 * l / (d * pi)"
  using prob_short_aux unfolding emeasure_eq_measure
  using l d by (subst (asm) ennreal_inj) auto

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: 
  shows "l * (1 - sqrt (1 - (d / l)2)) + arccos (d / l) * d  0"
  and   "emeasure lborel Buffon_set' =
           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)
  have "cos φ' = sqrt (1 - (d / l)^2)"
    unfolding φ'_def by (rule cos_arcsin) (insert d l l_ge_d, auto simp: field_simps)

  have "l * (1 - cos φ') + arccos (d / l) * d  0"
    using l d l_ge_d
    by (intro add_nonneg_nonneg mult_nonneg_nonneg arccos_lbound) (auto simp: field_simps)   
  thus "l * (1 - sqrt (1 - (d / l)2)) + arccos (d / l) * d  0"
    by (simp add: cos φ' = sqrt (1 - (d / l)^2))

  let ?f = "(λx. min (d / 2) (sin x * l / 2))"
  have "emeasure lborel Buffon_set' = 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_real_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 note cos φ' = sqrt (1 - (d / l)^2)
  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 "emeasure lborel Buffon_set' =
                  ennreal (l * (1 - sqrt (1 - (d / l)2)) + arccos (d / l) * d)" .
qed

lemma emeasure_set_long: "emeasure lborel Buffon_set = 
        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)

lemma prob_long_aux: 
  shows "2 / pi * ((l / d) - sqrt ((l / d)2 - 1) + arccos (d / l))  0"
  and   "Buffon {(x, φ). needle l x φ  {- d / 2, d / 2}  {}} = 
           ennreal (2 / pi * ((l / d) - sqrt ((l / d)2 - 1) + arccos (d / l)))"
  using emeasure_buffon_set'_long(1)
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 "l / d  sqrt ((l / d)2 - 1)"
    using l d l_ge_d by (intro real_le_lsqrt) (auto simp: field_simps)
  thus "2 / pi * ((l / d) - sqrt ((l / d)2 - 1) + arccos (d / l))  0"
    using d l l_ge_d
    by (intro mult_nonneg_nonneg add_nonneg_nonneg arccos_lbound) (auto simp: field_simps)

  have "emeasure Buffon {(x,φ). needle l x φ  {-d/2, d/2}  {}} = 
          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_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 "emeasure Buffon {(x,φ). needle l x φ  {-d/2, d/2}  {}} = 
                  ennreal (2 / pi * ((l / d) - sqrt ((l / d)2 - 1) + arccos (d / l)))" .
qed

lemma prob_long:
  "𝒫((x,φ) in Buffon. needle l x φ  {-d/2, d/2}  {}) =
     2 / pi * ((l / d) - sqrt ((l / d)2 - 1) + arccos (d / l))"
  using prob_long_aux unfolding emeasure_eq_measure
  by (subst (asm) ennreal_inj) simp_all

end

theorem prob_eq:
  defines "x  l / d"
  shows   "𝒫((x,φ) in Buffon. needle l x φ  {-d/2, d/2}  {}) =
             (if l  d then
                2 / pi * x
              else
                2 / pi * (x - sqrt (x2 - 1) + arccos (1 / x)))"
  using prob_short prob_long unfolding x_def by auto

end

end