The problem of characterizing those D for which the Diophantine equation
X
2
−
D
Y
2
=
−
1
{X^2}\, - \,D{Y^2}\, = \, - \,1
is solvable has been studied for two hundred years. This paper considers this problem from the viewpoint of determining the computational complexity of recognizing such D. For a given D, one can decide the solvability or unsolvability of
X
2
−
D
Y
2
=
−
1
{X^2}\, - \,D{Y^2}\, = \, - \,1
using the ordinary continued fraction expansion of
D
\sqrt D
, but for certain D this requires more than
1
3
D
(
log
D
)
−
1
\tfrac {1}{3}\,\sqrt D \,{(\log \,D)^{ - \,1}}
computational operations. This paper presents a new algorithm for answering this question and proves that this algorithm always runs to completion in
O
(
D
1
/
4
+
ε
)
O({D^{1/4\, + \,\varepsilon }})
bit operations. If the input to this algorithm includes a complete prime factorization of D and a quadratic nonresidue
n
i
{n_i}
for each prime
p
i
{p_i}
dividing D, then this algorithm is guaranteed to run to completion in
O
(
(
log
D
)
5
(
log
log
D
)
(
log
log
log
D
)
)
O({(\log \,D)^5}\,(\log \,\log \,D)(\log \,\log \,\log \,D))
bit operations. This algorithm is based on an algorithm that finds a basis of forms for the 2-Sylow subgroup of the class group of binary quadratic forms of determinant D.