Maple: RNG is not random

Question

i was "finding Pi" with Monte Carlo Method, but the answer was incorrect. The oryginal code was:

RandomTools[MersenneTwister]: with(Statistics):

tries := 10000:

s := 0; 
for i to tries do 
    if GenerateFloat()^2+GenerateFloat()^2 < 1 then s := s+1 end if; 
end do: 
evalf(4*s/tries)

It gives answer aroud 2.8-2.85

when I change the code to

s := 0; 
x := Array([seq(GenerateFloat(), i = 1 .. tries)]); 
y := Array([seq(GenerateFloat(), i = 1 .. tries)]); 
for i to tries do 
if x[i]^2+y[i]^2 < 1 then s := s+1 end if;
end do:
evalf(4*s/tries)

Then the answer is correct. I have no idea why i can't generate number in "for" loop.

I've founded that the mean of it is the same, but the variance is different. For:

tries := 100000; 
A := Array([seq(GenerateFloat(), i = 1 .. 2*tries)]); 
s1 := Array([seq(A[i]^2+A[tries+i]^2, i = 1 .. tries)]); 
Mean(s1); 
Variance(s1);
s2 := Array([seq(GenerateFloat()^2+GenerateFloat()^2, i = 1 .. tries)]);
Mean(s2); 
Variance(s2);

output is:

0.6702112097021581
0.17845439723457215
0.664707674135025
0.35463131700965245

What's wrong with it? GenerateFloat() should be as uniform as possible.

Answer 1

Automatic simplification is turning your,

GenerateFloat()^2+GenerateFloat()^2

into,

2*GenerateFloat()^2

before GenerateFloat() is evaluated.

One simple change to get it to work as you expected would be separate them. Eg,

restart:
with(RandomTools[MersenneTwister]):
tries := 10^4:
s := 0:
for i to tries do
  t1,t2 := GenerateFloat(),GenerateFloat();
  if t1^2+t2^2 < 1 then s := s+1 end if;
end do:
evalf(4*s/tries);

Another way is to use a slightly different construction which doesn't automatically simplify. Consider, single right quotes (uneval quotes) don't stop automatic simplification (which is a definition of the term if you want).

'f()^2 + f()^2';                                                
                                     2
                                2 f()

But the following does not automatically simplify,

a:=1:
'f()^2 + a*f()^2';
                                2        2
                             f()  + a f()

Therefore another easy workaround is,

restart:
with(RandomTools[MersenneTwister]):
tries := 10^4:
s := 0:
a := 1;
for i to tries do
  if GenerateFloat()^2 + a*GenerateFloat()^2 < 1 then s := s+1 end if;
end do:
evalf(4*s/tries);