function x=quadres(x0,m)
%first 30000 points of the quadratic residue generator
x=x0;
q=min(30000,m);
for i=1:q
   y=rem(x(i)^2,m);
   x=[x y];
end
%p=1:q;
%p=p(x==x(1));
%['period ' num2str(p)]