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

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

Section 1.3.3.2.2. The Good (or prime factor) algorithm

G. Bricognea

aGlobal Phasing Ltd, Sheraton House, Suites 14–16, Castle Park, Cambridge CB3 0AX, England, and LURE, Bâtiment 209D, Université Paris-Sud, 91405 Orsay, France

1.3.3.2.2. The Good (or prime factor) algorithm

| top | pdf |

1.3.3.2.2.1. Ring structure on [{\bb Z}/N{\bb Z}]

| top | pdf |

The set [{\bb Z}/N {\bb Z}] of congruence classes of integers modulo an integer N [see e.g. Apostol (1976)[link], Chapter 5] inherits from [{\bb Z}] not only the additive structure used in deriving the Cooley–Tukey factorization, but also a multiplicative structure in which the product of two congruence classes mod N is uniquely defined as the class of the ordinary product (in [{\bb Z}]) of representatives of each class. The multiplication can be distributed over addition in the usual way, endowing [{\bb Z}/N {\bb Z}] with the structure of a commutative ring.

If N is composite, the ring [{\bb Z}/N {\bb Z}] has zero divisors. For example, let [N = N_{1} N_{2}], let [n_{1} \equiv N_{1}] mod N, and let [n_{2} \equiv N_{2}] mod N: then [n_{1} n_{2} \equiv 0] mod N. In the general case, a product of nonzero elements will be zero whenever these elements collect together all the factors of N. These circumstances give rise to a fundamental theorem in the theory of commutative rings, the Chinese Remainder Theorem (CRT), which will now be stated and proved [see Apostol (1976[link]), Chapter 5; Schroeder (1986[link]), Chapter 16].

1.3.3.2.2.2. The Chinese remainder theorem

| top | pdf |

Let [N = N_{1} N_{2} \ldots N_{d}] be factored into a product of pairwise coprime integers, so that g.c.d. [(N_{i}, N_{j}) = 1] for [i \neq j]. Then the system of congruence equations[{\ell} \equiv {\ell}_{j} \hbox{ mod } N_{j},\qquad j = 1, \ldots, d,]has a unique solution [\ell] mod N. In other words, each [\ell \in {\bb Z}/N {\bb Z}] is associated in a one-to-one fashion to the d-tuple [(\ell_{1}, \ell_{2}, \ldots, \ell_{d})] of its residue classes in [{\bb Z}/N_{1} {\bb Z}, {\bb Z}/N_{2} {\bb Z}, \ldots, {\bb Z}/N_{d} {\bb Z}].

The proof of the CRT goes as follows. Let[Q_{j} = {N \over N_{j}} = {\prod\limits_{i \neq j}}\, N_{i}.]Since g.c.d. [(N_{j}, Q_{j}) = 1] there exist integers [n_{j}] and [q_{j}] such that[n_{j} N_{j} + q_{j} Q_{j} = 1,\qquad j = 1, \ldots, d,]then the integer[{\ell} = {\textstyle\sum\limits_{i = 1}^{d}}\, {\ell}_{i} q_{i} Q_{i} \hbox{ mod } N]is the solution. Indeed,[{\ell} \equiv {\ell}_{j} q_{j} Q_{j} \hbox{ mod } N_{j}]because all terms with [i \neq j] contain [N_{j}] as a factor; and[q_{j} Q_{j} \equiv 1 \hbox{ mod } N_{j}]by the defining relation for [q_{j}].

It may be noted that[\eqalign{(q_{i} Q_{i}) (q_{j} Q_{j}) &\equiv 0{\phantom{Q_j}}\quad\quad\hbox{ mod } N \hbox{ for } i \neq j,\cr (q_{j} Q_{j})^{2} &\equiv q_{j} Q_{j}\quad\quad\hbox{mod } N, \,\,j = 1, \ldots, d,}]so that the [q_{j} Q_{j}] are mutually orthogonal idempotents in the ring [{\bb Z}/N {\bb Z}], with properties formally similar to those of mutually orthogonal projectors onto subspaces in linear algebra. The analogy is exact, since by virtue of the CRT the ring [{\bb Z}/N {\bb Z}] may be considered as the direct product[{\bb Z}/N_{1} {\bb Z} \times {\bb Z}/N_{2} {\bb Z} \times \ldots \times {\bb Z}/N_{d} {\bb Z}]via the two mutually inverse mappings:

  • (i) [{\ell} \,\longmapsto\, (\ell_{1}, \ell_{2}, \ldots, \ell_{d})] by [\ell \equiv \ell_{j}] mod [N_{j}] for each j;

  • (ii) [(\ell_{1}, \ell_{2}, \ldots, \ell_{d}) \,\longmapsto\, \ell \hbox { by } \ell = {\textstyle\sum_{i = 1}^{d}} \ell_{i} q_{i} Q_{i}\hbox{ mod } N].

The mapping defined by (ii)[link] is sometimes called the `CRT reconstruction' of [\ell] from the [\ell_{j}].

These two mappings have the property of sending sums to sums and products to products, i.e:[\displaylines{\quad (\hbox{i})\quad\quad{\ell} + {\ell}' \,\longmapsto\, ({\ell}_{1} + {\ell}'_{1}, {\ell}_{2} + {\ell}'_{2}, \ldots, {\ell}_{d} + {\ell}'_{d}) \hfill\cr \quad\quad\quad\phantom{(\hbox{i})}{\ell} {\ell}' \,\longmapsto\, ({\ell}_{1} {\ell}'_{1}, {\ell}_{2} {\ell}'_{2}, \ldots, {\ell}_{d} {\ell}'_{d}) \quad\,\,\, \phantom{(\hbox{i})} \hfill\cr \quad (\hbox{ii}) \quad\quad\! ({\ell}_{1} + {\ell}'_{1}, {\ell}_{2} + {\ell}'_{2}, \ldots, {\ell}_{d} + {\ell}'_{d}) \,\longmapsto\, {\ell} + {\ell}' \,\,\hfill\cr \quad\quad\! \phantom{(\hbox{ii})}\quad({\ell}_{1} {\ell}'_{1}, {\ell}_{2} {\ell}'_{2}, \ldots, {\ell}_{d} {\ell}'_{d}) \,\longmapsto\, {\ell} {\ell}' \quad  \,\,\,\,\hfill}](the last proof requires using the properties of the idempotents [q_{j} Q_{j}]). This may be described formally by stating that the CRT establishes a ring isomorphism:[{\bb Z}/N {\bb Z} \cong ({\bb Z}/N_{1} {\bb Z}) \times \ldots \times ({\bb Z}/N_{d} {\bb Z}).]

1.3.3.2.2.3. The prime factor algorithm

| top | pdf |

The CRT will now be used to factor the N-point DFT into a tensor product of d transforms, the jth of length [N_{j}].

Let the indices k and [k^{*}] be subjected to the following mappings:

  • (i) [k \,\longmapsto\, (k_{1}, k_{2}, \ldots, k_{d}), k_{j} \in {\bb Z}/N_{j} {\bb Z}], by [k_{j} \equiv k] mod [N_{j}] for each j, with reconstruction formula[k = {\textstyle\sum\limits_{i = 1}^{d}} \,k_{i} q_{i} Q_{i} \hbox{ mod } N\hbox{\semi}]

  • (ii) [k^{*} \,\longmapsto\, (k_{1}^{*}, k_{2}^{*}, \ldots, k_{d}^{*}), k_{j}^{*} \in {\bb Z}/N_{j} {\bb Z}], by [k_{j}^{*} \equiv q_{j} k^{*}] mod [N_{j}] for each j, with reconstruction formula[k^{*} = {\textstyle\sum\limits_{i = 1}^{d}} \,k_{i}^{*} Q_{i} \hbox{ mod } N.]

Then[\eqalign{k^{*} k &= \left({\textstyle\sum\limits_{i = 1}^{d}}\, k_{i}^{*} Q_{i}\right) \left({\textstyle\sum\limits_{j = 1}^{d}} \,k_{j} q_{j} Q_{j}\right) \hbox{ mod } N\cr &= {\textstyle\sum\limits_{i, \, j = 1}^{d}} k_{i}^{*} k_{j} Q_{i} q_{j} Q_{j} \hbox{ mod } N.}]Cross terms with [i \neq j] vanish since they contain all the factors of N, hence[\eqalign{k^{*} k &= {\textstyle\sum\limits_{j = 1}^{d}}\, q_{j} Q_{j}^{2} k_{j}^{*} k_{j} \hbox{ mod } N\cr &= {\textstyle\sum\limits_{j = 1}^{d}} (1 - n_{j} N_{j}) Q_{j} k_{j}^{*} k_{j} \hbox{ mod } N.}]Dividing by N, which may be written as [N_{j} Q_{j}] for each j, yields[\eqalign{{k^{*} k \over N} &= {\sum\limits_{j = 1}^{d}} (1 - n_{j} N_{j}) {Q_{j} \over N_{j} Q_{j}} k_{j}^{*} k_{j} \hbox{ mod } 1\cr &= {\sum\limits_{j = 1}^{d}} \left({1 \over N_{j}} - n_{j}\right) k_{j}^{*} k_{j} \hbox{ mod } 1,}]and hence[{k^{*} k \over N} \equiv {\sum\limits_{j = 1}^{d}} {k_{j}^{*} k_{j} \over N_{j}} \hbox{ mod } 1.]Therefore, by the multiplicative property of [e(.)],[e \left({k^{*} k \over N}\right) \equiv \bigotimes\limits_{j = 1}^{d} e \left({k_{j}^{*} k_{j} \over N_{j}}\right).]

Let [{\bf X} \in L ({\bb Z}/N {\bb Z})] be described by a one-dimensional array [X(k)] indexed by k. The index mapping (i)[link] turns X into an element of [L ({\bb Z}/N_{1} {\bb Z} \times \ldots \times {\bb Z}/N_{d} {\bb Z})] described by a d-dimensional array [X (k_{1}, \ldots, k_{d})]; the latter may be transformed by [\bar{F} (N_{1}) \bigotimes \ldots \bigotimes \bar{F} (N_{d})] into a new array [X^{*} (k_{1}^{*}, k_{2}^{*}, \ldots, k_{d}^{*})]. Finally, the one-dimensional array of results [X^{*} (k^{*})] will be obtained by reconstructing [k^{*}] according to (ii)[link].

The prime factor algorithm, like the Cooley–Tukey algorithm, reindexes a 1D transform to turn it into d separate transforms, but the use of coprime factors and CRT index mapping leads to the further gain that no twiddle factors need to be applied between the successive transforms (see Good, 1971[link]). This makes up for the cost of the added complexity of the CRT index mapping.

The natural factorization of N for the prime factor algorithm is thus its factorization into prime powers: [\bar{F}(N)] is then the tensor product of separate transforms (one for each prime power factor [N_{j} = p_{j}^{\nu_{j}}]) whose results can be reassembled without twiddle factors. The separate factors [p_{j}] within each [N_{j}] must then be dealt with by another algorithm (e.g. Cooley–Tukey, which does require twiddle factors). Thus, the DFT on a prime number of points remains undecomposable.

References

Apostol, T. M. (1976). Introduction to Analytic Number Theory. New York: Springer-Verlag.
Good, I. J. (1971). The relationship between two fast Fourier transforms. IEEE Trans. C-20, 310–317.
Schroeder, M. R. (1986). Number Theory in Science and Communication, 2nd ed. Berlin: Springer-Verlag.








































to end of page
to top of page