What is the algebraic intuition behind Vieta jumping in IMO1988 Problem 6?

At the heart of these so-called "Vieta-jumping" techniques are certain symmetries (reflections) on conics. These symmetries govern descent in the group of integer points of the conic. If you wish to develop a deeper understanding of these proofs then I highly recommend that you study them from this more general perspective, where you will find much beauty and unification.

The group laws on conics can be viewed essentially as special cases of the group law on elliptic curves (e.g. see Franz Lemmermeyer's "poor man's" papers), which is a helpful perspective to know. See also Sam Northshield's expositions on associativity of the secant method (both linked here).

If memory serves correct, many of these contest problems are closely associated with so-called Richaud-Degert quadratic irrationals, which have short continued fraction expansions (or, equivalently, small fundamental units). Searching on "Richaud Degert" etc should locate pertinent literature (e.g. Lemmermeyer's Higher descent on Pell Conics 1). Many of the classical results are couched in the language of Pell equations, but it is usually not difficult to translate the results into more geometric language.

So, in summary, your query about a "natural or canonical way to see the answer to the problem" is given a beautiful answer when you study the group laws of conics (and closely related results such as the theory of Pell equations). Studying these results will provide much motivation and intuition for generalizations such as group laws on elliptic curves.

See also Aubry's beautiful reflective generation of primitive Pythagorean triples, which is a special case of modern general results of Wall, Vinberg, Scharlau et al. on reflective lattices, i.e. arithmetic groups of isometries generated by reflections in hyperplanes.

The idea has been known (at least) since Gauss' Disquisitiones Arithmeticae 200 years ago. "Vieta jumping" is a name used only in competition manuals, the accepted term in mathematics is "reduction theory of quadratic forms". The reduction-theory solution is the canonical solution, and I do not know any solutions using other methods, but some presentations of the method can make it seem artificial.

The reason more competitors did not solve the problem is that in the quaint old days, high school students were not learning heavy machinery before going to the IMO.

Even knowing the theory, it may not be easy to recognize within a few hours that the problem would fall to a straightforward application of the 'rotation' group of the integer points on the conic (showing that any positive integer solution can be moved to a smaller one with $ab=0$) by writing out the equation as $a^2 - kab + b^2=k$ and turning the handle. Stripped of computational details, that is what the solution does.

It was also known since the 1800's that this type of quadratic form has special properties and is easier to analyze.

For degree $(2,2,...,2)$ equations of higher total degree, such as Markov triples (solutions of $x^2+y^2+z^2=3xyz$), the "Vieta" transformation does not have a specific name, but has been the main tool for organizing the solutions since the first papers in the $19$th century where these equations were analyzed.

https://www.encyclopediaofmath.org/index.php/Hurwitz_equation

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}}} \\ $$

The same Wikipedia page quotes Arthur Engel about the history of the problem.

Nobody of the six members of the Australian problem committee could solve it … it was [then] sent to the four most renowned Australian number theorists. They were asked to work on it for six hours. None of them could solve it in this time … the jury finally had the courage to choose it as the last problem of the competition. Eleven students gave perfect solutions.

It is implied that the Australians had tried to solve the problem with more general tools and failed, and that the idea of Viète jumping was created specifically for this problem.

Another reason why Viète jumping is the canonical solution to the problem is that the problems of the IMO usually have only a few solutions, with special prizes given out for particularly ingenious ones, such as to Boreico Iurie for IMO 2005 Q3:

Let $x,y,z\in\Bbb R^+$ such that $xyz\ge1$. Prove that $$\frac{x^5-x^2}{x^5+y^2+z^2}+\frac{y^5-y^2}{y^5+z^2+x^2}+\frac{z^5-z^2}{z^5+x^2+y^2}\ge0.$$

(I got this from a book about the mathematical olympiads in China from 2003 to 2006, which had the IMO problems for completeness.) In particular, the Viète jumping proof is the only one with mathematics simple enough for a competitor to understand.

If you insist on a more "natural" solution, a geometric interpretation of the technique is also on the Wikipedia page, involving lattice points on a hyperbola. Dubuque's answer has more on this, including the possible sources for the problem.