Polynomial multiplication over Z/nZ for multiword n

These functions multiply polynomials over an mpn_mod ring. The dispatchers _mpn_mod_poly_mullow and _mpn_mod_poly_mulmid are the ring's _gr_poly_mullow and _gr_poly_mulmid methods, so every generic polynomial algorithm over the ring (gr_poly_mul, Newton iteration, half-gcd, composition, ...) multiplies through them. The other functions are the individual algorithms, useful for testing and profiling.

All of them multiply the residues as integers in $[0, n)$ and reduce each output coefficient once, at the end. The algorithms differ in how they compute the integer coefficients of the product: schoolbook multiplication and Karatsuba on fixed-size limb arrays, Kronecker substitution into one large integer product, or number-theoretic transforms modulo small primes.

Windows#

A middle product with window $[\mathit{nlo}, \mathit{nhi})$ computes only the coefficients $c_i$, $\mathit{nlo} \le i < \mathit{nhi}$, of $c = ab$, and stores $c_{\mathit{nlo}}$ at res[0]. Coefficients of $a$ and $b$ that cannot contribute to the window are skipped first: terms of index at least $\mathit{nhi}$, and, for $\mathit{nlo} > 0$, the lowest coefficients of an operand that is longer than the number $\mathit{len1} + \mathit{len2} - 1 - \mathit{nlo}$ of coefficients from $\mathit{nlo}$ to the end of the product (they only reach indices below $\mathit{nlo}$). After this, _mpn_mod_poly_mulmid chooses an algorithm for the remaining lengths.

Choice of algorithm#

Let $m = \min(\mathit{len1}, \mathit{len2})$ after the truncation.

1. If $m < 4$ or the window has at most 4 coefficients, the schoolbook algorithm (_mpn_mod_poly_mulmid_classical) is used.
2. A cutoff $c$ is read from one of five tables indexed by the bit size $b$ of $n$ in steps of 16 bits, for squaring, low squaring, unbalanced products ($\max(\mathit{len1}, \mathit{len2}) \ge 2m$), full products and other windows. The tables were generated by the profiling program p-poly_mullow; $c$ is between 25 and 173 and decreases as $b$ grows. For example, for full products of balanced operands $c$ is 91 at 80 bits, 95 at 256 bits, 54 at 512 bits and 25 at 1024 bits.
3. If $m < c$ and $\max(\mathit{len1}, \mathit{len2}) < 4m$, the product is computed without intermediate reductions by _mpn_mod_poly_mulmid_mpn, which leaves the choice between the schoolbook algorithm and Karatsuba to _flint_mpn_poly_mulmid of mpn_extras, with its own tuned cutoff.
4. Otherwise the small-prime FFT (_mpn_mod_poly_mulmid_fft_small) is tried, and if it declines (fft_small unavailable, or $n$ larger than about 190 bits), Kronecker substitution (_mpn_mod_poly_mulmid_KS) is used. The tables assume that the FFT always beats Kronecker substitution when it can be used, which is nearly but not always true.

Schoolbook and Karatsuba#

The schoolbook algorithm computes each $c_i = \sum_j a_j b_{i-j}$ as an exact sum of full products in $s$ limbs and reduces it with mpn_mod_set_mpn. A sum of at most $m$ products of residues below $n < 2^b$ is less than $m \, 2^{2b}$, so

limbs suffice, usually $2\ell$ or $2\ell + 1$. The sums are computed by the generated fixed-size kernels of mpn_extras (with the dot products inlined into the loop over the output coefficients) for $\ell \le 4$ and with flint_mpn_mul_n and mpn_add_n otherwise. The first and last coefficients, which are single products, are computed with mpn_mod_mul instead.

Karatsuba's algorithm [KO1962] splits each operand into a low and a high half and computes the three half-size products $a_0 b_0$, $a_1 b_1$ and $(a_0 + a_1)(b_0 + b_1)$ recursively. Over $\mathbb{Z}/n\mathbb{Z}$ with reduced coefficients the sums $a_0 + a_1$ would need reductions, and the middle term needs a subtraction; FLINT instead runs the whole recursion on the unreduced integer coefficients, which grow by about one bit per level for the sums and are allowed a sign, and reduces only the final coefficients. The output size bound of _flint_mpn_poly_mulmid_slimbs accounts for the growth: $2b + 2 \lfloor \log_2 m \rfloor + 3$ bits (with one extra limb when the operand coefficients have no spare high bit and the recursion must extend them). This is why Karatsuba pays off from short lengths: an addition of residues of $\ell$ limbs costs much less than a multiplication.

Kronecker substitution#

Kronecker substitution [Kro1882; Sch1982] evaluates both polynomials at $x = 2^k$ with $k = 2b + \lfloor \log_2 m \rfloor + 1$, which exceeds the bit size of every coefficient of the integer product, multiplies the two resulting integers with flint_mpn_mul (which uses GMP, or fft_small for large operands) and reads the coefficients off the product $k$ bits at a time, reducing each modulo $n$. Only the windowed coefficients are extracted. The packing wastes the padding bits but reduces the polynomial product to a single integer product, which is efficient for large $n$, where the small-prime FFT is not available.

Small-prime FFT#

_mpn_mod_poly_mulmid_fft_small reduces every coefficient modulo the 50-bit primes of fft_small, multiplies the polynomials modulo each prime with truncated number-theoretic transforms and reconstructs the integer coefficients by Chinese remaindering, reducing them modulo $n$ directly from the remainder representation. It needs a product of primes exceeding $\mathit{len2} \cdot 2^{2b + 1}$, and at most eight primes are available, which limits $b$ to about 190 (196 bits for very short products, 189 bits at length $10^5$); above that it returns GR_UNABLE. See Multiplication with several primes and CRT and FFTs modulo word-size primes in floating point for the transforms and the reconstruction.

Example#

The program below computes a full product and a middle product of two polynomials of length 3 modulo $2^{127} - 1$, the second time also with Kronecker substitution. Coefficients are stored in arrays of MPN_MOD_MAX_LIMBS limbs per slot to be independent of the word size; nl is the actual number of limbs.

#include <flint/fmpz.h>
#include <flint/gr.h>
#include <flint/mpn_mod.h>

#define MAXL MPN_MOD_MAX_LIMBS

static void print_poly(nn_srcptr a, slong len, gr_ctx_t ctx)
{
    slong i, nl = MPN_MOD_CTX_NLIMBS(ctx);
    fmpz_t c;
    fmpz_init(c);
    for (i = 0; i < len; i++)
    {
        GR_MUST_SUCCEED(mpn_mod_get_fmpz(c, a + i * nl, ctx));
        fmpz_print(c);
        printf(i + 1 < len ? " " : "\n");
    }
    fmpz_clear(c);
}

int main(void)
{
    gr_ctx_t ctx;
    fmpz_t p;
    ulong f[3 * MAXL], g[3 * MAXL], h[5 * MAXL];
    slong i, nl;
    int status = GR_SUCCESS;

    fmpz_init(p);
    fmpz_ui_pow_ui(p, 2, 127);
    fmpz_sub_ui(p, p, 1);
    GR_MUST_SUCCEED(gr_ctx_init_mpn_mod(ctx, p));
    nl = MPN_MOD_CTX_NLIMBS(ctx);

    /* f = 1 + 2x + 3x^2, g = 4 + 5x - 6x^2 */
    for (i = 0; i < 3; i++)
    {
        status |= mpn_mod_set_si(f + i * nl, i + 1, ctx);
        status |= mpn_mod_set_si(g + i * nl, i == 2 ? -6 : i + 4, ctx);
    }

    /* full product: coefficients 0 to 4 */
    status |= _mpn_mod_poly_mullow(h, f, 3, g, 3, 5, ctx);
    print_poly(h, 5, ctx);

    /* middle product: coefficients 1, 2 and 3 */
    status |= _mpn_mod_poly_mulmid(h, f, 3, g, 3, 1, 4, ctx);
    print_poly(h, 3, ctx);

    /* the same middle product by Kronecker substitution */
    status |= _mpn_mod_poly_mulmid_KS(h, f, 3, g, 3, 1, 4, ctx);
    print_poly(h, 3, ctx);

    gr_ctx_clear(ctx);
    fmpz_clear(p);
    printf("status %d\n", status);
    return 0;
}

4 13 16 3 170141183460469231731687303715884105709
13 16 3
13 16 3
status 0

The leading coefficient is $-18 \bmod (2^{127} - 1)$: residues are always printed in $[0, n)$.