Here is a very strong and impressive result of the Reduce command.
Reduce[a^2 + b^2 == 841*(a*b + 1), {a, b}, PositiveIntegers]
performs (a == 24389 && b == 29) and 6 infinite series of the solutions, one of these is (The results are too big for $\LaTeX$.)
(C[1] \[Element] Integers && C[1] >= 0 && a == (1/1414554)(-(24389/ 2) (-707277 (297410399 - 353640 Sqrt[707277])^C[1] + 841 Sqrt[707277] (297410399 - 353640 Sqrt[707277])^C[1] - 707277 (297410399 + 353640 Sqrt[707277])^C[1] - 841 Sqrt[707277] (297410399 + 353640 Sqrt[707277])^C[1]) - 20511033/ 2 (-841 (297410399 - 353640 Sqrt[707277])^C[1] + Sqrt[707277] (297410399 - 353640 Sqrt[707277])^C[1] - 841 (297410399 + 353640 Sqrt[707277])^C[1] - Sqrt[707277] (297410399 + 353640 Sqrt[707277])^C[1])) && b == -(1/1414554) 29 (-707277 (297410399 - 353640 Sqrt[707277])^C[1] + 841 Sqrt[707277] (297410399 - 353640 Sqrt[707277])^C[1] - 707277 (297410399 + 353640 Sqrt[707277])^C[1] - 841 Sqrt[707277] (297410399 + 353640 Sqrt[707277])^C[1]))
This is problem 16 from Kvant 2020, # 10, p. 42 (in Russian). Its solution is given in p. 60 ibid .
However, no explicit formulas are presented there, but only an algorithm to obtain all the solutions in PositiveIntegers.
My questions are: are all the solutions produced by Reduce? how to establish it? I hope the answer to the first question is yes.
Up to the cited algorithm, the solutions $a,b$ (provided $a\le b$) are those and only those pairs of numbers $(a,b)$ that for each $n \in \mathbb N $ are calculated by the formulas:$a_n=b_{n-1},\, b_n=k^2b_{n-1}-a_{n-1},\,a_0=0,\,b_0=k$, here $k=29$.
The
RSolvecommand of Mathematica cracks it:$$a(n)=-\frac{29 \left(2 \left(\sqrt{707277}+841\right)\right)^{-n} \left(4^n-\left(\sqrt{707277}+841\right)^{2 n}\right)}{\sqrt{707277}},$$ $$b(n)=29 U_n\left(\frac{841}{2}\right). $$ in the above $U_n$ stands for Chebyshev polynomials. It remains to compare this solution with the solution found by
Reduce.