A random quadratic with three uniform coefficients, Geometric Probability

Problem 8.1

Let $a, b, c$ be independent uniform random variables in $[0, 1]$ . What is the probability that the quadratic $a x^2 + b x + c = 0$ has real roots?

Solution. The quadratic has real roots iff its discriminant is non-negative: $b^2 \geq 4 a c, \qquad \text{i.e.}, \qquad a c \leq b^2 / 4.$

Step 1: $\mathbb{P}(a c \leq t)$ for fixed $t \in (0, 1)$ . With $a, c$ iid uniform in $[0, 1]$ , $\mathbb{P}(a c \leq t) = \int_0^1 \int_0^1 \mathbf{1}[a c \leq t] \, dc \, da.$ For $a \leq t$ the inner integral gives $1$ (since $a c \leq a \leq t$ always). For $a > t$ it gives $t/a$ . Hence $\mathbb{P}(a c \leq t) = t + \int_t^1 \tfrac{t}{a}\, da = t(1 - \ln t).$

Step 2: Integrate over $b$ . Set $t = b^2/4$ , which lies in $[0, 1/4]$ for $b \in [0, 1]$ : $\mathbb{P}(\text{real roots}) = \int_0^1 \tfrac{b^2}{4} \left( 1 - \ln \tfrac{b^2}{4} \right) db.$ Expand $\ln(b^2 / 4) = 2 \ln b - \ln 4 = 2 \ln b - 2 \ln 2$ , so $1 - \ln(b^2/4) = 1 + 2 \ln 2 - 2 \ln b$ . Hence $\mathbb{P}(\text{real roots}) = \frac{1}{4} \int_0^1 b^2 (1 + 2 \ln 2 - 2 \ln b) \, db.$ Using $\int_0^1 b^2 \, db = \tfrac13$ and $\int_0^1 b^2 \ln b \, db = -\tfrac19$ : $\begin{align*} \mathbb{P}(\text{real roots}) &= \tfrac{1}{4}\left[ (1 + 2 \ln 2) \cdot \tfrac13 - 2 \cdot (-\tfrac19) \right] \\ &= \tfrac{1}{4} \cdot \tfrac{3 + 6\ln 2 + 2}{9} = \frac{5 + 6 \ln 2}{36}. \end{align*}$ Therefore $\boxed{\mathbb{P}(\text{real roots}) = \frac{5 + 6 \ln 2}{36} \approx 0.2544.} \qedhere$

Compare with Chapter 4, where the parametrisation $x^2 + 2 p x + q$ with $p, q \in [-1, 1]$ gave the cleaner probability $\tfrac23$ . The present version has three random coefficients rather than two, and adds the rescaling $a$ on the leading coefficient; the resulting answer involves $\ln 2$ — an irrational constant with no simple expression in $\pi$ . This is the first sign of how adding an extra random parameter can break a closed-form expression into a transcendental one.

import numpy as np
a, b, c = np.random.rand(3, 10**7)
p = (b*b >= 4*a*c).mean()
print(f"sim: {p:.5f}   exact: {(5 + 6*np.log(2))/36:.5f}")
# sim: 0.25443   exact: 0.25441