Saturation seems buggy for quotient rings

[HechtiDerLachs]

Try the following example:

P, (a, b, x, y, z) = QQ[:a, :b, :x, :y, :z]
I = ideal(P, [x^3 + x*y^2 - z^2, a*y - b*z, b*x - y, a*x - z, -a*z + x^2 + y^2, -a^2 + b*y + x, -a^3 + b^2*z + z, -a^2*b + b^2*y + y])

A, _ = quo(P, I)

FA = free_module(A, 8)
rel_mat = [0 0 0 3*a*z - 2*y^2 2*x*y -2*z 0 0; 0 0 0 0 3*a*z - 2*y^2 0 2*x*y -2*z; -6*a*z^2 + 4*y^2*z -4*x*y*z 4*z^2 -6*x*z -4*y*z 3*a*z - 2*y^2 -2*x*z 2*x*y; 0 0 0 6*x*y^2 - 9*z^2 2*y^3 6*x*z 4*x*y^2 -4*y*z; 0 0 0 -3*a*z + 2*y^2 -2*x*y 2*z 0 0; 0 0 0 0 -2*x*y^2 + 3*z^2 0 2*b*z^2 - 2*y^3 -2*x*z; 4*x*y^2*z - 6*z^3 -4*b*z^3 + 4*y^3*z 4*x*z^2 -6*a*z^2 + 6*y^2*z -4*x*y*z -2*x*y^2 + 3*z^2 -2*a*z^2 + 2*y^2*z 2*b*z^2 - 2*y^3; 0 0 0 -6*x*y^2 + 9*z^2 -2*y^3 -6*x*z -4*x*y^2 4*y*z]
rel_mat = matrix_space(A, 8, 8)(rel_mat)

U, inc_U = sub(FA, rel_mat)

J = ideal(A, x)

# saturation(U, J) # throws an error, too!

U_sat = Oscar._saturation(U.sub, J)

U_sat, _ = sub(FA, gens(U_sat))

all(x in U_sat for x in gens(U))

In my opinion the last statement should be true, but it returns false for me. The _saturation method delegates directly to singular and its signature suggests that it is suitable for quotient rings. Apparently this is an example where we get wrong results and it seems to me that the problem is actually inside Singular.

Ping @hannes14 , @wdecker , @jankoboehm .

I will try to provide a temporary workaround, soon. But I thought that this might be something that Singular developers want to look at?

[wdecker]

saturation(U, J) # throws an error, too!

There is a problem here in

M_quo = isdefined(M, :quo) ? M.quo : SubModuleOfFreeModule(ambient_free_module(M), Vector{elem_type(ambient_free_module(M))}())
U = M.sub+M.quo

in SubquoModule.jl:1926. I will provide a fix, also for quotient.

[wdecker]

saturation(U, J) # throws an error, too!

There is a problem here in

M_quo = isdefined(M, :quo) ? M.quo : SubModuleOfFreeModule(ambient_free_module(M), Vector{elem_type(ambient_free_module(M))}()) U = M.sub+M.quo

in SubquoModule.jl:1926. I will provide a fix, also for quotient.

Done with PR #5498

[wdecker]

Here is Singular code to play with. I start with the example given by @HechtiDerLachs over the polynomial ring:

ring P =0, (a, b, x, y, z), dp;
vector s1 = [0,0,0,3*a*z - 2*y^2,2*x*y,-2*z,0,0];
vector s2 = [0,0,0,0,3*a*z - 2*y^2,0,2*x*y,-2*z];
vector s3 = [-6*a*z^2 + 4*y^2*z,-4*x*y*z,4*z^2,-6*x*z,-4*y*z,3*a*z - 2*y^2,-2*x*z, 2*x*y];
vector s4 = [0,0, 0,6*x*y^2 - 9*z^2,2*y^3,6*x*z,4*x*y^2,-4*y*z];
vector s5 = [0,0,0,-3*a*z + 2*y^2,-2*x*y,2*z,0,0];
vector s6 = [0,0,0,0,-2*x*y^2 + 3*z^2,0,2*b*z^2 - 2*y^3,-2*x*z];
vector s7 = [4*x*y^2*z - 6*z^3,-4*b*z^3 + 4*y^3*z,4*x*z^2,-6*a*z^2 + 6*y^2*z,-4*x*y*z,-2*x*y^2 + 3*z^2,-2*a*z^2 + 2*y^2*z,2*b*z^2 - 2*y^3];
vector s8 = [0,0,0,-6*x*y^2 + 9*z^2,-2*y^3,-6*x*z,-4*x*y^2,4*y*z];
module m = s1,s2,s3,s4,s5,s6,s7,s8;
ideal J = x;
module Q = quotient(m, J);
lift(Q, m[1]);

Result seems o.k.

Now do the same over the quotient ring: add two lines following the definition of P:

ideal I = x^3 + x*y^2 - z^2, a*y - b*z, b*x - y, a*x - z, -a*z + x^2 + y^2, -a^2 + b*y + x, -a^3 + b^2*z + z, -a^2*b + b^2*y + y;
qring A = std(I);

Here we are in trouble.

So it seems to me, as suggested by @HechtiDerLachs, that the problem w.r.t. quotient rings lies within Singular.

Ping @hannes14 , @ederc , @jankoboehm .