Solutions to Diophantine Equation $x^2 - D y^2 = m^2$ I'm looking for solutions to an equation of form $m =\sqrt{x^2 - Dy^2}$. I know that $m$ is a positive integer and so the inside of the square root has to be complete square. So I'm stuck with this diophantine equation $x^2 - Dy^2 = m^2$. 
In the case where $m = \pm1$ that's Pell Equation which I know how to find solutions to. I also know that I can find infinitely many solutions because solutions to $x^2 - Dy^2 = 1$ are also solutions to $(mx)^2 - D(my)^2 = m^2$ but that doesn't capture all solutions.
This site seems to derive the solutions by solving $s^2 - D = 0$ mod $m^2$ but I don't know why and I don't know what to do with it.
Edit: In the methods section, the hyperbolic case the author describes what seems to be this case, but I'm having a hard time following his steps.
Edit2: In fact, he makes the substitution $x = sy - m^2z$ and arrives at $\frac{s^2-D}{m^2}y^2 + (2z)y + (m^2z^2) = 1$ (hence the desire to find $s = D$ mod $m^2$). Fhen he jumps to the claim that a continuous fraction expansion of roots for $\frac{s^2 - D}{m^2}t + (2s)t + m^2 = 0$ will give answers to $z$ and $y$. I guess I'm having a brain-fart on why that might be.
Most papers I find on the topic say they have a solution if $\sqrt{D} > m^2$ but that's not sufficient here. 
For example: $m = 7, D = 8$ ->  $x^2 - 8y^2 = 49$ and I want to be able to arrive at $x = 11; y = 3$ and other solutions where $x$ and $y$ aren't $0$ mod $7$.
Any help is appreciated. Thank you.
 A: The main book you want is Binary Quadratic Forms: Classical Theory and Modern Computations by Duncan A. Buell. He does not actually quote Lagrange's theorem on bounds, for that Introduction to the Theory of Numbers by Leonard E. Dickson. On page 111 of Dickson, we have Theorem 85 due to Lagrange, any primitively represented value $n$ occurs as the first coefficient of a reduced form in the cycle as long as $|n| < \frac{1}{2} \sqrt \Delta.$ 
From your edit after Gerry's last suggestion, it seems you may expect to figure this out without looking at any books. Although I give a complete algorithm below, you still will need to do some reading.
Please see my Numbers representable as $x^2 + 2y^2$ , I do not wish to retype everything. 
Note: to actually find $\sqrt 8 \pmod p,$ see TONELLI or that other one
CIPOLLA. Or Berlekamp's factoring modulo primes, applied to $w^2 - 8.$
Other note: for an odd prime $q$ such that $(8|q)= -1,$ then whenever $x^2 - 8 y^2 \equiv 0 \pmod q,$ it follows that $x,y \equiv 0 \pmod q,$ so $x^2 - 8 y^2 \equiv 0 \pmod {q^2},$ and $(x/q)^2 - 8 (y/q)^2 $ is a smaller integer (in absolute value anyway).
Every odd prime $p$ such that $(8|p )= 1$ can be represented by either $x^2 - 8 y^2$ or $8 x^2 - y^2.$ These are in distinct genera. 
If $p \equiv 1 \pmod 8, $ we can solve $w^2 \equiv 8 \pmod p$ and $p t - w^2 = -8.$  The indefinite quadratic form $\langle p, 2 w, t \rangle$ can be reduced to $x^2 - 8 y^2$ by a matrix in $SL_2 \mathbb Z.$ The inverse of that matrix tells us how to write one solution to $\delta^2 - 8 \gamma^2 = p.$   Then
$$  (\delta^2 + 8 \gamma^2)^2 - 8 (2 \delta \gamma)^2 = p^2.    $$
If $p \equiv 7 \pmod 8, $ we can solve $w^2 \equiv 8 \pmod p$ and $p t - w^2 = -8.$  The indefinite quadratic form $\langle p, 2 w, t \rangle$ can be reduced to $8x^2 -  y^2$ by a matrix in $SL_2 \mathbb Z.$ The inverese of that matirx tells us how to write one solution to $8 \delta^2 -  \gamma^2 = p,$ or $  \gamma^2 - 8 \delta^2   = -p.$   Then
$$  (\gamma^2 + 8 \delta^2)^2 - 8 (2 \delta \gamma)^2 = (-p)^2 = p^2.    $$
Given any solution $x^2 - 8 y^2 = q, $ we get a new solution by negating either $x$ or $y.$ Other then that, infinitely many new solutions are created by multiplying by the automorph as in
$$
 \left(  \begin{array}{cc}
  3  &  8  \\
  1  &  3  
\end{array} 
  \right)  
 \left(  \begin{array}{c}
  x    \\
  y   
\end{array} 
  \right). 
$$
Written as a horizontal formula
$$ (3 x + 8 y)^2 - 8 (x + 3 y)^2 = x^2 - 8 y^2.  $$
You can then multiply by the automorph again and again. Note that, because there is no $xy$ term, we can be satisfied with negating entries as we like and do not need to introduce the inverse of the automorph, but it does no harm if you also use it.
I will need to think a bit about completeness for this process. But experiment with it for a while, see what you can do. 
POSTSCRIPT: Note that $\langle 1,0,-8 \rangle$ is not actually reduced in the sense of Gauss et al. However, it is one simple step away from $\langle 1,4,-4 \rangle,$ which is reduced.
=-=-=-=-=-=-=
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle
Input three coefficients a b c for indef f(x,y)= a x^2 + b x y + c y^2 
1  0 -8

  0  form              1           0          -8  delta      0
  1  form             -8           0           1  delta      2
  2  form              1           4          -4


          -1          -2
           0          -1

To Return  
          -1           2
           0          -1

0  form   1 4 -4   delta  -1
1  form   -4 4 1   delta  4
2  form   1 4 -4
minimum was   1rep 1 0 disc   32 dSqrt 5.6568542495  M_Ratio  32
Automorph, written on right of Gram matrix:  
-1  -4
-1  -5
 Trace:  -6   gcd(a21, a22 - a11, a12) : 1
=========================================
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ 

=-=-=-=-=-=-=
A: Various guises,
$$ s^2 + 4 s t + 4 t^2 = (s + 2t)^2 - 8 t^2,    $$
$$  -u^2 + 4 u v + 4 v^2 = 8 v^2 - (u - 2v)^2,   $$
and so on as needed.
