% maxprime must be „5
function  z=tplus3(m,n,maxprime)
format compact
tic
disp('   Screening for non-square QRs the numbers')
str1=num2str([m n]);str2=num2str(maxprime);
disp(['in the range: ',str1, '  (mod all the primes < ',str2,')'])
mp=mod(m,24);
x=(m-mp+1):24:n;
ArrayLength=length(x)
y=primes(maxprime); y(1:2)=[]; % primes 2,3  are removed from the list  y
disp(' PRIMES:')
for p=y
  disp(p);p1=(p-1)/2;
  if p==17, sqrtx=sqrt(x); x=x(round(sqrtx)~=sqrtx); end
  xxx=mod(x,p);
 %  Now producing the list  qr  of quadratic residues mod  p
  qr=mod((1:p1).^2,p);
  f=0*xxx;
  for kkk=1:p1
    f=f | (xxx==qr(kkk));
  end
  x=x(f);
  if isempty(x), z=x; toc, return, end
  if length(x)<10, x, disp(' PRIMES:'), end
end
z=x(x>=m);
z=sort(z);
toc

