diff --git a/Project.toml b/Project.toml
index 3e3755d..53116bd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "AlgebraicSolving"
 uuid = "66b61cbe-0446-4d5d-9090-1ff510639f9d"
 authors = ["ederc <ederc@mathematik.uni-kl.de>", "Mohab Safey El Din <Mohab.Safey@lip6.fr", "Rafael Mohr <rafael.mohr@lip6.fr>"]
-version = "0.4.5"
+version = "0.4.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -17,7 +17,7 @@ msolve_jll = "6d01cc9a-e8f6-580e-8c54-544227e08205"
 
 [compat]
 LoopVectorization = "0.12"
-Nemo = "0.35.1, 0.36, 0.37, 0.38"
+Nemo = "0.35.1, 0.36, 0.37, 0.38, 0.39"
 StaticArrays = "1"
 julia = "1.6"
 LinearAlgebra = "1.6"
diff --git a/docs/Project.toml b/docs/Project.toml
index b9eb78f..ee69314 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -9,4 +9,3 @@ Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 [compat]
 Documenter = "0.26"
 msolve_jll = "0.6.1"
-Nemo = "0.35.1"
diff --git a/src/algorithms/groebner-bases.jl b/src/algorithms/groebner-bases.jl
index 5bff26b..6fab9d2 100644
--- a/src/algorithms/groebner-bases.jl
+++ b/src/algorithms/groebner-bases.jl
@@ -27,13 +27,13 @@ At the moment the underlying algorithm is based on variants of Faugère's F4 Alg
 julia> using AlgebraicSolving
 
 julia> R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"], ordering=:degrevlex)
-(Multivariate polynomial ring in 3 variables over GF(101), fpMPolyRingElem[x, y, z])
+(Multivariate polynomial ring in 3 variables over GF(101), Nemo.FqMPolyRingElem[x, y, z])
 
 julia> I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
-fpMPolyRingElem[x + 2*y + 2*z + 100, x^2 + 2*y^2 + 2*z^2 + 100*x, 2*x*y + 2*y*z + 100*y]
+Nemo.FqMPolyRingElem[x + 2*y + 2*z + 100, x^2 + 2*y^2 + 2*z^2 + 100*x, 2*x*y + 2*y*z + 100*y]
 
 julia> eliminate(I, 2)
-1-element Vector{fpMPolyRingElem}:
+1-element Vector{Nemo.FqMPolyRingElem}:
  z^4 + 38*z^3 + 95*z^2 + 95*z
 ```
 """
@@ -81,20 +81,20 @@ At the moment the underlying algorithm is based on variants of Faugère's F4 Alg
 julia> using AlgebraicSolving
 
 julia> R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"], ordering=:degrevlex)
-(Multivariate polynomial ring in 3 variables over GF(101), fpMPolyRingElem[x, y, z])
+(Multivariate polynomial ring in 3 variables over GF(101), Nemo.FqMPolyRingElem[x, y, z])
 
 julia> I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
-fpMPolyRingElem[x + 2*y + 2*z + 100, x^2 + 2*y^2 + 2*z^2 + 100*x, 2*x*y + 2*y*z + 100*y]
+Nemo.FqMPolyRingElem[x + 2*y + 2*z + 100, x^2 + 2*y^2 + 2*z^2 + 100*x, 2*x*y + 2*y*z + 100*y]
 
 julia> groebner_basis(I)
-4-element Vector{fpMPolyRingElem}:
+4-element Vector{Nemo.FqMPolyRingElem}:
  x + 2*y + 2*z + 100
  y*z + 82*z^2 + 10*y + 40*z
  y^2 + 60*z^2 + 20*y + 81*z
  z^3 + 28*z^2 + 64*y + 13*z
 
 julia> groebner_basis(I, eliminate=2)
-1-element Vector{fpMPolyRingElem}:
+1-element Vector{Nemo.FqMPolyRingElem}:
  z^4 + 38*z^3 + 95*z^2 + 95*z
 ```
 """
diff --git a/src/algorithms/normal-forms.jl b/src/algorithms/normal-forms.jl
index 51e9e38..0e86389 100644
--- a/src/algorithms/normal-forms.jl
+++ b/src/algorithms/normal-forms.jl
@@ -27,10 +27,10 @@ Gröbner basis, then this one is first computed.
 julia> using AlgebraicSolving
 
 julia> R, (x,y) = polynomial_ring(GF(101),["x","y"])
-(Multivariate polynomial ring in 2 variables over GF(101), fpMPolyRingElem[x, y])
+(Multivariate polynomial ring in 2 variables over GF(101), Nemo.FqMPolyRingElem[x, y])
 
 julia> I = Ideal([y*x^3+12-y, x+y])
-fpMPolyRingElem[x^3*y + 100*y + 12, x + y]
+Nemo.FqMPolyRingElem[x^3*y + 100*y + 12, x + y]
 
 julia> f = 2*x^2+7*x*y
 2*x^2 + 7*x*y
@@ -74,18 +74,18 @@ Gröbner basis, then this one is first computed.
 julia> using AlgebraicSolving
 
 julia> R, (x,y) = polynomial_ring(GF(101),["x","y"])
-(Multivariate polynomial ring in 2 variables over GF(101), fpMPolyRingElem[x, y])
+(Multivariate polynomial ring in 2 variables over GF(101), Nemo.FqMPolyRingElem[x, y])
 
 julia> I = Ideal([y*x^3+12-y, x+y])
-fpMPolyRingElem[x^3*y + 100*y + 12, x + y]
+Nemo.FqMPolyRingElem[x^3*y + 100*y + 12, x + y]
 
 julia> F = [2*x^2+7*x*y, x+y]
-2-element Vector{fpMPolyRingElem}:
+2-element Vector{Nemo.FqMPolyRingElem}:
  2*x^2 + 7*x*y
  x + y
 
 julia> normal_form(F,I)
-2-element Vector{fpMPolyRingElem}:
+2-element Vector{Nemo.FqMPolyRingElem}:
  96*y^2
  0
 ```
diff --git a/src/siggb/siggb.jl b/src/siggb/siggb.jl
index 1c498cf..9b14fea 100644
--- a/src/siggb/siggb.jl
+++ b/src/siggb/siggb.jl
@@ -35,20 +35,20 @@ signature and the second the underlying polynomial.
 julia> using AlgebraicSolving
 
 julia> R, vars = polynomial_ring(GF(17), ["x$i" for i in 1:4])
-(Multivariate polynomial ring in 4 variables over GF(17), fpMPolyRingElem[x1, x2, x3, x4])
+(Multivariate polynomial ring in 4 variables over GF(17), Nemo.FqMPolyRingElem[x1, x2, x3, x4])
 
 julia> F = AlgebraicSolving.cyclic(R)
-fpMPolyRingElem[x1 + x2 + x3 + x4, x1*x2 + x1*x4 + x2*x3 + x3*x4, x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4, x1*x2*x3*x4 + 16]
+Nemo.FqMPolyRingElem[x1 + x2 + x3 + x4, x1*x2 + x1*x4 + x2*x3 + x3*x4, x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4, x1*x2*x3*x4 + 16]
 
 julia> Fhom = AlgebraicSolving._homogenize(F.gens)
-4-element Vector{fpMPolyRingElem}:
+4-element Vector{Nemo.FqMPolyRingElem}:
  x1 + x2 + x3 + x4
  x1*x2 + x2*x3 + x1*x4 + x3*x4
  x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4
  x1*x2*x3*x4 + 16*x5^4
 
 julia> sig_groebner_basis(Fhom)
-7-element Vector{Tuple{Tuple{Int64, fpMPolyRingElem}, fpMPolyRingElem}}:
+7-element Vector{Tuple{Tuple{Int64, Nemo.FqMPolyRingElem}, Nemo.FqMPolyRingElem}}:
  ((1, 1), x1 + x2 + x3 + x4)
  ((2, 1), x2^2 + 2*x2*x4 + x4^2)
  ((3, 1), x2*x3^2 + x3^2*x4 + 16*x2*x4^2 + 16*x4^3)
@@ -59,7 +59,7 @@ julia> sig_groebner_basis(Fhom)
 ```
 """
 function sig_groebner_basis(sys::Vector{T}; info_level::Int=0, degbound::Int=0) where {T <: MPolyRingElem}
-    R = first(sys).parent
+    R = parent(first(sys))
     Rchar = characteristic(R)
 
     # check if input is ok
@@ -131,9 +131,9 @@ function sig_groebner_basis(sys::Vector{T}; info_level::Int=0, degbound::Int=0)
             m = monomial(SVector{nv}((Exp).(exps[j])))
             eidx = insert_in_hash_table!(basis_ht, m)
             if isone(j)
-                inver = inv(Coeff(cfs[1].data), char)
+                inver = inv(Coeff(lift(ZZ, cfs[1])), char)
             end
-            cf = isone(j) ? one(Coeff) : mul(inver, Coeff(cfs[j].data), char)
+            cf = isone(j) ? one(Coeff) : mul(inver, Coeff(lift(ZZ, cfs[j])), char)
             mons[j] = eidx
             coeffs[j] = cf
         end
@@ -181,13 +181,13 @@ function sig_groebner_basis(sys::Vector{T}; info_level::Int=0, degbound::Int=0)
         exps = [basis_ht.exponents[m].exps for m in basis.monomials[i]]
         ctx = MPolyBuildCtx(R)
         for (e, c) in zip(exps, basis.coefficients[i])
-            push_term!(ctx, c, Vector{Int}(e))
+            push_term!(ctx, coefficient_ring(R)(c), Vector{Int}(e))
         end
         pol = finish(ctx)
 
         s = basis.sigs[i]
         ctx = MPolyBuildCtx(R)
-        push_term!(ctx, 1, Vector{Int}(monomial(s).exps))
+        push_term!(ctx, one(coefficient_ring(R)), Vector{Int}(monomial(s).exps))
         sig = (Int(index(s)), finish(ctx))
 
         push!(outp, (sig, pol))
diff --git a/test/interfaces/nemo.jl b/test/interfaces/nemo.jl
index 8f28415..d3ffdcc 100644
--- a/test/interfaces/nemo.jl
+++ b/test/interfaces/nemo.jl
@@ -3,12 +3,14 @@
     F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
     cmp = (Int32[2, 2, 2], BigInt[1, 1, -2, 1, 1, 1, -1, 1, 1, 1, -3, 1], Int32[2, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 2, 0, 2, 0])
     @test AlgebraicSolving._convert_to_msolve(F) == cmp
-    R, (x,y,z) = polynomial_ring(GF(2147483659),["x","y","z"], ordering=:degrevlex)
-    F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
-    # prime is bigger than 2^31, should throw an error
-    @test_throws ErrorException AlgebraicSolving._convert_to_msolve(F)
-    R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"], ordering=:degrevlex)
-    F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
-    res = AlgebraicSolving._convert_to_msolve(F)
-    @test AlgebraicSolving._convert_finite_field_array_to_abstract_algebra(Int32(3), res..., R) == F
+    for _GF in [GF, AlgebraicSolving.Nemo.Native.GF]
+      R, (x,y,z) = polynomial_ring(_GF(2147483659),["x","y","z"], ordering=:degrevlex)
+      F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
+      # prime is bigger than 2^31, should throw an error
+      @test_throws ErrorException AlgebraicSolving._convert_to_msolve(F)
+      R, (x,y,z) = polynomial_ring(_GF(101),["x","y","z"], ordering=:degrevlex)
+      F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
+      res = AlgebraicSolving._convert_to_msolve(F)
+      @test AlgebraicSolving._convert_finite_field_array_to_abstract_algebra(Int32(3), res..., R) == F
+    end
 end