Computing the double Integral using MonteCarlo techniques using Julia

I decided to try and learn Julia for doing scientific computing, and I decided to tackle the problem of finding

$$ \int_{D_{\frac{1}{4}}} x^4 + y^2 dA $$

where \$ D_{\frac{1}{4}} \$ is the part of the unit circle in the first cuadrant.

My code in Julia is the following:

using Distributions
e = 10.0^(-3);
p = 0.85;
variance = 4;

N = floor(Int, variance / ((1-p)*((e/2)^2))) + 1

u = Uniform(0,2);
x = rand(N);
y = rand(N);
z = rand(u, N);

result = sum((x.^2 + y.^2 .<= 1) & (z .<= x.^4 + y.^2))*2.0 / N

which gives the nice result \$ = 0.2945746303294543 \$

I kindly ask for how to improve my implementation, and reduce the footprint of memory (it uses almost 2 to 3gb in RAM).

Try swapping out the three arrays with calling a function iteratively---three Int64 arrays (each element 8 bytes) of ~100M elements is a lot. That is, switch out the body of the sum into a separate function:

function foo(u)
    x = rand()
    y = rand()
    z = rand(u)
    (x^2 + y^2 <= 1) & (z <= x^4 + y^2) ? 1 : 0
end

and replace the sum expression itself with just a loop:

count = 0
for i in 1:N
    count += foo(u)
end
count * 2.0 / N

This runs a bit faster, since the compiler can optimize the internal function, and with a much smaller memory footprint, since we sidestep the arrays entirely.

StackExchange Code Review Q#119998, answer score: 6