Products and divisions with a fixed operand

A product costs, per prime, two forward transforms, a pointwise product and one inverse transform. When one operand b is multiplied by many polynomials, its transforms can be computed once and kept, which saves about a third of the work of every later product. For middle products in which the fixed polynomial is the long operand, a transposed version keeps its images through the transposed inverse transform instead. nmod_poly uses both for evaluation, interpolation and extrapolation at geometric progressions, where the same polynomial multiplies several inputs. Division with remainder by a fixed divisor builds on it: the divisor and its power series inverse are both kept in transformed form.

Products#

1. Call _mul_precomp_init with b and a depth d \ge 4 such that 2^d is at least the lengths of the products (or windows) needed later. Pass btrunc = 2^d unless all later outputs are shorter, and len_bound = 2^d unless a smaller bound on the number of products per coefficient is known, such as \min(a_n, b_n). Pass a nonzero direct_len (at least 100) to allow a transform modulo n itself when n is a suitable prime.
2. Call _nmod_poly_mul_mid_precomp for each a. It returns 0, writing nothing, if the window does not fit in a cyclic convolution of length 2^d.
3. Free the data with _mul_precomp_clear.

A window [z_l, z_h) of the product of lengths a_n and b_n fits if a_n + b_n - 1 \le 2^d, or if a_n \le 2^d and z_h \le 2^d \le a_n + b_n - 1 \le z_l + 2^d: the coefficients of index i \ge 2^d then fold onto i - 2^d < z_l, outside the window. This is the wrap-around trick for middle products [BZ2011]: a middle product of lengths 2m and m fits in any transform length 2^d \ge 2m, where the full product would need 2^d \ge 3m - 1.

_nmod_poly_mul_mid_precomp only reads M: it works on a copy of the plan whose primes, Chinese remaindering data and scratch memory come from the context R. Several threads can use one M at the same time, each with its own context, normally get_default_mpn_ctx() of the thread. All default contexts have the same primes, which is required.

Transposed products#

The transposed product z_j = \sum_i r_i v_{i+j}, 0 \le j < z_n, is computed by _nmod_poly_mul_tr_mpn_ctx as F^T(\hat r \cdot I^T(v)), with \hat r the forward transform of r and F^T, I^T the transposes of the forward and inverse transforms (Transposed FFT). For a fixed v, the images I^T(v) can be stored:

1. Call _mul_precomp_init_tr with (v, v_n), a bound len_bound on the number of nonzero coefficients of the later r (at most v_n), and direct_len as above. The plan is that of the window [0, v_n) of a product of length v_n, whose truncated transform length is v_n rounded up to a multiple of 256 (for v_n > 128).
2. Call _nmod_poly_mul_tr_precomp for each r with r_n + z_n - 1 \le v_n (it returns 0 otherwise). Each call converts and transforms r, multiplies pointwise and applies F^T: two transforms instead of three.
3. Free the data with _mul_precomp_clear.

With r the reversal of a, preceded by k zeros, the output is the window [a_n + k - 1, a_n + k - 1 + z_n) of av. Storing the forward transform of v with _mul_precomp_init instead would need a cyclic length 2^d into which that window fits; the transposed version needs transforms of length about v_n only. As for _nmod_poly_mul_mid_precomp, M is only read.

A polynomial v of length 3000 prepared once, and the window [999, 2999) of its products with three polynomials of length 1000:

#include <flint/flint.h>
#include <flint/nmod.h>
#include <flint/nmod_vec.h>
#include <flint/nmod_poly.h>
#include <flint/fft_small.h>

int main(void)
{
    slong vn = 3000, an = 1000, zn = 2000, i, k;
    nmod_t mod;
    nn_ptr v, a, r, z, ref;
    mul_precomp_struct M[1];
    mpn_ctx_struct * R = get_default_mpn_ctx();

    nmod_init(&mod, 1000003);
    v = _nmod_vec_init(vn);
    a = _nmod_vec_init(an);
    r = _nmod_vec_init(an);
    z = _nmod_vec_init(zn + 2);
    ref = _nmod_vec_init(zn);
    for (i = 0; i < vn; i++)
        v[i] = nmod_pow_ui(3, (ulong) i * i, mod);

    /* store the images of v through the transposed inverse transform,
       for short operands of length at most an */
    if (!_mul_precomp_init_tr(M, v, vn, an, 0, mod, R))
        flint_abort();

    for (k = 1; k <= 3; k++)
    {
        for (i = 0; i < an; i++)
            a[i] = nmod_pow_ui(i + k, 3, mod);

        /* r = rev(a): the output is the window [an - 1, an - 1 + zn) of a*v */
        _nmod_poly_reverse(r, a, an, an);
        _nmod_poly_mulmid(ref, v, vn, a, an, an - 1, an - 1 + zn, mod);

        flint_printf("k = %wd: %s\n", k,
            (_nmod_poly_mul_tr_precomp(z, zn, r, an, M, mod, R)
              && _nmod_vec_equal(z, ref, zn)) ? "agrees" : "WRONG");
    }

    /* too long: an + zn - 1 > vn */
    flint_printf("window of length %wd returns %d\n", zn + 2,
        _nmod_poly_mul_tr_precomp(z, zn + 2, r, an, M, mod, R));

    _mul_precomp_clear(M);
    _nmod_vec_clear(v);
    _nmod_vec_clear(a);
    _nmod_vec_clear(r);
    _nmod_vec_clear(z);
    _nmod_vec_clear(ref);
    flint_cleanup_master();
    return 0;
}

k = 1: agrees
k = 2: agrees
k = 3: agrees
window of length 2002 returns 0

Division with remainder#

For a divisor b of length b_n \ge 2 with invertible leading coefficient, _nmod_poly_divrem_precomp_init computes the power series inverse B of the reversal of b to precision B_n and stores the transforms of B and of b. A dividend a of length a_n \le b_n + B_n - 1 is then divided as follows, with q_n = a_n - b_n + 1 \le B_n:

• the quotient q is the window [B_n - 1, B_n - 1 + q_n) of the product of the top q_n coefficients of a with B, which is the reversal of \operatorname{rev}(a) \cdot \operatorname{rev}(b)^{-1} \bmod x^{q_n} [GG2003, §9.1];
• the remainder is r = a - bq, computed modulo x^N - 1 for the power of two N \ge \max(b_n - 1, 256). Since \deg r < b_n - 1 \le N, the reduction modulo x^N - 1 loses nothing, and only the residue of a modulo x^N - 1 and the cyclic product bq are needed.

_nmod_poly_divrem_mpn_ctx performs the same steps without precomputation.

_nmod_poly_divrem_precomp returns 0 without computing anything when q_n > B_n, or when the quotient window does not fit in the prepared transform.

Examples#

A polynomial of length 1000 transformed once for length 2^{12} and multiplied by three others, with full products of length 3999 and with a middle product that uses the wrap-around:

#include <flint/flint.h>
#include <flint/nmod.h>
#include <flint/nmod_vec.h>
#include <flint/nmod_poly.h>
#include <flint/fft_small.h>

int main(void)
{
    slong an = 4000, bn = 1000, i, k;
    nmod_t mod;
    nn_ptr a, b, full, z;
    mul_precomp_struct M[1];
    mpn_ctx_struct * R = get_default_mpn_ctx();

    nmod_init(&mod, 1000003);
    a = _nmod_vec_init(an);
    b = _nmod_vec_init(bn);
    full = _nmod_vec_init(an + bn - 1);
    z = _nmod_vec_init(an + bn - 1);
    for (i = 0; i < bn; i++)
        b[i] = nmod_pow_ui(2, i, mod);

    /* transform b once, for cyclic convolutions of length 2^12 */
    if (!_mul_precomp_init(M, b, bn, 4096, 12, 4096, 0, mod, R))
        flint_abort();

    for (k = 1; k <= 3; k++)
    {
        for (i = 0; i < an; i++)
            a[i] = nmod_pow_ui(i + k, 3, mod);

        /* a product of lengths 3000 and 1000 fits into 4096 coefficients */
        _nmod_poly_mul(full, a, 3000, b, bn, mod);
        flint_printf("k = %wd: full product %s", k,
            (_nmod_poly_mul_mid_precomp(z, 0, 3999, a, 3000, M, mod, R)
              && _nmod_vec_equal(z, full, 3999)) ? "agrees" : "WRONG");

        /* lengths 4000 and 1000: only windows that wrap below 999 fit */
        _nmod_poly_mul(full, a, an, b, bn, mod);
        flint_printf(", full product of length 4999 returns %d",
            _nmod_poly_mul_mid_precomp(z, 0, 4999, a, an, M, mod, R));
        flint_printf(", middle product %s\n",
            (_nmod_poly_mul_mid_precomp(z, 999, 4000, a, an, M, mod, R)
              && _nmod_vec_equal(z, full + 999, 3001)) ? "agrees" : "WRONG");
    }

    _mul_precomp_clear(M);
    _nmod_vec_clear(a);
    _nmod_vec_clear(b);
    _nmod_vec_clear(full);
    _nmod_vec_clear(z);
    flint_cleanup_master();
    return 0;
}

k = 1: full product agrees, full product of length 4999 returns 0, middle product agrees
k = 2: full product agrees, full product of length 4999 returns 0, middle product agrees
k = 3: full product agrees, full product of length 4999 returns 0, middle product agrees

Division by a fixed divisor of length 500, prepared for quotients of length up to 1500:

#include <flint/flint.h>
#include <flint/nmod.h>
#include <flint/nmod_vec.h>
#include <flint/nmod_poly.h>
#include <flint/fft_small.h>

int main(void)
{
    ulong n = 998244353;
    slong bn = 500, Bn = 1500, an, i;
    nmod_t mod;
    nmod_poly_t A, B, Q, R;
    nn_ptr q, r;
    nmod_poly_divrem_precomp_struct M[1];
    mpn_ctx_struct * ctx = get_default_mpn_ctx();

    nmod_init(&mod, n);
    nmod_poly_init(A, n);
    nmod_poly_init(B, n);
    nmod_poly_init(Q, n);
    nmod_poly_init(R, n);
    for (i = 0; i < bn; i++)
        nmod_poly_set_coeff_ui(B, i, nmod_pow_ui(7, i + 1, mod));

    /* the divisor serves dividends of length up to bn + Bn - 1 = 1999 */
    _nmod_poly_divrem_precomp_init(M, B->coeffs, bn, Bn, mod, ctx);
    q = _nmod_vec_init(Bn);
    r = _nmod_vec_init(bn - 1);

    for (an = 1000; an <= 1999; an += 333)
    {
        nmod_poly_zero(A);
        for (i = 0; i < an; i++)
            nmod_poly_set_coeff_ui(A, i, nmod_pow_ui(i + 2, 5, mod));
        nmod_poly_divrem(Q, R, A, B);

        if (!_nmod_poly_divrem_precomp(q, r, A->coeffs, an, M, mod, ctx))
            flint_printf("length %wd: transform too short\n", an);
        else
            flint_printf("length %wd: quotient %s, remainder %s\n", an,
                _nmod_vec_equal(q, Q->coeffs, an - bn + 1) ? "agrees" : "WRONG",
                (R->length <= bn - 1
                  && _nmod_vec_equal(r, R->coeffs, R->length)
                  && _nmod_vec_is_zero(r + R->length, bn - 1 - R->length))
                    ? "agrees" : "WRONG");
    }

    _nmod_poly_divrem_precomp_clear(M);
    _nmod_vec_clear(q);
    _nmod_vec_clear(r);
    nmod_poly_clear(A);
    nmod_poly_clear(B);
    nmod_poly_clear(Q);
    nmod_poly_clear(R);
    flint_cleanup_master();
    return 0;
}

length 1000: quotient agrees, remainder agrees
length 1333: quotient agrees, remainder agrees
length 1666: quotient agrees, remainder agrees
length 1999: quotient agrees, remainder agrees