Continuing from @Batominovski's answer,
the problem asks us to find all the values taking $\pm 1$ in the three following recurrent linear sequences satisfying $a_n = 5a_{n-1} - 7a_{n-2}$,
- $1,2,3,1,-16,-87,-323,\ldots$
- $0,1,5,18,55,149,360,\ldots$
- $1,3,8,19,39,62,37,\ldots$
Looking at the recurrence relation mod $18$, we have $a_n \equiv 5a_{n-1}+11a_{n-2} \equiv 5(5a_{n-2}+11a_{n-3}) + 11a_{n-2} \equiv a_{n-3}$.
Looking at the first three terms of each sequence modulo $18$, we can already discard $2/3$ of all the terms and also remove the possibility of a $-1$ appearing, which reduces the problem to finding the $1$ values in the linear recurrent sequences satisfying $a_n = 20a_{n-1}-343_{n-2}$
- $1,1,-323,-6803,\ldots$
- $1,55,757,-3725,\ldots$
- $1,19,37,-5777,\ldots$
Now each sequence is a linear combination of the coefficient sequences of $(2-\omega)^{3n} = (1-18\omega)^n = (10 - 9 \sqrt{-3})^n$, that is, the two sequences
- $1,10,-143,-6290,\ldots$
- $0,-9,-180,-513, \ldots$
(with integer coefficents in all three cases)
With the binomial theorem, we can expand $(1+9(1-\sqrt{-3}))^n = 1 + 9(1-\sqrt{-3})n + 81(1-\sqrt{-3})^2n(n-1)/2 + \ldots$.
If we place ourselves in $\Bbb Z_3[\sqrt{-3}]$ and we notice that $v_3(9^k/k!) >= 3k/2 \to \infty$, we can reorder the summation and get $1 + [9(1-\sqrt{-3}) + O(81)]n + O(81)n^2 + \ldots$ i.e. each coefficient sequence can be extended as a function $\Bbb Z_3 \mapsto \Bbb Z_3$ which is a power series in $n$ with coefficients in $\Bbb Z_3$
Now, if a power series is of the form $a_0 + a_1 n + a_2n^2$ with $|a_1| > |a_2|,|a_3|,|a_4|,\ldots$, then it is a bijection from $\Bbb Z_3$ to $a_0+a_1\Bbb Z_3$.
The $a_1$ coefficients modulo $81$ in our three sequences are $9-9=0, 9-9 \times (-5) = 54 \neq 0, 9-9 \times (-1) = 18 \neq 0$.
In the last two cases we deduce that the sequence is injective so the initial occurence of $1$ is the only one.
As for the first sequence, we have that both $a_1$ and $a_2$ are of the order of $3^4$.
But then its derivative is of the right form to show that it has a single zero at some $n_0 \in \Bbb Z_3$.
By setting $m = n-n_0$, we get a new power series of the form $b_0 + b_2m^2 + b_3m^3 + \ldots$ where $|b_2| > |b_3|,|b_4|,\ldots$. This shows that the power series is $2$-to-$1$ (with the exception of $0 \mapsto b_0$), so the two values of $1$ we already know are the only ones.
$$ 8x^2 - 4 y^2 - 4 y - 1 = -1 $$
$$ (2y+1)^2 - 8 x^2 = 1 $$
$$ w^2 - 8 x^2 = 1 \; , $$
so that $w$ really is odd, then $y = (w-1)/2$
The first two are
$$ (w,x) = (1,0) $$
then
$$ (w,x) = (3,1) $$
After that, growth comes from
$$ (w,x) \mapsto (3w+8x, w + 3x ) $$
over and over
w: 1 x: 0
w: 3 x: 1
w: 17 x: 6
w: 99 x: 35
w: 577 x: 204
w: 3363 x: 1189
w: 19601 x: 6930
w: 114243 x: 40391
w: 665857 x: 235416
w: 3880899 x: 1372105
w: 22619537 x: 7997214
but still followed by $y = (w-1)/2$ for each
If preferred, there are separate recurrences
$$ w_{n+2} = 6 w_{n+1} - w_n \; , \; $$
$$ x_{n+2} = 6 x_{n+1} - x_n \; . \; $$
Best Answer
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.
=-=-=-=-=-=-=
=-=-=-=-=-=-=