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]
```

1Unfortunately, this is not possible, at least not without further effort. The function

`InverseErfc`

is not compilable. That means that the compiled function will contain calls to`MainEvaluate`

which will slow down the execution severely. You can speed it up a little by using the`Listable`

-attribute of`InverseCDF`

to process wholematricesat 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];] }`

. – Henrik Schumacher – 2018-09-02T12:48:33.340