Calculate z-score given probability using erfinv

Your method to compute $z$ is unstable. You suffer from catastrophic cancellation if you compute $2p-1$:

I do not know your algorithm, but let's assume that you are using an correctly rounded double precision function $\mathrm{erf}^{-1}$. Now let $p$ be smaller than half of the machine epsilon ($\approx 2.2 \cdot 10^{-16}$ for IEEE double), then $2p-1=-1$ and all precision is lost. For larger $p$ you are loosing 'only' some precision.

A solution is to use the inverse of $\mathrm{erfc}$

$$p = \frac{1}{2}\biggl(1 + \mathrm{erf}\left(\frac {z}{\sqrt2}\right)\biggr) = \frac{1}{2} \mathrm{erfc}\left(-\frac {z}{\sqrt2}\right)$$

This follows from the facts that $\mathrm{erfc}(x)=1-\mathrm{erf}(x)$ and $\mathrm{erf}(-x) = -\mathrm{erf}(x).$ Now solve for $z$:

$$z = -\sqrt{2} \; \mathrm{erfc}^{-1}(2p)$$

For a C++ implementation of erfc_inv see e.g. the Boost file erf_inv.hpp. My Pascal double routines give the following values:

    p             z with erf_inv          z with erfc_inv
 1.000E-03     -3.0902323061678136     -3.0902323061678136
 1.000E-04     -3.7190164854557088     -3.7190164854556809
 1.000E-05     -4.2648907939226017     -4.2648907939228247
 1.000E-06     -4.7534243088283059     -4.7534243088228987
 1.000E-07     -5.1993375821874714     -5.1993375821928174
 1.000E-08     -5.6120012442658487     -5.6120012441747891
 1.000E-09     -5.9978070105847330     -5.9978070150076874
 1.000E-10     -6.3613408896974226     -6.3613409024040566
 1.000E-11     -6.7060231434147486     -6.7060231554951368
 1.000E-12     -7.0344869100478356     -7.0344838253011313
 1.000E-13     -7.3488287482023118     -7.3487961028006774
 1.000E-14     -7.6507309051556440     -7.6506280929352704
 1.000E-15     -7.9414444874159802     -7.9413453261709970
 1.000E-16     -8.2095361516013874     -8.2220822161304366
 1.000E-17                     Nan     -8.4937932241095986
 1.000E-18                     Nan     -8.7572903487823162

$$ $$ Edit regarding algorithms: As already noted, there is http://www.boost.org/doc/libs/1_65_1/libs/math/doc/html/math_toolkit/sf_erf/error_inv.html which uses seven rational approximations.

The GSL routines gsl_cdf_gaussian_Pinv, gsl_cdf_gaussian_Qinv from https://www.gnu.org/software/gsl/doc/html/randist.html#the-gaussian-distribution use the method from M. J. Wichura, “Algorithm AS 241: the percentage points of the normal distribution,” Applied Statistics, vol. 37, no. 3, pp. 477–484, 1988.

The Cephes library function ndtri computes for small $p$ an expansion in $\sqrt{-2 \ln p}$.

An asymptotic expansion of erfc⁡_inv for small arguments is given at http://dlmf.nist.gov/7.17#iii.

The above results are from my Pascal port of the Boost routines.