Missing divrem(a::EuclideanRingResidueRingElem{QQPolyRingElem}, b::EuclideanRingResidueRingElem{QQPolyRingElem})

[YueRen]

My Durham Workshop group (@sfitz01,@DanielGreenTripp,@mlodyjesienin) are working on some deformation-theoretic computations, and we noticed that we were unable to work over the coefficient ring QQ[ε]/(ε^2):

Code to reproduce the error:

S,ε = polynomial_ring(QQ,"ε")
S,pi = quo(S,ε^2)
R,(x1,x2,x3) = S["x1","x2","x3"]
intersect(ideal([-ε*x1*x2 + x1*x3 - x2*x3 + x3^2]),ideal([x1^2*x2*x3 - x1*x2^2*x3 + x1*x2*x3^2]))

The error:

julia> intersect(ideal([-ε*x1*x2 + x1*x3 - x2*x3 + x3^2]),ideal([x1^2*x2*x3 - x1*x2^2*x3 + x1*x2*x3^2]))
ERROR: function divrem is not implemented for arguments
EuclideanRingResidueRingElem{QQPolyRingElem}: ε
EuclideanRingResidueRingElem{QQPolyRingElem}: -ε

Stacktrace:
 [1] divrem(a::EuclideanRingResidueRingElem{QQPolyRingElem}, b::EuclideanRingResidueRingElem{QQPolyRingElem})
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/AqzuZ/src/algorithms/GenericFunctions.jl:50
 [2] divrem(a::EuclideanRingResidueRingElem{QQPolyRingElem}, b::EuclideanRingResidueRingElem{QQPolyRingElem})
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/AqzuZ/src/AbstractAlgebra.jl:65
 [3] gcdx(a::EuclideanRingResidueRingElem{QQPolyRingElem}, b::EuclideanRingResidueRingElem{QQPolyRingElem})
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/AqzuZ/src/algorithms/GenericFunctions.jl:306
 [4] nemoRingExtGcd(a::Ptr{Nothing}, b::Ptr{Nothing}, s::Ptr{Ptr{Nothing}}, t::Ptr{Ptr{Nothing}}, cf::Ptr{Nothing})
   @ Singular.libSingular ~/.julia/packages/Singular/JyB5B/src/libsingular/nemo/Rings.jl:200
 [5] id_Intersection
   @ ~/.julia/packages/CxxWrap/eWADG/src/CxxWrap.jl:668 [inlined]
 [6] intersection(I::Singular.sideal{Singular.spoly{…}}, J::Singular.sideal{Singular.spoly{…}})
   @ Singular ~/.julia/packages/Singular/JyB5B/src/ideal/ideal.jl:508
 [7] intersect(I::MPolyIdeal{AbstractAlgebra.Generic.MPoly{…}}, Js::MPolyIdeal{AbstractAlgebra.Generic.MPoly{…}})
   @ Oscar ~/.julia/dev/Oscar/src/Rings/mpoly-ideals.jl:217
 [8] top-level scope
   @ REPL[265]:1
Some type information was truncated. Use `show(err)` to see complete types.

[lgoettgens]

I hoped that you could temporarily circumvent the problem by using a multivariate poly ring in one variable instead of a univariate (e.g. S,(ε,) = polynomial_ring(QQ,["ε"])), but unfortunately this just runs into a different error:

julia> intersect(ideal([-ε*x1*x2 + x1*x3 - x2*x3 + x3^2]),ideal([x1^2*x2*x3 - x1*x2^2*x3 + x1*x2*x3^2]))
ERROR: function is_zero_divisor_with_annihilator is not implemented for argument
MPolyQuoRingElem{QQMPolyRingElem}: -ε

Stacktrace:
 [1] is_zero_divisor_with_annihilator(a::MPolyQuoRingElem{QQMPolyRingElem})
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/aPVcM/src/algorithms/GenericFunctions.jl:441
 [2] nemoRingAnn(a::Ptr{Nothing}, cf::Ptr{Nothing})
   @ Singular.libSingular ~/.julia/packages/Singular/JyB5B/src/libsingular/nemo/Rings.jl:170
 [3] id_Intersection
   @ ~/.julia/packages/CxxWrap/eWADG/src/CxxWrap.jl:668 [inlined]
 [4] intersection(I::Singular.sideal{Singular.spoly{…}}, J::Singular.sideal{Singular.spoly{…}})
   @ Singular ~/.julia/packages/Singular/JyB5B/src/ideal/ideal.jl:508
 [5] intersect(I::MPolyIdeal{AbstractAlgebra.Generic.MPoly{…}}, Js::MPolyIdeal{AbstractAlgebra.Generic.MPoly{…}})
   @ Oscar ~/.julia/packages/Oscar/KByfV/src/Rings/mpoly-ideals.jl:217
 [6] top-level scope
   @ REPL[10]:1

Iirc @thofma has the most knowledge about the EuclideanRingResidueRingElem type

[YueRen]

A way to circumvent it is to add epsilon to the variables of your multivariate polynomial ring and epsilon^2 to the generators of your ideal.

[thofma]

Needs to be fixed upstream, but in the meantime you could do:

julia> function _gcdx(a, b, c)
         g, s, t = gcdx(a, b)
         h, u, v = gcdx(g, c)
         gg, ss, tt, uu = (h, s*u, t * u, v)
         @assert gg == ss * a + tt * b + uu * c
         return gg, ss, tt, uu
       end
_gcdx (generic function with 1 method)

julia> function Oscar.gcdx(a::EuclideanRingResidueRingElem{QQPolyRingElem}, b::EuclideanRingResidueRingElem{QQPolyRingElem})
         s = modulus(parent(a))
         g, s, t, u = _gcdx(lift(a), lift(b), s)
         return parent(a)(g), parent(a)(s), parent(b)(t)
       end

julia> begin S,ε = polynomial_ring(QQ,"ε")
       S,pi = quo(S,ε^2)
       R,(x1,x2,x3) = S["x1","x2","x3"]
       intersect(ideal([-ε*x1*x2 + x1*x3 - x2*x3 + x3^2]),ideal([x1^2*x2*x3 - x1*x2^2*x3 + x1*x2*x3^2]))
       end
Ideal generated by
  ε*x1^2*x2*x3 - ε*x1*x2^2*x3 + ε*x1*x2*x3^2

The trick with adding epsilon as a variable and epsilon^2 as a generator should probably be much faster. Would be nice if we could do this automatically.