Multiplication modulo 2^N + 1

Schönhage–Strassen multiplication reduces a large product to many products modulo $2^N + 1$ for a moderately large $N$ (Schönhage–Strassen multiplication). The fft module computes these pointwise products with the functions here when the ring has at least 16384 bits (FLINT_FFT_SD_FFT_MIN_RING), choosing $N$ with fft_adjust_limbs_sd_fft so that a good digit splitting exists. They compute a product modulo $2^N + 1$ as a negacyclic convolution of $b$-bit digits with the small-prime transforms.

Usage#

1. Create a context for $N$ with sd_fft_mpn_mulmod_2expp1_ctx_init, passing get_default_mpn_ctx() and $m = 0$ to let it choose the digit count. $N$ must be a positive multiple of FLINT_BITS; otherwise an exception is raised.
2. Create one scratch object per thread with sd_fft_mpn_mulmod_2expp1_scratch_init. The context is read-only during multiplications and can be shared.
3. Multiply with sd_fft_mpn_mulmod_2expp1. Operands and result have $n + 1$ limbs, $n = N / \mathtt{FLINT\_BITS}$; operands may have top limb 0 or 1 with arbitrary low limbs, and the result is fully reduced into $[0, 2^N]$. The result may alias an operand.
4. For many products with the same $y$, store its transform once with sd_fft_mpn_mulmod_2expp1_transform in an array of sd_fft_mpn_mulmod_2expp1_transformed_size(C) doubles, aligned to 32 bytes, and multiply with sd_fft_mpn_mulmod_2expp1_mul_cached.

Algorithm#

Write $x = \sum_{i<m} x_i 2^{bi}$ with digits $0 \le x_i < 2^b$, where $N = mb$ and a top limb 1 contributes $2^N \equiv -1$, which is subtracted from $x_0$. Since $2^{bm} \equiv -1$, the product is $\sum_k z_k 2^{bk}$ with the negacyclic convolution

where $m = 2^d$. Modulo each prime $p$, with $\theta$ a primitive $2m$-th root of unity ($\theta^m = -1$), the weighted sequences $\theta^i x_i$ and $\theta^j y_j$ have the cyclic convolution $\theta^k z_k$ [CF1994], which a transform of length $m$ computes. The residues are untwisted by $\theta^{-k}$, the factor $m^{-1}$ of the inverse transform and the Chinese remainder normalisation being folded into the same table. To use the unsigned reconstruction, $2^{2b + d}$ is added to every coefficient, which makes them lie in $[0, 2^{2b + d + 1})$, so the product of the primes must be at least $2^{2b + d + 1}$. After reconstructing $\sum_k (z_k + 2^{2b+d}) 2^{bk}$ the bias polynomial is subtracted, and the part of the signed result above $2^N$ is folded back by subtraction.

The number of primes is the least $k \ge 2$ with enough capacity for which digits of $b$ bits can be converted: digits of at most 50 bits are converted directly, larger ones need one of the vectorised conversions of the integer multiplication profiles. With $m = 0$ the context compares the admissible digit counts $m \ge 16$ (with $20 \le b \le 200$) by an estimate $w_k \, k \, m (d + 4 + b/16)$ with measured weights $w_k$ per prime count, which favours two primes slightly and penalises eight.

Example#

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

int main(void)
{
    slong N = 65536, n = N / FLINT_BITS, i;
    sd_fft_mpn_mulmod_2expp1_ctx_t C;
    sd_fft_mpn_mulmod_2expp1_scratch_t S;
    nn_ptr x = flint_calloc(n + 1, sizeof(ulong));
    nn_ptr y = flint_calloc(n + 1, sizeof(ulong));
    nn_ptr z = flint_calloc(n + 1, sizeof(ulong));
    fmpz_t X, Y, Z, M;

    sd_fft_mpn_mulmod_2expp1_ctx_init(C, get_default_mpn_ctx(), N, 0);
    sd_fft_mpn_mulmod_2expp1_scratch_init(S, C);
    flint_printf("N = %wd: %wd digits of %wd bits, %wd primes\n",
                 C->N, C->m, C->b, C->np);

    for (i = 0; i < n; i++)
    {
        x[i] = UWORD(0x0123456789abcdef) * (i + 1);
        y[i] = UWORD(0xfedcba9876543210) ^ i;
    }
    y[n] = 0;
    x[n] = 1;       /* x = 2^N + (low limbs), not reduced */

    sd_fft_mpn_mulmod_2expp1(C, z, x, y, S);

    fmpz_init(X); fmpz_init(Y); fmpz_init(Z); fmpz_init(M);
    fmpz_set_ui_array(X, x, n + 1);
    fmpz_set_ui_array(Y, y, n + 1);
    fmpz_one(M);
    fmpz_mul_2exp(M, M, N);
    fmpz_add_ui(M, M, 1);
    fmpz_mul(Z, X, Y);
    fmpz_mod(Z, Z, M);
    fmpz_set_ui_array(X, z, n + 1);
    flint_printf("x*y mod 2^N + 1: %s\n", fmpz_equal(X, Z) ? "agrees" : "WRONG");

    fmpz_clear(X); fmpz_clear(Y); fmpz_clear(Z); fmpz_clear(M);
    sd_fft_mpn_mulmod_2expp1_scratch_clear(S);
    sd_fft_mpn_mulmod_2expp1_ctx_clear(C);
    flint_free(x);
    flint_free(y);
    flint_free(z);
    flint_cleanup_master();
    return 0;
}

N = 65536: 1024 digits of 64 bits, 3 primes
x*y mod 2^N + 1: agrees