# How to compile inverse error function?

I have the following code:

function[u_] = InverseErf[2*u - 1]
comp = Hold@Compile[{{u, _Real}}, function[u]] /.
DownValues@function // ReleaseHold
<< CompiledFunctionTools
CompilePrint@comp


However, it returns MainEvaluate when calling function. Probably, it is because InverseErf is a special function. Is there any way to compile it?

This question appeared from the task of smearing some quantity by Gaussian distribution. This is how I do it:

distff[\[Sigma]_] =
ProbabilityDistribution[
1/(Sqrt[2*Pi] \[Sigma])
Exp[-((x)^2/(2 \[Sigma]^2))] , {x, -Infinity, Infinity}];
inversecdf[\[Sigma]_, u_] =
InverseCDF[distff[\[Sigma]], u][[1]][[2]];
quantitySmeared=RandomReal[{0.5,2},10^5]+0.5*inversecdf[1, RandomReal[{0,1},10^5]];


The problem is that I need to perform this smearing many times for large lists; this is why I need to compile.

Update

It looks like, for this particular case, I may generate the smearing using the so-called Box-Muller transform, which avoids using InverseErf. However, I wonder whether some procedure allows compiling the special functions.

You could use OpenCLLink if you need to run this on a huge amount of data and want to use your GPU:

kernelSource = "//based on https://stackoverflow.com/a/49743348
float erfinv(float a)
{
float p, r, t;
t = fma(a, 0.0f - a, 1.0f);
t = log(t);
if (fabs(t) > 6.125f) {
p = 3.03697567e-10f;
p = fma(p, t,  2.93243101e-8f); //  0x1.f7c9aep-26
p = fma(p, t,  1.22150334e-6f); //  0x1.47e512p-20
p = fma(p, t,  2.84108955e-5f); //  0x1.dca7dep-16
p = fma(p, t,  3.93552968e-4f); //  0x1.9cab92p-12
p = fma(p, t,  3.02698812e-3f); //  0x1.8cc0dep-9
p = fma(p, t,  4.83185798e-3f); //  0x1.3ca920p-8
p = fma(p, t, -2.64646143e-1f); // -0x1.0eff66p-2
p = fma(p, t,  8.40016484e-1f); //  0x1.ae16a4p-1
} else {
p = 5.43877832e-9f;
p = fma(p, t,  1.43285448e-7f); //  0x1.33b402p-23
p = fma(p, t,  1.22774793e-6f); //  0x1.499232p-20
p = fma(p, t,  1.12963626e-7f); //  0x1.e52cd2p-24
p = fma(p, t, -5.61530760e-5f); // -0x1.d70bd0p-15
p = fma(p, t, -1.47697632e-4f); // -0x1.35be90p-13
p = fma(p, t,  2.31468678e-3f); //  0x1.2f6400p-9
p = fma(p, t,  1.15392581e-2f); //  0x1.7a1e50p-7
p = fma(p, t, -2.32015476e-1f); // -0x1.db2aeep-3
p = fma(p, t,  8.86226892e-1f); //  0x1.c5bf88p-1
}
r = a * p;
return r;
}

__kernel void bulkinverf(__global float * input, __global float * output)
{
int id = get_global_id(0);
output[id] = erfinv(input[id]);
}
";

"bulkinverf", {{"Float", "Input"}, {"Float", "Output"}}, {16}];

input = RandomReal[{-0.9999, 0.9999}, 10^8];
result = First@fun[input, ConstantArray[0, Length[input]]];

• Thanks, that's very interesting! Perhaps by using GPU, I may speed up all the routines I have many times... Commented Feb 2 at 11:51
• @JohnTaylor Though be warned, it is very difficult to control in Mathematica's OpenCLLink except for simpler cases like this. I'd recommend trying things in this order first to see at what point you're satisfied with performance: Pure Mathematica < Compile < call out to a C library < OpenCLLink , because the further you get away from pure Mathematica the harder things get in my experience. Commented Feb 2 at 11:57

You can port one of the C algorithms for erfinv to Mathematica.

Here, I've used @njuffa's code for erfinvf, which is a single-precision variant.

inverseErf = Compile[{{x, _Real}}, Module[{p, r, t},
t = 1. - x^2;
t = Log[t];
If[Abs[t] > 6.125, p = 3.03697567*^-10;
p = p*t + 2.93243101*^-8;
p = p*t + 1.22150334*^-6;
p = p*t + 2.84108955*^-5;
p = p*t + 3.93552968*^-4;
p = p*t + 3.02698812*^-3;
p = p*t + 4.83185798*^-3;
p = p*t - 2.64646143*^-1;
p = p*t + 8.40016484*^-1;,
p = 5.43877832*^-9;
p = p*t + 1.43285448*^-7;
p = p*t + 1.22774793*^-6;
p = p*t + 1.12963626*^-7;
p = p*t - 5.61530760*^-5;
p = p*t - 1.47697632*^-4;
p = p*t + 2.31468678*^-3;
p = p*t + 1.15392581*^-2;
p = p*t - 2.32015476*^-1;
p = p*t + 8.86226892*^-1;
];
x*p
], RuntimeAttributes -> {Listable}, CompilationTarget -> "C"
]

Plot[InverseErf[x] - inverseErf[x], {x, -1, 1}, PlotRange -> All]