The Rader algorithm

Bricogne, G.

doi:10.1107/97809553602060000551

International
Tables for
Crystallography
Volume B
Reciprocal space
Edited by U. Shmueli

pdf | chapter contents | chapter index | related articles

International Tables for Crystallography (2006). Vol. B. ch. 1.3, pp. 53-54 | 1 | 2 |

Section 1.3.3.2.3. The Rader algorithm

G. Bricogne^a

^a MRC Laboratory of Molecular Biology, Hills Road, Cambridge CB2 2QH, England, and LURE, Bâtiment 209D, Université Paris-Sud, 91405 Orsay, France

1.3.3.2.3. The Rader algorithm

| top | pdf |

The previous two algorithms essentially reduce the calculation of the DFT on N points for N composite to the calculation of smaller DFTs on prime numbers of points, the latter remaining irreducible. However, Rader (1968) showed that the p-point DFT for p an odd prime can itself be factored by invoking some extra arithmetic structure present in $[{\bb Z} / p {\bb Z}]$ .

1.3.3.2.3.1. N an odd prime

| top | pdf |

The ring $[{\bb Z} / p {\bb Z} = \{0,1,2,\ldots,p - 1\}]$ has the property that its [p - 1] non-zero elements, called units, form a multiplicative group [U(p)] . In particular, all units $[r \in U(p)]$ have a unique multiplicative inverse in $[{\bb Z} / p {\bb Z}]$ , i.e. a unit $[s \in U(p)]$ such that $[rs \equiv 1\hbox { mod } p]$ . This endows $[{\bb Z} / p {\bb Z}]$ with the structure of a finite field.

Furthermore, [U(p)] is a cyclic group, i.e. consists of the successive powers $[g^{m}\hbox{ mod } p]$ of a generator g called a primitive root mod p (such a g may not be unique, but it always exists). For instance, for [p = 7] , $[U(7) = \{1,2,3,4,5,6\}]$ is generated by [g = 3] , whose successive powers mod 7 are: $[g^{0} = 1, \quad g^{1} = 3, \quad g^{2} = 2, \quad g^{3} = 6, \quad g^{4} = 4, \quad g^{5} = 5]$ [see Apostol (1976), Chapter 10].

The basis of Rader's algorithm is to bring to light a hidden regularity in the matrix [F(p)] by permuting the basis vectors $[{\bf u}_{k}]$ and $[{\bf v}_{k^{*}}]$ of $[L({\bb Z} / p {\bb Z})]$ as follows: $[\eqalign{{\bf u}'_{0} &= {\bf u}_{0} \cr {\bf u}'_{m} &= {\bf u}_{k} {\hbox to 12pt{}}\hbox{with } k = g^{m}, {\hbox to 15pt{}} m = 1, \ldots, p - 1\hbox{;} \cr {\bf v}'_{0} &= {\bf v}_{0} \cr {\bf v}'_{m^{*}} &= {\bf v}_{k^{*}} \quad \hbox{with } k^{*} = g^{m^{*}}, \quad m^{*} = 1, \ldots, p - 1\hbox{;}}]$ where g is a primitive root mod p.

With respect to these new bases, the matrix representing $[\bar{F}(p)]$ will have the following elements: $[\eqalign{\hbox{element } (0,0) &= 1 \cr \hbox{element } (0, m + 1) &= 1 \quad \hbox{for all } m = 0, \ldots p - 2, \cr \hbox{element } (m^{*} + 1,0) &= 1 \quad \hbox{for all } m^{*} = 0, \ldots, p - 2, \cr \hbox{element } (m^{*} + 1, m + 1) &= e \left({k^{*}k \over p}\right) \cr &= e(g^{(m^{*} + m)/p}) \cr &\qquad \quad \hbox{for all } m^{*} = 0, \ldots, p - 2.}]$ Thus the `core' $[\bar{C}(p)]$ of matrix $[\bar{F}(p)]$ , of size $[(p - 1) \times (p - 1)]$ , formed by the elements with two non-zero indices, has a so-called skew-circulant structure because element $[(m^{*}, m)]$ depends only on $[m^{*} + m]$ . Simplification may now occur because multiplication by $[\bar{C}(p)]$ is closely related to a cyclic convolution. Introducing the notation $[C(m) = e(g^{m/p})]$ we may write the relation $[{\bf Y}^{*} = \bar{F}(p){\bf Y}]$ in the permuted bases as $[\eqalign{Y^{*} (0) &= {\textstyle\sum\limits_{k}} Y(k) \cr Y^{*} (m^{*} + 1) &= Y(0) + {\textstyle\sum\limits_{m = 0}^{p - 2}} C(m^{*} + m) Y(m + 1) \cr &= Y(0) + {\textstyle\sum\limits_{m = 0}^{p - 2}} C(m^{*} - m) Z(m) \cr &= Y(0) + ({\bf C} * {\bf Z}) (m^{*}), \quad m^{*} = 0, \ldots, p - 2,}]$ where Z is defined by [Z(m) = Y(p - m - 2)] , $[m = 0, \ldots, p - 2]$ .

Thus $[{\bf Y}^{*}]$ may be obtained by cyclic convolution of C and Z, which may for instance be calculated by $[{\bf C} * {\bf Z} = F(p - 1) [\bar{F}(p - 1) [{\bf C}] \times \bar{F} (p - 1) [{\bf Z}]],]$ where × denotes the component-wise multiplication of vectors. Since p is odd, [p - 1] is always divisible by 2 and may even be highly composite. In that case, factoring $[\bar{F} (p - 1)]$ by means of the Cooley–Tukey or Good methods leads to an algorithm of complexity p log p rather than $[p^{2}]$ for $[\bar{F}(p)]$ . An added bonus is that, because $[g^{(p-1) / 2} = -1]$ , the elements of $[\bar{F} (p - 1) [{\bf C}]]$ can be shown to be either purely real or purely imaginary, which halves the number of real multiplications involved.

1.3.3.2.3.2. N a power of an odd prime

| top | pdf |

This idea was extended by Winograd (1976, 1978) to the treatment of prime powers $[N = p^{\nu}]$ , using the cyclic structure of the multiplicative group of units $[U(p^{\nu})]$ . The latter consists of all those elements of $[{\bb Z} / p^{\nu} {\bb Z}]$ which are not divisible by p, and thus has $[q_{\nu} = p^{\nu - 1} (p - 1)]$ elements. It is cyclic, and there exist primitive roots g modulo $[p^{\nu}]$ such that $[U(p^{\nu}) = \{1, g, g^{2}, g^{3}, \ldots, g^{q_{\nu} - 1}\}.]$ The $[p^{\nu - 1}]$ elements divisible by p, which are divisors of zero, have to be treated separately just as 0 had to be treated separately for [N = p] .

When $[k^{*} \not\in U(p^{\nu})]$ , then $[k^{*} = pk_{1}^{*}]$ with $[k_{1}^{*} \in {\bb Z} / p^{\nu - 1} {\bb Z}]$ . The results $[X^{*} (pk_{1}^{*})]$ are p-decimated, hence can be obtained via the $[p^{\nu - 1}]$ -point DFT of the $[p^{\nu - 1}]$ -periodized data Y: $[X^{*} (pk_{1}^{*}) = \bar{F} (p^{\nu - 1}) [{\bf Y}] (k_{1}^{*})]$ with $[Y(k_{1}) = {\textstyle\sum\limits_{k_{2} \in {\bb Z} / p {\bb Z}}} X(k_{1} + p^{\nu - 1} k_{2}).]$

When $[k^{*} \in U(p^{\nu})]$ , then we may write $[X^{*} (k^{*}) = X_{0}^{*} (k^{*}) + X_{1}^{*} (k^{*}),]$ where $[{\bf X}_{0}^{*}]$ contains the contributions from $[k\; \notin\; U(p^{\nu})]$ and $[{\bf X}_{1}^{*}]$ those from $[k \in U(p^{\nu})]$ . By a converse of the previous calculation, $[{\bf X}_{0}^{*}]$ arises from p-decimated data Z, hence is the $[p^{\nu - 1}]$ -periodization of the $[p^{\nu - 1}]$ -point DFT of these data: $[X_{0}^{*} (p^{\nu - 1} k_{1}^{*} + k_{2}^{*}) = \bar{F} (p^{\nu - 1}) [{\bf Z}] (k_{2}^{*})]$ with $[Z(k_{2}) = X(pk_{2}), \qquad k_{2} \in {\bb Z} / p^{\nu - 1} {\bb Z}]$ (the $[p^{\nu - 1}]$ -periodicity follows implicity from the fact that the transform on the right-hand side is independent of $[k_{1}^{*} \in {\bb Z} / p {\bb Z}]$ ).

Finally, the contribution $[X_{1}^{*}]$ from all $[k \in U(p^{\nu})]$ may be calculated by reindexing by the powers of a primitive root g modulo $[p^{\nu}]$ , i.e. by writing $[X_{1}^{*} (g^{m^{*}}) = {\textstyle\sum\limits_{m = 0}^{q_{\nu} - 1}} X(g^{m}) e(g^{(m^{*} + m) / p^{\nu}})]$ then carrying out the multiplication by the skew-circulant matrix core as a convolution.

Thus the DFT of size $[p^{\nu}]$ may be reduced to two DFTs of size $[p^{\nu - 1}]$ (dealing, respectively, with p-decimated results and p-decimated data) and a convolution of size $[q_{\nu} = p^{\nu - 1} (p - 1)]$ . The latter may be `diagonalized' into a multiplication by purely real or purely imaginary numbers (because $[g^{(q_{\nu} / 2)} = -1]$ ) by two DFTs, whose factoring in turn leads to DFTs of size $[p^{\nu - 1}]$ and [p - 1] . This method, applied recursively, allows the complete decomposition of the DFT on $[p^{\nu}]$ points into arbitrarily small DFTs.

1.3.3.2.3.3. N a power of 2

| top | pdf |

When $[N = 2^{\nu}]$ , the same method can be applied, except for a slight modification in the calculation of $[{\bf X}_{1}^{*}]$ . There is no primitive root modulo $[2^{\nu}]$ for $[\nu \gt 2]$ : the group $[U(2^{\nu})]$ is the direct product of two cyclic groups, the first (of order 2) generated by −1, the second (of order [N/4] ) generated by 3 or 5. One then uses a representation $[\eqalign{k &= (-1)^{m_{1}} 5^{m_{2}} \cr k^{*} &= (-1)^{m_{1}^{*}} 5^{m_{2}^{*}}}]$ and the reindexed core matrix gives rise to a two-dimensional convolution. The latter may be carried out by means of two 2D DFTs on $[2 \times (N/4)]$ points.

References

Apostol, T. M. (1976). Introduction to analytic number theory. New York: Springer-Verlag.Google Scholar

Rader, C. M. (1968). Discrete Fourier transforms when the number of data samples is prime. Proc. IEEE, 56, 1107–1108.Google Scholar

Winograd, S. (1976). On computing the discrete Fourier transform. Proc. Natl Acad. Sci. USA, 73, 1005–1006.Google Scholar

Winograd, S. (1978). On computing the discrete Fourier transform. Math. Comput. 32, 175–199.Google Scholar

International Tables for Crystallography (2006). Vol. B. ch. 1.3, pp. 53-54