Where FLINT uses fft_small

fft_small is FLINT's main engine for large products. Integer, polynomial and matrix arithmetic reach it through their own dispatchers, which compare the operand sizes with tuned thresholds; this guide collects those thresholds and explains what a call costs. The module exists only when FLINT_HAVE_FFT_SMALL is 1 in flint-config.h, which configure sets on 64-bit x86 with AVX2 and on 64-bit ARM with NEON. Elsewhere the same functions multiply large integers with GMP's mpn_mul and, from FLINT_FFT_MUL_THRESHOLD (32000) limbs, the Schönhage–Strassen code of the fft module.

Callers and thresholds#

Thresholds depend on the target and are defined in flint-mparam.h. The first value is that of the Intel targets from Skylake on, which an installation for a recent x86-64 processor uses; the others are given in parentheses.

The polynomial thresholds come from the table fft_mul_tab in nmod_poly/mullow_fft_small.c: a full product of lengths \ell_1 \ge \ell_2 uses fft_small when \min(\ell_1, 2 \ell_2) reaches the entry for the bit size of n, which is between 80 and 1326, for example 306 for 8 bits, 99 for 16 bits, 176 for 32 bits, 162 for 50 bits and 94 for 64 bits. The irregular entries reflect the Kronecker substitution variants that compete with it. Squares and low products have their own tables; moduli that are context primes use fft_small from \min(\ell_1, 2 \ell_2) = 24 (40 for low products), and moduli n \le 23 first pack two coefficients into one (from total lengths of about 370). Middle products use the threshold of the product they are the transpose of, and _nmod_poly_mulmid_fft_small computes them as that transposed product (_nmod_poly_mul_tr_default_mpn_ctx) when its transform is shorter than that of the window of the full product (Transposed products).

Several other modules build their own drivers on the plans, operands and the convolution engine of this module, or on the transforms directly:

• the transformed integer and polynomial rings of gr (gr_ctx_init_transformed_mpn and the transformed polynomial rings of nmod and mpn_mod), used by fmpz_powm and fmpz_mod_pow for moduli of at least 70000 bits, fmpz_mat_mul_fft_small, nmod_poly_mat_mulmid_fft_small, the matrix products of the polynomial half-gcd in gr_poly, and division with precomputed inverse over nmod from length 160;
• mpn_mod polynomial multiplication and the middle products of the decimal radix arithmetic, through _fft_small_conv (for radix, in its transposed mode for windows shorter than the longer operand, see radix_mulmid_fft_small);
• multipoint resultants in nmod_mpoly and fmpz_mpoly_factor, which transform with sd_fft_trunc modulo primes accepted by fft_small_prime_supports_depth, and fmpz_poly_gcd_modular, which tries the context primes first because products modulo them need a single transform.

Which pipeline and how many primes#

For integers, operands of up to 40000 limbs use two primes and chunks of 40 to 45 bits; larger ones use four to eight primes and chunks of 84 to 192 bits (Multiplying large integers). For polynomials over \mathbb{Z}/n\mathbb{Z} the number of primes grows with the bit size of n and, slowly, with the length (Number of primes): one prime for moduli of up to about 16 bits, two up to about 44 bits, three above. Moduli that are among the eight context primes, and primes n of at least 20 bits with a large power of two dividing n - 1, need a single transform without Chinese remaindering. The cost of a product is roughly proportional to the number of primes times the transform length.

The plan constructors make the same choices as the multiplication drivers (for one thread) and can be used to inspect them:

#include <flint/flint.h>
#include <flint/nmod.h>
#include <flint/fft_small.h>

int main(void)
{
    mpn_ctx_struct * R = get_default_mpn_ctx();
    slong sizes[5] = { 1000, 40000, 100000, 1000000, 10000000 };
    int bits[4] = { 16, 32, 50, 64 };
    int i;

    /* integer products of two n-limb operands */
    for (i = 0; i < 5; i++)
    {
        fft_small_plan_t P;
        if (fft_small_plan_init_mpn(P, R, sizes[i], sizes[i], 1, 0))
            flint_printf("%8wd limbs: %wu primes, %3wu-bit chunks, length 2^%wu\n",
                         sizes[i], P->np, P->bits, P->depth);
        fft_small_plan_clear(P);
    }

    /* polynomials of length 10^5 modulo n = 2^k - 1 */
    for (i = 0; i < 4; i++)
    {
        fft_small_plan_t P;
        nmod_t mod;
        ulong len = 100000;

        nmod_init(&mod, (bits[i] == 64) ? UWORD_MAX : (UWORD(1) << bits[i]) - 1);
        if (fft_small_plan_init_nmod(P, R, 0, 2 * len - 1, 2 * len - 1,
                n_round_up(len, BLK_SZ), len, 2 * NMOD_BITS(mod), mod, 0))
            flint_printf("%2d-bit modulus: %wu prime%s, length 2^%wu\n",
                         bits[i], P->np, P->np == 1 ? "" : "s", P->depth);
        fft_small_plan_clear(P);
    }

    flint_cleanup_master();
    return 0;
}

    1000 limbs: 2 primes,  44-bit chunks, length 2^12
   40000 limbs: 2 primes,  41-bit chunks, length 2^17
  100000 limbs: 5 primes, 116-bit chunks, length 2^17
 1000000 limbs: 6 primes, 136-bit chunks, length 2^20
10000000 limbs: 7 primes, 160-bit chunks, length 2^23
16-bit modulus: 1 prime, length 2^18
32-bit modulus: 2 primes, length 2^18
50-bit modulus: 3 primes, length 2^18
64-bit modulus: 3 primes, length 2^18

Memory and threads#

Each thread that multiplies has its own default context, created on first use: about 170 KB of tables, plus twiddle factors for transforms longer than 2^{12}, which take half a double per point and prime (4 MB per prime at length 2^{20}), plus a scratch buffer that grows to the largest product computed so far and is kept for reuse. flint_cleanup frees the default context of the calling thread, and flint_cleanup_master also stops the thread pool; long-running programs that once multiplied very large numbers can call them to return the memory.

The multiplication functions use the threads of FLINT's pool, as set by flint_set_num_threads. Integer products above 40000 limbs use up to eight threads, polynomial products use one thread per prime from lengths of about 1500 and up to eight threads for the reconstruction. The two-prime integer pipeline, and therefore every product below 40000 limbs, is single-threaded. Threads are only borrowed for the duration of a call, so the functions can be called from several threads at once, each using its own default context.

Calling the module directly#

Calling fft_small directly is useful to

• multiply without the size dispatch, for example to measure;
• keep the transform of an operand that multiplies many others (Products and divisions with a fixed operand for polynomials, Integer products with a common operand for integers);
• compute a middle product as a transposed product (Transposed products);
• add several products before transforming back, as in a dot product of polynomials or integers or a matrix product, which saves one inverse transform and one reconstruction per product (Arithmetic in transform space);
• compute products modulo 2^N + 1 (Multiplication modulo 2^N + 1) or with the transforms themselves (Transforms modulo one prime).

Code that must also build without fft_small should test FLINT_HAVE_FFT_SMALL and fall back to flint_mpn_mul, _nmod_poly_mul or _fmpz_poly_mul.