TAOCP 4.7 Exercise 13
Let W(z)=\sum_{n\ge0}W_n z^n, and define the truncated polynomial
Section 4.7: Manipulation of Power Series
Exercise 13. [M27] (Rational function approximation.) It is occasionally desirable to find polynomials whose quotient has the same initial terms as a given power series. For example, if $W(z) = 1 + z + 3z^2 + 7z^3 + \cdots$, there are essentially four different ways to express $W(z)$ as $w_1(z)/w_2(z) + O(z^k)$ where $w_1(z)$ and $w_2(z)$ are polynomials with $\deg(w_1) + \deg(w_2) < 4$:
$$ (1 + z + 3z^2 + 7z^3)/1 = 1 + z + 3z^2 + 7z^3 + 0z^4 + \cdots, $$ $$ (3 - 4z + 2z^2)/(3 - 7z) = 1 + z + 3z^2 + 7z^3 + \tfrac{49}{3}z^4 + \cdots, $$ $$ (1 - z)/(1 - 2z - z^2) = 1 + z + 3z^2 + 7z^3 + 17z^4 + \cdots, $$ $$ 1/(1 - z - 2z^2 - 2z^3) = 1 + z + 3z^2 + 7z^3 + 15z^4 + \cdots. $$
Rational functions of this kind are commonly called Padé approximations, since they were studied extensively by H. E. Padé [Annales Scient. de l'École Normale Supérieure (3) 9 (1892), S1–S93; (3) 16 (1899), 395–426].
Show that all Padé approximations $W(z) \approx w_1(z)/w_2(z) + O(z^N)$ with $\deg(w_1) + \deg(w_2) < N$ can be obtained by applying an extended Euclidean algorithm to the polynomials $z^N$ and $W_0 + W_1 z + \cdots + W_{N-1} z^{N-1}$, and design an all-integer algorithm for the case that each $W_i$ is an integer. [Hint: See exercise 4.6.1–26.]
Verified: yes
Solve time: 11m50s
Let
$$ W(z)=\sum_{n\ge0}W_n z^n, $$
and define the truncated polynomial
$$ P(z)=W_0+W_1z+\cdots+W_{N-1}z^{N-1}. $$
We seek all rational approximations
$$ W(z)=\frac{w_1(z)}{w_2(z)}+O(z^N) $$
such that
$$ \deg w_1+\deg w_2<N. $$
Since
$$ W(z)-P(z)=O(z^N), $$
the Padé condition is equivalent to
$$ \frac{w_1(z)}{w_2(z)}-P(z)=O(z^N). $$
Multiplying by $w_2(z)$,
$$ w_1(z)-P(z)w_2(z)\equiv0\pmod{z^N}. $$
Hence there exists a polynomial $c(z)$ such that
$$ w_1(z)-P(z)w_2(z)=c(z)z^N. \eqno(1) $$
Thus Padé approximants are exactly the polynomial solutions of (1) satisfying
$$ \deg w_1+\deg w_2<N. $$
We now derive all such solutions from the extended Euclidean algorithm applied to $z^N$ and $P(z)$.
Set
$$ r_{-1}(z)=z^N,\qquad r_0(z)=P(z), $$
and perform the Euclidean algorithm:
$$ r_{k-2}=q_k r_{k-1}+r_k, \qquad \deg r_k<\deg r_{k-1}. \eqno(2) $$
Simultaneously define polynomials $a_k,b_k$ by
$$ r_k=a_k z^N+b_k P. \eqno(3) $$
Initially,
$$ a_{-1}=1,\quad b_{-1}=0, \qquad a_0=0,\quad b_0=1. $$
Substituting (3) into (2), we obtain
$$ a_k=a_{k-2}-q_k a_{k-1}, \qquad b_k=b_{k-2}-q_k b_{k-1}. \eqno(4) $$
Equation (3) may be rewritten as
$$ r_k-b_kP=a_k z^N. \eqno(5) $$
Therefore
$$ r_k(z)-P(z)b_k(z)\equiv0\pmod{z^N}, $$
hence every pair
$$ w_1=r_k,\qquad w_2=b_k \eqno(6) $$
yields a Padé approximation.
We now establish the degree relations.
For every $k\ge0$,
$$ \deg b_k=N-\deg r_{k-1}. \eqno(7) $$
Proof is by induction on $k$.
For $k=0$,
$$ b_0=1, \qquad \deg b_0=0, $$
while
$$ \deg r_{-1}=N, $$
so (7) holds.
For $k=1$,
$$ b_1=b_{-1}-q_1b_0=-q_1. $$
Since
$$ \deg q_1=\deg r_{-1}-\deg r_0=N-\deg r_0, $$
we obtain
$$ \deg b_1=N-\deg r_0. $$
Now assume $k\ge2$, and that (7) holds for $b_{k-1}$ and $b_{k-2}$. From (4),
$$ b_k=b_{k-2}-q_k b_{k-1}. $$
Since
$$ \deg q_k=\deg r_{k-2}-\deg r_{k-1}, $$
we have
$$ \deg(q_k b_{k-1}) = \deg r_{k-2}-\deg r_{k-1} + N-\deg r_{k-2} = N-\deg r_{k-1}. $$
Also,
$$ \deg b_{k-2}=N-\deg r_{k-3} < N-\deg r_{k-1}, $$
because
$$ \deg r_{k-3}>\deg r_{k-1}. $$
Hence the leading term of $b_k$ comes from $-q_k b_{k-1}$, and therefore
$$ \deg b_k=N-\deg r_{k-1}. $$
This proves (7).
Consequently,
$$ \deg r_k+\deg b_k = \deg r_k+N-\deg r_{k-1} < N, $$
since
$$ \deg r_k<\deg r_{k-1}. $$
Thus every pair (6) satisfies the required degree bound.
It remains to prove that every Padé approximation arises in this way.
Let $(w_1,w_2)$ satisfy
$$ w_1-Pw_2=c z^N \eqno(8) $$
and
$$ \deg w_1+\deg w_2<N. \eqno(9) $$
If $w_1=w_2=0$, the pair is trivial. Assume otherwise.
Choose $k$ so that
$$ \deg r_k\le \deg w_1<\deg r_{k-1}. \eqno(10) $$
Divide $w_1$ by $r_k$:
$$ w_1=s r_k+t, \qquad \deg t<\deg r_k. \eqno(11) $$
Using (5),
$$ r_k=a_k z^N+b_k P, $$
therefore
$$ w_1 = s a_k z^N+s b_k P+t. $$
Substitute into (8):
$$ t-P(w_2-sb_k)=(c-sa_k)z^N. \eqno(12) $$
Now
$$ \deg t<\deg r_k, $$
and from (9) and (10),
$$ \deg w_2 < N-\deg w_1 \le N-\deg r_k. $$
Also, by (7),
$$ \deg b_k = N-\deg r_{k-1} < N-\deg r_k. $$
Hence
$$ \deg(w_2-sb_k)<N-\deg r_k. $$
Since $\deg P<N$,
$$ \deg P(w_2-sb_k) < N. $$
Therefore the left side of (12) has degree $<N$, while the right side is divisible by $z^N$. Hence both sides vanish identically:
$$ t=P(w_2-sb_k). \eqno(13) $$
Now
$$ \deg t<\deg r_k, $$
while
$$ \deg(w_2-sb_k)<N-\deg r_k. $$
Suppose $u=w_2-sb_k\neq0$. Since $P\neq0$,
$$ \deg t=\deg(Pu)=\deg P+\deg u. $$
Because $\deg P=N-1$,
$$ \deg t\ge N-1. $$
But
$$ \deg t<\deg r_k. $$
Hence
$$ \deg r_k>N-1. $$
This is impossible, because $\deg r_k<\deg r_{k-1}\le N$. Therefore $u=0$. Thus
$$ w_2=sb_k. $$
Equation (13) now gives $t=0$, so from (11),
$$ w_1=sr_k. $$
Hence every Padé approximation is a scalar multiple of one of the pairs
$$ (r_k,b_k) $$
produced by the Euclidean algorithm.
Therefore all Padé approximations satisfying
$$ \deg w_1+\deg w_2<N $$
are obtained from the extended Euclidean algorithm applied to $z^N$ and $P(z)$.
We now describe an all-integer algorithm when all $W_i$ are integers.
Ordinary polynomial division over $\mathbf Z[z]$ may introduce fractions. Instead we use pseudo-division.
Suppose
$$ A(z),B(z)\in\mathbf Z[z], \qquad \deg A=m, \quad \deg B=n, \quad m\ge n, $$
and let $b$ be the leading coefficient of $B$. Then there exist polynomials $Q,R\in\mathbf Z[z]$ such that
$$ b^{,m-n+1}A=QB+R, \qquad \deg R<n. \eqno(14) $$
Apply this at every Euclidean step. Starting with
$$ r_{-1}=z^N, \qquad r_0=P, $$
compute
$$ \lambda_k r_{k-2} = q_k r_{k-1}+r_k, \qquad \deg r_k<\deg r_{k-1}, \eqno(15) $$
where $\lambda_k$ is a suitable power of the leading coefficient of $r_{k-1}$. All coefficients remain integral.
Simultaneously maintain
$$ r_k=a_k z^N+b_k P $$
by the recurrences
$$ a_k=\lambda_k a_{k-2}-q_k a_{k-1}, $$
$$ b_k=\lambda_k b_{k-2}-q_k b_{k-1}. \eqno(16) $$
Since all operations involve only integer multiplication and subtraction, every polynomial occurring in the computation has integer coefficients.
Each pair
$$ (w_1,w_2)=(r_k,b_k) $$
then yields a Padé approximation. Common integer factors may be removed after each step if desired.
Finally, consider the example
$$ W(z)=1+z+3z^2+7z^3+\cdots, \qquad N=4, $$
so
$$ P(z)=1+z+3z^2+7z^3. $$
The Euclidean algorithm gives
$$ z^4 = \left(\frac17 z-\frac{3}{49}\right)P(z) + \left( \frac{2}{49}z^2-\frac{4}{49}z+\frac{3}{49} \right). $$
Thus
$$ r_1(z)=\frac{2}{49}z^2-\frac{4}{49}z+\frac{3}{49}, $$
and
$$ b_1(z)=-\frac17 z+\frac{3}{49}. $$
Multiplying numerator and denominator by $49$,
$$ \frac{r_1(z)}{b_1(z)} = \frac{2z^2-4z+3}{3-7z} = \frac{3-4z+2z^2}{3-7z}, $$
which is the Padé approximation listed in the statement.
Continuing the Euclidean algorithm yields
$$ \frac{1-z}{1-2z-z^2}, \qquad \frac1{1-z-2z^2-2z^3}, $$
together with the truncation
$$ 1+z+3z^2+7z^3. $$
These are exactly the Padé approximations satisfying
$$ \deg w_1+\deg w_2<4. $$
$\square$