Convolutions of large integers

fft_convolution multiplies two polynomials whose coefficients are residues modulo $p = 2^N + 1$, $N = \mathit{limbs} \cdot \mathtt{FLINT\_BITS}$. When the coefficients of the true integer product are known to be smaller than $p$ (or than $p/2$ in absolute value), it computes the product of two polynomials with integer coefficients exactly. This is how fmpz_poly implements Schönhage–Strassen multiplication, and with integer coefficients obtained by cutting integers into chunks, it is also the core of mul_mfa_truncate_sqrt2. The guide Multiplying integers and polynomials explains the whole procedure.

Parameters and layout#

The transform has length $4n$, $n = 2^{\mathit{depth}}$, and only trunc of its points are computed. To multiply polynomials of lengths $\ell_1$ and $\ell_2$, take $\mathit{trunc} = \ell_1 + \ell_2 - 1$ and $\mathit{depth} = \lceil \log_2 \mathit{trunc} \rceil - 2$, so that $2n < \mathit{trunc} \le 4n$ as required (for $\mathit{trunc} \ge 3$). The ring size $N$ must be a multiple of $n$, so that the transform has $w = N/n \ge 1$ bits per point, and if $\mathit{limbs} > 128$ it must be a value returned by fft_adjust_limbs (or fft_adjust_limbs_sd_fft) because the pointwise products use fft_mulmod_2expp1.

• ii and jj are arrays of $4n$ pointers to blocks of $\mathit{limbs} + 1$ limbs. The first $\ell_1$, respectively $\ell_2$, entries hold the input coefficients as $(\mathit{limbs} + 1)$-limb two's complement integers (negative values are allowed); all other entries must be zero. Only the entries below trunc, rounded up to an even number or, for $\mathit{depth} > 6$, to a multiple of $2^{\lfloor \mathit{depth}/2 \rfloor + 1}$, are read, but zeroing all of them is simplest.
• On return, ii[k] holds the coefficient $c_k$ of the product, reduced to $[0, p)$, for $k < \mathit{trunc}$. If the true coefficients satisfy $|c_k| < p/2$, then $c_k$ is the residue if it is less than $p/2$ and the residue minus $p$ otherwise.
• jj is overwritten by its transform. Passing the same array twice (ii == jj) computes a square with one forward transform less.
• t1, t2 and s1 are arrays of at least flint_get_num_threads() pointers to $\mathit{limbs} + 1$ limbs each, and tt an array of as many pointers to $2(\mathit{limbs} + 1)$ limbs each. For $\mathit{depth} \le 6$ only the first entries are used.

The bound on $N$ follows from the sizes of the inputs: if $|a_i| < 2^{b_1}$ and $|b_j| < 2^{b_2}$, then $|c_k| < \min(\ell_1, \ell_2) \, 2^{b_1 + b_2}$, so

suffices, with $s = 1$ if any coefficient is negative and $s = 0$ otherwise. _fmpz_poly_mulmid_SS uses exactly this bound, rounds it up to a multiple of $n$ and of $\mathtt{FLINT\_BITS}$, and adjusts the limb count with fft_adjust_limbs_sd_fft.

If the product has more than trunc coefficients, the result is meaningless, except for $\mathit{trunc} = 4n$, where fft_convolution computes the cyclic product modulo $x^{4n} - 1$ (each coefficient $c_k$ with $k \ge 4n$ is added to $c_{k - 4n}$). fft_convolution_choose_mid exploits this.

Algorithm#

For $\mathit{depth} \le 6$, fft_convolution uses fft_truncate_sqrt2 on both operands, multiplies the transformed coefficients with fft_mulmod_2expp1, applies ifft_truncate_sqrt2 and divides by $4n$. For larger depths it uses the matrix Fourier algorithm with $n_1 = 2^{\lfloor \mathit{depth}/2 \rfloor}$ columns, in three threaded stages: fft_mfa_truncate_sqrt2_outer (sqrt2 layer and column transforms) on each operand, then fft_mfa_truncate_sqrt2_inner, which handles one row at a time (row transforms of both operands, pointwise products, inverse row transform) while the row is in cache, and finally ifft_mfa_truncate_sqrt2_outer, which also divides by $4n$. The total cost is that of three truncated transforms of trunc points plus trunc multiplications modulo $2^N + 1$, each a Schönhage–Strassen product of its own for large $N$.

fft_convolution_sd_fft computes the pointwise products with sd_fft_mpn_mulmod_2expp1 from fft_small instead, a negacyclic convolution modulo several word-size primes (Multiplication modulo 2^N + 1), which is about twice as fast on a recent x86-64 processor for rings of $16384$ bits and more. Its context must be initialised with sd_fft_mpn_mulmod_2expp1_ctx_init for $N = \mathit{limbs} \cdot \mathtt{FLINT\_BITS}$ and the digit count $m$ that fft_adjust_limbs_sd_fft returned together with limbs. For large depths it keeps the fused inner stage only when there are at least eight work units (columns plus truncated rows) per thread; otherwise it performs the transforms in separate single-threaded passes and only the pointwise products, which dominate, in threads.

Partial products#

When only the coefficients with indices in $[\mathit{nlo}, \mathit{nhi})$ of a product of length $\ell = \ell_1 + \ell_2 - 1$ are wanted, a cyclic transform can be shorter than the truncated one. In the product modulo $x^M - 1$, coefficient $k$ receives $c_k + c_{k+M}$. If $M \ge \mathit{nhi}$ and $M \ge \ell - \mathit{nlo}$, every wrapped term $c_{k+M}$ with $k \ge \mathit{nlo}$ has $k + M \ge \ell$ and is zero, so the window is unaffected. fft_convolution_choose_mid takes the least power of two $M = 4 \cdot 2^{\mathit{depth}}$ with this property and returns it as trunc if it is smaller than the number of points the truncated transform of length $\ell$ would process; otherwise it returns $\mathit{depth} = \lceil \log_2 \ell \rceil - 2$ (at least 0) with $\mathit{trunc} = \ell$. For the middle product of $2m$ by $m$ coefficients, where the wanted window has length $m$, this gives $M = 2m$ points instead of about $3m$ when $2m$ is a power of two.

Examples#

The product of two polynomials of length 300 with signed 31-bit coefficients. The coefficients of the product are below $300 \cdot 2^{60} < 2^{69}$ in absolute value, so $N = 128$ would do, but $N$ must be a multiple of $n = 256$, so the example uses four limbs.

#include <gmp.h>
#include <flint/flint.h>
#include <flint/fft.h>
#include <flint/fmpz.h>
#include <flint/fmpz_poly.h>
#include <flint/mpn_extras.h>

/* one coefficient c with |c| < 2^63, as a (limbs + 1)-limb two's complement value */
static void set_si(mp_limb_t * x, slong c, slong limbs)
{
    flint_mpn_store(x, limbs + 1, c < 0 ? ~UWORD(0) : 0);
    x[0] = (mp_limb_t) c;
}

/* a reduced residue modulo 2^N + 1, read in the symmetric range */
static void get_fmpz(fmpz_t c, mp_limb_t * x, slong limbs, const mpz_t p)
{
    mpz_t t;
    mpz_init(t);
    fermat_to_mpz(t, x, limbs);
    if (mpz_sizeinbase(t, 2) >= limbs * FLINT_BITS)   /* t >= 2^(N-1) */
        mpz_sub(t, t, p);
    fmpz_set_mpz(c, t);
    mpz_clear(t);
}

int main(void)
{
    slong len = 300, len_out = 2 * len - 1;
    slong depth = FLINT_CLOG2(len_out) - 2;     /* 4n >= len_out > 2n */
    slong n = WORD(1) << depth;                  /* 256 */
    slong limbs = 4, size = limbs + 1;           /* N = 256 >= 70 bits, n | N */
    slong nt = flint_get_num_threads(), i;
    mp_limb_t ** ii, ** jj, ** t1, ** t2, ** s1, ** tt, * mem;
    fmpz_poly_t f, g, h, h2;
    mpz_t p;

    fmpz_poly_init(f); fmpz_poly_init(g); fmpz_poly_init(h); fmpz_poly_init(h2);
    for (i = 0; i < len; i++)
    {
        fmpz_poly_set_coeff_si(f, i, (i * 1103515245 + 12345) % (WORD(1) << 30) - (WORD(1) << 29) * (i % 3));
        fmpz_poly_set_coeff_si(g, i, (i * 69069 + 1) % (WORD(1) << 30) - (WORD(1) << 30) * (i % 2));
    }

    ii = flint_malloc((8 * n + 4 * nt) * sizeof(mp_limb_t *));
    jj = ii + 4 * n;
    t1 = jj + 4 * n; t2 = t1 + nt; s1 = t2 + nt; tt = s1 + nt;
    mem = flint_malloc((8 * n + 5 * nt) * size * sizeof(mp_limb_t));
    for (i = 0; i < 8 * n; i++)
        ii[i] = mem + i * size;
    for (i = 0; i < nt; i++)
    {
        t1[i] = mem + (8 * n + 5 * i) * size;
        t2[i] = t1[i] + size;
        s1[i] = t2[i] + size;
        tt[i] = s1[i] + size;               /* 2 (limbs + 1) limbs */
    }

    for (i = 0; i < 4 * n; i++)
    {
        set_si(ii[i], i < len ? fmpz_get_si(f->coeffs + i) : 0, limbs);
        set_si(jj[i], i < len ? fmpz_get_si(g->coeffs + i) : 0, limbs);
    }

    fft_convolution(ii, jj, depth, limbs, len_out, t1, t2, s1, tt);

    mpz_init(p);
    mpz_ui_pow_ui(p, 2, limbs * FLINT_BITS);
    mpz_add_ui(p, p, 1);
    fmpz_poly_fit_length(h, len_out);
    for (i = 0; i < len_out; i++)
        get_fmpz(h->coeffs + i, ii[i], limbs, p);
    _fmpz_poly_set_length(h, len_out);
    _fmpz_poly_normalise(h);

    fmpz_poly_mul(h2, f, g);
    flint_printf("depth = %wd, trunc = %wd, equal: %d\n", depth, len_out, fmpz_poly_equal(h, h2));
    flint_printf("coefficient 299: ");
    fmpz_print(h->coeffs + 299);
    flint_printf("\n");

    mpz_clear(p);
    flint_free(mem); flint_free(ii);
    fmpz_poly_clear(f); fmpz_poly_clear(g); fmpz_poly_clear(h); fmpz_poly_clear(h2);
    return 0;
}

depth = 8, trunc = 599, equal: 1
coefficient 299: -889640045340373524

A product with 8000-bit coefficients, where the ring has 16384 bits and the pointwise products can use the negacyclic engine of fft_small. Here fft_adjust_limbs_sd_fft rounds 251 limbs up to 256 and returns $m = 256$ digits of 64 bits. In builds without fft_small it returns $m = 0$ and the program uses fft_convolution.

#include <gmp.h>
#include <flint/flint.h>
#include <flint/fft.h>
#include <flint/fmpz.h>
#include <flint/fmpz_poly.h>
#include <flint/mpn_extras.h>
#if FLINT_HAVE_FFT_SMALL
#include <flint/fft_small.h>
#endif

int main(void)
{
    slong len = 40, len_out = 2 * len - 1, i, m;
    slong depth = FLINT_CLOG2(len_out) - 2, n = WORD(1) << depth;   /* n = 32 */
    slong bits = 8000, res_bits, limbs, size;
    slong nt = flint_get_num_threads();
    mp_limb_t ** ii, ** jj, ** t1, ** t2, ** s1, ** tt, * mem;
    fmpz_poly_t f, g, h, h2;
    flint_rand_t state;

    flint_rand_init(state);
    fmpz_poly_init(f); fmpz_poly_init(g); fmpz_poly_init(h); fmpz_poly_init(h2);
    fmpz_poly_randtest_unsigned(f, state, len, bits);
    fmpz_poly_randtest_unsigned(g, state, len, bits);
    fmpz_poly_set_coeff_ui(f, len - 1, 1);
    fmpz_poly_set_coeff_ui(g, len - 1, 1);

    /* product coefficients < len 2^(2 bits): round up to a multiple of n */
    res_bits = 2 * bits + FLINT_CLOG2(len);
    res_bits = ((res_bits + n - 1) / n) * n;
    limbs = fft_adjust_limbs_sd_fft((res_bits - 1) / FLINT_BITS + 1, n, &m);
    size = limbs + 1;

    ii = flint_malloc((8 * n + 4 * nt) * sizeof(mp_limb_t *));
    jj = ii + 4 * n;
    t1 = jj + 4 * n; t2 = t1 + nt; s1 = t2 + nt; tt = s1 + nt;
    mem = flint_calloc((8 * n + 5 * nt) * size, sizeof(mp_limb_t));
    for (i = 0; i < 8 * n; i++)
        ii[i] = mem + i * size;
    for (i = 0; i < nt; i++)
    {
        t1[i] = mem + (8 * n + 5 * i) * size;
        t2[i] = t1[i] + size; s1[i] = t2[i] + size; tt[i] = s1[i] + size;
    }

    for (i = 0; i < len; i++)
    {
        fmpz_get_ui_array(ii[i], limbs, f->coeffs + i);
        fmpz_get_ui_array(jj[i], limbs, g->coeffs + i);
    }

#if FLINT_HAVE_FFT_SMALL
    if (m != 0)
    {
        sd_fft_mpn_mulmod_2expp1_ctx_t ctx;
        sd_fft_mpn_mulmod_2expp1_ctx_init(ctx, get_default_mpn_ctx(),
                                          limbs * FLINT_BITS, m);
        fft_convolution_sd_fft(ctx, ii, jj, depth, limbs, len_out, t1, t2, s1, tt);
        sd_fft_mpn_mulmod_2expp1_ctx_clear(ctx);
    }
    else
#endif
        fft_convolution(ii, jj, depth, limbs, len_out, t1, t2, s1, tt);

    fmpz_poly_fit_length(h, len_out);
    for (i = 0; i < len_out; i++)
        fmpz_set_ui_array(h->coeffs + i, ii[i], limbs);
    _fmpz_poly_set_length(h, len_out);
    fmpz_poly_mul(h2, f, g);
    flint_printf("equal: %d\n", fmpz_poly_equal(h, h2));

    flint_free(mem); flint_free(ii);
    fmpz_poly_clear(f); fmpz_poly_clear(g); fmpz_poly_clear(h); fmpz_poly_clear(h2);
    flint_rand_clear(state);
    return 0;
}

equal: 1