Greatest common divisors over generic rings

gr_poly_gcd and gr_poly_xgcd compute a greatest common divisor of two polynomials, and with it the Bézout cofactors, over any ring whose context says that the GCD makes sense: a field (or a ring that pretends to be one, see gr_ctx_set_is_pretend_field) or a unique factorisation domain. The algorithms are compared in Polynomial GCDs and resultants.

Normalisation of the result#

A GCD is unique only up to a unit factor. The non-underscore functions return the canonical associate (gr_poly_canonical_associate): the monic GCD over a field, and over $\mathbb{Z}$ the GCD with positive leading coefficient. $\gcd(0, 0) = 0$, and $\gcd(a, 0)$ is the canonical associate of $a$. For the extended GCD, $G = SA + TB$ holds exactly; over a field $G$ is the monic GCD, $\deg S < \deg B - \deg G$ and $\deg T < \deg A - \deg G$. Over a ring that is not a field, there may be no cofactors for the GCD itself (in $\mathbb{Z}[x]$, $1 = \gcd(2, x)$ is not a combination $2S + xT$), and gr_poly_xgcd then returns a constant multiple of the GCD: the subresultant pseudo-remainder sequence of the primitive parts of $A$ and $B$, which carries the cofactors along, ends with a multiple of their GCD, and this is multiplied by the least common multiple of the contents of $A$ and $B$.

Both functions return GR_UNABLE if a leading coefficient of an input is not provably nonzero, which can happen over inexact rings (see Inexact coefficients and weak normalisation), and if the ring is neither a field nor a unique factorisation domain. An algorithm that meets a non-invertible leading coefficient returns GR_DOMAIN, or GR_UNABLE when the ring pretends to be a field.

Algorithm selection#

gr_poly_gcd and gr_poly_xgcd call the ring methods GR_METHOD_POLY_GCD and GR_METHOD_POLY_XGCD. The defaults, _gr_poly_gcd_generic and _gr_poly_xgcd_generic, choose:

• over a field (or a pretended field) that is finite, the half-gcd algorithm (Half-GCD) when both inputs have length at least 200, with the Euclidean algorithm for remainders shorter than 200 and inside the half-gcd recursion below length 100;
• over other fields, such as $\mathbb{Q}$, number fields, real and complex balls, the Euclidean algorithm, since the half-gcd algorithm suffers from coefficient growth over exact infinite fields and from numerical instability over approximate fields;
• over unique factorisation domains that are not fields, such as $\mathbb{Z}$ and $\mathbb{Z}[y]$, the subresultant algorithm: the contents are removed, the primitive parts are run through the subresultant pseudo-remainder sequence, and the GCD of the contents is multiplied back;
• otherwise GR_UNABLE.

Several rings override the methods with tuned algorithms:

There is no override for $\mathbb{Q}$: over fmpq the generic Euclidean algorithm is used, whose intermediate coefficients can grow large; for large inputs fmpq_poly_gcd is much faster.

The underscore functions _gr_poly_gcd_euclidean, _gr_poly_gcd_hgcd and _gr_poly_gcd_subresultant, and the corresponding xgcd functions, run one algorithm unconditionally; _gr_poly_gcd_hgcd and _gr_poly_xgcd_hgcd take the half-gcd recursion cutoff inner_cutoff (respectively hgcd_cutoff) and the length cutoff below which the outer loop switches to the Euclidean algorithm.

Conventions of the underscore functions#

_gr_poly_gcd and the algorithm-specific underscore GCD functions require $\mathit{lenA} \ge \mathit{lenB} \ge 1$ and nonzero leading coefficients, write at most $\mathit{lenB}$ coefficients to $G$ and its length to *lenG, and need not canonicalise the result.

_gr_poly_xgcd and its variants require $\mathit{lenA} \ge \mathit{lenB} \ge 2$ and do not allow aliasing. The output arrays must have room for $\mathit{lenB}$ coefficients ($G$), $\mathit{lenB} - 1$ ($S$) and $\mathit{lenA} - 1$ ($T$), even if the results are shorter. The length of $G$ is written to *lenG, and coefficients of $G$ beyond it may contain garbage; on success the unused high coefficients of $S$ and $T$ are zeroed. A length $\mathit{lenG} = 0$ is possible if the algorithm fails, for example at a non-invertible leading coefficient.

To wrap a different underscore implementation with the memory management, aliasing rules and normalisation of the non-underscore functions, use gr_poly_gcd_wrapper and gr_poly_xgcd_wrapper.

Example#

The program computes an extended GCD of $A = (x + 1)(x - 2)$ and $B = (x + 1)(x + 3)$ over $\mathbb{Z}$ and over $\mathbb{F}_7$. The GCD over $\mathbb{Z}$ is $x + 1$, but it is not a $\mathbb{Z}[x]$-combination of $A$ and $B$, since $B - A = 5(x + 1)$; the subresultant algorithm returns $G = 5(x + 1)$ with $S = -1$, $T = 1$. The factor $5$ is the resultant of the cofactors $x - 2$ and $x + 3$. Over $\mathbb{F}_7$ the monic GCD is returned.

#include <flint/flint.h>
#include <flint/gr.h>
#include <flint/gr_poly.h>

int main(void)
{
    gr_ctx_t ZZ, F;
    gr_poly_t a, b, g, s, t;

    gr_ctx_init_fmpz(ZZ);
    gr_poly_init(a, ZZ);
    gr_poly_init(b, ZZ);
    gr_poly_init(g, ZZ);
    gr_poly_init(s, ZZ);
    gr_poly_init(t, ZZ);

    /* a = (x + 1)(x - 2), b = (x + 1)(x + 3) */
    GR_MUST_SUCCEED(gr_poly_set_str(a, "x^2 - x - 2", "x", ZZ));
    GR_MUST_SUCCEED(gr_poly_set_str(b, "x^2 + 4*x + 3", "x", ZZ));

    GR_MUST_SUCCEED(gr_poly_xgcd(g, s, t, a, b, ZZ));
    flint_printf("over Z:    g = "); gr_poly_print(g, ZZ);
    flint_printf(", s = "); gr_poly_print(s, ZZ);
    flint_printf(", t = "); gr_poly_print(t, ZZ); flint_printf("\n");

    /* the same inputs over GF(7) */
    gr_ctx_init_nmod(F, 7);
    {
        gr_poly_t a7, b7, g7, s7, t7;
        gr_poly_init(a7, F); gr_poly_init(b7, F); gr_poly_init(g7, F);
        gr_poly_init(s7, F); gr_poly_init(t7, F);
        GR_MUST_SUCCEED(gr_poly_set_gr_poly_other(a7, a, ZZ, F));
        GR_MUST_SUCCEED(gr_poly_set_gr_poly_other(b7, b, ZZ, F));
        GR_MUST_SUCCEED(gr_poly_xgcd(g7, s7, t7, a7, b7, F));
        flint_printf("over GF(7): g = "); gr_poly_print(g7, F);
        flint_printf(", s = "); gr_poly_print(s7, F);
        flint_printf(", t = "); gr_poly_print(t7, F); flint_printf("\n");
        gr_poly_clear(a7, F); gr_poly_clear(b7, F); gr_poly_clear(g7, F);
        gr_poly_clear(s7, F); gr_poly_clear(t7, F);
    }
    gr_ctx_clear(F);

    gr_poly_clear(a, ZZ);
    gr_poly_clear(b, ZZ);
    gr_poly_clear(g, ZZ);
    gr_poly_clear(s, ZZ);
    gr_poly_clear(t, ZZ);
    gr_ctx_clear(ZZ);
    return 0;
}

over Z:    g = 5 + 5*x, s = -1, t = 1
over GF(7): g = 1 + x, s = 4, t = 3

The GCD over $\mathbb{F}_7$ satisfies $4A + 3B = 7x^2 + 8x + 1 \equiv x + 1$.