scipy norm.sf(z) error grows like z^2 ulp; don't assert tails at rtol=1e-15
SciPy 1.18.1, checked against mpmath at 50 digits by int64-view ulp distance: st.norm.sf(z) is 4, 12, 44, 76, 224 and 376 ulp off at z = 2, 5, 8, 12, 20, 30. scipy.special.erfc evaluated at the already-rounded argument z/np.sqrt(2) is within 0-56 ulp of the exact erfc at that same float, so most of the sf error comes from rounding the argument. The tail's relative condition number is about 2x^2 for erfc(x), which amplifies the 1-ulp rounding of x = z/sqrt 2. 0.5*erfc(z/np.sqrt(2)) is similarly off (46 ulp at z=8, 225 at z=20), and in the opposite direction from sf, so the two differ by 90 and 449 ulp (relative 1.4e-14 and 5.7e-14). That is irrelevant for a p-value, but assert_allclose(rtol=1e-15) between two tail implementations fails at large z, and so does rtol = z^2 * 1e-16 at z=20. z^2 * 5e-16 covered both measured points. Comparing logsf instead turns a relative error of r into an absolute error of r on a value of about -z^2/2.