% maxprime must be „11,  m¾n
function  z=tplus7(m,n,maxprime)
format compact
if maxprime<11,error('maxprime must be >10'), end
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,840);
x=[(m-mp+1):840:n, (m-mp+121):840:n  (m-mp+361):840:n];
x2=[(m-mp+169):840:n  (m-mp+289):840:n  (m-mp+529):840:n];
x=[x x2]; clear x2
% ArrayLength=length(x)
y=primes(maxprime); y(1:4)=[]; % primes 2,3,5,7  are removed from the list  y
% disp(' PRIMES:')
for p=y
  p1=(p-1)/2; %disp(p);
  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; return, end
  %if length(x)<6, x, disp(' PRIMES:'), end
end
z=x(x>=m);
z=sort(z);
