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

Question

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

Answer 1

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 CompiledFunctionTools`CompilePrint[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;