Let me demonstrate two approaches in this answer.
P.J. Acklam (WayBack archive of his page) devised a not-too-long method to approximate the quantile of the Gaussian distribution, with absolute errors on the order of $10^{-9}$. I have tweaked the Mathematica implementation of Acklam's approximation by William Shaw so that one can obtain a compilable function:
AcklamQuantile = Block[{u},
Compile[{{u, _Real}}, #, RuntimeAttributes -> {Listable}] & @@
Hold[With[{a = {-39.69683028665376, 220.9460984245205, -275.9285104469687,
138.3577518672690, -30.66479806614716, 2.506628277459239},
b = {-54.47609879822406, 161.5858368580409, -155.6989798598866,
66.80131188771972, -13.28068155288572, 1},
c = {-0.007784894002430293, -0.3223964580411365, -2.400758277161838,
-2.549732539343734, 4.374664141464968, 2.938163982698783},
d = {0.007784695709041462, 0.3224671290700398, 2.445134137142996,
3.754408661907416, 1.}},
Which[0.02435 <= u <= 0.97575,
With[{v = u - 1/2},
v Fold[(#1 v^2 + #2) &, 0, a]/
Fold[(#1 v^2 + #2) &, 0, b]] // Evaluate,
u > 0.97575,
With[{q = Sqrt[-2 Log[1 - u]]},
-Fold[(#1 q + #2) &, 0, c]/
Fold[(#1 q + #2) &, 0, d]] // Evaluate,
True,
With[{q = Sqrt[-2 Log[u]]},
Fold[(#1 q + #2) &, 0, c]/
Fold[(#1 q + #2) &, 0, d]] // Evaluate]]]];
You can use CompiledFunctionTools`CompilePrint[]
to check that the function was compiled properly.
Here is a plot of the absolute error over $(0,1)$:
nq[p_] = InverseCDF[NormalDistribution[], p];
Plot[AcklamQuantile[p] - nq[p], {p, 0, 1}, PlotRange -> All]

Of course, one can polish this further with a few steps of Newton-Raphson or Halley, if seen fit.
As it turns out, however, some spelunking shows that Mathematica does provide a compilable implementation of the normal distribution quantile. The undocumented function SpecialFunctions`Probit[]
is the routine of interest, but there is a minor botch in its implementation that I will show how to fix.
Spelunking the code of SpecialFunctions`Probit[]
shows that the function relies on a compiled function, System`StatisticalFunctionsDump`CompiledProbit[]
for numerical implementation. Using CompilePrint[]
to inspect its code, however, one sees a number of unsightly MainEvaluate
calls in the code. Thus, I offer a fix to this routine that yields a fully compiled function:
SpecialFunctions`Probit; (* force autoload *)
probit = With[{d = N[17/40],
cf1 = System`StatisticalFunctionsDump`CompiledProbitCentralMinimax,
cf2 = System`StatisticalFunctionsDump`CompiledProbitAsymptotic},
Compile[{{u, _Real}}, If[Abs[u - 0.5] <= d, cf1[u], cf2[u]],
CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True},
RuntimeAttributes -> {Listable}]];
where the two compiled sub-functions System`StatisticalFunctionsDump`CompiledProbitCentralMinimax
and System`StatisticalFunctionsDump`CompiledProbitAsymptotic
implement different approximations; for hopefully obvious reasons, I will not be reproducing their code here.
(Thanks to Henrik Schumacher for suggesting a better reformulation.)
CompilePrint[probit]
will show that this version has compiled properly. Here is a plot of the absolute error:
Plot[probit[p] - nq[p], {p, 0, 1}, PlotRange -> All]

InverseErfc
is not compilable. That means that the compiled function will contain calls toMainEvaluate
which will slow down the execution severely. You can speed it up a little by using theListable
-attribute ofInverseCDF
to process whole matrices at once:T = RandomReal[{0., 1.}, {10000, 64}]; { First@AbsoluteTiming[ P1 = Table[InverseCDF[NormalDistribution[0, 1], t], {t, T}];], First@AbsoluteTiming[P2 = InverseCDF[NormalDistribution[0, 1], T];] }
. $\endgroup$