Constructing plans

Every transform in this module is applied through a plan, which is constructed once for a length n, a ring and a root of unity w, and can then be used for any number of forward and inverse transforms. Plan construction chooses the algorithm, splits the length into the sub-plans that the algorithm needs, and computes the table of powers w^0, \ldots, w^{n-1} (or the parts of it that the algorithm reads). Its cost is comparable to that of a few transforms, so a plan should be reused whenever several vectors of the same length are transformed.

Which constructor#

gr_dft_default_root computes the root used by the last two constructors. It evaluates e^{\pi i q} with q = -2/n by gr_exp_pi_i, or, if the ring lacks that method, \exp(-2\pi i/n) from gr_pi, gr_i and gr_exp. Over a ring without these operations, or without the value (real numbers for n \ge 3), the constructors return its status, GR_UNABLE or GR_DOMAIN. For n \le 2 the default root \pm 1 exists in every ring.

With the default root, the constructors do not compute the table from powers of w: over acb each entry is rounded from a running product at a precision raised by about \log_2 n + 16 bits, completed by the symmetries of the unit circle, and over other rings it is computed entry by entry with gr_exp_pi_i (over the fixed-point contexts, see Plans). The entries therefore have about the accuracy of single evaluations, and over acb those equal to \pm 1 and \pm i are exact.

Roots of unity#

gr_dft_precomp_init_root requires a principal n-th root of unity: w^n = 1, and w^{n/p} - 1 is not a zero divisor for every prime p \mid n. Over a field this means that w has order n. Over \mathbb{Z}/m\mathbb{Z} with composite m it means that w has order n modulo every prime factor of m. Roots of unity in a ring explains why the fast algorithms need this. The constructor checks two necessary conditions, using gr_is_neg_one and gr_is_one on the table entries:

• if n is even, w^{n/2} = -1;
• w^{n/p} \ne 1 for every odd prime p \mid n.

A definite failure gives GR_DOMAIN. When the ring cannot decide (an inexact ring returns T_UNKNOWN), the root is accepted. The condition w^n = 1 itself is not checked; the table is simply w^0, \ldots, w^{n-1}. For n a power of two the checks are complete (and only w^{n/2} = -1 is needed); in a ring with zero divisors they are not, as the second example shows: the algorithms then return wrong results without any error. The constructors with default roots perform no checks.

Over \mathbb{Z}/p\mathbb{Z} with p prime, a root of order n exists if and only if n \mid p - 1, and g^{(p-1)/n} is one for any primitive root g (n_primitive_root_prime computes one). Primes p = c 2^k + 1 with large k, such as 998244353 = 119 \cdot 2^{23} + 1, allow long power-of-two transforms.

Over the complex fixed-point contexts (gr_dft_ctx_init_nfixed_complex), the root table is always that of e^{-2\pi i/n}, computed to the full precision of the context, and w only selects this root: it is accepted if it lies within 1/n of e^{-2\pi i/n}, so a value rounded to a lower precision is fine. gr_dft_precomp_init_root returns GR_UNABLE if w approximates another primitive n-th root of unity, and GR_DOMAIN if it approximates none.

Complex mode#

gr_dft_precomp_init_karatsuba with a real ring real_ctx builds a plan in complex mode. The ring ctx must store each element as a real part followed by an imaginary part, both elements of real_ctx; only the element sizes are checked: GR_UNABLE is returned unless an element of ctx is exactly twice as large as an element of real_ctx. With real_ctx equal to NULL the function is gr_dft_precomp_init. In complex mode:

• multiplications by \pm 1 and \pm i (the powers w^{jn/4}) are done by swapping and negating parts, with no multiplication;
• multiplications by w^{jn/8} for odd j, which have the form c(1 \pm i), take two real multiplications;
• multiplications by other roots take three real multiplications with precomputed tables (the complex Karatsuba product), but only over exact rings and over rings whose gr_ctx_get_real_prec is at least 1024 bits, and never over the fixed-point contexts. Otherwise they take one multiplication in ctx;
• GR_DFT_ALG_AUTO chooses split-radix rather than radix-2 for powers of two, unless GR_DFT_SCRAMBLED is set.

The special values are recognised from their exponents, which is valid because complex mode always uses w = e^{-2\pi i/n}; over balls a comparison could not prove that a computed value is exactly -i. Since the formulas for c(1 \pm i) use only the enclosure of c, they still give enclosures of the exact products. The gain is largest for small powers of two: a split-radix transform of length 16 in complex mode needs 20 real multiplications.

Failure and lifetime#

On any return value other than GR_SUCCESS, the plan is left in a cleared state: it must not be used for transforms, and calling gr_dft_precomp_clear on it is harmless. The plan stores pointers to ctx and real_ctx (sub-plans included), so the contexts must outlive the plan, and transforms must be given the same ctx. A plan is not modified by transforms.

The status values of the constructors are:

Example#

The number 2 is a primitive root modulo 13, so it is a principal 12th root of unity, while 4 has order 6. The automatic plan of length 12 = 4 \cdot 3 uses the prime factor algorithm.

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

int main(void)
{
    gr_ctx_t ctx;
    gr_dft_pre_t P;
    gr_ptr x, y, w;
    slong j, n = 12;
    int status;

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

    /* 4 has order 6 modulo 13: not a principal 12th root of unity */
    GR_MUST_SUCCEED(gr_set_ui(w, 4, ctx));
    status = gr_dft_precomp_init_root(P, w, n, GR_DFT_ALG_AUTO, 0, ctx);
    flint_printf("w = 4: %s\n", status == GR_DOMAIN ? "GR_DOMAIN" : "?");
    gr_dft_precomp_clear(P);   /* harmless after a failure */

    /* 2 is a primitive root modulo 13 */
    GR_MUST_SUCCEED(gr_set_ui(w, 2, ctx));
    GR_MUST_SUCCEED(gr_dft_precomp_init_root(P, w, n, GR_DFT_ALG_AUTO, 0, ctx));

    for (j = 0; j < n; j++)
        GR_MUST_SUCCEED(gr_set_si(GR_ENTRY(x, j, ctx->sizeof_elem), j, ctx));

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

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

    gr_dft_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;
}

w = 4: GR_DOMAIN
[1, 12, 4, 11, 6, 5, 7, 9, 8, 3, 10, 2]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]

Example: a root that is not principal#

Modulo 91 = 7 \cdot 13, the number w = 79 satisfies w^3 = 1 and w \ne 1, so it passes the checks for n = 3. But w \equiv 1 \pmod{13}, so w - 1 is a zero divisor and 1 + w + w^2 = 42 \ne 0. The naive algorithm evaluates the definition and is correct for any w with w^n = 1, while the mixed-radix kernel chosen by GR_DFT_ALG_AUTO uses 1 + w + w^2 = 0 and returns a different vector.

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

/* w = 79 satisfies w^3 = 1 and w != 1 modulo 91 = 7 * 13,
   but w - 1 = 78 is a zero divisor */
int main(void)
{
    gr_ctx_t ctx;
    gr_dft_pre_t P1, P2;
    gr_ptr x, y, w;
    slong j, n = 3;

    gr_ctx_init_nmod(ctx, 91);
    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, 79, ctx));
    for (j = 0; j < n; j++)
        GR_MUST_SUCCEED(gr_set_si(GR_ENTRY(x, j, ctx->sizeof_elem), j + 1, ctx));

    /* both constructors accept w */
    GR_MUST_SUCCEED(gr_dft_precomp_init_root(P1, w, n, GR_DFT_ALG_NAIVE, 0, ctx));
    GR_MUST_SUCCEED(gr_dft_precomp_init_root(P2, w, n, GR_DFT_ALG_AUTO, 0, ctx));

    GR_MUST_SUCCEED(gr_dft_precomp(y, x, P1, ctx));
    flint_printf("naive: ");
    _gr_vec_print(y, n, ctx);
    GR_MUST_SUCCEED(gr_dft_precomp(y, x, P2, ctx));
    flint_printf("\nauto:  ");
    _gr_vec_print(y, n, ctx);
    flint_printf("\n");

    gr_dft_precomp_clear(P1);
    gr_dft_precomp_clear(P2);
    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;
}

naive: [6, 45, 71]
auto:  [6, 10, 78]

The naive result is the transform as defined, and the automatic one differs from it by -(1 + w + w^2)(0, x_2, x_1) = (0, -126, -84) modulo 91. Neither is useful for convolutions, because for this w the inverse transform does not undo the forward one.