Discrete Fourier transforms over a generic ring

The module gr_dft computes the discrete Fourier transform

of a vector x \in R^n over an arbitrary commutative ring R, given an element w \in R with w^n = 1. Unlike the transforms of fft_small and fft, which are specialised to particular rings and lengths, it works with any length and only uses the ring operations of gr, so the same code serves complex balls, exact algebraic numbers, \mathbb{Z}/m\mathbb{Z}, finite fields and fixed-point numbers. The classical references are the survey of Duhamel and Vetterli [DV1990] and the book of Nussbaumer [Nus1982]; transforms over general rings are treated by von zur Gathen and Gerhard [GG2003, §8.2], Pollard [Pol1971] and Bernstein [Ber2001].

Roots of unity in a ring#

Over \mathbb{C} the fast algorithms use identities of the roots of unity such as w^{n/2} = -1 and 1 + v + \cdots + v^{p-1} = 0 for a root v of order p. In a general ring these identities do not follow from w having order n, because the ring can have zero divisors.

Over \mathbb{Z}/m\mathbb{Z} the lemma says that w must have order exactly n modulo every prime factor of m. Over \mathbb{Z}/p\mathbb{Z} with p prime, such roots exist exactly when n \mid p - 1.

The inverse transform of gr_dft is therefore the transform with w^{-1} followed by a division by n, and it fails with GR_DOMAIN when the ring cannot divide by n. The same computation shows that the transform maps cyclic convolution to pointwise multiplication.

The algorithms below need less than principality to compute the forward transform as defined. Their derivations use only the following identities:

The prime factor algorithm passes the roots w^{n_2} and w^{n_1} to its sub-plans, which need their own identities. For n = 2^d, the condition w^{n/2} = -1 already implies orthogonality: for k = 2^s k' with k' odd and s < d, the terms j and j + n/2^{s+1} of \sum_j w^{jk} cancel, because w^{k n / 2^{s+1}} = (w^{n/2})^{k'} = -1. So for powers of two w^{n/2} = -1 is all that is needed.

gr_dft_precomp_init_root checks w^{n/2} = -1 for even n and w^{n/p} \ne 1 for odd primes p \mid n, as far as the ring can decide them with gr_is_neg_one and gr_is_one. The second check is necessary but not sufficient: over a ring with zero divisors, w^{n/p} \ne 1 does not imply 1 + v + \cdots + v^{p-1} = 0, and the mixed-radix kernel then returns a wrong result without an error. Example: a root that is not principal gives an example modulo 91.

The naive algorithm#

GR_DFT_ALG_NAIVE evaluates the definition, accumulating x_j w^{jk \bmod n} with the exponent updated incrementally. It performs n(n-1) calls to the root multiplication (those with exponent 0 are copies), is correct for any w with w^n = 1, and serves as the reference in the tests. For powers of two it also supports the bit-reversed order of GR_DFT_SCRAMBLED.

Radix-2 Cooley–Tukey#

For n = 2^d, splitting the input into halves (decimation in frequency, [GS1966]) gives

because w^{(j + n/2) \cdot 2k} = w^{2jk} and w^{(j + n/2)(2k + 1)} = -w^{j(2k+1)} by w^n = 1 and w^{n/2} = -1. Since (w^2)^{n/4} = w^{n/2} = -1, the two halves are transforms of the same kind with root w^2, and induction gives the algorithm of Cooley and Tukey [CT1965].

The forward transform of gr_dft uses Algorithm 1 followed by a bit-reversal permutation, which GR_DFT_SCRAMBLED skips. The inverse transform runs the transposed algorithm (decimation in time, (a, b) \mapsto (a + w^{-jr} b, a - w^{-jr} b) for s = 1, \ldots, d) on bit-reversed input, so a scrambled forward transform followed by a scrambled inverse transform gives n times the identity without any permutation. There are \frac{n}{2} d butterflies, each with one multiplication by a root; those with j = 0 need none, and in complex mode those with jr \in \{n/8, n/4, 3n/8\} are cheap (Complex mode). The integer FFTs of fft (The ring) and fft_small (FFTs modulo word-size primes in floating point) are built on the same recursion, specialised to their rings, and truncate it so that only the needed values are computed (Truncated FFTs).

Split-radix#

The split-radix algorithm of Yavne [Yav1968] and Duhamel and Hollmann [DH1984] splits a transform of length m with root v into one transform U of length m/2 (the even entries) and two transforms Z, Z' of length m/4 (the entries \equiv 1 and \equiv 3 \pmod 4). With t = v^k Z_k, t' = v^{3k} Z'_k and \iota = v^{m/4},

for 0 \le k < m/4, which again uses only v^{m/2} = -1. Each step costs two multiplications by general roots per k plus one by \iota, which plays the role of -i and is free in complex mode. The total is about \frac{n}{3} \log_2 n general multiplications instead of \frac{n}{2} \log_2 n. In complex mode, with three real multiplications per general root and two per eighth root, it uses n \log_2 n - 3n + 4 real multiplications, which is minimal for n \le 16 [Win1978] and close to the lower bound of Heideman and Burrus [HB1986] for larger n. Later variants reduce the number of additions further [JF2007]; they are not used. gr_dft computes split-radix recursively, out of place, in natural order.

The four-step algorithm#

For n = n_1 n_2 with n_1 = 2^{\lfloor d/2 \rfloor}, Bailey's four-step algorithm [Bai1990] writes j = j_1 + n_1 j_2, k = k_2 + n_2 k_1 and

which gives n_1 transforms of length n_2, n twiddle multiplications by w^{j_1 k_2}, n_2 transforms of length n_1 (all radix-2 sub-plans) and a transposition, which GR_DFT_SCRAMBLED skips. The method and its variants are described in The matrix Fourier algorithm. In gr_dft it is available as GR_DFT_ALG_BAILEY but not chosen automatically.

Mixed-radix Cooley–Tukey#

For a transform of length \ell = pm with p prime and root u, split the input into its p residue classes modulo p and let Y_b be the transform of length m, with root u^p, of the class b. Then

so after twiddling by u^{bq}, each q requires a transform of prime length p with root v. GR_DFT_ALG_MIXED applies this recursively, one prime factor of n per level, with all exponents taken in the root table of the top-level plan. The transforms of prime length use the following kernel.

The output for r = 0 is \sum_b z_b, so the kernel takes (p-1)(p-2) multiplications instead of (p-1)^2. For p = 3, for example, X_1 = (x_0 - x_2) + v(x_1 - x_2) and X_2 = (x_0 - x_1) + v(x_2 - x_1). With the (p - 1) twiddles per kernel, a transform of length p^e costs at most e\, p^{e-1} (p-1)^2 multiplications, which is O(n \log n) only for bounded p. For n = p^e with p > 19, the transforms of length p use Bluestein's algorithm instead when it is available.

The prime factor algorithm#

If n = n_1 n_2 with \gcd(n_1, n_2) = 1, the Chinese remainder theorem removes the twiddle factors altogether. This is the algorithm of Good [Goo1958] and Thomas [Tho1963].

The right side is a two-dimensional transform: n_1 transforms of length n_2 with root w^{n_1} followed by n_2 transforms of length n_1 with root w^{n_2}, with a gather before and a scatter after but no multiplications in between. If w is principal, so are w^{n_2} and w^{n_1}, since (w^{n_2})^{n_1/p} - 1 = w^{n/p} - 1. GR_DFT_ALG_PFA takes n_1 to be the full power of the smallest prime factor of n and plans n_2 automatically, so a general length is reduced to its prime power factors. With the default root, its sub-plans build their own tables and the plan itself stores none.

Bluestein's algorithm#

Bluestein's algorithm [Blu1970], also known as the chirp z-transform [RSR1969], reduces a transform of any odd length n to a cyclic convolution, which can then be computed with transforms of a convenient length.

The cyclic convolution of length n is embedded in one of length N = 2^{\lceil \log_2(2n - 1) \rceil}: the vector (c_j x_j) is padded with zeros, and the kernel is h'_j = h_j for 0 \le j < n, h'_{N - j} = h_j for 1 \le j < n and 0 otherwise. Because N \ge 2n - 1 the two parts do not overlap, and for 0 \le k < n the cyclic convolution of length N agrees with that of length n by the evenness of h.

The chirps c_j and h_j are entries of the root table of w, and \hat h is in the same scrambled order as Z, so no permutation is needed. The inverse transform of length n is computed as the forward transform of x_{(n - j) \bmod n}, since \sum_j x_j w^{-jk} = \sum_j x_{(n-j) \bmod n} w^{jk}. The cost is 2n chirp multiplications, N pointwise products and two transforms of length N < 4n; the precomputation costs one more transform. The convolution needs a root of order N, unrelated to w, which gr_dft obtains from gr_dft_default_root. Over rings without default roots, such as \mathbb{Z}/m\mathbb{Z} (even when m has roots of order N), the algorithm is therefore not available, and mixed-radix plans fall back to the direct kernel. Over the fixed-point contexts the kernel is scaled by 1/(2N) instead of 1/N to keep the intermediate values below 1, and the outputs are doubled at the end.

Products of cyclic groups and Dirichlet characters#

The DFT on \mathbb{Z}/c_0\mathbb{Z} \times \cdots \times \mathbb{Z}/c_{r-1}\mathbb{Z}, X_k = \sum_j x_j \prod_a w_a^{j_a k_a} with w_a a principal c_a-th root of unity, factors over the axes. It is computed by transforming all lines along each axis in turn with a cyclic plan of length c_a, with no twiddle factors, in \sum_a (n / c_a) C(c_a) operations for n = \prod c_a, where C(m) is the cost of a transform of length m. Plans are shared between axes of equal length, and the inverse divides by n once. By Theorem 6, a cyclic transform of coprime length n_1 n_2 is a product transform up to permutations of the input and the output.

The group of Dirichlet characters modulo q is isomorphic to (\mathbb{Z}/q\mathbb{Z})^\times, which the Conrey decomposition writes as a product of cyclic groups with generators of orders \varphi_k. If x has logarithms b_k and the character \chi_a is indexed by the logarithms a_k of a, then \chi_a(x) = \prod_k e^{2\pi i a_k b_k / \varphi_k}, and the transform \hat v_a = \sum_x \overline{\chi_a(x)} v_x is the product transform over the \varphi_k in the lexicographic order of the logarithms. gr_dft_dirichlet adds the permutation between this order and the order of the residues x.

Complex mode#

When the ring is C = R[i] with elements stored as pairs of elements of R and w = e^{-2\pi i/n}, multiplications by roots can be specialised from the exponent alone, without inspecting values (which over balls could not prove that a computed value is exactly -i):

• w^{jn/4} \in \{1, -i, -1, i\}: swaps and negations of the parts;
• w^{jn/8} for odd j equals c(1 \pm i) with c = \pm\sqrt{2}/2, so that, for example, c(1 + i)(x + iy) = c(x - y) + i\, c(x + y) takes two real multiplications;
• for a general root c + di, the product with x + iy is k_1 - y(d + c) + i\,(k_1 + x(d - c)) with k_1 = c(x + y): three real multiplications and three additions with the precomputed c, d - c and d + c, instead of four multiplications and two additions. This form is attributed to Gauss [Knu1997, §4.6.4].

The last trick pays when multiplications in R are expensive compared to additions, so gr_dft uses it over exact rings and when the real precision is at least 1024 bits; below that, a native complex multiplication is faster. In fixed point the entries d \pm c can reach \sqrt 2 and do not fit the representation, so complex products use the schoolbook formula below 22 limbs and a carry-extended three-product variant from 22 limbs.

Error handling#

The algorithms only call gr operations and combine their status codes with |. Over exact rings the forward transforms cannot fail, but the division by n of the inverse can, with GR_DOMAIN, as can the division by N in the Bluestein precomputation. Over inexact rings the root checks may be inconclusive (T_UNKNOWN), in which case the root is accepted. Over the fixed-point contexts, an overflow during a transform returns GR_UNABLE; Rigorous numerical DFTs in fixed point shows how to scale the input so that this cannot happen and bounds the rounding errors.

Summary#