TAOCP 4.7 Exercise 14
Let U(z)=z+U_k z^k+\cdots,\qquad k\ge 2, and define the iterates of $U$ by
Section 4.7: Manipulation of Power Series
Exercise 14. ▶ [HM30] Fill in the details of Brent and Traub's method for calculating $U^{[n]}(z)$ when $U(z) = z + U_k z^k + \cdots$, using (27) and (28).
Verified: yes
Solve time: 6m57s
Setup
Let
$$ U(z)=z+U_k z^k+\cdots,\qquad k\ge 2, $$
and define the iterates of $U$ by
$$ U^{[0]}(z)=z,\qquad U^{[m+1]}(z)=U!\left(U^{[m]}(z)\right). $$
The problem is to make explicit Brent and Traub's procedure for computing the first $N$ coefficients of $U^{[n]}(z)$, using equations (27) and (28).
Equation (27) gives the semigroup law
$$ U^{[m+n]}=U^{[m]}\circ U^{[n]}, $$
and equation (28) is the special case
$$ U^{[2m]}=U^{[m]}\circ U^{[m]}. $$
The task is therefore to compute $U^{[n]}$ by repeated composition, exactly as ordinary exponentiation is computed by repeated squaring.
Solution
Write the binary expansion of $n$ as
$$ n=\sum_{j=0}^{r} \varepsilon_j 2^j, \qquad \varepsilon_j\in{0,1}, \qquad \varepsilon_r=1 . $$
Define power series
$$ A_0(z)=U(z), $$
and recursively
$$ A_{j+1}(z)=A_j(A_j(z)), \qquad 0\le j<r. $$
By (28),
$$ A_j(z)=U^{[2^j]}(z) $$
for every $j$. Indeed, $A_0=U^{[1]}$; if
$A_j=U^{[2^j]}$, then
$$ A_{j+1} =A_j\circ A_j =U^{[2^j]}\circ U^{[2^j]} =U^{[2^{j+1}]} $$
by (27).
Now set
$$ B(z)=z. $$
Scanning the binary digits of $n$, replace $B$ by
$$ B \leftarrow A_j(B) $$
whenever $\varepsilon_j=1$.
After all digits have been processed,
$$ B =U^{[\varepsilon_0]} \circ U^{[\varepsilon_1 2]} \circ \cdots \circ U^{[\varepsilon_r 2^r]} . $$
Since iterates of the same series commute,
$$ U^{[a]}\circ U^{[b]} =U^{[a+b]} $$
by repeated application of (27). Hence
$$ B =U^{\left[\sum_{j=0}^{r}\varepsilon_j2^j\right]} =U^{[n]}. $$
Therefore $B(z)$ is the desired iterate.
To obtain only the first $N$ coefficients, every composition is performed modulo $z^N$. Exercise 11 provides an algorithm for composing two power series to precision $N$, so each step above can be carried out with truncated series.
The complete procedure is therefore:
- Compute
$$ A_0=U,\qquad A_{j+1}=A_j\circ A_j \pmod{z^N} $$
for $0\le j<r$. 2. Initialize $B=z$. 3. For each $j$ with $\varepsilon_j=1$, replace
$$ B\leftarrow A_j(B)\pmod{z^N}. $$ 4. Output $B$.
The result is $U^{[n]}(z)$ through terms of degree $N-1$.
Verification
The construction of the series $A_j$ is correct because induction gives
$$ A_j=U^{[2^j]}. $$
The update step is correct because after processing digits
$\varepsilon_0,\ldots,\varepsilon_s$,
$$ B = U^{\left[\sum_{j=0}^{s}\varepsilon_j2^j\right]}. $$
This statement is proved by induction on $s$. When
$\varepsilon_s=0$, $B$ is unchanged. When $\varepsilon_s=1$,
$$ B \leftarrow U^{[2^s]}\circ U^{\left[\sum_{j=0}^{s-1}\varepsilon_j2^j\right]} = U^{\left[\sum_{j=0}^{s}\varepsilon_j2^j\right]} $$
by (27).
After the final digit is processed,
$$ B=U^{[n]}. $$
This completes the proof.
∎
Notes
The method is the exact analogue of binary exponentiation. Equation (28) supplies the repeated-squaring step, while equation (27) supplies the multiplication step. The only additional ingredient is a routine for composing power series modulo $z^N$, such as the algorithm of Exercise 11.