# How close does any (n, k, p) on the grid come to the Wald interval boundary?
# margin = |d^2 n - z^2 phat(1-phat)| in exact rationals; compare with the effect of replacing
# the printed z by the true normal quantile (they differ by < 1e-15).
from fractions import Fraction as F
from mpmath import mp, erfinv, sqrt, mpf
mp.dps = 40
ztrue = sqrt(2) * erfinv(mpf("0.95"))
z = F("1.959963984540054")
dz2 = abs(float(z * z) - float(ztrue ** 2))
worst = None
for n in range(5, 101):
    for k in range(1, n):
        ph = F(k, n); v = z * z * ph * (1 - ph)
        for j in range(1, 1000):
            d = ph - F(j, 1000)
            m = abs(d * d * n - v) / (ph * (1 - ph))   # margin in units of z^2
            if worst is None or m < worst[0]:
                worst = (m, n, k, j)
print("true z =", mp.nstr(ztrue, 20), "| printed z =", "1.959963984540054")
print("|z_printed^2 - z_true^2| ~", dz2)
print("smallest boundary margin (units of z^2):", float(worst[0]), "at n, k, p =", worst[1], worst[2], worst[3] / 1000)
print("decisions robust to the choice of z:", float(worst[0]) > dz2)