import mpmath as mp
mp.mp.dps = 130
p = mp.mpf(19) / 20
mp.findroot(lambda x: mp.jtheta(4, 0, mp.exp(-2*x*x)) - p,
(mp.mpf('0.25'), mp.mpf(4)), solver='bisect',
maxsteps=600, tol=mp.mpf('1e-120'))# Machine precision, using the lower-tail probability.
from scipy.stats import kstwobign
kstwobign.ppf(0.95)The generator bisects exact rational brackets with 461-bit Arb arithmetic, increasing the precision if a sign is undecidable. Every cumulative probability includes the alternating-series remainder bound $2\exp(-2(N+1)^2x^2)$ after $N$ terms.
The returned ball encloses the whole final bracket, whose width is less than $10^{-112}$; every returned ball retained at least 110 decimal digits before rounding to 100. Independently, 560-bit MPFI interval arithmetic evaluates the positive series in (1) with a geometric bound on its omitted terms. It certifies opposite signs of $F(x)-p$ at the two endpoints represented by every stored decimal. The checks also test strict monotonicity, the upper-tail conversion, a special-value control, and rejection of shifted intervals. All 14 values agree with mpmath at the stored precision and with SciPy at machine precision.