Multidimensional transforms

The DFT on a finite abelian group $G = \mathbb{Z}/c_0\mathbb{Z} \times \cdots \times \mathbb{Z}/c_{r-1}\mathbb{Z}$ is the $r$-dimensional DFT of an array with $c_0 \cdots c_{r-1}$ entries. These functions compute it with one cyclic plan per distinct component length. They are the basis of the transforms on groups of Dirichlet characters (gr_dft_dirichlet), whose components are given by the Conrey decomposition and often repeat.

Layout and algorithm#

The array is stored in row-major order: the entry with index $(j_0, \ldots, j_{r-1})$ is at position $j_{r-1} + c_{r-1}(j_{r-2} + c_{r-2}(\cdots + c_1 j_0))$, so the last index varies fastest. The transform is computed axis by axis: for each axis $a$, every line of $c_a$ entries along that axis is transformed with the cyclic plan of length $c_a$, in place for the last axis and through a buffer for the others, which are strided. Since

factors over the axes, there are no twiddle factors between the axes, in contrast to the Cooley–Tukey decomposition of a cyclic transform of the same total length. The inverse transform applies the inverse cyclic transforms without their divisions, and divides by $n = \prod c_a$ once at the end. Axes of length 1 are skipped.

The order of the axes does not matter mathematically, and the product transform is not the cyclic transform of length $n$, even when the $c_a$ are pairwise coprime: the Good–Thomas algorithm reduces the cyclic transform to a product transform only after permuting the input and the output (The prime factor algorithm).

Roots#

gr_dft_prod_precomp_init and the one-shot functions use the default roots $e^{-2\pi i/c_a}$, so they need a ring such as acb. Over other rings, gr_dft_prod_precomp_init_root takes one principal root $w$ of some order $N$ divisible by every $c_a$, for example $N = \mathrm{lcm}(c_0, \ldots, c_{r-1})$, and uses $w^{N/c_a}$ for the axis of length $c_a$. Each of these is then a principal $c_a$-th root of unity, and the component plans check them as described in Roots of unity.

Threads#

The lines of one axis are independent, and when the ring is thread-safe and flint_get_num_threads is greater than 1, they are distributed over at most $n / \mathit{serial\_block}$ chunks. Threads are requested from the global pool for each axis separately, and only for an axis with at least two lines; an axis with a single line, as in a product with one long axis, is transformed by its component plan, which can then use the threads of the pool itself. The granularity is set with gr_dft_prod_precomp_set_serial_block (default GR_DFT_SERIAL_BLOCK_DEFAULT).

Example#

The transform on $\mathbb{Z}/2\mathbb{Z} \times \mathbb{Z}/3\mathbb{Z}$ over $\mathbb{Z}/7\mathbb{Z}$, with the root 3 of order 6: the axes use $3^3 = -1$ and $3^2 = 2$. The input is the $2 \times 3$ array with rows $(1, 2, 3)$ and $(4, 5, 6)$.

#include <flint/flint.h>
#include <flint/gr.h>
#include <flint/gr_vec.h>
#include <flint/gr_dft.h>

/* DFT on Z/2 x Z/3 over Z/7Z, using w = 3 of order 6 */
int main(void)
{
    gr_ctx_t ctx;
    gr_dft_prod_pre_t P;
    gr_ptr x, y, w;
    ulong cyc[2] = { 2, 3 };
    slong j, n = 6;

    gr_ctx_init_nmod(ctx, 7);
    x = gr_heap_init_vec(n, ctx);
    y = gr_heap_init_vec(n, ctx);
    w = gr_heap_init(ctx);

    GR_MUST_SUCCEED(gr_set_ui(w, 3, ctx));
    GR_MUST_SUCCEED(gr_dft_prod_precomp_init_root(P, w, 6, cyc, 2, 0, ctx));

    /* entry 3 j0 + j1 holds x[j0][j1] */
    for (j = 0; j < n; j++)
        GR_MUST_SUCCEED(gr_set_si(GR_ENTRY(x, j, ctx->sizeof_elem), j + 1, ctx));

    GR_MUST_SUCCEED(gr_dft_prod_precomp(y, x, P, ctx));
    _gr_vec_print(y, n, ctx);
    flint_printf("\n");

    GR_MUST_SUCCEED(gr_dft_prod_inverse_precomp(y, y, P, ctx));
    _gr_vec_print(y, n, ctx);
    flint_printf("\n");

    gr_dft_prod_precomp_clear(P);
    gr_heap_clear(w, ctx);
    gr_heap_clear_vec(x, n, ctx);
    gr_heap_clear_vec(y, n, ctx);
    gr_ctx_clear(ctx);
    return 0;
}

[0, 6, 2, 5, 0, 0]
[1, 2, 3, 4, 5, 6]

For instance, the entry $(1, 0)$ is $(1 + 2 + 3) - (4 + 5 + 6) = -9 \equiv 5$.