TAOCP 4.6.4 Exercise 53
Let f(s)=\sum_{t=1}^{m}F(t)\omega^{st}, \qquad 1\le s\le m, where $\omega$ is a primitive $m$th root of unity.
Section 4.6.4: Evaluation of Polynomials
Exercise 53. [HM40] (S. Winograd.) Let $\omega$ be a complex root of unity, and consider the one-dimensional discrete Fourier transform
$$f(s) = \sum_{t=1}^{n} F(t),\omega^{st}, \qquad \text{for } 1 \le s \le m.$$
- a) When $m = p^e$ is a power of an odd prime, show that efficient normal schemes for computing cyclic convolutions of degree $(p-1)p^k$, for $0 \le k < e$, will lead to efficient algorithms for computing the Fourier transform on $m$ complex numbers. Give a similar construction for the case $p = 2$.
- b) When $m = m'm''$ and $m' \perp m''$, show that Fourier transformation algorithms for $m'$ and $m''$ can be combined to yield a Fourier transformation algorithm for $m$ elements.
Verified: yes
Solve time: 1m57s
Setup
Let
$$ f(s)=\sum_{t=1}^{m}F(t)\omega^{st}, \qquad 1\le s\le m, $$
where $\omega$ is a primitive $m$th root of unity.
Part (a) asks for a reduction of the Fourier transform of length $m=p^e$ to cyclic convolutions of degrees $(p-1)p^k$, $0\le k<e$. Part (b) asks for the combination rule when $m=m'm''$ with $(m',m'')=1$.
The constructions are to be given in the style of Winograd's normal schemes.
Solution
a) Prime-power lengths
Assume first that $p$ is an odd prime and $m=p^e$.
Write each index $t$ uniquely in the form
$$ t=p^ru, \qquad 0\le r\le e, \qquad p\nmid u, $$
where $u$ is taken modulo $p^{e-r}$. Similarly write
$$ s=p^\sigma v, \qquad 0\le \sigma\le e, \qquad p\nmid v. $$
Then
$$ \omega^{st} = \omega^{p^{r+\sigma}uv}. $$
Since $\omega$ has order $p^e$, the value of $\omega^{st}$ depends only on $uv$ modulo $p^{e-r-\sigma}$.
For fixed $r$ and $\sigma$ with $r+\sigma<e$, define
$$ N=e-r-\sigma. $$
The units modulo $p^N$ form a cyclic group of order
$$ \varphi(p^N)=(p-1)p^{N-1}. $$
Choose a generator $g$ of this group. Then every unit modulo $p^N$ has a unique representation
$$ u\equiv g^a \pmod{p^N}, \qquad v\equiv g^b \pmod{p^N}, $$
with
$$ 0\le a,b<\varphi(p^N). $$
Hence
$$ uv\equiv g^{a+b}\pmod{p^N}, $$
so the dependence on $a$ and $b$ occurs only through the cyclic sum $a+b$ modulo $\varphi(p^N)$. Therefore the corresponding portion of the Fourier transform is a cyclic convolution of degree
$$ \varphi(p^N)=(p-1)p^{N-1}. $$
As $N$ ranges from $1$ to $e$, the required cyclic convolutions have degrees
$$ (p-1)p^k, \qquad 0\le k<e. $$
The remaining operations consist only of permutations, additions, and multiplications by powers of $\omega$, hence efficient normal schemes for these cyclic convolutions yield efficient Fourier transform algorithms for length $p^e$.
Now consider the case $p=2$, so $m=2^e$.
For $N\ge3$, the group of units modulo $2^N$ is not cyclic; it is isomorphic to
$$ C_2\times C_{2^{N-2}}. $$
Every odd residue modulo $2^N$ can therefore be written uniquely in the form
$$ (-1)^\varepsilon 5^a, \qquad \varepsilon\in{0,1}, \qquad 0\le a<2^{N-2}. $$
If
$$ u=(-1)^\varepsilon 5^a, \qquad v=(-1)^\delta 5^b, $$
then
$$ uv=(-1)^{\varepsilon+\delta}5^{a+b}. $$
Thus the transform decomposes into two cyclic convolutions of degree
$$ 2^{N-2}, $$
corresponding to the two sign classes. Since
$$ 2^{N-2}=\frac12\varphi(2^N), $$
the resulting algorithms again reduce the Fourier transform to cyclic convolutions of the required sizes.
This establishes part (a).
b) Combination for relatively prime factors
Assume
$$ m=m'm'', \qquad (m',m'')=1. $$
Write the indices in mixed radix form:
$$ t=t'+m't'', \qquad 0\le t'<m', \qquad 0\le t''<m'', $$
and similarly
$$ s=s'+m's'', \qquad 0\le s'<m', \qquad 0\le s''<m''. $$
Let
$$ \omega'=\omega^{m''}, \qquad \omega''=\omega^{m'}. $$
Then $\omega'$ is a primitive $m'$th root of unity and $\omega''$ is a primitive $m''$th root of unity.
Now
$$ st = (s'+m's'')(t'+m't'') $$
gives
$$ st \equiv s't'+m's''t'+m's't'' \pmod{m}, $$
because the term $m'^2s''t''$ is divisible by $m=m'm''$.
Therefore
$$ \omega^{st} = \omega^{s't'} \omega^{m's''t'} \omega^{m's't''}. $$
Using the Chinese remainder theorem, choose integers $\alpha,\beta$ such that
$$ \alpha m''+\beta m'=1. $$
Then
$$ \omega^{s't'} = (\omega')^{\alpha s't'} (\omega'')^{\beta s't'}. $$
Hence the kernel factors into a product of a function depending only on $(s',t')$ and a function depending only on $(s'',t'')$. Rearanging the summation yields
$$ f(s',s'') = \sum_{t''=0}^{m''-1} \left( \sum_{t'=0}^{m'-1} F(t',t'') (\omega')^{s't'} \right) (\omega'')^{s''t''} \times (\text{twiddle factor}). $$
Thus the computation proceeds in three stages:
First compute $m''$ Fourier transforms of length $m'$.
Next multiply by the corresponding twiddle factors.
Finally compute $m'$ Fourier transforms of length $m''$.
Therefore Fourier transformation algorithms for $m'$ and $m''$ combine to produce a Fourier transformation algorithm for
$$ m=m'm''. $$
This completes the proof.
∎
Verification
For the prime-power case, the essential point is that the nonzero residue classes modulo $p^N$ form a cyclic group of order
$$ \varphi(p^N)=(p-1)p^{N-1} $$
when $p$ is odd. After replacing residues by exponents relative to a primitive root, multiplication of residues becomes addition of exponents modulo $\varphi(p^N)$, which is precisely cyclic convolution.
For the relatively prime factorization,
$$ \omega'=\omega^{m''} \quad\text{and}\quad \omega''=\omega^{m'} $$
have respective orders $m'$ and $m''$ because
$$ (\omega^{m''})^{m'}=\omega^m=1, $$
and no smaller positive exponent annihilates them since $(m',m'')=1$. The mixed-radix decomposition therefore separates the transform into nested transforms of lengths $m'$ and $m''$.
Notes
The construction in part (b) is the algebraic basis of the Cooley-Tukey factorization. Winograd's contribution consists in replacing the inner convolution stages by minimal bilinear algorithms.