intrinsic HasSolution(a::RngIntElt, b::RngIntElt, yrange::SeqEnum) -> RngIntElt
{}
   if a eq 0 or a eq 1 or b eq 0 or b eq 1 or GCD(a,b) ne 1 then 
      return 0,0;
   end if;
   for y in [z : z in yrange |z ne 0 and GCD(a,z) eq 1] do
      disc := b^2*y^2 - 4*(b/a*y^2-1/a);
      val, s := IsSquare(disc);
      if val then
         x := (-b*y + s) / 2;
         if x ne 0 and x in Integers() then
            return x, y;
         end if;
      end if;
   end for;
   return 0,0;
end intrinsic;

intrinsic DoIt(aRange::SeqEnum, range::RngIntElt)
{}
   for a in aRange do 
      t:=Cputime();
      "Trying a=",a;
      for b in [-range..range] do 
         x,y := HasSolution(a,b,[1..range]); 
         if x ne 0 then     
            printf "%o\t%o\t%o\t%o\n", a,b,x,y;
         end if;
      end for; 
      "time =", Cputime(t);
   end for;
end intrinsic;
