Simulation with your sample-size of $n=10^7$ gave you an observed proportion $\hat p=0.4285833,$ so an approximate confidence interval for $p:=P(a^2<bc)$, with, say a $95\%$ "confidence level", is
$$\hat p\pm 1.96\,\sqrt{\hat p(1-\hat p)\over n}=(0.4283, 0.4289)$$
Although $3/7=0.42857...$ does lie in this interval, the interval only esimates $p$ to about three digits of precision. Simulations using a larger sample-size suggest, with a very high degree of confidence/credibility, that $p\ne {3\over 7},$ the difference showing up in the fourth decimal place. Here's a picture of what I find with sample-size $10^{11}$:
The posterior distribution is $\text{Beta}(a,b)$, with $(a,b)=(n\hat p+{1\over 2}, n(1-\hat p)+{1\over 2}).$
The $100(1-\alpha)\%$ credibility interval is therefore $$(B_{\alpha\over 2},B_{1-{\alpha\over 2}})$$
where $B_q$ is the $q$-th quantile of the $\text{Beta}(a,b)$ distribution.
The $100(1-\alpha)\%$ confidence interval is $$\hat p\pm z_{\alpha\over 2}\,\sqrt{\hat p(1-\hat p)\over n}$$
where $z_q$ is the $q$-th quantile of the standard normal distribution.
With the sample-size this large the various approaches to confidence/credibility intervals for $p$ should all give practically the same results (assuming a non-informative prior distribution in the Bayesian methods); e.g., they agree to about six decimal places for a $99.9\%$ confidence/credibility interval, namely
$$0.428790 \pm 0.000005 = (0.428785, 0.428795).$$
(I found no OEIS entries for decimal expansions in that interval, so I can't guess at the exact value.)
Best Answer
The answer I propose is:
$$Pr(a)=\dfrac{b+c-a}{a+b+c}, \ \ \ \ Pr(b)=\dfrac{c+a-b}{a+b+c}, \ \ \ \ Pr(c)=\dfrac{a+b-c}{a+b+c} \tag{0}$$
Why that ? This question deserves to be formulated at first in terms of "measure theory", probabilities being considered in a second step, in the framework of the domain called "Integral Geometry" or "Geometric Probability" (see Historical note below).
As all plane lines can be given a unique standardized form:
$x \cos \theta+y \sin \theta=p, \tag{1}$
Therefore, $(p,\theta)\in \mathbb{R^+} \times [0,2\pi)$ is a parameterization of the set of plane lines.
With this parametrization, a natural measure (invariant for the group of direct isometries) is
$$\mu(S)=\int_S dpd\theta$$
(see first reference at the bottom of this text).
For example the measure of the set of straight lines intersecting the disk centered in $0$ with radius $R$ is:
$$ \int_{(p,\theta)\in [0,R] \times [0,2\pi)}dpd\theta=\int_{p=0}^R\int_{\theta=0}^{2 \pi}dpd\theta=2 \pi R$$
(i.e., its perimeter).
Let $S_{A|B,C}$ be the set of lines that separate vertex $A$ from vertices $B,C$, otherwise said, that intersect $AB$ (length $c$) and $AC$ (length $b$) but not $BC$ (length $a$).
One can compute the measure of set $S_{A|B,C}$ to be:
$$ \mu(S_{A|B,C}) \ = \ \int_{S_{A|B,C}} dpd\theta \ = \ b+c-a \tag{2}$$
and similar formulas for the other cases.
It is remarkable that the positivity of measure (2) is equivalent to triangle inequality...
One has to work with conditional probabilities, i.e., by assuming that we limit our attention to lines $(L)$ intersecting triangle $ABC$.
Due to the fact that a line that intersects a triangle intersects two and only two of its sides (setting apart lines passing through one of the vertices, which have measure $0$, therefore are neglectable from the point of view of measure theory):
$$(L) \in S_{A|B,C}, \ \ (L) \in S_{B|C,A}, \ \ (L) \in S_{C|A,B}$$
are disjoint and their union is the whole set of lines (complete set of events in probabilistic terms). Therefore, we have a total measure:
$$(b+c-a) \ + \ (c+a-b) \ + \ (a+b-c) \ = \ a+b+c$$
i.e, the perimeter of the triangle (as was found for the disk above) !
As a consequence, the (conditional) probabilities you are looking for can "logicaly" been taken as given by (0).
Historical corner: "Integral Geometry" (aka "Geometric probability") was initiated by Buffon in the mid-18th century, then by Crofton in the $1880$s, followed by Poincaré (cinematic formulas, see reference below) around 1900, but especially developed by Blaschke and his students circa $1930-1950$ (see Santalo's book "Integral Geometry and Geometric Probability" for a rather clear exposition).
Remark: A cousin domain of Integral Geometry is "Stereology"
A nice presentation of Integral Geometry can be found here.
A rather accessible reference here.
For a generalisation of the perimeter result, see Cauchy-Crofton formula
For cinematic formulas see this (more difficult!) reference.
Another (very theoretical) reference.
Edit: For 3D, use Crofton formula (one of the many formulas of stereology dealing with a number of intercepts $N(...)$ of a variable plane $P$ with a fixed triangle $T$):
$$\int N(P \cap T) d\mu \ = \ \pi L(T)$$
where $L(T)$ is the length (= perimeter) of triangle $T$.
($N=0,2$ in the case of a triangle, $N=1$ being a case with measure $0$).
Of course, practically speaking, the integral in the LHS will be approximated using a mean value...