=== Part A: p by Simpson quadrature of log W (standard library only) ===
-- BHT critical rate k = ceil(d/2 + gamma sqrt d), 2 gamma = sqrt(e pi); predicted p -> 1/2
d=1e3  k=547             m=93        y=-0.3631  p_quad=0.929392  Phi(-y/rt2)=0.601312
d=1e4  k=5147            m=293       y=-0.2012  p_quad=0.885675  Phi(-y/rt2)=0.556563
d=1e5  k=50463           m=925       y=-0.1105  p_quad=0.814606  Phi(-y/rt2)=0.531136
d=1e6  k=501462          m=2923      y=-0.0656  p_quad=0.736209  Phi(-y/rt2)=0.518498
d=1e8  k=50014612        m=29223     y=-0.0239  p_quad=0.611331  Phi(-y/rt2)=0.506754
d=1e10 k=5000146115      m=292229    y=-0.0065  p_quad=0.545519  Phi(-y/rt2)=0.501830
d=1e12 k=500001461142    m=2922283   y=-0.0021  p_quad=0.517685  Phi(-y/rt2)=0.500600
-- window m = round(m_c + y sqrt(m_c)), y in {-1,0,1}
y0=-1 k=1e4  d=19606           p_quad=0.951293  Phi(-y/rt2)=0.759426
y0=-1 k=1e6  d=1995931         p_quad=0.890556  Phi(-y/rt2)=0.761761
y0=-1 k=1e8  d=199958875       p_quad=0.826470  Phi(-y/rt2)=0.760274
y0=-1 k=1e10 d=19999587369     p_quad=0.789004  Phi(-y/rt2)=0.760343
y0=-1 k=1e12 d=1999995869301   p_quad=0.771668  Phi(-y/rt2)=0.760298
y0=+0 k=1e4  d=19586           p_quad=0.843836  Phi(-y/rt2)=0.503503
y0=+0 k=1e6  d=1995866         p_quad=0.699198  Phi(-y/rt2)=0.498812
y0=+0 k=1e8  d=199958672       p_quad=0.592878  Phi(-y/rt2)=0.500435
y0=+0 k=1e10 d=19999586726     p_quad=0.538121  Phi(-y/rt2)=0.500059
y0=+0 k=1e12 d=1999995867268   p_quad=0.514827  Phi(-y/rt2)=0.500049
y0=+1 k=1e4  d=19565           p_quad=0.634715  Phi(-y/rt2)=0.235246
y0=+1 k=1e6  d=1995802         p_quad=0.432230  Phi(-y/rt2)=0.239803
y0=+1 k=1e8  d=199958468       p_quad=0.318211  Phi(-y/rt2)=0.239323
y0=+1 k=1e10 d=19999586083     p_quad=0.270533  Phi(-y/rt2)=0.239750
y0=+1 k=1e12 d=1999995865235   p_quad=0.251467  Phi(-y/rt2)=0.239779
-- shift k = d/2 + (sqrt(e pi)/2) sqrt d + c d^(1/4): predicted p -> Phi(-sqrt2 c/(pi e)^(1/4))
c=-1 d=1e6  p_quad=0.927286  predicted limit=0.795962
c=-1 d=1e8  p_quad=0.866239  predicted limit=0.795962
c=-1 d=1e10 p_quad=0.827328  predicted limit=0.795962
c=-1 d=1e12 p_quad=0.808309  predicted limit=0.795962
c=+1 d=1e6  p_quad=0.436142  predicted limit=0.204038
c=+1 d=1e8  p_quad=0.294145  predicted limit=0.204038
c=+1 d=1e10 p_quad=0.238309  predicted limit=0.204038
c=+1 d=1e12 p_quad=0.216870  predicted limit=0.204038
=== Part B: explicit inequalities (own digamma/polygamma code) ===
checked 8593 admissible pairs, d from 4 to 6309573
  (N a/b - F^2)/F^2 (must be <= 0)                 worst=-3.1673e-07  violations=0
  (F^2 - (N a/b)(1+2/b))/F^2 (<= 0)                worst=-3.1674e-07  violations=0
  |A-T| / ((17/24) N/b) (must be <= 1)             worst=+8.9673e-01  violations=0
  kappa4/kappa2^2 - bound (<= 0)                   worst=-6.4219e-07  violations=0
  (1-p) - sqrt(e rho), N>=pi, d<=2e4 (<= 0)        worst=-5.6419e-03  violations=0
  F * sup density of N log W (<= 1)                worst=+1.0000e+00  violations=0
  F * sup density / (2e^2) (<= 1)                  worst=+6.7668e-02  violations=0
