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$
=-=-=-=-=-=-=
There is something wrong here. The given set of formulas does give infinitely many solutions, that part is fine.
Example:
$$ x^2 + 8 xy + y^2 = z^2 $$
Set of solutions that does not fit those formulas:
$$ x = 16 u^2 - 14 u v + 3 v^2, \; \; y = -2 u^2 + 2 u v, \; \; z = 2 u^2 + 6 u v - 3 v^2 $$
I know this is new because the indefinite binary form $2 u^2 + 6 u v - 3 v^2$ is neither the principal form nor its negative: it does not represent $+1$ or $-1$ over the integers.
I have requested two books by Andreescu through my city library, a local college has them. I cannot imagine that they claim all solutions come up through one set of formulas.
pari
? x = 16 * u^2 - 14 * u * v + 3 * v^2
%1 = 16*u^2 - 14*v*u + 3*v^2
? y = -2 * u^2 + 2 * u * v
%2 = -2*u^2 + 2*v*u
? z = 2 * u^2 + 6 * u * v - 3 * v^2
%3 = 2*u^2 + 6*v*u - 3*v^2
? x^2 + 8 * x * y + y^2 - z^2
%4 = 0
==========================================================
jagy@phobeusjunior:~$ ./Conway_Positive_Primes 1 8 1 1000 5
1 8 1 original form
1 6 -6 Lagrange-Gauss reduced
Represented (positive) primes up to 1000
61 109 181 229 241 349 409 421 541 601
661 709 769 829
==========================================================
jagy@phobeusjunior:~$ ./Conway_Positive_Primes 2 6 -3 1000 5
2 6 -3 original form
2 6 -3 Lagrange-Gauss reduced
Represented (positive) primes up to 1000
2 5 17 53 113 137 173 197 233 257
293 317 353 557 593 617 653 677 773 797
857 953 977
Best Answer
Let $u$ be an element of $\mathbb{Z}[\sqrt 5]$ of norm 1, i.e. $u = r + s \sqrt 5$ with $r^2-5s^2 = 1$.
The multiplication by $u$ in $\mathbb{Z}[\sqrt 5]$ turns any element $y$ of norm $44$ into another element $uy$ of norm $44$. View this multiplication operation on $\mathbb{Z}[\sqrt 5]$ as the transformation of the plane $f : (p,q) \rightarrow (pr+5qs,ps+qr)$, and look for its eigenvalues :
$f(\sqrt5,1) = (r\sqrt5+5s,r+\sqrt5s) = u(\sqrt5,1)$, and we have $f(- \sqrt5,1) = \frac 1u (- \sqrt5,1)$ as well.
If $u>1$ this means that $f$ is an operation that, when iterated, takes elements near the line $(p = - \sqrt5 q)$ and moves them over to the line $(p = \sqrt5 q)$ Now you want to find a sector of the plane so that you can reach the whole plane by taking its images by the iterates of $f$ and $f^{-1}$
Define $g(p,q) = \frac {p + \sqrt5 q}{p - \sqrt5 q}$, which is the ratio of the coordinates of $(p,q)$ in the eigenbasis of $f$. $g(f(p,q)) = \frac {pr+5qs + \sqrt5 (ps+qr)}{pr+5qs - \sqrt5 (ps+qr)} = \frac{(r+\sqrt5 s)(p + \sqrt5 q)}{(r-\sqrt5 s)(p - \sqrt5 q)} = (r+\sqrt5 s)^2 g(p,q)$.
Or alternately, define $g(y) = y/\overline{y}$, so that $g(uy) = uy/\overline{uy} = u^2 g(y)$.
Thus if you look at any point $(p,q)$, you know you can apply $f$ or $f^{-1}$ to turn it into $(p',q')$ such that $g(p',q') \in [1 ; u^2[$
Thus, a suitable sector of the plane is the set of points $(p,q)$ such that $g(p,q) \in [1 ; u^2[$ : if you find all the elements $y$ of norm $44$ such that $g(y) \in [1 ; u^2[$, then this means that the $u^ky$ will cover all the elements of norm $44$
Finally, the good thing is that $ \{y \in \mathbb{Z}[ \sqrt 5] / g(y) \in [1; u^2[, y\overline{y} \in [0; M] \}$ is a finite set, so a finite computation can give you all the elements of norm $44$ you need.
In the case of $p²-10q²=9$, a fundamental unit is $u = 19+6\sqrt{10}$, so replace $\sqrt 5,r,s$ with $ \sqrt {10},19,6$ in everything I wrote above.
In order to find all the solutions, you only need to check potential solutions in the sector of the plane between the lines $g(p,q) = 1$ and $g(p,q) = u^2$.
You can look at the intersection of the line $g(p,q)=1$ with the curve $p^2-10q^2 = 9$. $g(p,q)=1$ implies that $p+\sqrt{10}q = p- \sqrt{10}q$, so $q=0$, and then the second equations has two solutions $p=3$ and $p= -3$. It so happens that the intersection points have integer coordinates so they give solutions to the original equation.
Next, the intersection of the line $g(p,q) = u^2$ with the curve will be $u \times (3,0) = f(3,0) = (19*3+60*0, 6*3+19*0) = (57,18)$ and $u \times (-3,0) = (-57,-18)$.
So you only have to look for points on the curve $p^2-10q^2=9$ with integers coordinates in the section of the curve between $(3,0)$ and $(57,18)$ (and the one between $(-3,0)$ and $(-57,-18)$ but it is essentially the same thing).
You can write a naïve program :
for q = 0 to 17 do :
let square_of_p = 9+10*q*q.
if square_of_p is a square, then add (sqrt(square_of_p),q) to the list of solutions.
Which will give you the list $\{(3,0) ; (7,2) ; (13,4)\}$. These three solutions, together with their opposite, will generate, using the forward and backward interations of the function $f$, all the solution in $\mathbb{Z}^2$.
If you only want solution with positive coordinates, the forward iteration of $f$ on those three solutions are enough.
Also, as Gerry points out, the conjugate of $(7,2)$ generates $(13,4)$ because $f(7,-2) = (13,4)$. Had we picked a sector of the plane symmetric around the $x$-axis, we could have halved the search space thanks to that symmetry, and we would have obtained $\{(7,-2),(3,0),(7,2)\}$ instead.
One loop of this hypnotic animation represents one application of the function $f$. Each dot corresponds to one point of the plane with integer coordinates, and is moved to its image by $f$ in the course of the loop. The points are colored according to their norm (and as you can see, each of them stay on their hyperbolic branch of points sharing their norm), and I've made the yellow-ish points of norm 9 (the solutions of $x^2-10y^2 = 9$) a bit bigger. For example, the point at (3,0) is sent outside the graph, and the point at (-7,2) is sent on (13,4) (almost vanishing).
You can see that there are three points going through (3,0) during the course of one loop. They correspond to three representants of the three fundamental solutions of the equation. For each yellowish point on the curve $x^2-10y^2=9$, no matter how far along the asymptote it may be, there is an iterate of $f$ or $f^{-1}$ that sends it to one of those three fundamental solutions.
In order to find all fundamental solutions, it is enough to explore only a fundamental portion of the curve (a portion whose iterates by $f$ covers the curve), for example the fundamental portion of the curve between (-7,2) and its image by $f$, (13,4). To find the solutions on that portion, you set $y=-2,-1,0,1,2,3$ and look if there is an integer $x$ that makes a solution for each of those $y$.
Whichever fundamental portion of the curve you choose, you will find 3 solutions inside it, whose images by $f$ are sent to the next three solutions in the next portion of the curve, and so on.
Now there is a better procedure than the "brute search" I did to get all the solutions. It is an adaptation of the procedure to obtain a fundamental unit :
Start with the equation $x^2-10y^2 = 9$, and suppose we want all the positive solutions.
We observe that we must have $x > 3y$, or else $-y^2 \ge 9$, which is clearly impossible.
So, replace $x$ with $x_1 + 3y$.
We get the equation $x_1^2 + 6x_1 y - y^2 = 9$.
We observe that we must have $y > 6x_1$, or else $x_1^2 \le 9$.
In this case we quickly get the three small solutions $(x_1,y) = (1,2),(1,4),(3,0)$ which correspond to the solutions $(x,y) = (7,2),(13,4),(3,0)$.
Otherwise, continue and replace $y$ with $y_1 + 6x_1$.
We get the equation $x_1^2 - 6x_1y_1 - y_1^2 = 9$.
We observe that we must have $x_1 > 6y_1$, or else $-y_1^2 \ge 9$, which is clearly impossible.
So, replace $x_1$ with $x_2 + 6y_1$.
We get the equation $x_2^2 + 6x_2y_1 - y_1^2 = 9$.
But we already encountered that equation so we know how to solve it.