I know you were just copying what is here, but I'd be surprised if you can't use rtol equal to the precision of the reference p-value, since it is always much coarser than the precision of the floating point dtype.

I'd also be willing to review a refactor that parameterizes over case, where case includes the method and reference p-value (and preferably statistic, since it's looks like it's currently untested). Testing Stoufer weights could be moved to a separate test, and testing the cases in which the p-values are equal could either be a separate test (parameterized over all methods) or removed if it's not actually important.

Tightening the tolerance causes a failure on all backends

Them the reference p-value is probably wrong. Fine to leave it alone for this PR, but maybe when this is refactored I'll look for a better reference that we can cite in the test.

Well this is indeed what R gives.

library(metap)
x = c(0.01, 0.2, 0.3)
sumz(x)

sumz =  2.131791 p =  0.01651203

I'll see if there are any hints in the code or documentation as to why there is such a discrepancy.

No such hints. I think It's worth follow-up in a refactoring PR. I don't think an issue is needed if you were willing to submit such a PR soon. LMK either way.