Earlier i posted this reply which doesn't make use of Vieta jumping and proves $ k = \gcd(a,b)^2 $. The proof below gives a complete solution for the equation.
First note that:
$$ \gcd(a,1+ab) = \gcd(b,1+ab) =1
$$So the common prime factors in the sum of $a^2$ and $b^2$ will remain unchanged and indivisible by $ab+1$.
That suggests that we should take $\gcd(a,b)$ out like this:
$$ k={\gcd(a,b)^2({({a\over{\gcd(a,b)}})^2+({b\over{\gcd(a,b)}})^2})\over{1+ab}} \tag{1}
$$
Observe that also:
$$ \gcd(\gcd(a,b)^2,1+ab)) = 1 \tag{2}$$
Because no prime factor in $a$ or $b$ can divide $ab+1$.
Therefore:
$$ k'={k\over{\gcd(a,b)^2}} = {{({a\over{\gcd(a,b)}})^2+({b\over{\gcd(a,b)}})^2}\over{1+ab}} $$ also must be an integer because $\gcd(a,b)^2$ divides the numerator of $(1)$ while having no prime factors in common with the denominator because of $(2)$.
For clarity write: $$ k'= {{{x_{n-1}}^2+{x_{n}}^2}\over{1+g^2{x_{n-1}}{x_{n}}}} \tag{3} $$ with: $$ g=\gcd(a,b) \text{ and } {x_{n-1}}= {a\over{g}} \text{ and } {x_{n}}= {b\over{g}} \text{ and } \gcd({x_{n-1}},{x_{n}}) = 1 $$
all integers of course.
Any prime divisor of $k'$ must divide ${x_{n-1}}^2 +{x_{n}}^2$. Because $ {x_{n-1}} $ and $ {x_{n}} $ are coprime it cannot divide $ {x_{n-1}} $ and $ {x_{n}} $ both, also it can't divide only one of $ {x_{n-1}} $ and $ {x_{n}} $ so it divides neither:
$$
\gcd(k',{x_{n-1}})=\gcd(k',{x_{n}}) =\gcd({x_{n-1}},{x_{n}}) =1
$$
If we can prove $k' = 1$ then we are done because that would mean:
$$ k = \gcd(a,b)^2
$$
Now the proof of $ k'= 1 $ (with the assumption : $ 0<{x_{n-1}}<{x_{n}}$ )
Write $(3)$ as: $ {x_{n-1}}^2+{x_{n}}^2=k'g^2{x_{n-1}}{x_{n}}+k' $ . Now if we assume $ k' \ge {x_{n-1}}^2 +{x_{n}} $ then:
$$
1={{x_{n-1}}^2+{x_{n}}^2\over{k'g^2{x_{n-1}}{x_{n}}+k'}} \le {{x_{n-1}}^2+{x_{n}}^2\over{g^2{x_{n-1}}^3{x_{n}} + g^2{x_{n-1}}{x_{n}}^2+{x_{n-1}}^2+{x_{n}}}} < 1 \enspace \implies contradiction \\
$$
So: $ 0 \lt k' \lt {x_{n-1}}^2 + {x_{n}} $ .
Or also: $ -{x_{n}} \lt {x_{n-1}}^2 -k' \tag{*}$
From this it follows that ${x_{n-1}}^2 -k' \ge 0 $ :
We know $\enspace {x_{n}} \mid {x_{n-1}}^2 -k' $. And $(*)$ implies that if $ {x_{n-1}}^2 -k' < 0 $ then $ {x_{n-1}}^2 -k' $ would have to be in $ \left\langle -{x_{n}} , 0 \right\rangle $ . In that interval it would not be divisible by $ {x_{n}} $.
So we are left with the following two possibilities:
Case 1: $ {x_{n-1}}^2 -k' =0 $
Now because $ \gcd(k',{x_{n-1}}) = 1 $ we have $ {x_{n-1}} =1 $ and $ k'=1 $.
Case 2: $ {x_{n-1}}^2 -k' \gt 0 $
In this case the requirement $ {x_{n}} \mid {x_{n-1}}^2 -k' $ tells us that there must be an integer $ {x_{n-2}} $ with: $ 0 < {x_{n-2}} < {x_{n-1}} $ and: $ {x_{n-2}}{x_{n}} = {x_{n-1}}^2 -k' $.
With $ (3) $ this means: $\enspace {x_{n}}=k'g^2{x_{n-1}} -{x_{n-2}} $. Substituting this again in $ (3) $ gives: $ {x_{n-1}}^2+ {x_{n-2}}^2 = k'g^2{x_{n-1}}{x_{n-2}}+k' $. So we have the same equation back with $ {x_{n}} $ replaced by a smaller term.
From: $\enspace {x_{n}}=k'g^2{x_{n-1}} -{x_{n-2}} \enspace $ it's also evident that $\enspace \gcd(k',{x_{n-2}})=\gcd({x_{n-1}},{x_{n-2}}) =1 \enspace $.
From the above we conclude that we can repeat taking smaller terms using: $\enspace {x_{i}}=k'g^2{x_{i-1}} -{x_{i-2}} \enspace $ until we reach a pair $ 0,{x_{0}} $ for which case $ 1 $ applies. We see that $ {x_{0}}=1 $ and conclude that $ k'=1 $
$ \enspace \square $
We see that the general solution can be generated using the following recursive formula:
$$
x_0=1, x_1 = g^2\\
{x_{n}}= g^2{x_{n-1}} -{x_{n-2}} \tag{4}\\
$$
This gives the following form for the $x_n$ :
$$
\text{if $n$ is even:}\enspace x_n= { \sum\limits_{i=0}^{{n\over{2}}} {(-1)^{i+{n\over{2}}}{{n\over{2}}+i\choose {n\over{2}}-i}g^{4i}}} \tag{5}\\
\text{if $n$ is odd:}\enspace x_n= \sum\limits_{i=0}^{{n-1\over{2}}} (-1)^{i+{n-1\over{2}}}{1+{n-1\over{2}}+i \choose {n-1\over{2}}-i}g^{4i+2} \\
$$
Proof of this is by induction:
1) $n$ is even:
$g^2x_{n-1} = \sum\limits_{i=0}^{{n-2\over{2}}} (-1)^{i+{n-2\over{2}}}{1+{n-2\over{2}}+i \choose {n-2\over{2}}-i}g^{4i+4} = $
$\sum\limits_{i=1}^{{n\over{2}}} (-1)^{i +{n\over{2}}}{{n\over{2}}+i-1 \choose {n\over{2}}-i}g^{4i} $
$ = \sum\limits_{i=1}^{{n\over{2}}-1} (-1)^{i +{n\over{2}}}{{n\over{2}}+i-1 \choose {n\over{2}}-i}g^{4i} +
(-1)^{n}{n-1 \choose 0}g^{2n}
$
$-x_{n-2}= { \sum\limits_{i=0}^{{n-2\over{2}}} {-(-1)^{i+{n-2\over{2}}}{{n-2\over{2}}+i\choose {n-2\over{2}}-i}g^{4i}}} = $
${ \sum\limits_{i=0}^{{n\over{2}}-1} {(-1)^{i+{n\over{2}}}{{n\over{2}}+i-1\choose {n\over{2}}-i-1}g^{4i}}}$
${= \sum\limits_{i=1}^{{n\over{2}}-1} {(-1)^{i+{n\over{2}}}{{n\over{2}}+i-1\choose {n\over{2}}-i-1}g^{4i}}} + {(-1)^{{n\over{2}}}{{n\over{2}}-1\choose {n\over{2}}-1}g^{0}} $
${x_{n}}= g^2{x_{n-1}} -{x_{n-2}} \implies $
$ {x_{n}}= {(-1)^{{n\over{2}}}{{n\over{2}}-1\choose {n\over{2}}-1}g^{0} + \sum\limits_{i=1}^{{n\over{2}}-1} (-1)^{i +{n\over{2}}} \bigg[ {{n\over{2}}+i-1 \choose {n\over{2}}-i} +{{n\over{2}}+i-1\choose {n\over{2}}-i-1} \bigg] g^{4i}
+
(-1)^{n}{n-1 \choose 0}g^{2n}} ={ \sum\limits_{i=0}^{{n\over{2}}} {(-1)^{i+{n\over{2}}}{{n\over{2}}+i\choose {n\over{2}}-i}g^{4i}}} $
2) $n$ is odd:
$x_n={ \sum\limits_{i=0}^{{n-1\over{2}}} {(-1)^{i+{n-1\over{2}}}{{n-1\over{2}}+i\choose {n-1\over{2}}-i}g^{4i+2}}} -
\sum\limits_{i=0}^{{n-3\over{2}}} (-1)^{i+{n-3\over{2}}}{1+{n-3\over{2}}+i \choose {n-3\over{2}}-i}g^{4i+2}$
$={ \sum\limits_{i=0}^{{n-3\over{2}}} {(-1)^{i+{n-1\over{2}}}\bigg[{{n-1\over{2}}+i\choose {n-1\over{2}}-i} + {{n-1\over{2}}+i \choose {n-1\over{2}}-i-1} \bigg] g^{4i+2}}} + {(-1)^{{n-1} }{{n-1}\choose 0}g^{2n} } = \sum\limits_{i=0}^{{n-1\over{2}}} (-1)^{i+{n-1\over{2}}}{1+{n-1\over{2}}+i \choose {n-1\over{2}}-i}g^{4i+2} $
Some values:
$$
n=1:\enspace g^2 \\
n=2:\enspace g^4-1 \\
n=3:\enspace g^6-2g^2 \\
n=4:\enspace g^8-3g^4+1\\
n=5:\enspace g^{10}-4g^6+3g^2 \\
n=6:\enspace g^{12}-5g^8+6g^4 -1 \\
n=7:\enspace g^{14}-6g^{10}+10g^6 -4g^2 \\
n=8:\enspace g^{16}-7g^{12}+15g^8 -10g^4 +1 \\
n=9:\enspace g^{18}-8g^{14}+21g^{10} -20g^{6} +5g^{2} \\
$$
A solution for the original equation with $ a $ and $ b $:
$$
\text{if $n$ is even:}\enspace a_n= { \sum\limits_{i=0}^{{n\over{2}}} {(-1)^{i+{n\over{2}}}{{n\over{2}}+i\choose {n\over{2}}-i}g^{4i+1}}} \tag{6}\\
\text{if $n$ is odd:}\enspace a_n= \sum\limits_{i=0}^{{n-1\over{2}}} (-1)^{i+{n-1\over{2}}}{1+{n-1\over{2}}+i \choose {n-1\over{2}}-i}g^{4i+3} \implies \\
{{a_n}^2+{a_{n+1}}^2\over{a_na_{n+1}+1}}=g^2\\
$$
To answer the question: I started out trying to prove the problem in a more logical way by taking out the $ \gcd(a,b) $ first. But found out I couldn't get any futher without using essentially the same arguments as used in Vieta jumping ( the fact that the equation has two integer roots and the sequence has to reach $ 0 $ eventually ).
But it does seem to me to give more insight if you do it this way. Also this gives a more complete answer to the problem. The general and more abstract Vieta jumping approach in Wikipedia gives you only that $ k $ must be square.
Update 12/14/16. Added a bit of Python code to check $(6)$ . With this you can also generate numbers yourself.
import math
N = 100 # fill in number of iterations
G = 25 # fill in gcd(a,b)
def binom(x, y):
if y == x:
return 1
if y == 1:
return x
if y > x:
return 0
a = math.factorial(x)
b = math.factorial(y)
c = math.factorial(x-y)
return a // (b*c)
def nextterm(a, gv):
b = 0
if a % 2 == 0: # even case
for n in range(a//2 + 1):
b += ((-1)**(n + (a // 2)) *
binom(a//2 + n, a//2 - n) * gv**(4*n))
else: # odd case
for n in range((a-1)//2 + 1):
b += ((-1)**(n + (a-1)//2) *
binom((a-1)//2 + n+1, (a-1)//2 - n) * gv**(2 + 4*n))
return b
def quotient(a, b):
return (a**2 + b**2) // (a*b + 1)
for n in range(1, N+1):
u = nextterm(n, G)
v = nextterm(n-1, G)
g = G
j = (((math.sqrt(g**4 - 4) + g**2) / 2)**n *
(g / (math.sqrt(g**4 - 4))) -
(((-1 * (math.sqrt(g**4 - 4))) + g**2) / 2)**n *
(g / (math.sqrt(g**4 - 4))))
print("------------- iteration %s -------------" % n)
print("a_%d = %d" % (n-1, G*v))
print("a_%d = %d" % (n, G*u))
print("a*_%d = %d" % (n-1, j))
print("k = %d" % quotient(G*u, G*v))
12/20/16: One last update to give also a closed-form expression in the style of Binet's formula:
$$
\text {The recursive formula : }\\
x_n=g^2x_{n-1}-x_{n-2}\\
\text {Can be written in matrix form : }\\
\begin{pmatrix} x_{n} \\ x_{n-1} \end{pmatrix} = \begin{pmatrix} g^2 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} x_{n-1} \\ x_{n-2} \end{pmatrix}\\
\text {Or like this : }\\
\vec{F}_{n}=A\vec{F}_{n-1}\\
\text {Which means after applying matrix $A$ several times from start value : }\\
$$
$$
\implies \vec{F}_{n}=A^n\vec{F}_{0} \tag{7}\\
$$
$$
\text{eigenvalues of $A$ are : } \lambda_1={g^2-\sqrt{g^4-4}\over{2}} \text{ and: } \lambda_2={g^2+\sqrt{g^4-4}\over{2}} \\
\text{eigenvectors of $A$ are : }\vec{\mu}=\begin{pmatrix} \lambda_1 \\ 1 \end{pmatrix} \text{ and: } \vec{\nu}=\begin{pmatrix} \lambda_2 \\ 1 \end{pmatrix} \\
\vec{F}_{0}=\begin{pmatrix} 1 \\ 0 \end{pmatrix} = {1\over{\sqrt{g^4-4}}}\vec{\nu} - {1\over{\sqrt{g^4-4}}} \vec{\mu}\\
\text{With $(7)$ this means : }\\
\vec{F}_{n}={{\bigg({g^2+\sqrt{g^4-4}\over{2}}\bigg)}^n\over{\sqrt{g^4-4}}} \begin{pmatrix} 1 \\ 0 \end{pmatrix} - {{\bigg({g^2-\sqrt{g^4-4}\over{2}} \bigg)}^n\over{\sqrt{g^4-4}}} \begin{pmatrix} 1 \\ 0 \end{pmatrix}\\
\text{We can read off values for $x$ and $a$ as : }\\
$$
$$
x_n={{\bigg({g^2+\sqrt{g^4-4}\over{2}}\bigg)}^n\over{\sqrt{g^4-4}}} - {{\bigg({g^2-\sqrt{g^4-4}\over{2}} \bigg)}^n\over{\sqrt{g^4-4}}} \tag{8}\\
a_n={g{\bigg({g^2+\sqrt{g^4-4}\over{2}}\bigg)}^n\over{\sqrt{g^4-4}}} - {g{\bigg({g^2-\sqrt{g^4-4}\over{2}} \bigg)}^n\over{\sqrt{g^4-4}}} \\
$$
