Question: How do I solve 2nd degree Diophantine equation?

eq := 61*x^2 + 1 = y^2;
sols := (isolve(eq) assuming (x::posint, y::posint)):
for sol in sols do
  print(subs(_Z1 = 1, sol));
end do;

isolve returns the solution in terms of _Z1, indicating there may be infinite number of solutions?
I think {x = 226153980, y = 1766319049} is the only solution with x::posint, y::posint.

How do I restrict this to only positive integers in the solution?

Please Wait...