What is this change of variable in a polynomial? The autocovariance generating function, $\gamma(z)$, is defined as the $z$-transform of the autocovariance function, $\gamma_\tau$:
$$
\gamma(z) = \gamma_0 + \gamma_1(z+z^{-1}) + \gamma_2(z^2+z^{-2}) + \gamma_3(z^3+z^{-3}) + \cdots = 
\displaystyle\sum_{\tau=-\infty}^\infty \gamma_\tau z^\tau \,, 
\qquad\quad (1)
$$
where $z$ is a complex variable.
I am studying the code in this software package [1] and have found a function that 
performs the following operation on the equation given above for $\gamma(z)$. According to the documentation:

This function transforms a polynomial in the variables S(n)=z^n + z^(-n)
  into a polynomial in the variable U=z + z^(-1).

To do so, the following recursions are implemented in the code:
$$
\begin{eqnarray}
\begin{array}{ll}
S(0) = 2, \quad S(1) = 1, \quad \hbox{initialisation} \\
S(n) = U \times S(n-1) - S(n-2), \quad   n >= 2
\end{array}
\end{eqnarray}
$$
where $S(n)$ is the input polynomial of the form:
$$
c(S) = c_1 + c_2 S(1) + \dots + c_{n-1} S(n-2) + c_n S(n-1)
$$
and the output contains the coefficients of the polynomial:
$$
p(U) = p_1 U(n-1) + p_2 U(n-2) + \dots + p_{n-1}U + p_n \,.
$$
What may be the purpose of this transformation? For what purposes can this transformation be convenient?

Is it equivalent to work with the original polynomial or the other? One of the reasons why I don't understand this transformation is the following. If $z$ is replaced in equation (1) by $e^{-iw}$, with $\omega \in[0, \pi]$, then the Fourier transform of the autocovariances $\gamma_\tau$ is obtained (multiplied by $2\pi$). I think there should be a mapping for the variable $U$ that also returns the Fourier transform of $\gamma_\tau$, but I don't know which values should $U$ take in order to get it. 

Hints: Maybe the algorithms employed in the code require the polynomials in that format, so here is some context:
The transformation of the polynomial is performed before running an algorithm that computes the greatest common divisor of polynomials and partial fraction decomposition of a polynomial.

The code is available in the file SSMMATLAB/sn2u.m of the source code.

[1] SSMMATLAB, a Set of MATLAB Programs for the Statistical Analysis of State Space Models by Víctor Gómez.
 A: I have received valuable feedback from other sources and have come up with 
an understanding of the transformation. For $e^{-iw}$ we have that the reciprocal $1/z$ is the same as the conjugate of $z$, so that $(e^{-iw})^n + (e^{iw})^{-n} = \cos(nw) - i\sin(nw) + \cos(nw) + i\sin(nw) = 2\cos(n w)$.
Thus, for a given frequency $w$ the original polynomial is defined as:
$$
A(z) = a_0 + a_1 (2\cos(w)) + a_2 (2\cos(2w)) + a_3 (2\cos(3w)) + ...
$$
The function sn2u returns the coefficients $b_i$ of the following polynomial:
$$
B(u) = b_0 + b_1 (2cos(w)) + b_2 (2cos(w))^2 + b_3 (2cos(w))^3 + ...
$$
rewriting $x = 2cos(w)$:
$$
B(x) = b_0 + b_1 x + b_2 x^2 + b_3 x^3 + ...
$$
The last representation is the usual expression of a polynomial. This makes 
the quotient of these kinds of polynomials liable to be decomposed according 
to standard techniques for partial fraction decomposition.
A: If $f(z)$ is an analytic function, then
$$f(z)=f(z_0)+f'(z_0)(z-z_0)+\frac{f''(z_0)}2(z-z_0)^2+\dots$$
Set $z_0=0:$
$$f(z)=\sum_{n=0}^\infty\frac{f^{(n)}(0)}{n!}z^n$$
Set $z\to\frac1z:$
$$f(1/z)=\sum_{n=0}^\infty\frac{f^{(n)}(0)}{n!}z^{-n}$$
Add the two together and you get
$$f(z)+f(1/z)=\sum_{n=0}^\infty\frac{f^{(n)}(0)}{n!}(z^n+z^{-n})$$
So that $\gamma(z)=f(z)+f(1/z)$ and $\gamma_\tau=\frac{f^{(\tau)}(0)}{\tau!}$.  One can easily calculate $\gamma(z)$ from $f(z)$ and the series expansion, or given the series expansion, calculate $\gamma(z)$ from $f(z)$.
