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 single-word exponents and 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).