Finding the square root of a polynomial with radical I have two positive integers, $a$ and $b$, such that
$$a^2=28b\sqrt{8b^2+1}+80b^2+5, \tag{$\star$}$$
and I’d like to find $a$ in terms of $b$ (including $\sqrt{8b^2+1}$, if necessary/appropriate), with a range like $$6b + \sqrt{8b^2+1} ≤ a ≤ 2b^2+ 3\sqrt{8b^2+1}$$ being an acceptable compromise.
I recently learned how to use Laurent series, evaluated at infinity, to approximate the square root of a polynomial, but can't figure out exactly how to apply that technique here (Wolfram gives somewhat unhelpful results) — or even figure out if it’s applicable [because of the radical].
Obviously, $b^2$ is a square triangular number (because $8b^2+1$ is a square), but I haven’t figured out how to use that detail to help me.
What’s the best approximation to the square root of ($\star$)?
 A: Suppose that $a$ and $b$ are integers such that
$$a^2=28b\sqrt{8b^2+1}+80b^2+5.$$
Then in particular $\sqrt{8b^2+1}$ is an integer, and so $8b^2+1=c^2$ for some integer $c$. Then
$$a^2=28bc+10c^2-5,$$
and a bit of rearranging shows that
$$(10c^2-(a^2+5))^2=(28bc)^2=28^2b^2c^2=98(8b^2)c^2=98(c^2-1)c^2.$$
Expanding the left hand side and rearranging, we see that $c^2$ is a root of the quadratic
$$2X^2-(20a^2+2)X+(a^2+5)^2=0.$$
As $c^2$ is an integer the discriminant $\Delta$ of this quadratic is a perfect square, where
$$\Delta=(20a^2+2)^2-4\cdot2\cdot(a^2+5)^2=14^2(2a^4-1).$$
It follows that also $2a^4-1$ is a perfect square, say $d^2$, and so
$$d^2-2a^4=-1.$$
This shows that $(X,Y)=(d,a^2)$ is a solution to the Pell equation $X^2-2Y^2=-1$, and hence that $a^2$ is a Pell number that is a perfect square. That $a=13$ is the only such number is shown in this article.
A: Wolfram Alpha tells me that integer solutions to $8b^2+1 = n^2$ are of the form
$b = \frac{(3 + 2 \sqrt{2})^k - (3 - 2 \sqrt{2})^k} {4 \sqrt{2}} $
$n =  (\frac{1}{2} ((3 - 2 \sqrt{2})^k + (3 + 2 \sqrt{2})^k))$
for positive integer $k$.
Alternatively, solutions to $8b^2+1=n^2$ can be obtained by treating it as a Pell equation and using  Brahmagupta's fomula to generate solutions.
Trying values of $k$ up to $25000$, where $b$ has 19139 digits,  results in only one solution where $a$ is an integer, for $k=1$, where $b=1$,  $a=13$.
My suspicion is that this could be the only solution.
Code for the search:

from gmpy2 import isqrt

def add( p, q ):
  return ( p[0]*q[0]+8*p[1]*q[1], p[0]*q[1]+q[0]*p[1])

X=(1,0)
for k in range(1,100000):
  X = add ( X,(3,1))
  (n,b) = X
  assert 8*b*b+1==n*n

  a2 = 28*b*n + 80*b*b + 5

  a = isqrt(a2)

  assert a*a <= a2 and a2 <= (a+1)**2
  if a*a==a2:
    print ("solution", b,n, a)

  if k%1000==0:
    print ("tried k=",k,", b has ",len(str(b)),"digits")


