Multiplication modulo 2^N + 1

The pointwise products of a Schönhage–Strassen convolution are products modulo 2^N + 1. fft_mulmod_2expp1 computes them, by a full product followed by a reduction for N up to 8192 bits and by a negacyclic convolution, the original idea of Schönhage and Strassen [SS1971], for larger N. Since the negacyclic convolution needs N to have a special form, the outer convolution must choose its coefficient ring with fft_adjust_limbs (or fft_adjust_limbs_sd_fft, which also considers the faster pointwise engine of fft_small).

Usage#

The operands are reduced residues of \mathit{limbs} + 1 limbs, where \mathit{limbs} = nw/\mathtt{FLINT\_BITS}: only the product nw matters, and fft_convolution passes the transform parameters n and w of the outer transform. The function accepts any N up to FFT_MULMOD_2EXPP1_CUTOFF limbs, which is 128 on all targets. Beyond that, N must be a value produced by fft_adjust_limbs; for other values the result is silently wrong. For example, fft_adjust_limbs maps 129 to 144, 200 to 208, 1000 to 1024 and 12345 to 12544, and products modulo 2^{64 \cdot 129} + 1 or 2^{64 \cdot 1000} + 1 come out wrong. Every power of two and every multiple of three by a power of two above 128 is returned unchanged.

tt must provide 2 \cdot \mathit{limbs} limbs of scratch space for the full product of the basecase; FLINT's callers pass 2(\mathit{limbs} + 1) limbs. The output may alias either input.

Algorithm#

If one of the inputs is 2^N \equiv -1, the result is the negation of the other. Up to the cutoff, flint_mpn_mulmod_2expp1_basecase computes the full 2N-bit product with flint_mpn_mul_n (or flint_mpn_sqr) and subtracts the high half from the low half. Beyond the cutoff, _fft_mulmod_2expp1 cuts both operands into 2n pieces of b = N/(2n) bits, with n = 2^{\mathit{depth}} and \mathit{depth} = \lfloor \lceil \log_2 N \rceil / 2 \rfloor - k, where k comes from the tuning table MULMOD_TAB (between 1 and 4). Since 2^{2nb} = 2^N \equiv -1, the product modulo 2^N + 1 is the negacyclic convolution of the pieces, evaluated at 2^b. The convolution is computed with fft_negacyclic and ifft_negacyclic in the ring \mathbb{Z}/(2^{2b} + 1)\mathbb{Z}, with flint_mpn_mulmod_2expp1_basecase for its pointwise products, and, to recover coefficients that exceed 2^{2b}, combined by the Chinese remainder theorem with the same convolution modulo B, computed by fft_naive_convolution_1 from the lowest limb of each piece. The proof that this gives the exact coefficients is in Multiplication modulo 2^N + 1. Some values:

(generic tuning; the inner products are handled by the basecase, so there is a single level of recursion).

fft_adjust_limbs rounds limbs up so that the depth d_1 chosen for limbs and for the next power of two divides it correctly: the limb count becomes a multiple of 2^{d_1 + 1} and the bit count a multiple of 2^{2 d_1}, where d_1 is the larger of the two depths.

The negacyclic engine of fft_small#

In builds with fft_small, sd_fft_mpn_mulmod_2expp1 computes the same products with a negacyclic convolution modulo three to eight word-size primes in floating point, for rings N' = mb with m \ge 16 a power of two and b one of 18 instantiated digit sizes between 64 and 192 bits. fft_adjust_limbs_sd_fft(tight, n_outer, m_out) looks for such a ring for a convolution with transform parameter n = n_outer and coefficients of at least tight limbs:

• m is the least power of two \ge 16 with mb \ge \mathtt{FLINT\_BITS} \cdot \mathit{tight};
• mb must be a multiple of \mathtt{FLINT\_BITS}, at least FLINT_FFT_SD_FFT_MIN_RING = 16384 bits, and a multiple of 2n, so that the outer transform has an even number w of bits per point;
• the number of primes needed for 2b + \log_2 m + 1 bits must be available, and the ring may exceed fft_adjust_limbs(tight) by at most 35, 25, 15, 12, 5 or 2 percent for 3, 4, …, 8 primes;
• among the admissible rings, the one with the least estimated cost \mathit{np} \cdot m \cdot (\log_2 m + 4 + b/16) is chosen.

It returns the limb count mb/\mathtt{FLINT\_BITS} and stores m in *m_out, or returns fft_adjust_limbs(tight) with *m_out = 0 if no ring qualifies or fft_small is not available. The caps and weights are tuning parameters. fmpz_poly and fmpzi select their rings this way and create a context with sd_fft_mpn_mulmod_2expp1_ctx_init when m \ne 0; see the second example in Examples.

Example#

A product modulo 2^N + 1 with N near 64000 bits, checked against GMP.

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

int main(void)
{
    /* multiply modulo 2^N + 1 with N a multiple of 64 near 64000 bits */
    slong limbs = fft_adjust_limbs(1000);
    mp_limb_t * a, * b, * r, * tt;
    mpz_t A, B, R, P;

    flint_printf("limbs: 1000 -> %wd\n", limbs);

    a = flint_calloc(limbs + 1, sizeof(mp_limb_t));
    b = flint_calloc(limbs + 1, sizeof(mp_limb_t));
    r = flint_calloc(limbs + 1, sizeof(mp_limb_t));
    tt = flint_malloc(2 * (limbs + 1) * sizeof(mp_limb_t));

    mpz_init(A); mpz_init(B); mpz_init(R); mpz_init(P);
    mpz_ui_pow_ui(P, 2, limbs * FLINT_BITS);
    mpz_add_ui(P, P, 1);
    mpz_ui_pow_ui(A, 3, 40000);
    mpz_ui_pow_ui(B, 5, 27000);
    mpz_mod(A, A, P);
    mpz_mod(B, B, P);

    /* reduced operands: limbs + 1 limbs, top limb 0 (or 1 for 2^N) */
    mpz_export(a, NULL, -1, sizeof(mp_limb_t), 0, 0, A);
    mpz_export(b, NULL, -1, sizeof(mp_limb_t), 0, 0, B);

    /* n * w = N; here n = limbs and w = FLINT_BITS */
    fft_mulmod_2expp1(r, a, b, limbs, FLINT_BITS, tt);

    fermat_to_mpz(R, r, limbs);
    mpz_mul(A, A, B);
    mpz_mod(A, A, P);
    flint_printf("equal to mpz: %d\n", mpz_cmp(A, R) == 0);

    flint_free(a); flint_free(b); flint_free(r); flint_free(tt);
    mpz_clear(A); mpz_clear(B); mpz_clear(R); mpz_clear(P);
    return 0;
}

limbs: 1000 -> 1024
equal to mpz: 1