function out = elkies(n, m)
% generate 4 random numbers in interval [1,n] and
% determine prob(gcd(a,b)== gcd(c,d)) by m trials
succ = 0;
for j = 1: m
    a = ceil(rand()*n);
    b = ceil(rand()*n);
    c = ceil(rand()*n);
    d = ceil(rand()*n);
    if gcd(a,b)==gcd(c,d) succ = succ+1;
    end;
end;
out = succ/m;