Validity of generalized Euclidean algorithm without primitive part

[blegat]

With integers, it works:

julia> using TypedPolynomials

julia> @polyvar x y z
(x, y, z)

julia> p = y^6 + x^3*z^3 - x^4*y^2 - x^5*z
y⁶ + x³z³ - x⁴y² - x⁵z

julia> q = z^3 + 2y^3 + x^2*z
z³ + 2y³ + x²z

julia> gcd(p, q)
1

With floating points, it fails:

julia> p = 1.0y^6 + x^3*z^3 - x^4*y^2 - x^5*z
y⁶ + x³z³ - x⁴y² - x⁵z

julia> q = 1.0z^3 + 2y^3 + x^2*z
z³ + 2.0y³ + x²z

julia> gcd(p, q)
z³ + 2.0y³ + x²z

Internally, it's computing this one:

julia> a = -y^6 - 4.0x^3*y^3 + x^4*y^2 - 4x^6
-y⁶ - 4.0x³y³ + x⁴y² - 4.0x⁶

julia> b = -2.0y^6 - 6x^3*y^3 - 4x^6
-2.0y⁶ - 6.0x³y³ - 4.0x⁶

julia> gcd(a, b)
-9.2908732547072e13y - 9.2908732547072e13x

[blegat]

With Int:

julia> p = -2y^5*z - x^3*z^3 - 2x^6
-2y⁵z - x³z³ - 2x⁶

julia> q = 2y*z^2 - 2y^3 - 2x*z^2 - x*y^2 + x^3
2yz² - 2y³ - 2xz² - xy² + x³

julia> gcd(p, q, GeneralizedEuclideanAlgorithm(false, false))
2yz⁵ - 8y²z⁴ - 2y³z³ + 16y⁴z² - 2y⁵z - 8y⁶ - 2xz⁵ + 16xyz⁴ - xy²z³ - 8xy³z² - 8xy⁵ - 8x²z⁴ - 8x²y²z² - 2x²y⁴

It is failing because internally it is calling

julia> a = (62*z^10) + (24*z^9)*y + (z^8)*y^2 + (-36*z^7)*y^3 + (131*z^6)*y^4 + (-16*z^5)*y^5 + (-159*z^4)*y^6 + (20*z^3)*y^7 + (-60*z^2)*y^8 + (8*z)*y^9 + (28)*y^10
62z¹⁰ + 24yz⁹ + y²z⁸ - 36y³z⁷ + 131y⁴z⁶ - 16y⁵z⁵ - 159y⁶z⁴ + 20y⁷z³ - 60y⁸z² + 8y⁹z + 28y¹⁰

julia> b = (-124*z^10) + (-48*z^9)*y + (-124*z^8)*y^2 + (-511*z^6)*y^4 + (20*z^5)*y^5 + (-304*z^4)*y^6 + (-4*z^3)*y^7 + (56*z^2)*y^8 + (-4*z)*y^9 + (44)*y^10
-124z¹⁰ - 48yz⁹ - 124y²z⁸ - 511y⁴z⁶ + 20y⁵z⁵ - 304y⁶z⁴ - 4y⁷z³ + 56y⁸z² - 4y⁹z + 44y¹⁰

julia> gcd(a, b, GeneralizedEuclideanAlgorithm(false, false))
8108232830836140562z² + 4899916394579099648yz + 5207037920024917257y²

The euclidean gcd algorithm makes the integer grows until overflow.