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

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.

• On develop.wolframcloud.com still takes around: {9.85082,True} seconds. Nov 14, 2018 at 19:33
• Could you please explain, are you interested in a general analytic formula for any a and b or a general analytic formula for fixed a and b, or a random sample of solutions produced quickly? Nov 14, 2018 at 20:11
• @VitaliyKaurov Given a fixed pair of integers, say $a=99$ and $b=100$, I want to know whether the equation is solvable or not; the actual solution being irrelevant. And I want to know this as fast as possible: ideally, it should take it less than a second for pairs of order $a,b\sim 10^4$ or higher, to decide whether the equation is solvable or not. I believe there is some algorithm (Gauss? Legendre?) for this kind of polynomial equations, although I'm not sure. Nov 14, 2018 at 20:14
• @AccidentalFourierTransform understood, thank you ! Nov 14, 2018 at 20:16
• Strange it's so slow. There's a general approach for quadratic Diophantine equations. See this solver: alpertron.com.ar/QUAD.HTM Nov 14, 2018 at 22:24

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

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]]] =!= {}
]]
• Why not use the substitution x -> (w-b y)/a instead of x -> (w - 2 b y)/(2 a)? It seems that the resulting equation can be solved more quickly by Mathematica. Nov 15, 2018 at 0:29