Bivariate resultants by evaluation and interpolation

Let A, B \in \mathbb{F}_p[x][y] have degrees n_A = \deg_y A \ge 1 and n_B = \deg_y B \ge 1 in y and degrees d_A, d_B in x. Their resultant with respect to y is the determinant of the Sylvester matrix, a polynomial R = \mathrm{Res}_y(A, B) \in \mathbb{F}_p[x]. The subresultant algorithm computes it by pseudo-remainder sequences in \mathbb{F}_p[x][y], whose intermediate coefficients have growing degree in x. When \mathbb{F}_p is large enough, it is faster to evaluate x at many points, compute univariate resultants over \mathbb{F}_p, and interpolate, as in the classical modular resultant algorithm [GG2003, §6.9]. FLINT implements this method for nmod_mpoly with fast evaluation and interpolation at structured points: a geometric progression [BS2005], or roots of unity when p is an FFT prime.

Degree bound and specialisation#

Hence R is determined by its values at N = n_B d_A + n_A d_B + 1 distinct points. The value at x = a must be computed from the specialisations f = A(a, y) and g = B(a, y), whose degrees may be lower than n_A and n_B when a is a root of a leading coefficient.

Algorithm#

The evaluation of one coefficient of length at most \max(d_A, d_B) + 1 at N points of a geometric progression costs O(M(N)) operations, so step 3 costs O((n_A + n_B) M(N)), step 5 costs O(M(N)), and step 4 costs N univariate resultants of degree n_A, each O(M(n_A) \log n_A) with the half-gcd algorithm (Half-GCD) or O(n_A n_B) with Euclid's algorithm for small degrees. In the FFT case the points are split into 2^{\mathrm{depth} - m} cosets of the 2^m-th roots of unity, each evaluated by a substitution x \mapsto w^t x and one transform of length 2^m; the values come out in the order that the truncated inverse transform needs, so interpolation is a single inverse transform (or a geometric interpolation with r = \sqrt{w} when that is cheaper). Both loops over the points run in parallel on the threads set with flint_set_num_threads when the work is large enough.

Example#

Two random polynomials of degree 8 in each variable modulo the prime 10^6 + 3; nmod_mpoly_resultant takes the evaluation and interpolation path, and the result agrees with the subresultant algorithm applied to the univariate forms. The resultant reaches the degree bound of Lemma 1:

#include <flint/flint.h>
#include <flint/nmod_mpoly.h>

int main(void)
{
    nmod_mpoly_ctx_t ctx;
    nmod_mpoly_t a, b, r1, r2;
    nmod_mpoly_univar_t ua, ub;
    ulong bounds[2] = {9, 9};     /* degrees < 9 in x and in y */
    flint_rand_t state;

    flint_rand_init(state);
    nmod_mpoly_ctx_init(ctx, 2, ORD_LEX, UWORD(1000003));
    nmod_mpoly_init(a, ctx); nmod_mpoly_init(b, ctx);
    nmod_mpoly_init(r1, ctx); nmod_mpoly_init(r2, ctx);
    nmod_mpoly_univar_init(ua, ctx); nmod_mpoly_univar_init(ub, ctx);

    nmod_mpoly_randtest_bounds(a, state, 81, bounds, ctx);
    nmod_mpoly_randtest_bounds(b, state, 81, bounds, ctx);

    /* eliminate y (variable 1): evaluation/interpolation in x */
    nmod_mpoly_resultant(r1, a, b, 1, ctx);

    /* the same resultant by subresultants over F_p[x] */
    nmod_mpoly_to_univar(ua, a, 1, ctx);
    nmod_mpoly_to_univar(ub, b, 1, ctx);
    nmod_mpoly_univar_resultant(r2, ua, ub, ctx);

    flint_printf("deg_y a = %wd, deg_x a = %wd\n", nmod_mpoly_degree_si(a, 1, ctx),
                 nmod_mpoly_degree_si(a, 0, ctx));
    flint_printf("deg_x res = %wd (bound %wd), equal: %d\n",
                 nmod_mpoly_degree_si(r1, 0, ctx),
                 nmod_mpoly_degree_si(b, 1, ctx) * nmod_mpoly_degree_si(a, 0, ctx)
                 + nmod_mpoly_degree_si(a, 1, ctx) * nmod_mpoly_degree_si(b, 0, ctx),
                 nmod_mpoly_equal(r1, r2, ctx));

    nmod_mpoly_univar_clear(ua, ctx); nmod_mpoly_univar_clear(ub, ctx);
    nmod_mpoly_clear(a, ctx); nmod_mpoly_clear(b, ctx);
    nmod_mpoly_clear(r1, ctx); nmod_mpoly_clear(r2, ctx);
    nmod_mpoly_ctx_clear(ctx);
    flint_rand_clear(state);
    return 0;
}

deg_y a = 8, deg_x a = 8
deg_x res = 128 (bound 128), equal: 1

Use in FLINT#

• nmod_mpoly_resultant uses the algorithm for polynomials in exactly two variables with positive degrees in the eliminated variable, after converting to dense bivariate form (nmod_mpoly_get_bpoly), when the dense boxes are not too large (at most 2^{32} entries and points, and not more than 64 times the number of terms beyond 2^{26} entries), the inputs are not tiny, and n_bpoly_mod_resultant_cutoff predicts a gain: (n_A + 1)(n_B + 1) \ge 36, or \ge 16 with at least 4096 points. Otherwise, or if \mathbb{F}_p is too small (Algorithm 1 fails), it uses the subresultant algorithm of Ducos [Duc2000] on the univariate forms (nmod_mpoly_univar_resultant).
• n_bpoly_mod_resultant and _n_bpoly_mod_resultant implement Algorithm 1 on n_bpoly_t.
• fmpz_mpoly_resultant uses _n_bpoly_mod_resultant modulo many primes and combines the results by the Chinese remainder theorem (see fmpz_mpoly.h).
• gr_poly_resultant_multipoint, for polynomials over \mathbb{F}_p[x] in the generic rings framework, passes its arguments to _n_bpoly_mod_resultant (the memory layouts of gr_poly over nmod and of n_bpoly coincide).