What is the correct behavior of gcd for polynomials over an inexact base ring?

The following behavior is observed (compare #16190):

sage: P.<x> = PolynomialRing(RR)
sage: gcd((x-1)^2*(x-2)*(x-3), (x-1)*(x+1))
x - 1.00000000000000
sage: gcd((x-1)^2*(x-2)*(x-3), (x-1)*(x+1)^2)
1.00000000000000

Since the base ring RR is inexact, whether or not there is a common factor is not mathematically well-defined (i.e., it is not uniform over the possible real polynomials represented by a given floating-point approximation). However, I would contend that the most useful behavior would be to return the largest possible common factor that is consistent with the approximation (i.e., modulo which the two polynomials each reduce to machine zero).

This comes up in particular when trying to tell whether a polynomial is squarefree. Over RR, this is needed in order to carry out root-counting algorithms; see #20243.

sage: P.<x> = PolynomialRing(RR)
sage: pol = (x-1)^2*(x-2)^2*(x-3)
sage: pol.is_squarefree() # computed using GCD
True
sage: pari(pol).issquarefree # should return 0
1
sage: pari(pol).polsturm()
...
PariError: domain error in polsturm: issquarefree(pol) = 0

Note that PARI isn't getting it right either; I've reported that upstream.

Upstream: Reported upstream. No feedback yet.

Component: algebra

Keywords: polynomials, inexact base rings, gcd

Issue created by migration from https://trac.sagemath.org/ticket/20282

sagemath / sage

What is the correct behavior of gcd for polynomials over an inexact base ring? #20282