%generate all sums of two squares up to 10^8 = (10^4)^2.
t = 10^4;
a(1,1:(t+1)) = (0:t).^2;
temp = repmat(a,(t+1),1);
sumOfSquares = temp + temp';
u2 = sort(unique(sumOfSquares));
b1 = find(u2>t^2);
u1 = u2(1:(b1(1)-1));

%generate nice prime moduli, from 1 to 10^8, with (roughly) constant spacings on a logarithmic scale.
q1 = primes(t^2)';
q2 = zeros(floor(log(size(q1,1))/log(1.1)),1);
for j=1:floor(log(size(q1,1))/log(1.1))
    q2(j,1) = q1(round((1.1)^j));
end
q = unique(q2);
qcount = zeros(size(q,1),1);

%compute variance in arithmetic progressions.
parfor k = 1:size(q,1)
    buckets = zeros(q(k,1),1);
    for a = 1:size(u1,1)
        buckets(mod(u1(a,1), q(k,1))+1,1) = buckets(mod(u1(a,1), q(k,1))+1,1) + 1;
    end
    b1 = buckets(2:end,1);
    normq = t^2/((sqrt(log(t^2)))*(q(k,1)-1));
    qcount(k,1) = var(b1)/normq;
end
x_axis = log(q)/(log(t^2));

save('results.mat', 'x_axis', 'qcount')