Resultants and discriminants

The resultant \operatorname{res}(f, g) is the determinant of the Sylvester matrix of f and g; it vanishes exactly when f and g have a common factor, and for f = a \prod (x - \alpha_i) it equals a^{\deg g} \prod_i g(\alpha_i). The discriminant of f is, up to sign and a power of the leading coefficient, the resultant of f and f', and it vanishes exactly when f has a repeated root. Resultants over \mathbb{Z}/p\mathbb{Z} for many primes p are the usual way to compute resultants of integer polynomials by Chinese remaindering. These functions require a prime modulus.

Algorithms#

The resultant is computed along the Euclidean remainder sequence: for r = f \bmod g,

and \operatorname{res}(f, c) = c^{\deg f} for a constant c [GG2003, §6]. nmod_poly_resultant_euclidean follows this recursion in O(\deg f \deg g) operations. nmod_poly_resultant_hgcd performs the same divisions with the half-GCD algorithm in O(M(n) \log n) operations; because a recursive half-GCD call only knows the top coefficients of the remainders, the factors contributed by its last division are recorded and applied after the remainders have been completed. nmod_poly_resultant handles small degrees by explicit formulas, and otherwise chooses between the two at the same cutoff as nmod_poly_gcd, using Montgomery arithmetic for the Euclidean algorithm under the same conditions.

The discriminant is computed as

which is (-1)^{d(d-1)/2} \operatorname{res}(f, f') / \operatorname{lc}(f) when \deg f' = d - 1, that is, when the characteristic does not divide d.

Examples#

For the cubic x^3 + ax + b the discriminant is -4a^3 - 27b^2. Modulo p = 101, the program checks this formula, finds the values of b for which x^3 + 3x + b has a repeated root, and compares the two resultant algorithms on long random input:

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

int main(void)
{
    nmod_poly_t f, g, h;
    flint_rand_t state;
    nmod_t mod;
    ulong a = 3, b, d;

    nmod_init(&mod, 101);
    flint_rand_init(state);
    nmod_poly_init_mod(f, mod);
    nmod_poly_init_mod(g, mod);
    nmod_poly_init_mod(h, mod);

    nmod_poly_set_coeff_ui(f, 3, 1);
    nmod_poly_set_coeff_ui(f, 1, a);
    nmod_poly_set_coeff_ui(f, 0, 5);
    d = nmod_sub(nmod_neg(nmod_mul(4, nmod_pow_ui(a, 3, mod), mod), mod),
                 nmod_mul(27, 25, mod), mod);
    flint_printf("disc = %wu, -4a^3 - 27b^2 = %wu\n",
                 nmod_poly_discriminant(f), d);

    flint_printf("repeated root for b =");
    for (b = 0; b < mod.n; b++)
    {
        nmod_poly_set_coeff_ui(f, 0, b);
        if (nmod_poly_discriminant(f) == 0)
            flint_printf(" %wu", b);
    }
    flint_printf("\n");

    /* res(f, g) = 0 exactly when gcd(f, g) != 1 */
    nmod_poly_set_coeff_ui(f, 0, 5);
    nmod_poly_set_str(g, "2 101  1 1");               /* x + 1 */
    flint_printf("res(f, x + 1) = %wu = -f(-1) = %wu\n",
                 nmod_poly_resultant(f, g),
                 nmod_neg(nmod_poly_evaluate_nmod(f, mod.n - 1), mod));

    nmod_poly_rand(g, state, 3000);
    nmod_poly_rand(h, state, 2500);
    flint_printf("euclidean == hgcd: %d\n",
                 nmod_poly_resultant_euclidean(g, h) ==
                 nmod_poly_resultant_hgcd(g, h));

    nmod_poly_clear(f);
    nmod_poly_clear(g);
    nmod_poly_clear(h);
    flint_rand_clear(state);
    return 0;
}

disc = 25, -4a^3 - 27b^2 = 25
repeated root for b = 20 81
res(f, x + 1) = 100 = -f(-1) = 100
euclidean == hgcd: 1