acb_dft.h – Discrete Fourier transform

Warning: the interfaces in this module are experimental and may change without notice.

All functions support aliasing.

Let G be a finite abelian group, and \(\chi\) a character of G. For any map \(f:G\to\mathbb C\), the discrete fourier transform \(\hat f:\hat G\to \mathbb C\) is defined by

\[\hat f(\chi) = \sum_{x\in G}\overline{\chi(x)}f(x)\]

Note that by the inversion formula

\[\widehat{\hat f}(\chi) = \#G\times f(\chi^{-1})\]

it is straightforward to recover \(f\) from its DFT \(\hat f\).

Main DFT functions

If \(G=\mathbb Z/n\mathbb Z\), we compute the DFT according to the usual convention

\[w_x = \sum_{y\bmod n} v_y e^{-\frac{2i \pi}nxy}\]
void acb_dft(acb_ptr w, acb_srcptr v, slong n, slong prec)

Set w to the DFT of v of length len, using an automatic choice of algorithm.

void acb_dft_inverse(acb_ptr w, acb_srcptr v, slong n, slong prec)

Compute the inverse DFT of v into w.

If several computations are to be done on the same group, the FFT scheme should be reused.

type acb_dft_pre_struct
type acb_dft_pre_t

Stores a fast DFT scheme on \(\mathbb Z/n\mathbb Z\) as a recursive decomposition into simpler DFT with some tables of roots of unity.

An acb_dft_pre_t is defined as an array of acb_dft_pre_struct of length 1, permitting it to be passed by reference.

void acb_dft_precomp_init(acb_dft_pre_t pre, slong len, slong prec)

Initializes the fast DFT scheme of length len, using an automatic choice of algorithms depending on the factorization of len.

The length len is stored as pre->n.

void acb_dft_precomp_clear(acb_dft_pre_t pre)

Clears pre.

void acb_dft_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec)

Computes the DFT of the sequence v into w by applying the precomputed scheme pre. Both v and w must have length pre->n.

void acb_dft_inverse_precomp(acb_ptr w, acb_srcptr v, const acb_dft_pre_t pre, slong prec)

Compute the inverse DFT of v into w.

DFT on products

A finite abelian group is isomorphic to a product of cyclic components

\[G = \bigoplus_{i=1}^r \mathbb Z/n_i\mathbb Z\]

Characters are product of component characters and the DFT reads

\[\hat f(x_1,\dots x_r) = \sum_{y_1\dots y_r} f(y_1,\dots y_r) e^{-2i \pi \sum\frac{x_i y_i}{n_i}}\]

We assume that \(f\) is given by a vector of length \(\prod n_i\) corresponding to a lexicographic ordering of the values \(y_1,\dots y_r\), and the computation returns the same indexing for values of \(\hat f\).

void acb_dirichlet_dft_prod(acb_ptr w, acb_srcptr v, slong *cyc, slong num, slong prec)

Computes the DFT on the group product of num cyclic components of sizes cyc. Assume the entries of v are indexed according to lexicographic ordering of the cyclic components.

type acb_dft_prod_struct
type acb_dft_prod_t

Stores a fast DFT scheme on a product of cyclic groups.

An acb_dft_prod_t is defined as an array of acb_dft_prod_struct of length 1, permitting it to be passed by reference.

void acb_dft_prod_init(acb_dft_prod_t t, slong *cyc, slong num, slong prec)

Stores in t a DFT scheme for the product of num cyclic components whose sizes are given in the array cyc.

void acb_dft_prod_clear(acb_dft_prod_t t)

Clears t.

void acb_dirichlet_dft_prod_precomp(acb_ptr w, acb_srcptr v, const acb_dft_prod_t prod, slong prec)

Sets w to the DFT of v. Assume the entries are lexicographically ordered according to the product of cyclic groups initialized in t.

Convolution

For functions \(f\) and \(g\) on \(G\) we consider the convolution

\[(f \star g)(x) = \sum_{y\in G} f(x-y)g(y)\]
void acb_dft_convol_naive(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)
void acb_dft_convol_rad2(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)
void acb_dft_convol(acb_ptr w, acb_srcptr f, acb_srcptr g, slong len, slong prec)

Sets w to the convolution of f and g of length len.

The naive version simply uses the definition.

The rad2 version embeds the sequence into a power of 2 length and uses the formula

\[\widehat{f \star g}(\chi) = \hat f(\chi)\hat g(\chi)\]

to compute it using three radix 2 FFT.

The default version uses radix 2 FFT unless len is a product of small primes where a non padded FFT is faster.

FFT algorithms

Fast Fourier transform techniques allow to compute efficiently all values \(\hat f(\chi)\) by reusing common computations.

Specifically, if \(H\triangleleft G\) is a subgroup of size \(M\) and index \([G:H]=m\), then writing \(f_x(h)=f(xh)\) the translate of \(f\) by representatives \(x\) of \(G/H\), one has a decomposition

\[\hat f(\chi) = \sum_{x\in G/H} \overline{\chi(x)} \hat{f_x}(\chi_{H})\]

so that the DFT on \(G\) can be computed using \(m\) DFT on \(H\) (of appropriate translates of \(f\)), then \(M\) DFT on \(G/H\), one for each restriction \(\chi_{H}\).

This decomposition can be done recursively.

Naive algorithm

void acb_dft_naive(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of v into w, where v and w have size n, using the naive \(O(n^2)\) algorithm.

type acb_dft_naive_struct
type acb_dft_naive_t
void acb_dft_naive_init(acb_dft_naive_t t, slong len, slong prec)
void acb_dft_naive_clear(acb_dft_naive_t t)

Stores a table of roots of unity in t. The length len is stored as t->n.

void acb_dft_naive_precomp(acb_ptr w, acb_srcptr v, const acb_dft_naive_t t, slong prec)

Sets w to the DFT of v of size t->n, using the naive algorithm data t.

CRT decomposition

void acb_dft_crt(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of v into w, where v and w have size len, using CRT to express \(\mathbb Z/n\mathbb Z\) as a product of cyclic groups.

type acb_dft_crt_struct
type acb_dft_crt_t
void acb_dft_crt_init(acb_dft_crt_t t, slong len, slong prec)
void acb_dft_crt_clear(acb_dft_crt_t t)

Initialize a CRT decomposition of \(\mathbb Z/n\mathbb Z\) as a direct product of cyclic groups. The length len is stored as t->n.

void acb_dft_crt_precomp(acb_ptr w, acb_srcptr v, const acb_dft_crt_t t, slong prec)

Sets w to the DFT of v of size t->n, using the CRT decomposition scheme t.

Cooley-Tukey decomposition

void acb_dft_cyc(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of v into w, where v and w have size n, using each prime factor of \(m\) of \(n\) to decompose with the subgroup \(H=m\mathbb Z/n\mathbb Z\).

type acb_dft_cyc_struct
type acb_dft_cyc_t
void acb_dft_cyc_init(acb_dft_cyc_t t, slong len, slong prec)
void acb_dft_cyc_clear(acb_dft_cyc_t t)

Initialize a decomposition of \(\mathbb Z/n\mathbb Z\) into cyclic subgroups. The length len is stored as t->n.

void acb_dft_cyc_precomp(acb_ptr w, acb_srcptr v, const acb_dft_cyc_t t, slong prec)

Sets w to the DFT of v of size t->n, using the cyclic decomposition scheme t.

Radix 2 decomposition

void acb_dft_rad2(acb_ptr w, acb_srcptr v, int e, slong prec)

Computes the DFT of v into w, where v and w have size \(2^e\), using a radix 2 FFT.

void acb_dft_inverse_rad2(acb_ptr w, acb_srcptr v, int e, slong prec)

Computes the inverse DFT of v into w, where v and w have size \(2^e\), using a radix 2 FFT.

type acb_dft_rad2_struct
type acb_dft_rad2_t
void acb_dft_rad2_init(acb_dft_rad2_t t, int e, slong prec)
void acb_dft_rad2_clear(acb_dft_rad2_t t)

Initialize and clear a radix 2 FFT of size \(2^e\), stored as t->n.

void acb_dft_rad2_precomp(acb_ptr w, acb_srcptr v, const acb_dft_rad2_t t, slong prec)

Sets w to the DFT of v of size t->n, using the precomputed radix 2 scheme t.

Bluestein transform

void acb_dft_bluestein(acb_ptr w, acb_srcptr v, slong n, slong prec)

Computes the DFT of v into w, where v and w have size n, by conversion to a radix 2 one using Bluestein’s convolution trick.

type acb_dft_bluestein_struct
type acb_dft_bluestein_t

Stores a Bluestein scheme for some length n : that is a acb_dft_rad2_t of size \(2^e \geq 2n-1\) and a size n array of convolution factors.

void acb_dft_bluestein_init(acb_dft_bluestein_t t, slong len, slong prec)
void acb_dft_bluestein_clear(acb_dft_bluestein_t t)

Initialize and clear a Bluestein scheme to compute DFT of size len.

void acb_dft_bluestein_precomp(acb_ptr w, acb_srcptr v, const acb_dft_bluestein_t t, slong prec)

Sets w to the DFT of v of size t->n, using the precomputed Bluestein scheme t.