import mpmath as mp
mp.mp.dps = 130
for n in range(1, 21):
seed = mp.sqrt(2) * (n - mp.mpf(3)/8) * mp.pi
root = mp.findroot(lambda x: mp.ber(0, x),
(seed-mp.mpf('0.25'), seed+mp.mpf('0.25')))
print(n, mp.nstr(root, 100))
The attached generator uses 800-bit Arb evaluations of $I_0(xe^{\pi i/4})$ to prove opposite signs at rational endpoints enclosing each root. Every returned interval has positive width and lies inside a counted root interval. Its decimal representation contains the root within one unit of the last written place.
The index certificate covers the whole positive axis segment through the twentieth root. For $F=I_0(xe^{\pi i/4})=u+iv$ and $W=uv'-vu'$, the differential equation gives $(xW)'=x|F|^2$. The origin expansion implies $W>0$ for $x>0$, so the continuous phase increases strictly [1] [3]. The origin series rules out a zero in $[0,1]$. A Gronwall bound proves that the phase changes by less than $\pi/2$ in each consecutive interval of width $1/8$ thereafter. Endpoint signs therefore count every zero without omissions. The full derivation and the numerical bounds are attached in certificate.md and checks-evidence.json.
An independent 900-bit real power-series computation, with an explicit geometric tail bound, proves opposite signs at both endpoints of every written hundred-digit interval. It also rejects intervals displaced by ten units of the last place. Nonzero origin and small-argument controls precede the checks. All roots agree with independent hypergeometric computations at 160 and 210 decimal digits, and the first five agree with the printed values [5] within half a unit in their fifth decimal place. The rotated Bessel representations and their derivatives agree at every root midpoint. The Programs snippet gives numerical approximations; the attached generator supplies the certificates.