On the solvability of the Diophantine equation $a(x^2-y^2)+2bxy=1$

Question

Take two integers $a,b$ with $\mathrm{GCD}(a,b)=1$ and $a$ odd. Consider the Diophantine equation $$ a(x^2-y^2)+2bxy=1,\qquad x,y\in\mathbb Z $$

Is there any efficient way to determine whether such an equation is solvable for a given $(a,b)$? For example, for $a=9$ and $b=10$, FindInstance is unacceptably slow:

FindInstance[9 (x^2 - y^2) + 2 10 x y == 1, {x, y}, Integers] != {} // AbsoluteTiming
(* around 9 seconds on my laptop *)

I am interested in pairs $(a,b)$ much larger than $\sim 10$. Surely there must be a better, much more specific way to solve this.

Answer 1

To determine if a solution exists, it suffices to look at convergents of a quadratic irrational.

SolutionQ[-1|1, _] = True;
SolutionQ[0, _] = False;
SolutionQ[_, 0] = False;

SolutionQ[a_, b_] /; !CoprimeQ[a, 2b] = False;

SolutionQ[a_, b_] := 
  Or[
    AnyTrue[Most[Convergents[(-b + Sqrt[b^2+a^2])/a]], solQ[a, b]],
    AnyTrue[Most[Convergents[(-b - Sqrt[b^2+a^2])/a]], solQ[a, b]]
  ]

solQ[a_, b_][r_] := 
  With[{x = Numerator[r], y = Denominator[r]}, 
    a(x^2 - y^2) + 2 b x y == 1
  ]

Some comparisons:

SolutionQ[9, 10] // AbsoluteTiming
SolutionQ[13, 9] // AbsoluteTiming
SolutionQ[99, 100] // AbsoluteTiming
SolutionQ[997, 1100] // AbsoluteTiming
SolutionQ[99711, 11011] // AbsoluteTiming

{0.00199, True}
{0.002114, False}
{0.009289, True}
{0.03192, True}
{0.743921, True}

FindInstance[9 (x^2 - y^2) + 2*10 x y == 1, {x, y}, Integers] // AbsoluteTiming
FindInstance[13 (x^2 - y^2) + 2*9 x y == 1, {x, y}, Integers] // AbsoluteTiming
TimeConstrained[FindInstance[99 (x^2 - y^2) + 2*100 x y == 1, {x, y}, Integers], 300]
TimeConstrained[FindInstance[997 (x^2 - y^2) + 2*1100 x y == 1, {x, y}, Integers], 300]

{7.4759, {{x -> -8485, y -> 3256}}}
{0.324242, {}}
$Aborted
$Aborted

Answer 2

The following code relies on Reduce, which is the rate-determining step and can be slow. Nevertheless, it solves $a=99$ and $b=100$ in about half a second. It solves $a=997$ and $b=1100$ in about 0.8 s. Remove "=!= {}" to see a few of the smallest solutions, if any.

AccidentalSolutionQ[a_, b_] :=
   Block[{r, s, w, y},
      r = Flatten[Table[
        Simplify[{ToRules[
          Reduce[{w^2 - 4 (a^2 + b^2) y^2 == 4 a}, {w, y}, Integers] /. C[1] -> k]}],
        {k, -1, 1}], 1];
      If[r == {},
         False,
         s = {w, y} /. r;
         Map[{(#[[1]] - 2 b #[[2]])/(2 a), #[[2]]} &, 
             Pick[s, Map[IntegerQ[(#[[1]] - (2 b)*#[[2]])/(2 a)] &, s]]] =!= {}
      ]]