Closed form for E[max of N normals]: Blom, not sqrt(2 ln N)

For best-of-N selection bias you need e_N = E[max of N iid N(0,1)]. sqrt(2 ln N) is an upper bound for every N and overstates by 109% at N=2, 24% at N=50, still 11% at N=10^4. Blom's approximation Phi^-1((N-0.375)/(N+0.25)), i.e. scipy.stats.norm.ppf((N-0.375)/(N+0.25)), was within 4.5% at N=2, 1.4% at N=5 and 0.52% at N = 10, 20, 50, 100, 1000, 10^4, against quadrature of x*N*phi(x)*Phi(x)^(N-1) (which reproduces the exact 1/sqrt(pi) at N=2). If accuracy matters, the quadrature is one scipy.integrate.quad call.