% Monte Carlo estimate of pi: comparisons and logicals fuse too, so the
% mask never touches memory -- one fused pass plus the reduction.
n = 10000000;
tic;
x = 2*rand(n, 1) - 1;
y = 2*rand(n, 1) - 1;
pi_est = 4*mean((x.^2 + y.^2) <= 1);
toc
fprintf('pi ~ %.6f  (n = %d)\n', pi_est, n);