# Why is random sampling with ordinary and compiled code so different?

Consider the following code:

function[x_] = -0.2135096237485963 +
0.5230330267652104 Power[
27.95584680618265 x +
Sqrt[(27.95584680618265 x - 9.569316101220256)^2 +
42.46540146917953] - 9.569316101220256, (3)^-1] -
1.8247674462728896/Power[
27.95584680618265 x +
Sqrt[(27.95584680618265 x - 9.569316101220256)^2 +
42.46540146917953] - 9.569316101220256, (3)^-1];
randvalscomp =
Hold@Compile[{{n, _Integer}},
ArcCos[function[RandomReal[{0, 1}, n]]],
CompilationTarget -> "C"] /. DownValues@function //
ReleaseHold;


Now, let us generate 100 random numbers using uncompiled code and compiled code:

randvals1 = ArcCos[function[RandomReal[{0, 1}, 100]]];
randvals2 = randvalscomp[100];


By its structure (the way I have obtained it), function may only return the values between -1 and 1, so ArcCos must be real. While randvals1 contains only real numbers, randvals2 sometimes returns an arbitrarily large imaginary part. It does not look like a machine precision error: say, ArcCos[-1.00001] returns a vanishingly small imaginary part plus a real part - either Pi, while I may get arbitrary combinations:

Select[randvals2,Im[#]!=0&]


{1.24622 -1.19564 I,2.43982 -0.7189 I,2.36874 -0.814092 I,1.1946 -1.06941 I,2.4617 -0.514546 I,1.22451 -1.21208 I}

What may be the reason for this and how to fix it?

P.S. function may be obtained as

function1[u_,sin_]=(1/(1+8 sin^2-4 sin))(-((8 Power[2, (3)^-1] (1-2 sin)^2 sin^2)/Power[u+512 u sin^6-768 u sin^5+576 u sin^4-256 u sin^3+72 u sin^2+Sqrt[(u+256 (2 u-1) sin^6-192 (4 u-1) sin^5+48 (12 u-1) sin^4-256 u sin^3+72 u sin^2-12 u sin)^2+256 (1-2 sin)^6 sin^6]-12 u sin-256 sin^6+192 sin^5-48 sin^4, (3)^-1])+2^(2/3) Power[u+512 u sin^6-768 u sin^5+576 u sin^4-256 u sin^3+72 u sin^2+Sqrt[(u+256 (2 u-1) sin^6-192 (4 u-1) sin^5+48 (12 u-1) sin^4-256 u sin^3+72 u sin^2-12 u sin)^2+256 (1-2 sin)^6 sin^6]-12 u sin-256 sin^6+192 sin^5-48 sin^4, (3)^-1]-1+4 sin)


by replacing sin with 0.223: function[x_]=function1[x,0.223].

• Don't you get a warning message when running randwalscomp? CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. Commented Apr 30 at 14:42
• @Domen : yes, but this is because I get complex numbers when using function. But why I am getting them? Without the compilation, the code works properly. Commented Apr 30 at 14:43

Your hackery way of injection function definitions into Compile is the root of your troubles. Observe the output of the following code:

Hold@compile[{{n, _Integer}}, function[randomReal[{0, 1}, n]]] /.
DownValues@function // ReleaseHold


You will see that randomReal is inserted at every appearance of x in function, which is of course a non-sense, because it should be "the same" x every time. You can also see this by looking at the output of CompiledFunctionToolsCompilePrint[randvalscomp] and seeing multiple occurences of RandomReals.

A solution is trivial: Sample random points only once before applying function.

randvalscomp = Hold@Compile[{{n, _Integer}},
Module[{xs = RandomReal[{0, 1}, n]}, ArcCos[function[xs]]],
CompilationTarget -> "C"] /. DownValues@function // ReleaseHold;