Tail p-values: 1 - norm.cdf(x) is 0 from x=9, use sf/logsf
Measured with SciPy 1.18.1. 1 - st.norm.cdf(x) cancels catastrophically: at x=8 it gives 6.661e-16 vs the true 6.221e-16 (7% off), and from x=9 it returns exactly 0 (true 1.129e-19). st.norm.sf(x) matches 0.5*erfc(x/sqrt 2) to the last digit down to 2.754e-89 at x=20. At x=40, sf also underflows to 0, but logsf gives -804.61. So compute p-values, likelihood tails and z-score thresholds with dist.sf / dist.logsf (also isf instead of ppf(1-p)), never 1 - cdf. Heavy tails make it less visible: t(df=3) at x=1e3 still agrees (1.103e-9). The same overflow trap holds for determinants: np.linalg.det of a 1000x1000 standard-normal matrix is inf while slogdet gives sign +1, log|det| 2956.0. Rescaling does not rescue it, because det(A/10) is still 6e283 and can flip to 0 or inf for other n.
Self-audit, SciPy 1.18.1: the substance reproduces (x=8: 6.661e-16 vs 6.221e-16; 0 from x=9; logsf(40) = -804.61; t3 at 1e3 1.10265e-9 both ways). One overstatement: sf does not match 0.5*erfc(x/sqrt 2) 'to the last digit'. Measured against mpmath at 50 digits via the int64 view, sf is 44 ulp off at x=8 and 224 ulp at x=20. 0.5*erfc(x/np.sqrt(2)) is 46 and 225 ulp off in the other direction, so the two differ by 90 and 449 ulp. Relative 1e-14 is irrelevant for a p-value, but neither is last-digit exact. The det example is seed-dependent. For default_rng(0..2), 1000x1000, det was +inf/-inf/-inf and log|det| 2952.6/2949.3/2951.6. det(A/10) was 1.9e282, -7.2e280 and -7.0e281. So the point holds, but my 2956.0 and 6e283 were one draw.