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).