Probability that a quadratic equation with random coefficients has real roots
If $B$ is uniformly distributed over $[0,1]$ and $X=B^2$, the pdf of $X$ can be computed through:
$$\mathbb{P}[X\leq t] = \mathbb{P}[B\leq\sqrt{t}], $$ from which $f_X(x)=\frac{\mathbb{1}_{(0,1)}(x)}{2\sqrt{x}}$. In a similar way, if $A,C$ are uniformly distributed over $(0,1)$, independent, and $Y=AC$, $$\mathbb{P}[Y\leq t]=\int_{0}^{1}\mathbb{P}\left[C\leq\frac{t}{u}\right]\,du=\int_{0}^{1}\min\left(1,\frac{t}{u}\right)\,du=t-t\log(t)$$ hence $f_Y(y) = -\mathbb{1}_{(0,1)}(y)\cdot\log(y)$. It follows that:
$$\mathbb{P}[B^2\geq 4AC] = \int_{0}^{1}\frac{1}{2\sqrt{x}}\int_{0}^{x/4}-\log(y)\,dy\,dx =\color{red}{\frac{5+6\log 2}{36}\approx 25,44\%}.$$
Sketch of analytic solution. An analytic solution is based on noting that the density of $Q = B^2$ is $f(q) = \frac{1}{2\sqrt{q}},$ for $q \in (0,1),$ the density of $X = 4AC$ is $g(x) = \frac{-\log(x/4)}{4},$ for $x \in (0,4).$
An appropriate double integration of the joint density $h(q, x) = f(q)g(x)$ gives $P(\text{Real Roots}) = P(Q > X) = \frac{5 + \log(64)}{36} = 0.254413.$ [Note: Details of the integration are shown in Horton (2015) referenced in the Question, and in a subsequent Answer.]
Simulation. A simulation based on a million simulated equations is shown below.
m = 10^6
a = runif(m); b = runif(m); c = runif(m)
q = b^2; x = 4*a*c
d = q - x # discriminant
real = (d > 0)
mean(real)
## 0.254302 # approximates analytic result
The following figure is based on 30,000 simulated equations. Histograms of simulated values of $Q$ and $X$ show the theoretical density curves. The brown lines in the scatterplot and histogram of values of the discriminant $D\,$ separate equations with real (blue) and complex solutions.
Alternate discrete version. The following simulation of the discrete version of the problem shows that a little over 20% of the computer generated quadratic equations have real roots.
m = 10^6; val=seq(.1, 1.0, by=.1)
a = sample(val,m,rep=T)
b = sample(val, m, rep=T)
c = sample(val, m, rep=T)
q = b^2; x = 4*a*c
d = q - x # discriminant
real = (d > 0)
mean(real)
## 0.206176
single = (d == 0)
mean(single)
## 0.007964
The final result suggests that eight of the 1000 possible equations have only a single root. It is not difficult to see that the discriminant can be zero only if $B = .2, .4, .8,$ or $1.0$. From there, simple arithmetic shows that there are indeed eight possible combinations of values of $A$ and $C$ that produce $D = 0.$
A complete analytic solution to the discrete version would seem to involve some tedious bookkeeping, beginning with the ten possible values of $B^2.$ Perhaps there is a clever way to get an exact analytic solution using convolutions of discrete distributions in Matlab.
Someone asked almost this a few days ago. I did want to point out that a cube is not the most natural shape to consider for this problem, although it was the one chosen. Better, in some ways, to consider the ball $A^2 + B^2 + C^2 \leq R^2.$ In that case we take rotated coordinates $$ u = B; v = (A - C)/ \sqrt 2; w = (A + C)/ \sqrt 2. $$
Then the condition $B^2 \geq 4AC$ becomes $u^2 + 2 v^2 \geq 2 w^2.$ It makes no difference here whether we use volume or surface area, so we are asking for the total surface area of the two peculiar elliptical patches $u^2 + 2 v^2 \geq 2 w^2$ on the sphere $u^2 + v^2 + w^2 = 1,$ divided by $4 \pi.$ On second thought, the figure we want is one minus this, the funny annular region.
Not sure I know how to calculate this.