% Monte Carlo estimate of the number of integer solutions <x,y,z>
% of 0 <= x,y,z <= N, x^{2}+y^{2}+z^{2} is prime.
function C = MCPrimes(N,T)
Count = 0;
for (i=1:T)
  x=floor(rand*N);
  y=floor(rand*N);
  z=floor(rand*N);
  w=x^2+x+y^2+y+z^2+z+1;
  Count=Count+IsPrime(w);
end
Count
C= Count*N^3/T;
end

function P = IsPrime(Q)
P = 1;
for (i=2:sqrt(Q))
  if (abs(Q-round(Q/i)*i) < 1/Q) 
    P=0;
    break;
  end
end
end
