From c336a38f0fe8a8190aac4cdadea938fa986d25ba Mon Sep 17 00:00:00 2001
From: Tommy Hofmann <thofma@gmail.com>
Date: Mon, 19 Feb 2024 08:18:17 +0100
Subject: [PATCH] fix: multivariate gcd over number fields (#1407)

- handle the case of non-nice defining polynomial better
---
 src/NumField/NfAbs/MPolyGcd.jl     | 6 ++++--
 test/NfAbs/MPolyGcd.jl             | 9 +++++++--
 test/NumField/NfAbs/MPolyFactor.jl | 9 +++++++++
 3 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/NumField/NfAbs/MPolyGcd.jl b/src/NumField/NfAbs/MPolyGcd.jl
index b4d924662a..08ff7935d2 100644
--- a/src/NumField/NfAbs/MPolyGcd.jl
+++ b/src/NumField/NfAbs/MPolyGcd.jl
@@ -267,10 +267,12 @@ function _gcd(f::Hecke.Generic.MPoly{AbsSimpleNumFieldElem}, g::Hecke.Generic.MP
   Zx = Hecke.Globals.Zx
   R = RecoCtx(K)
 
-  de = lcm(lcm(map(denominator, coefficients(f))), lcm(map(denominator, coefficients(g))))
+  E = any_order(K)
+  de = lcm(lcm(map(c -> denominator(c, E), coefficients(f))),
+           lcm(map(c -> denominator(c, E), coefficients(g))))
+
   f*=de
   g*=de
-  E = any_order(K)
   lI = E*E(leading_coefficient(f)) + E*E(leading_coefficient(g))
   gl = Hecke.short_elem(lI)
   gl *= evaluate(derivative(K.pol), gen(K))  # use Kronnecker basis
diff --git a/test/NfAbs/MPolyGcd.jl b/test/NfAbs/MPolyGcd.jl
index 06857f53e9..0904c547bb 100644
--- a/test/NfAbs/MPolyGcd.jl
+++ b/test/NfAbs/MPolyGcd.jl
@@ -15,6 +15,11 @@
     @test isone(gcd(f+1,g))
 
   end
-end
-
 
+  # test non-integral defining equation
+  Qx, x = QQ["x"]
+  K, a = number_field(x^2 - 1//3*x + 1);
+  R, (u, v) = polynomial_ring(K, ["u", "v"])
+  @test gcd(a*u, a*u) == u
+  @test is_one(gcd(a*u, a^2*v))
+end
diff --git a/test/NumField/NfAbs/MPolyFactor.jl b/test/NumField/NfAbs/MPolyFactor.jl
index 3b9380a799..fac5142a74 100644
--- a/test/NumField/NfAbs/MPolyFactor.jl
+++ b/test/NumField/NfAbs/MPolyFactor.jl
@@ -65,4 +65,13 @@
    f = (x^3+5*y^3)*(x^2+2*y^2)
    # not a real test, just check that it does not crash
    @test length(factor_absolute(f)) == 3
+
+   # non-integral defining equation
+   Qx, x = QQ["x"]
+   K, a = number_field(x^2 - 1//3*x + 1);
+   R, (u, v) = polynomial_ring(K, ["u", "v"])
+   f = (u + a)*(v^2 + a)
+   fa = factor(f)
+   @test unit(fa) * prod(g^e for (g, e) in fa) == f
+
 end