diff --git a/src/GroebnerWalk/Examples.jl b/src/GroebnerWalk/Examples.jl new file mode 100644 index 000000000..34103e9bc --- /dev/null +++ b/src/GroebnerWalk/Examples.jl @@ -0,0 +1,791 @@ +function katsura4() + dim = 5 + ve = [1, 1, 1, 1, 1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x0, x1, x2, x3, x4) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x0", "x1", "x2", "x3", "x4"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = 2 * x4^2 + 2 * x3^2 + 2 * x2^2 + 2 * x1^2 + x0^2 - x0 + f2 = 2 * x3 * x4 + 2 * x3 * x2 + 2 * x2 * x1 + 2 * x1 * x0 - x1 + f3 = 2 * x4 * x2 + 2 * x3 * x1 + 2 * x1^2 + 2 * x2 * x0 - x2 + f4 = 2 * x4 * x1 + 2 * x2 * x1 + 2 * x3 * x0 - x3 + f5 = 2 * x4 + 2 * x3 + 2 * x2 + 2 * x1 + x0 - 1 + return Singular.Ideal(R, [f1, f2, f3, f4, f5]) +end + + +#Katsura6 +function katsura5() + dim = 6 + ve = [1, 1, 1, 1, 1, 1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x0, x1, x2, x3, x4, x5) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x0", "x1", "x2", "x3", "x4", "x5"], + ordering = Singular.ordering_M(StartOrd), + ) + f1 = 2 * x5^2 + 2 * x4^2 + 2 * x3^2 + 2 * x2^2 + 2 * x1^2 + x0^2 - x0 + f2 = 2 * x5 * x4 + 2 * x4 * x3 + 2 * x3 * x2 + 2 * x2 * x1 + 2 * x1 * x0 - x1 + f3 = 2 * x5 * x3 + 2 * x4 * x2 + 2 * x3 * x1 + 2 * x1^2 + 2 * x2 * x0 - x2 + f4 = 2 * x5 * x2 + 2 * x4 * x1 + 2 * x2 * x1 + 2 * x3 * x0 - x3 + f5 = x2^2 + 2 * x5 * x0 + 2 * x4 * x0 + 2 * x3 * x0 - x4 + f6 = 2 * x5 + 2 * x4 + 2 * x3 + 2 * x2 + 2 * x1 + x0 - 1 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6]) +end + + +function katsura6() + dim = 7 + ve = [1, 1, 1, 1, 1, 1,1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x0,x1,x2,x3,x4,x5,x6) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x0","x1","x2","x3","x4","x5","x6"], + ordering = Singular.ordering_M(StartOrd), + ) + f1 = x0+2*x1+2*x2+2*x3+2*x4+2*x5+2*x6-1 + f2 = 2*x2*x3+2*x1*x4+2*x0*x5+2*x1*x6-x5 + f3 = x2^2+2*x1*x3+2*x0*x4+2*x1*x5+2*x2*x6-x4 + f4 = 2*x1*x2+2*x0*x3+2*x1*x4+2*x2*x5+2*x3*x6-x3 + f5 = x1^2+2*x0*x2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6-x2 + f6 = 2*x0*x1+2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6-x1 + f7 = x0^2+2*x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2-x0 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7]) +end + + + +function katsura7() + dim = 8 + ve = [1, 1, 1, 1, 1, 1, 1,1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (u0, u1, u2, u3, u4, u5, u6, u7) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["u0","u1","u2","u3", "u4","u5","u6", "u7"], + ordering = Singular.ordering_M(StartOrd), + ) +f1= u0+2*u1+2*u2+2*u3+2*u4+2*u5+2*u6+2*u7-1 +f2= u3^2+2*u2*u4+2*u1*u5+2*u0*u6+2*u1*u7-u6 +f3= 2*u2*u3+2*u1*u4+2*u0*u5+2*u1*u6+2*u2*u7-u5 +f4= u2^2+2*u1*u3+2*u0*u4+2*u1*u5+2*u2*u6+2*u3*u7-u4 +f5= 2*u1*u2+2*u0*u3+2*u1*u4+2*u2*u5+2*u3*u6+2*u4*u7-u3 +f6= u1^2+2*u0*u2+2*u1*u3+2*u2*u4+2*u3*u5+2*u4*u6+2*u5*u7-u2 +f7= 2*u0*u1+2*u1*u2+2*u2*u3+2*u3*u4+2*u4*u5+2*u5*u6+2*u6*u7-u1 +f8= u0^2+2*u1^2+2*u2^2+2*u3^2+2*u4^2+2*u5^2+2*u6^2+2*u7^2-u0 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7, f8]) +end + +function katsura8() + dim = 9 + ve = [1, 1, 1, 1, 1, 1,1,1,1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (u0,u1,u2,u3,u4,u5,u6,u7,u8) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["u0","u1","u2","u3","u4","u5","u6","u7","u8"], + ordering = Singular.ordering_M(StartOrd), + ) + f1 = u0+2*u1+2*u2+2*u3+2*u4+2*u5+2*u6+2*u7+2*u8-1 + f2 = 2*u3*u4+2*u2*u5+2*u1*u6+2*u0*u7+2*u1*u8-u7 + f3 = u3^2+2*u2*u4+2*u1*u5+2*u0*u6+2*u1*u7+2*u2*u8-u6 + f4 = 2*u2*u3+2*u1*u4+2*u0*u5+2*u1*u6+2*u2*u7+2*u3*u8-u5 + f5 = u2^2+2*u1*u3+2*u0*u4+2*u1*u5+2*u2*u6+2*u3*u7+2*u4*u8-u4 + f6 = 2*u1*u2+2*u0*u3+2*u1*u4+2*u2*u5+2*u3*u6+2*u4*u7+2*u5*u8-u3 + f7 = u1^2+2*u0*u2+2*u1*u3+2*u2*u4+2*u3*u5+2*u4*u6+2*u5*u7+2*u6*u8-u2 + f8 = 2*u0*u1+2*u1*u2+2*u2*u3+2*u3*u4+2*u4*u5+2*u5*u6+2*u6*u7+2*u7*u8-u1 + f9 = u0^2+2*u1^2+2*u2^2+2*u3^2+2*u4^2+2*u5^2+2*u6^2+2*u7^2+2*u8^2-u0 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7, f8, f9]) +end + + +function cyclic7() + dim = 7 + ve = [1, 1, 1, 1, 1, 1, 1] + example = "Cyclic7" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5, x6, x7) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5", "x6", "x7"], + ordering = Singular.ordering_M(StartOrd), + ) + f1 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + f2 = x1 * x2 + x2 * x3 + x3 * x4 + x4 * x5 + x5 * x6 + x1 * x7 + x6 * x7 + f3 = + x1 * x2 * x3 + + x2 * x3 * x4 + + x3 * x4 * x5 + + x4 * x5 * x6 + + x1 * x2 * x7 + + x1 * x6 * x7 + + x5 * x6 * x7 + f4 = + x1 * x2 * x3 * x4 + + x2 * x3 * x4 * x5 + + x3 * x4 * x5 * x6 + + x1 * x2 * x3 * x7 + + x1 * x2 * x6 * x7 + + x1 * x5 * x6 * x7 + + x4 * x5 * x6 * x7 + f5 = + x1 * x2 * x3 * x4 * x5 + + x2 * x3 * x4 * x5 * x6 + + x1 * x2 * x3 * x4 * x7 + + x1 * x2 * x3 * x6 * x7 + + x1 * x2 * x5 * x6 * x7 + + x1 * x4 * x5 * x6 * x7 + + x3 * x4 * x5 * x6 * x7 + f6 = + x1 * x2 * x3 * x4 * x5 * x6 + + x1 * x2 * x3 * x4 * x5 * x7 + + x1 * x2 * x3 * x4 * x6 * x7 + + x1 * x2 * x3 * x5 * x6 * x7 + + x1 * x2 * x4 * x5 * x6 * x7 + + x1 * x3 * x4 * x5 * x6 * x7 + + x2 * x3 * x4 * x5 * x6 * x7 + f7 = x1 * x2 * x3 * x4 * x5 * x6 * x7 - 1 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7]) +end + + +function cyclic6() + dim = 6 + ve = [1, 1, 1, 1, 1, 1] + example = "Cyclic6" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5", "x6"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = x1 + x2 + x3 + x4 + x5 + x6 + f2 = x1 * x2 + x2 * x3 + x3 * x4 + x4 * x5 + x1 * x6 + x5 * x6 + f3 = + x1 * x2 * x3 + + x2 * x3 * x4 + + x3 * x4 * x5 + + x1 * x2 * x6 + + x1 * x5 * x6 + + x4 * x5 * x6 + f4 = + x1 * x2 * x3 * x4 + + x2 * x3 * x4 * x5 + + x1 * x2 * x3 * x6 + + x1 * x2 * x5 * x6 + + x1 * x4 * x5 * x6 + + x3 * x4 * x5 * x6 + f5 = + x1 * x2 * x3 * x4 * x5 + + x1 * x2 * x3 * x4 * x6 + + x1 * x2 * x3 * x5 * x6 + + x1 * x2 * x4 * x5 * x6 + + x1 * x3 * x4 * x5 * x6 + + x2 * x3 * x4 * x5 * x6 + f6 = x1 * x2 * x3 * x4 * x5 * x6 - 1 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6]) +end + +function cyclic5() + dim = 5 + ve = [1, 1, 1, 1, 1] + example = "Cyclic5" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (v, w, x, y, z) = Singular.PolynomialRing( + Singular.n_ZpField(32003), + ["v", "w", "x", "y", "z"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = v + w + x + y + z + f2 = v * w + w * x + x * y + v * z + y * z + f3 = v * w * x + w * x * y + v * w * z + v * y * z + x * y * z + f4 = + v * w * x * y + + v * w * x * z + + v * w * y * z + + v * x * y * z + + w * x * y * z + f5 = v * w * x * y * z - 1 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5]) +end + +function cyclic4() + dim = 4 + ve = [1, 1, 1, 1] + example = "Cyclic4" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (w, x, y, z) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["w", "x", "y", "z"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = w + x + y + z + f2 = w * x + x * y + w * z + y * z + f3 = w * x * y + w * x * z + w * y * z + x * y * z + f4 = w * x * y * z - 1 + return Singular.Ideal(R, [f1, f2, f3, f4]) +end + +function cyclic8() + dim = 8 + ve = [1, 1, 1, 1, 1 ,1 ,1 ,1] + example = "Cyclic8" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1,x2,x3,x4,x5,x6,x7,x8) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1","x2", "x3","x4","x5","x6","x7","x8"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = x1+x2+x3+x4+x5+x6+x7+x8 + f2 = x1*x2+x2*x3+x3*x4+x4*x5+x5*x6+x6*x7+x1*x8+x7*x8 + f3 = x1*x2*x3+x2*x3*x4+x3*x4*x5+x4*x5*x6+x5*x6*x7+x1*x2*x8+x1*x7*x8+x6*x7*x8 + f4 = x1*x2*x3*x4+x2*x3*x4*x5+x3*x4*x5*x6+x4*x5*x6*x7+x1*x2*x3*x8+x1*x2*x7*x8+x1*x6*x7*x8+x5*x6*x7*x8 + f5 = x1*x2*x3*x4*x5+x2*x3*x4*x5*x6+x3*x4*x5*x6*x7+x1*x2*x3*x4*x8+x1*x2*x3*x7*x8+x1*x2*x6*x7*x8+x1*x5*x6*x7*x8+x4*x5*x6*x7*x8 + f6 = x1*x2*x3*x4*x5*x6+x2*x3*x4*x5*x6*x7+x1*x2*x3*x4*x5*x8+x1*x2*x3*x4*x7*x8+x1*x2*x3*x6*x7*x8+x1*x2*x5*x6*x7*x8+x1*x4*x5*x6*x7*x8+x3*x4*x5*x6*x7*x8 + f7 = x1*x2*x3*x4*x5*x6*x7+x1*x2*x3*x4*x5*x6*x8+x1*x2*x3*x4*x5*x7*x8+x1*x2*x3*x4*x6*x7*x8+x1*x2*x3*x5*x6*x7*x8+x1*x2*x4*x5*x6*x7*x8+x1*x3*x4*x5*x6*x7*x8+x2*x3*x4*x5*x6*x7*x8 + f8 = x1*x2*x3*x4*x5*x6*x7*x8-1 + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7, f8]) +end + + +function eco6() + dim = 6 + ve = [1, 1, 1, 1, 1, 1] + example = "eco6" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5", "x6"], + ordering = Singular.ordering_M(StartOrd), + ) + + + f1 = x1 + x2 + x3 + x4 + x5 + 1 + f2 = x5 * x6 - 5 + f3 = x1 * x5 * x6 + x4 * x6 - 4 + f4 = x1 * x4 * x6 + x2 * x5 * x6 + x3 * x6 - 3 + f5 = x1 * x3 * x6 + x2 * x4 * x6 + x3 * x5 * x6 + x2 * x6 - 2 + f6 = x1 * x2 * x6 + x2 * x3 * x6 + x3 * x4 * x6 + x4 * x5 * x6 + x1 * x6 - 1 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6]) +end + +function eco7() + dim = 6 + ve = [1, 1, 1, 1, 1, 1] + example = "eco7" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5", "x6"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = x1 + x2 + x3 + x4 + x5 + x6 + 1 + f2 = x6 * x7 - 6 + f3 = x1 * x6 * x7 + x5 * x7 - 5 + f4 = x1 * x5 * x7 + x2 * x6 * x7 + x4 * x7 - 4 + f5 = x1 * x4 * x7 + x2 * x5 * x7 + x3 * x6 * x7 + x3 * x7 - 3 + f6 = x1 * x3 * x7 + x2 * x4 * x7 + x3 * x5 * x7 + x4 * x6 * x7 + x2 * x7 - 2 + f7 = + x1 * x2 * x7 + + x2 * x3 * x7 + + x3 * x4 * x7 + + x4 * x5 * x7 + + x5 * x6 * x7 + + x1 * x7 - 1 + + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7]) +end + +function noon5() + dim = 5 + ve = [1, 1, 1, 1, 1] + example = "noon5" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = + 10 * x1^2 * x5 + 10 * x2^2 * x5 + 10 * x3^2 * x5 + 10 * x4^2 * x5 - + 11 * x5 + 10 + f2 = + 10 * x1^2 * x4 + 10 * x2^2 * x4 + 10 * x3^2 * x4 + 10 * x4 * x5^2 - + 11 * x4 + 10 + f3 = + 10 * x1^2 * x3 + 10 * x2^2 * x3 + 10 * x3 * x4^2 + 10 * x3 * x5^2 - + 11 * x3 + 10 + f4 = + 10 * x1 * x2^2 + 10 * x1 * x3^2 + 10 * x1 * x4^2 + 10 * x1 * x5^2 - + 11 * x1 + 10 + f5 = + 10 * x1^2 * x2 + 10 * x2 * x3^2 + 10 * x2 * x4^2 + 10 * x2 * x5^2 - + 11 * x2 + 10 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5]) +end + +function noon6() + dim = 6 + ve = [1, 1, 1, 1, 1, 1] + example = "noon6" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5", "x6"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = + 10 * x1^2 * x6 + + 10 * x2^2 * x6 + + 10 * x3^2 * x6 + + 10 * x4^2 * x6 + + 10 * x5^2 * x6 - 11 * x6 + 10 + f2 = + 10 * x1^2 * x5 + + 10 * x2^2 * x5 + + 10 * x3^2 * x5 + + 10 * x4^2 * x5 + + 10 * x5 * x6^2 - 11 * x5 + 10 + f3 = + 10 * x1^2 * x4 + + 10 * x2^2 * x4 + + 10 * x3^2 * x4 + + 10 * x4 * x5^2 + + 10 * x4 * x6^2 - 11 * x4 + 10 + f4 = + 10 * x1^2 * x3 + + 10 * x2^2 * x3 + + 10 * x3 * x4^2 + + 10 * x3 * x5^2 + + 10 * x3 * x6^2 - 11 * x3 + 10 + f5 = + 10 * x1 * x2^2 + + 10 * x1 * x3^2 + + 10 * x1 * x4^2 + + 10 * x1 * x5^2 + + 10 * x1 * x6^2 - 11 * x1 + 10 + f6 = + 10 * x1^2 * x2 + + 10 * x2 * x3^2 + + 10 * x2 * x4^2 + + 10 * x2 * x5^2 + + 10 * x2 * x6^2 - 11 * x2 + 10 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6]) +end + +function noon7() + dim = 7 + ve = [1, 1, 1, 1, 1, 1, 1] + example = "noon7" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5, x6, x7) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5", "x6", "x7"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = + 10 * x1^2 * x7 + + 10 * x2^2 * x7 + + 10 * x3^2 * x7 + + 10 * x4^2 * x7 + + 10 * x5^2 * x7 + + 10 * x6^2 * x7 - 11 * x7 + 10 + f2 = + 10 * x1^2 * x6 + + 10 * x2^2 * x6 + + 10 * x3^2 * x6 + + 10 * x4^2 * x6 + + 10 * x5^2 * x6 + + 10 * x6 * x7^2 - 11 * x6 + 10 + f3 = + 10 * x1^2 * x5 + + 10 * x2^2 * x5 + + 10 * x3^2 * x5 + + 10 * x4^2 * x5 + + 10 * x5 * x6^2 + + 10 * x5 * x7^2 - 11 * x5 + 10 + f4 = + 10 * x1^2 * x4 + + 10 * x2^2 * x4 + + 10 * x3^2 * x4 + + 10 * x4 * x5^2 + + 10 * x4 * x6^2 + + 10 * x4 * x7^2 - 11 * x4 + 10 + f5 = + 10 * x1^2 * x3 + + 10 * x2^2 * x3 + + 10 * x3 * x4^2 + + 10 * x3 * x5^2 + + 10 * x3 * x6^2 + + 10 * x3 * x7^2 - 11 * x3 + 10 + f6 = + 10 * x1 * x2^2 + + 10 * x1 * x3^2 + + 10 * x1 * x4^2 + + 10 * x1 * x5^2 + + 10 * x1 * x6^2 + + 10 * x1 * x7^2 - 11 * x1 + 10 + f7 = + 10 * x1^2 * x2 + + 10 * x2 * x3^2 + + 10 * x2 * x4^2 + + 10 * x2 * x5^2 + + 10 * x2 * x6^2 + + 10 * x2 * x7^2 - 11 * x2 + 10 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7]) +end + +function noon8() + dim = 8 + ve = [1, 1, 1, 1, 1, 1, 1, 1] + example = "noon8" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, x5, x6, x7, x8) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = + 10 * x1^2 * x8 + + 10 * x2^2 * x8 + + 10 * x3^2 * x8 + + 10 * x4^2 * x8 + + 10 * x5^2 * x8 + + 10 * x6^2 * x8 + + 10 * x7^2 * x8 - 11 * x8 + 10 + f2 = + 10 * x1^2 * x7 + + 10 * x2^2 * x7 + + 10 * x3^2 * x7 + + 10 * x4^2 * x7 + + 10 * x5^2 * x7 + + 10 * x6^2 * x7 + + 10 * x7 * x8^2 - 11 * x7 + 10 + f3 = + 10 * x1^2 * x6 + + 10 * x2^2 * x6 + + 10 * x3^2 * x6 + + 10 * x4^2 * x6 + + 10 * x5^2 * x6 + + 10 * x6 * x7^2 + + 10 * x6 * x8^2 - 11 * x6 + 10 + f4 = + 10 * x1^2 * x5 + + 10 * x2^2 * x5 + + 10 * x3^2 * x5 + + 10 * x4^2 * x5 + + 10 * x5 * x6^2 + + 10 * x5 * x7^2 + + 10 * x5 * x8^2 - 11 * x5 + 10 + f5 = + 10 * x1^2 * x4 + + 10 * x2^2 * x4 + + 10 * x3^2 * x4 + + 10 * x4 * x5^2 + + 10 * x4 * x6^2 + + 10 * x4 * x7^2 + + 10 * x4 * x8^2 - 11 * x4 + 10 + f6 = + 10 * x1^2 * x3 + + 10 * x2^2 * x3 + + 10 * x3 * x4^2 + + 10 * x3 * x5^2 + + 10 * x3 * x6^2 + + 10 * x3 * x7^2 + + 10 * x3 * x8^2 - 11 * x3 + 10 + f7 = + 10 * x1 * x2^2 + + 10 * x1 * x3^2 + + 10 * x1 * x4^2 + + 10 * x1 * x5^2 + + 10 * x1 * x6^2 + + 10 * x1 * x7^2 + + 10 * x1 * x8^2 - 11 * x1 + 10 + f8 = + 10 * x1^2 * x2 + + 10 * x2 * x3^2 + + 10 * x2 * x4^2 + + 10 * x2 * x5^2 + + 10 * x2 * x6^2 + + 10 * x2 * x7^2 + + 10 * x2 * x8^2 - 11 * x2 + 10 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7, f8]) +end + +function redeco7() + dim = 7 + ve = [1, 1, 1, 1, 1, 1, 1] + example = "redeco7" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, u7, x5, x6) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "u7", "x5", "x6"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = -6 * u7 + x6 + f2 = x1 + x2 + x3 + x4 + x5 + x6 + 1 + f3 = x1 * x6 - 5 * u7 + x5 + f4 = x1 * x5 + x2 * x6 + x4 - 4 * u7 + f5 = x1 * x4 + x2 * x5 + x3 * x6 + x3 - 3 * u7 + f6 = x1 * x3 + x2 * x4 + x3 * x5 + x4 * x6 + x2 - 2 * u7 + f7 = x1 * x2 + x2 * x3 + x3 * x4 + x4 * x5 + x5 * x6 + x1 - u7 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7]) +end + +function redeco8() + dim = 8 + ve = [1, 1, 1, 1, 1, 1, 1, 1] + example = "redeco8" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3, x4, u8, x5, x6, x7) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3", "x4", "u8", "x5", "x6", "x7"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = -7 * u8 + x7 + f2 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + 1 + f3 = x1 * x7 - 6 * u8 + x6 + f4 = x1 * x6 + x2 * x7 + x5 - 5 * u8 + f5 = x1 * x5 + x2 * x6 + x3 * x7 + x4 - 4 * u8 + f6 = x1 * x4 + x2 * x5 + x3 * x6 + x4 * x7 + x3 - 3 * u8 + f7 = x1 * x3 + x2 * x4 + x3 * x5 + x4 * x6 + x5 * x7 + x2 - 2 * u8 + f8 = x1 * x2 + x2 * x3 + x3 * x4 + x4 * x5 + x5 * x6 + x6 * x7 + x1 - u8 + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6, f7, f8]) +end + +function wang91() + dim = 6 + ve = [1, 1, 1, 1, 1, 1] + example = "Wang-91" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x3, x2, x1, x0, b, a) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x3", "x2", "x1", "x0", "b", "a"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = 3 * x2 * x1 * a + 3 * x0^2 + f2 = 3 * x2 * x1 * b + 3 * x3^2 + f3 = 3 * x3 * x1 * b + 3 * x1 * x0 * a + 3 * x2^2 + f4 = 3 * x3 * x2 * b + 3 * x2 * x0 * a + 3 * x1^2 + + return Singular.Ideal(R, [f1, f2, f3, f4]) +end + +function cohn4() + dim = 4 + ve = [1, 1, 1, 1] + example = "cohn4" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x, y, z, t) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x", "y", "z", "t"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = + -x^3 * y^2 + 2 * x^2 * y^2 * z - x^2 * y * z^2 - 144 * x^2 * y^2 - + 207 * x^2 * y * z + + 288 * x * y^2 * z + + 78 * x * y * z^2 + + x * z^3 - 3456 * x^2 * y - 5184 * x * y^2 - 9504 * x * y * z - + 432 * x * z^2 - 248832 * x * y + 62208 * x * z - 2985984 * x + f2 = + y^3 * t^3 - y^2 * z * t^3 + 4 * y^3 * t^2 - 2 * y^2 * z * t^2 + + 72 * y^2 * t^3 + + 71 * y * z * t^3 + + 288 * y^2 * t^2 + + 360 * y * z * t^2 + + 6 * z^2 * t^2 + + 1728 * y * t^3 - 464 * z * t^3 + + 432 * y * z * t + + 8 * z^2 * t + + 6912 * y * t^2 - 4320 * z * t^2 + + 13824 * t^3 + + z^2 - 13824 * z * t + 55296 * t^2 - 13824 * z + f3 = + x^2 * y * t^3 - 2 * x * y^2 * t^3 + y^3 * t^3 + 8 * x^2 * y * t^2 - + 12 * x * y^2 * t^2 + 4 * y^3 * t^2 - 24 * x * y * t^3 + + 24 * y^2 * t^3 + + 20 * x^2 * y * t - 20 * x * y^2 * t - 160 * x * y * t^2 + + 96 * y^2 * t^2 + + 128 * x * t^3 + + 16 * x^2 * y + + 96 * x * y * t + + 2304 * x * t^2 + + 1152 * x * y + + 13824 * x * t + + 27648 * x + f4 = + -x^3 * z * t^2 + x^2 * z^2 * t^2 - 6 * x^3 * z * t + + 4 * x^2 * z^2 * t + + 32 * x^3 * t^2 - 72 * x^2 * z * t^2 - 87 * x * z^2 * t^2 - z^3 * t^2 - + 8 * x^3 * z - 432 * x^2 * z * t - 414 * x * z^2 * t + + 2592 * x * z * t^2 + + 864 * z^2 * t^2 - 1728 * x^2 * z - 20736 * x * z * t + 3456 * z^2 * t - + 186624 * z * t^2 - 124416 * x * z - 1492992 * z * t - 2985984 * z + return Singular.Ideal(R, [f1, f2, f3, f4]) +end + +function oberfr() + dim = 3 + ve = [9, 9, 8] + example = "Oberfranz" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x1, x2, x3) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x1", "x2", "x3"], + ordering = Singular.ordering_M(StartOrd), + ) + f1 = x2^3 + x1 * x2 * x3 + x2^2 * x3 + x1 * x3^3 + f2 = 3 + x1 * x2 + x1^2 * x2 + x2^2 * x3 + + return Singular.Ideal(R, [f1, f2]) +end + +function trinks1() + dim = 6 + ve = [1, 1, 1, 1, 1, 1] + example = "trinks1" + StartOrd = ordering_as_matrix(:degrevlex, dim) + R, (x, y, z, t, u, v) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["x", "y", "z", "t", "u", "v"], + ordering = Singular.ordering_M(StartOrd), + ) + + f1 = 45 * y + 35 * u - 165 * v - 36 + f2 = 36 * y + 25 * z + 40 * t - 27 * u + f3 = 25 * y * u - 165 * v^2 + 15 * x - 18 * z + 30t + f4 = 15 * y * z + 20 * t * u - 9 * x + f5 = -11 * v^3 + x * y + 2 * z * t + f6 = -11 * u * v + 3 * v^2 + 99 * x + + + return Singular.Ideal(R, [f1, f2, f3, f4, f5, f6]) +end + +function ex1() + dim = 4 + ve = [1, 1, 1, 1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + TarOrd = ordering_as_matrix(:lex, dim) + R, (a, b, c, d) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["a", "b", "c", "d"], + ordering = Singular.ordering_M(StartOrd), + ) + S = change_order(R, TarOrd) + + return Singular.Ideal( + R, + [ + d^2 + + 3 * a^2 * d + + 3 * c * d^2 + + 5 * a^3 * c + + 5 * a^2 * c^2 + + 4 * a^2 * c * d, + 3 + 2 * b + 3 * a * b + 3 * c^2 * d + 5 * a^4 + a * b * c^2, + ], + ) +end + +function ex2() + dim = 4 + ve = [1, 1, 1, 1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + TarOrd = ordering_as_matrix(:lex, dim) + R, (a, b, c, d) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["a", "b", "c", "d"], + ordering = Singular.ordering_M(StartOrd), + ) + + S = change_order(R, TarOrd) + return Singular.Ideal( + R, + [ + 2 * a^2 * b + + 3 * a * b^2 + + 3 * b^3 + + 4 * c^3 + + 4 * a * b * d + + c^2 * d + + 2 * b * d^2 + + 2 * d^3 + + 4 * c^2 + + 2 * c * d + + 2 * c, + 2 * a^2 * b + + 5 * a^2 * c + + 2 * b * c^2 + + c^2 * d + + a * c + + 2 * b * d + + 2 * c * d + + d^2 + + a + + 4 * d, + ], + ) +end + +function v3g7d50() + + dim = 3 + ve = [1, 1, 1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + TarOrd = ordering_as_matrix(:lex, dim) + R, (a, b, c) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["a", "b", "c"], + ordering = Singular.ordering_M(StartOrd), + ) + + S = change_order(R, TarOrd) + return Singular.Ideal(R,[5*a^7+4*a^5*b^2+2*a^3*b^4+5*a^2*b^5+4*b^7+5*a^3*b^3*c+3*a^2*b^4*c+4*a^5*c^2+5*a^2*b^3*c^2+5*a^4*c^3+a^3*b*c^3+2*a^2*b^2*c^3+3*a*b^3*c^3+4*a^2*c^5+5*a*b*c^5+b^2*c^5+a^6+3*a^4*b^2+4*a*b^5+3*a^5*c+3*a*b^4*c+5*a^3*b*c^2+2*a^3*c^3+5*a^2*b*c^3+2*a^2*c^4+2*b^2*c^4+a*c^5+2*b*c^5+2*c^6+3*a^3*b^2+a*b^4+5*a*b^3*c+4*a^2*b*c^2+a*b^2*c^2+a^2*c^3+5*a*c^4+4*c^5+2*a^3*b+2*b^4+4*a*b^2*c+3*a*b*c^2+5*b*c^3+2*a*b^2+5*b^2*c+5*a*c^2+5*c^3+3*a^2+a*c+b*c+2*c^2+a+2*c,4*a^6*b+4*a^5*b^2+3*a^4*b^3+4*a^3*b^4+5*a*b^6+4*a^6*c+4*a^4*b^2*c+5*a^4*b*c^2+a^3*b^2*c^2+4*a^2*b^3*c^2+a^4*c^3+5*a^3*b*c^3+4*b^4*c^3+3*a^2*b*c^4+a*b^2*c^4+3*b^3*c^4+4*a^2*c^5+5*a*b*c^5+3*a*c^6+c^7+3*a^5*b+2*a^3*b^3+5*a^2*b^4+3*b^6+a^5*c+5*a^4*b*c+4*a^2*b^3*c+5*b^5*c+a^4*c^2+4*a^3*b*c^2+4*a^2*b^2*c^2+4*a^3*c^3+4*a^2*b*c^3+4*b^3*c^3+3*a^2*c^4+a^5+5*a^2*b^3+4*a*b^4+5*b^5+5*a^4*c+5*a^2*b^2*c+4*b^4*c+4*a^3*c^2+3*a^2*b*c^2+5*b^3*c^2+3*a^2*c^3+4*b^2*c^3+3*c^5+5*a^4+2*a^3*b+b^4+2*a^3*c+3*a^2*b*c+4*a*b^2*c+2*b^3*c+4*a^2*c^2+3*b^2*c^2+2*c^4+3*a^3+3*b^3+a*b*c+2*b^2*c+2*a*c^2+c^3+2*b^2+4*a*c+5*a+5]) +end + +function v4g3() + + dim = 4 + ve = [1, 1, 1, 1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + TarOrd = ordering_as_matrix(:lex, dim) + R, (a, b, c, d) = Singular.PolynomialRing( + Singular.N_ZpField(32003), + ["a", "b", "c", "d"], + ordering = Singular.ordering_M(StartOrd), + ) + + S = change_order(R, TarOrd) + return Singular.Ideal(R,[3*a^3+3*a^2*b+5*a*b^2+2*b^3+a^2*c+2*b^2*c+3*a*c^2+3*a^2*d+5*b^2*d+3*a*c*d+a*d^2+b*d^2+3*c*d^2+5*d^3+4*b^2+3*a*c+5*b*c+3*c^2+4*a*d+2*b*d+c*d+d^2+4*a+4*b+c,3*a^2*b+5*a*b^2+2*a^2*c+5*a*b*c+4*b^2*c+4*a*c^2+5*c^3+5*a^2*d+a*b*d+4*b^2*d+2*a*c*d+3*b*c*d+4*c^2*d+5*b*d^2+4*a*b+3*b^2+2*a*c+5*c^2+5*a*d+2*d^2+4*a+5*c +]) +end diff --git a/src/GroebnerWalk/FractalWalkUtilities.jl b/src/GroebnerWalk/FractalWalkUtilities.jl new file mode 100644 index 000000000..84414492f --- /dev/null +++ b/src/GroebnerWalk/FractalWalkUtilities.jl @@ -0,0 +1,131 @@ +################################################################# +# Procedures of the fractal walk. +# The fractal walk is proposed by Amrhein & Gloor (1998). +################################################################# + +# lifts the Groebner basis G to the Groebner basis w.r.t. in the Fractal Walk like it´s done in Fukuda et. al (2005). +function lift_fractal_walk( + G::Singular.sideal, + H::Singular.sideal, + Rn::Singular.PolyRing, +) + R = base_ring(G) + G.isGB = true + G = Singular.Ideal( + Rn, + [ + change_ring(gen, Rn) - + change_ring(Singular.reduce(change_ring(gen, R), G), Rn) for + gen in Singular.gens(H) + ], + ) + G.isGB = true + return G +end + +# returns ´true´ if all polynomials of the given array are monomials. +function ismonomial(Gw::Vector{spoly{L}}) where {L<:Nemo.RingElem} + for g in Gw + if length(Singular.coefficients(g)) > 1 + return false + end + end + return true +end + +# returns ´true´ if all polynomials of the given array are binomials or less. +function isbinomial(Gw::Vector{spoly{L}}) where {L<:Nemo.RingElem} + for g in Gw + if length(Singular.coefficients(g)) > 2 + return false + end + end + return true +end + +# returns the next t to compute the next weight vector w(t) = w + t * (tw - w) like it´s done in Amrhein & Gloor (1998). This Method is NOT tested sufficiently. +function nextT( + G::Singular.sideal, + w::Array{T,1}, + tw::Array{K,1}, +) where {T<:Number,K<:Number} + if (w == tw) + return [0] + end + tmin = 2 + t = 0 + for g in gens(G) + a = Singular.leading_exponent_vector(g) + d = Singular.exponent_vectors(tail(g)) + for v in d + frac = (dot(w, a) - dot(w, v) + dot(tw, v) - dot(tw, a)) + if frac != 0 + t = (dot(w, a) - dot(w, v)) // frac + end + if t > 0 && t < tmin + tmin = t + end + end + end + if tmin <= 1 + return tmin + else + return [0] + end +end + +# returns the next t to compute the next weight vector w(t) = w + t * (tw - w) like it´s done in Cox, Little & O'Sheao (2005). +function next_weightfr( + G::Singular.sideal, + cweight::Array{T,1}, + tweight::Array{K,1}, +) where {T<:Number,K<:Number} + if (cweight == tweight) + return [0] + end + tmin = 1 + for v in difference_lead_tail(G) + tw = dot(tweight, v) + if tw < 0 + cw = dot(cweight, v) + t = cw // (cw - tw) + if t < tmin + tmin = t + end + end + end + + # BigInt is needed to prevent overflows in the conversion of the weight vectors. + return BigInt(numerator(tmin)) // BigInt(denominator(tmin)) +end + +# returns 'true' if the leading terms of G w.r.t the matrixordering T are the same as the leading terms of G w.r.t the weighted monomial ordering with weight vector of pvecs[p] (pvecs[p-1]) and the matrixordering T. +function inCone( + G::Singular.sideal, + T::Matrix{Int}, + pvecs::Vector{Vector{Int}}, + p::Int, +) + if p == 1 + return true + end + R = change_order(G.base_ring, T) + cvzip = zip( + Singular.gens(G), + initials(R, Singular.gens(G), pvecs[p-1]), + initials(R, Singular.gens(G), pvecs[p]), + ) + for (g, in, in2) in cvzip + if !isequal( + Singular.leading_exponent_vector(change_ring(g, R)), + Singular.leading_exponent_vector(in), + ) || + !isequal( + Singular.leading_exponent_vector(change_ring(g, R)), + Singular.leading_exponent_vector(in2), + ) + return false + end + end + return true +end diff --git a/src/GroebnerWalk/GenericWalkPertubed.jl b/src/GroebnerWalk/GenericWalkPertubed.jl new file mode 100644 index 000000000..06617b948 --- /dev/null +++ b/src/GroebnerWalk/GenericWalkPertubed.jl @@ -0,0 +1,312 @@ +############################################################### +# Not tested version of the Generic Walk with different degree +# of pertubation. This implementation is not verified. +############################################################### + +#= +############################################################### +# Generic Walk with choosable dergree of pertubation. +# This version is not tested yet. +############################################################### + +function pgeneric_walk( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + p::Int, + infoLevel::Int, +) + R = base_ring(G) + Rn = change_order(G.base_ring, T) + v = next_gamma(G, [0], S, T, p) + Lm = [change_ring(Singular.leading_term(g), Rn) for g in gens(G)] + G = Singular.Ideal(Rn, [change_ring(x, Rn) for x in gens(G)]) + + println("generic_walk results") + println("Crossed Cones with facetNormal: ") + while !isempty(v) + global counter = getCounter() + 1 + println(v) + G, Lm = generic_step(G, Lm, v, T, R) + v = next_gamma(G, Lm, v, S, T, p) + end + return Singular.interreduce(G) +end + +function generic_step( + G::Singular.sideal, + Lm::Vector{Singular.spoly{L}}, + v::Vector{Int}, + T::Matrix{Int}, + R::Singular.PolyRing, +) where {L<:Nemo.RingElem} + + Rn = Singular.base_ring(G) + facet_Generators = facet_initials(G, Lm, v) + H = Singular.std( + Singular.Ideal(Rn, facet_Generators), + complete_reduction = true, + ) + H, Lm = lift_generic(G, Lm, H) + G = interreduce(H, Lm) + G = Singular.Ideal(Rn, G) + G.isGB = true + return G, Lm +end +=# + +############################################################### +# Procedures +############################################################### + +#= +function facet_initials( + G::Singular.sideal, + lm::Vector{spoly{L}}, + v::Vector{Int}, +) where {L<:Nemo.RingElem} + Rn = base_ring(G) + initials = Array{Singular.elem_type(Rn),1}(undef, 0) + count = 1 + for g in Singular.gens(G) + inw = Singular.MPolyBuildCtx(Rn) + el = first(Singular.exponent_vectors(lm[count])) + for (e, c) in + zip(Singular.exponent_vectors(g), Singular.coefficients(g)) + if el == e || isParallel(el - e, v) + Singular.push_term!(inw, c, e) + end + end + h = finish(inw) + push!(initials, h) + count += 1 + end + return initials +end + +function difference_lead_tail( + I::Singular.sideal, + Lm::Vector{spoly{L}}, +) where {L<:Nemo.RingElem} + v = Vector{Int}[] + for i = 1:ngens(I) + ltu = Singular.leading_exponent_vector(Lm[i]) + for e in Singular.exponent_vectors(gens(I)[i]) + if ltu != e + push!(v, ltu .- e) + end + end + end + return unique!(v) +end + +function isParallel(u::Vector{Int}, v::Vector{Int}) + count = 1 + x = 0 + for i = 1:length(u) + if u[i] == 0 + if v[count] == 0 + count += +1 + else + return false + end + else + x = v[count] // u[i] + count += 1 + break + end + end + if count > length(v) + return true + end + for i = count:length(v) + @inbounds if v[i] != x * u[i] + return false + end + end + return true +end + +function lift_generic( + G::Singular.sideal, + Lm::Vector{Singular.spoly{L}}, + H::Singular.sideal, +) where {L<:Nemo.RingElem} + Rn = base_ring(G) + Newlm = Array{Singular.elem_type(Rn),1}(undef, 0) + liftPolys = Array{Singular.elem_type(Rn),1}(undef, 0) + for g in Singular.gens(H) + r = modulo(g, gens(G), Lm) + diff = g - r + if diff != 0 + push!(Newlm, Singular.leading_term(g)) + push!(liftPolys, diff) + end + end + return liftPolys, Newlm +end + +function filter_btz(S::Matrix{Int}, V::Vector{Vector{Int}}) + btz = Set{Vector{Int}}() + for v in V + if bigger_than_zero(S, v) + push!(btz, v) + end + end + return btz +end + +function filter_ltz(S::Matrix{Int}, V::Set{Vector{Int}}) + btz = Set{Vector{Int}}() + for v in V + if less_than_zero(S, v) + push!(btz, v) + end + end + return btz +end +function filter_lf( + w::Vector{Int}, + S::Matrix{Int}, + T::Matrix{Int}, + V::Set{Vector{Int}}, +) + btz = Set{Vector{Int}}() + for v in V + if less_facet(w, v, S, T) + push!(btz, v) + end + end + return btz +end + +function next_gamma( + G::Singular.sideal, + Lm::Vector{spoly{L}}, + w::Vector{Int}, + S::Matrix{Int}, + T::Matrix{Int}, + p::Int, +) where {L<:Nemo.RingElem} + V = filter_btz(S, difference_lead_tail(G, Lm)) + V = filter_ltz(T, V) + if (w != [0]) + V = filter_lf(w, S, T, V) + end + if isempty(V) + return V + end + minV = first(V) + for v in V + if less_facet(v, minV, S, T, p) + minV = v + end + end + return minV +end + +function next_gamma( + G::Singular.sideal, + w::Vector{Int}, + S::Matrix{Int}, + T::Matrix{Int}, + p::Int, +) + V = filter_btz(S, difference_lead_tail(G)) + V = filter_ltz(T, V) + if (w != [0]) + V = filter_lf(w, S, T, V) + end + if isempty(V) + return V + end + minV = first(V) + for v in V + if less_facet(v, minV, S, T, p) + minV = v + end + end + return minV +end + +function bigger_than_zero(M::Matrix{Int}, v::Vector{Int}) + nrows, ncols = size(M) + for i = 1:nrows + d = 0 + for j = 1:ncols + @inbounds d += M[i, j] * v[j] + end + if d != 0 + return d > 0 + end + end + return false +end + +function less_than_zero(M::Matrix{Int}, v::Vector{Int}) + nrows, ncols = size(M) + for i = 1:nrows + d = 0 + for j = 1:ncols + @inbounds d += M[i, j] * v[j] + end + if d != 0 + return d < 0 + end + end + return false +end + +function less_facet( + u::Vector{Int}, + v::Vector{Int}, + S::Matrix{Int}, + T::Matrix{Int}, + d::Int, +) + for i = 1:d + for j = 1:d + @inbounds Tuv = dot(T[i, :], u) * dot(S[j, :], v) + @inbounds Tvu = dot(T[i, :], v) * dot(S[j, :], u) + if Tuv != Tvu + return Tuv < Tvu + end + end + end + return false +end + +function modulo( + p::Singular.spoly, + G::Vector{spoly{L}}, + Lm::Vector{spoly{L}}, +) where {L<:Nemo.RingElem} + for i = 1:length(G) + (q, b) = dividesGW(p, Lm[i], parent(p)) + + if b + return modulo(p - (q * G[i]), G, Lm) + end + end + return p +end + +function interreduce( + G::Vector{spoly{L}}, + Lm::Vector{spoly{L}}, +) where {L<:Nemo.RingElem} + Rn = parent(first(G)) + for i = 1:Singular.length(G) + gensrest = Array{Singular.elem_type(Rn),1}(undef, 0) + Lmrest = Array{Singular.elem_type(Rn),1}(undef, 0) + for j = 1:length(G) + if i != j + push!(gensrest, G[j]) + push!(Lmrest, Lm[j]) + end + end + G[i] = modulo(G[i], gensrest, Lmrest) + end + return G +end +=# diff --git a/src/GroebnerWalk/GenericWalkUtilities.jl b/src/GroebnerWalk/GenericWalkUtilities.jl new file mode 100644 index 000000000..418cd5835 --- /dev/null +++ b/src/GroebnerWalk/GenericWalkUtilities.jl @@ -0,0 +1,313 @@ +################################################################# +# Procedures of the generic walk. +# The generic walk is proposed by Fukuda, Lauritzen & Thomas (2005). +################################################################# + +# returns the initials of the polynomials w.r.t. the vector v. +function facet_initials( + G::Vector{Singular.spoly{L}}, + lm::Vector{spoly{L}}, + v::Vector{Int}, +) where {L<:Nemo.RingElem} + Rn = parent(first(G)) + initials = Array{Singular.elem_type(Rn),1}(undef, 0) + count = 1 + for g in G + inw = Singular.MPolyBuildCtx(Rn) + el = first(Singular.exponent_vectors(lm[count])) + for (e, c) in + zip(Singular.exponent_vectors(g), Singular.coefficients(g)) + if el == e || isparallel(el - e, v) + Singular.push_term!(inw, c, e) + end + end + h = finish(inw) + push!(initials, h) + count += 1 + end + return initials +end + +# returns the differences of the exponent vectors of the leading terms and the polynomials of the generators of I. +function difference_lead_tail( + G::Vector{spoly{L}}, + Lm::Vector{spoly{L}}, +) where {L<:Nemo.RingElem} + v = Vector{Int}[] + for i = 1:length(G) + ltu = Singular.leading_exponent_vector(Lm[i]) + for e in Singular.exponent_vectors(G[i]) + if ltu != e + push!(v, ltu .- e) + end + end + end + return unique!(v) +end + +# returns true if the vector u is parallel to the vector v. +function isparallel(u::Vector{Int}, v::Vector{Int}) + count = 1 + x = 0 + for i = 1:length(u) + if u[i] == 0 + if v[count] == 0 + count += +1 + else + return false + end + else + x = v[count] // u[i] + count += 1 + break + end + end + if count > length(v) + return true + end + for i = count:length(v) + @inbounds if v[i] != x * u[i] + return false + end + end + return true +end + +# performs the lifting in the generic Walk like it´s proposed by Fukuda et al. (2005). +function lift_generic( + G::Vector{spoly{L}}, + Lm::Vector{Singular.spoly{L}}, + H::Singular.sideal, +) where {L<:Nemo.RingElem} + Rn = parent(first(G)) + Newlm = Array{Singular.elem_type(Rn),1}(undef, 0) + liftPolys = Array{Singular.elem_type(Rn),1}(undef, 0) + for g in Singular.gens(H) + push!(Newlm, Singular.leading_term(g)) + push!(liftPolys, g - reduce_walk(g, G, Lm)) + end + return liftPolys, Newlm +end + +# returns all v \in V if v<0 w.r.t. the ordering represented by T and v>0 w.r.t the ordering represented by S. +function filter_by_ordering(S::Matrix{Int},T::Matrix{Int}, V::Vector{Vector{Int}}) + btz = Set{Vector{Int}}() + for v in V + if less_than_zero(T, v) && bigger_than_zero(S, v) + push!(btz, v) + end + end + return btz +end + +# returns all v \in V if w0 w.r.t. the ordering M. +function bigger_than_zero(M::Matrix{Int}, v::Vector{Int}) + nrows, ncols = size(M) + for i = 1:nrows + d = 0 + for j = 1:ncols + @inbounds d += M[i, j] * v[j] + end + if d != 0 + return d > 0 + end + end + return false +end + +# tests if v<0 w.r.t. the ordering M. +function less_than_zero(M::Matrix{Int}, v::Vector{Int}) + nrows, ncols = size(M) + for i = 1:nrows + d = 0 + for j = 1:ncols + @inbounds d += M[i, j] * v[j] + end + if d != 0 + return d < 0 + end + end + return false +end + +# tests if u standard_walk(x, S, T, infoLevel) + elseif walktype == :generic + walk = (x) -> generic_walk(x, S, T, infoLevel) + elseif walktype == :pertubed + walk = (x) -> pertubed_walk(x, S, T, p, infoLevel) + elseif walktype == :fractal + walk = (x) -> fractal_walk(x, S, T, infoLevel) + elseif walktype == :fractal_start_order + walk = (x) -> fractal_walk_start_order(x, S, T, infoLevel) + elseif walktype == :fractal_lex + walk = (x) -> fractal_walk_lex(x, S, T, infoLevel) + elseif walktype == :fractal_look_ahead + walk = (x) -> fractal_walk_look_ahead(x, S, T, infoLevel) + elseif walktype == :tran + walk = (x) -> tran_walk(x, S, T, infoLevel) + elseif walktype == :fractal_combined + walk = (x) -> fractal_walk_combined(x, S, T, infoLevel) + end + + delete_counter() + !check_order_M(S, T, G) && throw( + error( + "The matrices representing the monomial order have to be nxn-matrices with full rank.", + ), + ) + + R = base_ring(G) + Gb = walk(Singular.Ideal(R, [R(x) for x in gens(G)])) + + if infoLevel >= 1 + println("Cones crossed: ", delete_counter()) + end + + S = change_order(R, T) + return Singular.Ideal(S, [change_ring(gen, S) for gen in gens(Gb)]) +end + +########################################### +# Counter for the steps in the Fractal Walk. +########################################### +counter = 0 +function delete_counter() + global counter + temp = counter + counter = 0 + return temp +end +function getcounter() + global counter + return counter +end +function raise_counter() + global counter = getcounter() + 1 +end +############################################################### +# Implementation of the standard walk. +############################################################### + +function standard_walk( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + if infoLevel >= 1 + println("standard_walk results") + println("Crossed Cones in: ") + end + + Gb = standard_walk(G, S, T, S[1, :], T[1, :], infoLevel) + + return Gb +end + +function standard_walk( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + currweight::Vector{Int}, + tarweight::Vector{Int}, + infoLevel::Int, +) + while true + G = standard_step(G, currweight, T) + if infoLevel >= 1 + println(currweight) + if infoLevel == 2 + println(G) + end + end + raise_counter() + if currweight == tarweight + return G + else + currweight = next_weight(G, currweight, tarweight) + end + end +end + +############################################################### +# The standard step is used for the strategies standard and pertubed. +############################################################### + +function standard_step(G::Singular.sideal, w::Vector{Int}, T::Matrix{Int}) + R = base_ring(G) + Rn = 0 + Gw = 0 + + #check if no entry of w is bigger than Int32. If it´s bigger multiply it by 0.1 and round. + if !checkInt32(w) + Gw = initials(R, gens(G), w) + w, b = truncw(G, w, Gw) + if !b + throw( + error( + "Some entries of the intermediate weight-vector $w are bigger than int32", + ), + ) + end + Rn = change_order(R, w, T) + Gw = [change_ring(x, Rn) for x in Gw] + else + Rn = change_order(R, w, T) + Gw = initials(Rn, gens(G), w) + end + + H = Singular.std(Singular.Ideal(Rn, Gw), complete_reduction = true) + #H = liftGW2(G, R, Gw, H, Rn) + H = lift(G, R, H, Rn) + return interreduce_walk(H) +end + +############################################################### +# Generic-version of the Groebner Walk. +############################################################### + +function generic_walk( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + R = base_ring(G) + Rn = change_order(G.base_ring, T) + Lm = [change_ring(Singular.leading_term(g), Rn) for g in gens(G)] + G = [change_ring(x, Rn) for x in gens(G)] + v = next_gamma(G, Lm, [0], S, T) + + if infoLevel >= 1 + println("generic_walk results") + println("Crossed Cones with: ") + end + while !isempty(v) + G, Lm = generic_step(G, Lm, v, Rn) + raise_counter() + + if infoLevel >= 1 + println(v) + if infoLevel == 2 + println(G) + end + end + + v = next_gamma(G, Lm, v, S, T) + end + G = Singular.Ideal(Rn, G) + G.isGB = true + return G +end + +function generic_step( + G::Vector{Singular.spoly{L}}, + Lm::Vector{Singular.spoly{L}}, + v::Vector{Int}, + Rn::Singular.PolyRing, +) where {L<:Nemo.RingElem} + facet_Generators = facet_initials(G, Lm, v) + H = Singular.std( + Singular.Ideal(Rn, facet_Generators), + complete_reduction = true, + ) + H, Lm = lift_generic(G, Lm, H) + G = interreduce(H, Lm) + return G, Lm +end + +############################################################### +#Pertubed-version of the Groebner Walk. +############################################################### + +function pertubed_walk( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + p::Int, + infoLevel::Int, +) + if infoLevel >= 1 + println("pertubed_walk results") + println("Crossed Cones in: ") + end + + currweight = pertubed_vector(G, S, p) + + while true + tarweight = pertubed_vector(G, T, p) + Tn = add_weight_vector(tarweight, T) + G = standard_walk(G, S, Tn, currweight, tarweight, infoLevel) + if same_cone(G, T) + return G + else + p = p - 1 + currweight = tarweight + S = Tn + end + end +end + +############################################################### +# The Fractal Walk +############################################################### + +########################################## +# global weightvectors +########################################## +pTargetWeights = [] +pStartWeights = [] +firstStepMode = false + +############################################################### +# Combined version of the extensions of the Fractal Walk. +# This version +# - checks if the starting weight vector represents the monomial order and pertubes it if necessary. +# - analyses the Groebner basis Gw of the initialforms and uses the Buchberger-algorithm if the generators of Gw are binomial or less. +# - skips a step in top level in the last step. +# - checks if an entry of an intermediate weight vector is bigger than int32. In case of that the Buchberger-Algorithm is used to compute the Groebner basis of the ideal of the initialforms. +############################################################### +function fractal_walk_combined( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + if infoLevel >= 1 + println("fractal_walk_combined results") + println("Crossed Cones in: ") + end + + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(Singular.base_ring(G))] + return fractal_walk_combined(G, S, T, S[1, :], pTargetWeights, 1, infoLevel) +end + +function fractal_walk_combined( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + currweight::Vector{Int}, + pTargetWeights::Vector{Vector{Int}}, + p::Int, + infoLevel::Int, +) + R = Singular.base_ring(G) + G.isGB = true + + # Handling the weight of the start order. + if (p == 1) + if !ismonomial(initials(R, Singular.gens(G), currweight)) + global pStartWeights = [pertubed_vector(G, S, i) for i = 1:nvars(R)] + global firstStepMode = true + end + end + if firstStepMode + w = pStartWeights[p] + else + w = currweight + end + + # main loop + while true + t = next_weightfr(G, w, pTargetWeights[p]) + + # Handling the final step in the current depth. + # Next_weightfr may return 0 if the target vector does not lie in the cone of T while G already defines the Groebner basis w.r.t. T. + # -> Checking if G is already a Groebner basis w.r.t. T solves this problem and reduces computational effort since next_weightfr returns 1 in the last step on every local path. + if t == 1 && p != 1 + if same_cone(G, T) + if infoLevel >= 1 + println("depth $p: in cone ", currweight, ".") + end + + # Check if a target weight of pTargetWeights[p] and pTargetWeights[p-1] lies in the wrong cone. + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + elseif t == [0] # The Groebner basis w.r.t. the target weight and T is already computed. + if inCone(G, T, pTargetWeights, p) + if infoLevel >= 1 + println("depth $p: in cone ", pTargetWeights[p], ".") + end + return G + end + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + continue + end + + # skip a step for target monomial order lex. + if t == 1 && p == 1 + if infoLevel >= 1 + println("depth $p: recursive call in ", pTargetWeights[p]) + end + return fractal_walk_combined( + G, + S, + T, + w, + pTargetWeights, + p + 1, + infoLevel, + ) + else + w = w + t * (pTargetWeights[p] - w) + w = convert_bounding_vector(w) + Gw = initials(R, gens(G), w) + + # handling the current weight with regards to Int32-entries. If an entry of w is bigger than Int32 use the Buchberger-algorithm. + if !checkInt32(w) + w, b = truncw(G, w, Gw) + if !b + Rn = change_order(R, T) + w = T[1, :] + G = Singular.std( + Singular.Ideal( + Rn, + [change_ring(x, Rn) for x in gens(G)], + ), + complete_reduction = true, + ) + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println( + "depth $p: not in cone ", + pTargetWeights[p], + ".", + ) + end + end + return G + end + end + Rn = change_order(R, w, T) + + # converting the Groebner basis + if (p == Singular.nvars(R) || isbinomial(Gw)) + H = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in Gw]), + complete_reduction = true, + ) + if infoLevel >= 1 + println("depth $p: conversion in ", w, ".") + end + raise_counter() + else + if infoLevel >= 1 + println("depth $p: recursive call in $w.") + end + H = fractal_walk_combined( + Singular.Ideal(R, Gw), + S, + T, + deepcopy(currweight), + pTargetWeights, + p + 1, + infoLevel, + ) + global firstStepMode = false + end + end + #H = liftGW2(G, R, Gw, H, Rn) + H = lift_fractal_walk(G, H, Rn) + G = interreduce_walk(H) + R = Rn + currweight = w + end +end + +############################################################### +# Plain version of the Fractal Walk. +# This version checks if an entry of an intermediate weight vector is bigger than int32. In case of that the Buchberger-Algorithm is used to compute the Groebner basis of the ideal of the initialforms. +############################################################### + +function fractal_walk( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + if infoLevel >= 1 + println("FractalWalk_standard results") + println("Crossed Cones in: ") + end + + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(base_ring(G))] + return fractal_recursiv(G, S, T, S[1, :], pTargetWeights, 1, infoLevel) +end + +function fractal_recursiv( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + currweight::Vector{Int}, + pTargetWeights::Vector{Vector{Int}}, + p::Int, + infoLevel::Int, +) + R = base_ring(G) + G.isGB = true + w = currweight + + while true + t = next_weightfr(G, w, pTargetWeights[p]) + + # Handling the final step in the current depth. + # Next_weightfr may return 0 if the target vector does not lie in the cone of T while G already defines the Groebner basis w.r.t. T. + # -> Checking if G is already a Groebner basis w.r.t. T solves this problem and reduces computational effort since next_weightfr returns 1 in the last step on every local path. if t == 1 && p != 1 + if t == 1 && p != 1 + if same_cone(G, T) + if infoLevel >= 1 + println("depth $p: in cone ", currweight, ".") + end + + # Check if a target weight of pTargetWeights[p] and pTargetWeights[p-1] lies in the wrong cone. + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + elseif t == [0] # The Groebner basis w.r.t. the target weight and T is already computed. + if inCone(G, T, pTargetWeights, p) + if infoLevel >= 1 + println("depth $p: in cone ",pTargetWeights[p], ".") + end + return G + end + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + continue + end + + w = w + t * (pTargetWeights[p] - w) + w = convert_bounding_vector(w) + Gw = initials(R, Singular.gens(G), w) + + # Handling the current weight with regards to Int32-entries. If an entry of w is bigger than Int32 use the Buchberger-algorithm. + if !checkInt32(w) + w, b = truncw(G, w, Gw) + if !b + Rn = change_order(R, T) + w = T[1, :] + G = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in gens(G)]), + complete_reduction = true, + ) + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + end + + # Converting the Groebner basis + Rn = change_order(R, w, T) + if p == nvars(R) + H = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in Gw]), + complete_reduction = true, + ) + if infoLevel >= 1 + println("depth $p: conversion in ", w, ".") + end + raise_counter() + else + if infoLevel >= 1 + println("depth $p: recursive call in $w.") + end + H = fractal_recursiv( + Singular.Ideal(R, Gw), + S, + T, + deepcopy(currweight), + pTargetWeights, + p + 1, + infoLevel, + ) + end + #H = liftGW2(G, R, Gw, H, Rn) + H = lift_fractal_walk(G, H, Rn) + G = interreduce_walk(H) + R = Rn + currweight = w + end +end + +############################################################### +# Extends the plain Fractal Walk by checking the start order. +############################################################### + +function fractal_walk_start_order( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + if infoLevel >= 1 + println("fractal_walk_withStartorder results") + println("Crossed Cones in: ") + end + + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(Singular.base_ring(G))] + return fractal_walk_recursiv_startorder( + G, + S, + T, + S[1, :], + pTargetWeights, + 1, + infoLevel, + ) +end + +function fractal_walk_recursiv_startorder( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + currweight::Vector{Int}, + pTargetWeights::Vector{Vector{Int}}, + p::Int, + infoLevel::Int, +) + R = Singular.base_ring(G) + G.isGB = true + + # Handling the starting weight. + if (p == 1) + if !ismonomial(initials(R, Singular.gens(G), currweight)) + global pStartWeights = [pertubed_vector(G, S, i) for i = 1:nvars(R)] + global firstStepMode = true + end + end + if firstStepMode + w = pStartWeights[p] + else + w = currweight + end + + while true + t = next_weightfr(G, w, pTargetWeights[p]) + + # Handling the final step in the current depth. + # Next_weightfr may return 0 if the target vector does not lie in the cone of T while G already defines the Groebner basis w.r.t. T. + # -> Checking if G is already a Groebner basis w.r.t. T solves this problem and reduces computational effort since next_weightfr returns 1 in the last step on every local path. if t == 1 && p != 1 + if t == 1 && p != 1 + if same_cone(G, T) + if infoLevel >= 1 + println("depth $p: in cone ", currweight, ".") + end + + # Check if a target weight of pTargetWeights[p] and pTargetWeights[p-1] lies in the wrong cone. + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + elseif t == [0] # The Groebner basis w.r.t. the target weight and T is already computed. + if inCone(G, T, pTargetWeights, p) + if infoLevel >= 1 + println("depth $p: in cone ",pTargetWeights[p], ".") + end + return G + end + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + continue + end + + w = w + t * (pTargetWeights[p] - w) + w = convert_bounding_vector(w) + Gw = initials(R, gens(G), w) + + # Handling the current weight with regards to Int32-entries. If an entry of w is bigger than Int32 use the Buchberger-algorithm. + if !checkInt32(w) + w, b = truncw(G, w, Gw) + if !b + Rn = change_order(R, T) + w = T[1, :] + G = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in gens(G)]), + complete_reduction = true, + ) + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + end + Rn = change_order(R, w, T) + + # Converting the Groebner basis + if p == Singular.nvars(R) + H = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in Gw]), + complete_reduction = true, + ) + if infoLevel >= 1 + println("depth $p: conversion in ", w, ".") + end + raise_counter() + else + if infoLevel >= 1 + println("depth $p: recursive call in $w.") + end + H = fractal_walk_recursiv_startorder( + Singular.Ideal(R, Gw), + S, + T, + deepcopy(currweight), + pTargetWeights, + p + 1, + infoLevel, + ) + global firstStepMode = false + end + #H = liftGW2(G, R, Gw, H, Rn) + H = lift_fractal_walk(G, H, Rn) + G = interreduce_walk(H) + R = Rn + currweight = w + end +end + +############################################################### +# Plain version of the Fractal Walk in case of a lexicographic target order. +############################################################### + +function fractal_walk_lex( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + if infoLevel >= 1 + println("fractal_walk_lex results") + println("Crossed Cones in: ") + end + + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(base_ring(G))] + return fractal_walk_recursive_lex( + G, + S, + T, + S[1, :], + pTargetWeights, + 1, + infoLevel, + ) +end + +function fractal_walk_recursive_lex( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + currweight::Vector{Int}, + pTargetWeights::Vector{Vector{Int}}, + p::Int, + infoLevel::Int, +) + R = Singular.base_ring(G) + G.isGB = true + w = currweight + + while true + t = next_weightfr(G, w, pTargetWeights[p]) + + # Handling the final step in the current depth. + # Next_weightfr may return 0 if the target vector does not lie in the cone of T while G already defines the Groebner basis w.r.t. T. + # -> Checking if G is already a Groebner basis w.r.t. T solves this problem and reduces computational effort since next_weightfr returns 1 in the last step on every local path. if t == 1 && p != 1 + if t == 1 && p != 1 + if same_cone(G, T) + if infoLevel >= 1 + println("depth $p: in cone ", currweight, ".") + end + + # Check if a target weight of pTargetWeights[p] and pTargetWeights[p-1] lies in the wrong cone. + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + elseif t == [0] # The Groebner basis w.r.t. the target weight and T is already computed. + if inCone(G, T, pTargetWeights, p) + if infoLevel >= 1 + println("depth $p: in cone ",pTargetWeights[p], ".") + end + return G + end + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + continue + end + + # Skipping a step in lex. + if t == 1 && p == 1 + return fractal_walk_recursive_lex( + G, + S, + T, + w, + pTargetWeights, + p + 1, + infoLevel, + ) + else + w = w + t * (pTargetWeights[p] - w) + w = convert_bounding_vector(w) + Gw = initials(R, Singular.gens(G), w) + + # Handling the current weight with regards to Int32-entries. If an entry of w is bigger than Int32 use the Buchberger-algorithm. + if !checkInt32(w) + w, b = truncw(G, w, Gw) + if !b + Rn = change_order(R, T) + w = T[1, :] + G = Singular.std( + Singular.Ideal( + Rn, + [change_ring(x, Rn) for x in gens(G)], + ), + complete_reduction = true, + ) + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println( + "depth $p: not in cone ", + pTargetWeights[p], + ".", + ) + end + end + return G + end + end + Rn = change_order(R, w, T) + + # Converting the Groebner basis + if p == Singular.nvars(R) + H = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in Gw]), + complete_reduction = true, + ) + if infoLevel >= 1 + println("depth $p: conversion in ", w, ".") + end + raise_counter() + else + if infoLevel >= 1 + println("depth $p: recursive call in $w.") + end + H = fractal_walk_recursive_lex( + Singular.Ideal(R, Gw), + S, + T, + deepcopy(currweight), + pTargetWeights, + p + 1, + infoLevel, + ) + global firstStepMode = false + end + end + #H = liftGW2(G, R, Gw, H, Rn) + H = lift_fractal_walk(G, H, Rn) + G = interreduce_walk(H) + R = Rn + currweight = w + end +end + +############################################################### +# Plain version of the Fractal Walk with look-ahead extension. +############################################################### + +function fractal_walk_look_ahead( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + if infoLevel >= 1 + println("fractal_walk_look_ahead results") + println("Crossed Cones in: ") + end + + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(base_ring(G))] + return fractal_walk_look_ahead_recursiv( + G, + S, + T, + S[1, :], + pTargetWeights, + 1, + infoLevel, + ) + return Gb +end + +function fractal_walk_look_ahead_recursiv( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + currweight::Vector{Int}, + pTargetWeights::Vector{Vector{Int}}, + p::Int, + infoLevel, +) + R = Singular.base_ring(G) + G.isGB = true + w = currweight + + while true + t = next_weightfr(G, w, pTargetWeights[p]) + + # Handling the final step in the current depth. + # Next_weightfr may return 0 if the target vector does not lie in the cone of T while G already defines the Groebner basis w.r.t. T. + # -> Checking if G is already a Groebner basis w.r.t. T solves this problem and reduces computational effort since next_weightfr returns 1 in the last step on every local path. if t == 1 && p != 1 + if t == 1 && p != 1 + if same_cone(G, T) + if infoLevel >= 1 + println("depth $p: in cone ", currweight, ".") + end + + # Check if a target weight of pTargetWeights[p] and pTargetWeights[p-1] lies in the wrong cone. + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + elseif t == [0] # The Groebner basis w.r.t. the target weight and T is already computed. + if inCone(G, T, pTargetWeights, p) + if infoLevel >= 1 + println("depth $p: in cone ",pTargetWeights[p], ".") + end + return G + end + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + continue + end + + w = w + t * (pTargetWeights[p] - w) + w = convert_bounding_vector(w) + Gw = initials(R, Singular.gens(G), w) + + # Handling the current weight with regards to Int32-entries. If an entry of w is bigger than Int32 use the Buchberger-algorithm. + if !checkInt32(w) + w, b = truncw(G, w, Gw) + if !b + Rn = change_order(R, T) + w = T[1, :] + G = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in gens(G)]), + complete_reduction = true, + ) + if !inCone(G, T, pTargetWeights, p) + global pTargetWeights = + [pertubed_vector(G, T, i) for i = 1:nvars(R)] + if infoLevel >= 1 + println("depth $p: not in cone ",pTargetWeights[p], ".") + end + end + return G + end + end + Rn = change_order(R, w, T) + + # Converting the Groebner basis + if (p == Singular.nvars(R) || isbinomial(Gw)) + H = Singular.std( + Singular.Ideal(Rn, [change_ring(x, Rn) for x in Gw]), + complete_reduction = true, + ) + if infoLevel >= 1 + println("depth $p: conversion in ", w, ".") + end + raise_counter() + else + if infoLevel >= 1 + println("depth $p: recursive call in $w.") + end + H = fractal_walk_look_ahead_recursiv( + Singular.Ideal(R, Gw), + S, + T, + deepcopy(currweight), + pTargetWeights, + p + 1, + infoLevel, + ) + end + #H = liftGW2(G, R, Gw, H, Rn) + H = lift_fractal_walk(G, H, Rn) + G = interreduce_walk(H) + R = Rn + currweight = w + end +end + +############################################################### +# Tran´s version of the Groebner Walk. +# Returns the intermediate Groebner basis if an entry of an intermediate weight vector is bigger than int32. +############################################################### + +function tran_walk( + G::Singular.sideal, + S::Matrix{Int}, + T::Matrix{Int}, + infoLevel::Int, +) + if infoLevel >= 1 + println("tran_walk results") + println("Crossed Cones in: ") + end + + currweight = S[1, :] + tarweight = T[1, :] + R = base_ring(G) + if !ismonomial(initials(R, Singular.gens(G), currweight)) + currweight = pertubed_vector(G, S, nvars(R)) + end + + while true + w = next_weight(G, currweight, tarweight) + + # return the Groebner basis if an entry of w is bigger than int32. + if !checkInt32(w) + w, b = truncw(G, w, initials(R, gens(G), w)) + if !b + return G + end + end + Rn = change_order(R, w, T) + if w == tarweight + if same_cone(G, T) + if infoLevel >= 1 + println("Cones crossed: ", counter) + end + return G + elseif inSeveralCones(initials(base_ring(G), gens(G), tarweight)) + tarweight = representation_vector(G, T) + continue + end + end + G = standard_step_without_int32_check(G, w, T) + if infoLevel >= 1 + println(w) + if infoLevel == 2 + println(G) + end + end + R = Rn + currweight = w + raise_counter() + end +end + +############################################################### +# Standard step without checking of the entries of a given weight vector. +############################################################### + +function standard_step_without_int32_check( + G::Singular.sideal, + w::Vector{Int}, + T::Matrix{Int}, +) + R = base_ring(G) + Rn = change_order(R, w, T) + Gw = initials(Rn, gens(G), w) + H = Singular.std(Singular.Ideal(Rn, Gw), complete_reduction = true) + #H = liftGW2(G, R, Gw, H, Rn) + H = lift(G, R, H, Rn) + return interreduce_walk(H) +end diff --git a/src/GroebnerWalk/GroebnerWalkUtilities.jl b/src/GroebnerWalk/GroebnerWalkUtilities.jl new file mode 100644 index 000000000..ba7e226c9 --- /dev/null +++ b/src/GroebnerWalk/GroebnerWalkUtilities.jl @@ -0,0 +1,350 @@ +############################################################### +# Several Procedures for the Groebner Walk +############################################################### + +# returns the next intermediate weight vector. +function next_weight( + G::Singular.sideal, + cw::Vector{Int}, + tw::Vector{Int}, +) where {K<:Number} + tmin = 1 + for v in difference_lead_tail(G) + tdotw = dot(tw, v) + if tdotw < 0 + cdotw = dot(cw, v) + t = cdotw // (cdotw - tdotw) + if t < tmin + tmin = t + end + end + end + return convert_bounding_vector(cw + BigInt(numerator(tmin))//BigInt(denominator(tmin)) * (tw - cw)) +end + +# multiplys every entry of the given weight w with 0.1 as long as it stays on the same halfspace as w. +function truncw( + G::Singular.sideal, + w::Vector{Int}, + inw::Vector{L}, +) where {L<:Nemo.RingElem} + while !checkInt32(w) + R = base_ring(G) + for i = 1:length(w) + w[i] = round(w[i] * 0.10) + end + w = convert_bounding_vector(w) + if inw != initials(R, gens(G), w) + # initials are different - return unrounded weight + return w, false + end + end + # converted to Vector w of the same face + return w, true +end + +# returns the initialform of G w.r.t. the given weight vector. +function initials( + R::Singular.PolyRing, + G::Vector{L}, + w::Vector{Int}, +) where {L<:Nemo.RingElem} + inits = spoly{elem_type(base_ring(R))}[] + indexw = Tuple{Vector{Int},elem_type(base_ring(R))}[] + for g in G + empty!(indexw) + maxw = 0 + eczip = zip(Singular.exponent_vectors(g), Singular.coefficients(g)) + for (e, c) in eczip + tmpw = dot(w, e) + if maxw == tmpw + push!(indexw, (e, c)) + elseif maxw < tmpw + empty!(indexw) + push!(indexw, (e, c)) + maxw = tmpw + end + end + inw = MPolyBuildCtx(R) + for (e, c) in indexw + Singular.push_term!(inw, c, e) + end + h = finish(inw) + push!(inits, h) + end + return inits +end + +function difference_lead_tail(I::Singular.sideal) + v = Vector{Int}[] + for g in gens(I) + ltu = Singular.leading_exponent_vector(g) + for e in Singular.exponent_vectors(tail(g)) + push!(v, ltu .- e) + end + end + return unique!(v) +end + +# computes a p-pertubed vector from the matrix M. +function pertubed_vector(G::Singular.sideal, M::Matrix{Int}, p::Integer) + m = Int[] + n = size(M, 1) + for i = 1:p + max = M[i, 1] + for j = 1:n + temp = abs(M[i, j]) + if temp > max + max = temp + end + end + push!(m, max) + end + msum = 0 + for i = 2:p + msum += m[i] + end + maxdeg = 0 + for g in gens(G) + td = deg(g, n) + if (td > maxdeg) + maxdeg = td + end + end + e = maxdeg * msum + 1 + w = M[1, :] * e^(p - 1) + for i = 2:p + w += e^(p - i) * M[i, :] + end + return convert_bounding_vector(w) +end + +# returns 'true' if the leading terms of G w.r.t the matrixorder T are the same as the leading terms of G w.r.t the weighted monomial order with weight vector t and matrix T. +function inCone(G::Singular.sideal, T::Matrix{Int}, t::Vector{Int}) + R = change_order(G.base_ring, T) + I = Singular.Ideal(R, [change_ring(x, R) for x in gens(G)]) + cvzip = zip(Singular.gens(I), initials(R, Singular.gens(I), t)) + for (g, ing) in cvzip + if !isequal(Singular.leading_exponent_vector(g), Singular.leading_exponent_vector(ing)) + return false + end + end + return true +end + +# returns 'true' if the leading terms of G w.r.t the matrixordering T are the same as the leading terms of G w.r.t the weighted monomial order with weight vector t and the matrix order T. +function inCone(G::Singular.sideal, t::Vector{Int}) + cvzip = zip(Singular.gens(G), initials(base_ring(G), Singular.gens(G), t)) + for (g, ing) in cvzip + if !isequal(Singular.leading_exponent_vector(g), Singular.leading_exponent_vector(ing)) + return false + end + end + return true +end + +# returns 'true' if the leading terms of G w.r.t the matrixordering T are the same as the leading terms of G with the current ordering. +function same_cone(G::Singular.sideal, T::Matrix{Int}) + R = change_order(G.base_ring, T) + for g in gens(G) + if !isequal(Singular.leading_exponent_vector(change_ring(g,R)), Singular.leading_exponent_vector(g)) + return false + end + end + return true +end + +# lifts the Groebner basis G to the Groebner basis w.r.t. the Ring Rn like it´s presented in Fukuda et al. (2005). +function lift( + G::Singular.sideal, + R::Singular.PolyRing, + H::Singular.sideal, + Rn::Singular.PolyRing, +) + G.isGB = true + G = Singular.Ideal( + Rn, + [ + gen - change_ring(Singular.reduce(change_ring(gen, R), G), Rn) + for gen in gens(H) + ], + ) + G.isGB = true + return G +end + +# lifts the Groebner basis G to the Groebner basis w.r.t. the Ring Rn like it´s done in Collart et al. (1997). +function liftGW2( + G::Singular.sideal, + R::Singular.PolyRing, + inG::Vector{spoly{L}}, + H::Singular.sideal, + Rn::Singular.PolyRing, +) where {L<:Nemo.RingElem} + + gH = collect(gens(H)) + gG = collect(gens(G)) + inG = [change_ring(x, R) for x in inG] + for i = 1:length(gH) + q = division_algorithm(change_ring(gH[i], R), inG, R) + gH[i] = Rn(0) + for j = 1:length(inG) + gH[i] = + change_ring(gH[i], Rn) + + change_ring(q[j], Rn) * change_ring(gG[j], Rn) + end + end + G = Singular.Ideal(Rn, [x for x in gH]) + G.isGB = true + return G +end + +# divisionalgorithm that returns q with q_1*f_1 + ... + q_s *f_s=p. +function division_algorithm( + p::spoly{L}, + f::Vector{spoly{L}}, + R::Singular.PolyRing, +) where {L<:Nemo.RingElem} + q = Array{Singular.elem_type(R),1}(undef, length(f)) + for i = 1:length(f) + q[i] = R(0) + end + while !isequal(p, R(0)) + i = 1 + div = false + while (div == false && i <= length(f)) + b, m = divides(leading_term(p), leading_term(f[i])) + if b + q[i] = q[i] + m + p = p - (m * f[i]) + div = true + else + i = i + 1 + end + end + if div == false + p = p - leading_term(p) + end + end + return q +end + +# converts a vector wtemp by dividing the entries with gcd(wtemp). +function convert_bounding_vector(wtemp::Vector{T}) where {T<:Rational{BigInt}} + w = Vector{Int}() + g = gcd(wtemp) + for i = 1:length(wtemp) + push!(w, float(divexact(wtemp[i], g))) + end + return w +end + +# converts a vector wtemp by dividing the entries with gcd(wtemp). +function convert_bounding_vector(wtemp::Vector{T}) where {T<:Number} + w = Vector{Int}() + g = gcd(wtemp) + for i = 1:length(wtemp) + push!(w, float(divexact(wtemp[i], g))) + end + return w +end + +# returns a copy of the PolynomialRing I, equipped with the ordering a(cw)*ordering_M(T). +function change_order( + R::Singular.PolyRing, + cw::Array{L,1}, + T::Matrix{Int}, +) where {L<:Number,K<:Number} + G = Singular.gens(R) + Gstrich = string.(G) + s = size(T) + if s[1] == s[2] + S, H = Singular.PolynomialRing( + R.base_ring, + Gstrich, + ordering = Singular.ordering_a(cw) * Singular.ordering_M(T), + cached = false, + ) + elseif s[1] - s[2] == 1 + S, H = Singular.PolynomialRing( + R.base_ring, + Gstrich, + ordering = Singular.ordering_a(cw) * + Singular.ordering_a(T[1, :]) * + Singular.ordering_M(T[2:end, :]), + cached = false, + ) + else + S, H = Singular.PolynomialRing( + R.base_ring, + Gstrich, + ordering = Singular.ordering_a(cw) * + Singular.ordering_a(T[1, :]) * + Singular.ordering_a(T[2, :]) * + Singular.ordering_M(T[3:end, :]), + cached = false, + ) + end + return S +end + +# returns a copy of the PolynomialRing I, equipped with the ordering a(w)*a(t)*ordering_M(T). +function change_order( + R::Singular.PolyRing, + w::Vector{Int}, + t::Vector{Int}, + T::Matrix{Int}, +) where {} + G = Singular.gens(R) + Gstrich = string.(G) + s = size(T) + S, H = Singular.PolynomialRing( + R.base_ring, + Gstrich, + ordering = Singular.ordering_a(w) * + Singular.ordering_a(t) * + Singular.ordering_M(T), + cached = false, + ) + return S +end + +# returns a copy of the PolynomialRing I, equipped with the ordering ordering_M(T). +function change_order( + R::Singular.PolyRing, + M::Matrix{Int}, +) where {T<:Number,K<:Number} + G = Singular.gens(R) + Gstrich = string.(G) + S, H = Singular.PolynomialRing( + R.base_ring, + Gstrich, + ordering = Singular.ordering_M(M), + cached = false, + ) + return S +end + +# recreates the polynomials p equipped with ring R. +function change_ring(p::Singular.spoly, R::Singular.PolyRing) + cvzip = zip(Singular.coefficients(p), Singular.exponent_vectors(p)) + M = MPolyBuildCtx(R) + for (c, v) in cvzip + Singular.push_term!(M, c, v) + end + return finish(M) +end + +# interreduces the Groebner basis G. +function interreduce_walk(G::Singular.sideal) where {L<:Nemo.RingElem} + Rn = base_ring(G) + Generator = collect(gens(G)) + I = 0 + for i = 1:ngens(G) + I = Singular.Ideal(Rn, Generator[1:end.!=i]) + I.isGB = true + Generator[i] = reduce(Generator[i], I) + end + G = Singular.Ideal(Rn, Generator) + return G +end diff --git a/src/GroebnerWalk/Helper.jl b/src/GroebnerWalk/Helper.jl new file mode 100644 index 000000000..bf33c5914 --- /dev/null +++ b/src/GroebnerWalk/Helper.jl @@ -0,0 +1,165 @@ +############################################# +# unspecific help functions +############################################# + +function ident_matrix(n::Int) + M = zeros(Int, n, n) + for i = 1:n + M[i, i] = 1 + end + return M +end + +function anti_diagonal_matrix(n::Int) + M = zeros(Int, n, n) + for i = 1:n + M[i, n+1-i] = -1 + end + return M +end + +# Singular.isequal depends on order of generators +function equalitytest(G::Singular.sideal, K::Singular.sideal) + if Singular.ngens(G) != Singular.ngens(K) + return false + end + generators = Singular.gens(G) + count = 0 + for gen in generators + for r in Singular.gens(K) + if gen * leading_coefficient(r) - r * leading_coefficient(gen) == 0 + count += 1 + break + end + end + end + if count == Singular.ngens(G) + return true + end + return false +end + +function dot(v::Vector{Int}, w::Vector{Int}) + sum = 0 + for i = 1:length(v) + @inbounds sum += v[i] * w[i] + end + return sum +end + +function ordering_as_matrix(w::Vector{Int}, ord::Symbol) + if length(w) > 2 + if ord == :lex + return [ + w' + ident_matrix(length(w))[1:length(w)-1, :] + ] + end + if ord == :deglex + return [ + w' + ones(Int, length(w))' + ident_matrix(length(w))[1:length(w)-2, :] + ] + end + if ord == :degrevlex + return [ + w' + ones(Int, length(w))' + anti_diagonal_matrix(length(w))[1:length(w)-2, :] + ] + end + if ord == :revlex + return [ + w' + anti_diagonal_matrix(length(w))[1:length(w)-1, :] + ] + end + else + error("not implemented") + end +end + +function change_weight_vector(w::Vector{Int}, M::Matrix{Int}) + return [ + w' + M[2:length(w), :] + ] +end + +function insert_weight_vector(w::Vector{Int}, M::Matrix{Int}) + return [ + w' + M[1:length(w)-1, :] + ] +end + +function add_weight_vector(w::Vector{Int}, M::Matrix{Int}) + return [ + w' + M + ] +end + +function ordering_as_matrix(ord::Symbol, nvars::Int) + if ord == :lex + return ident_matrix(nvars) + end + if ord == :deglex + return [ + ones(Int, nvars)' + ident_matrix(nvars)[1:nvars-1, :] + ] + end + if ord == :degrevlex + return [ + ones(Int, nvars)' + anti_diagonal_matrix(nvars)[1:nvars-1, :] + ] + end + if ord == :revlex + return [ + w' + anti_diagonal_matrix(length(w))[1:length(w)-1, :] + ] + else + error("not implemented") + end +end + +function deg(p::Singular.spoly, n::Int) + max = 0 + for mon in Singular.monomials(p) + ev = Singular.exponent_vectors(mon) + sum = 0 + for e in ev + for i = 1:n + sum += e[i] + end + if (max < sum) + max = sum + end + end + end + return max +end + +function check_order_M(S::Matrix{Int}, T::Matrix{Int}, G::Singular.sideal) + (nrows, ncols) = size(T) + return ( + (nrows, ncols) == size(S) && + nrows == ncols && + nrows == nvars(base_ring(G)) && + rank(T) == nvars(base_ring(G)) && + rank(S) == nvars(base_ring(G)) + ) +end + +function checkInt32(w::Vector{Int}) + for i = 1:length(w) + if tryparse(Int32, string(w[i])) == nothing + return false + end + end + return true +end diff --git a/src/GroebnerWalk/README.md b/src/GroebnerWalk/README.md new file mode 100644 index 000000000..2d228a88b --- /dev/null +++ b/src/GroebnerWalk/README.md @@ -0,0 +1,42 @@ +# GroebnerWalk + +### Usage + +The following function calls perform different versions of the Gröbner Walk by computing the Gröbner basis w.r.t. the monomial order ```:degrevlex``` and converting it to the Groebner basis w.r.t. the monomial order ```:lex```. + +```julia +using Singular +using Nemo +include("Examples.jl") +include("GroebnerWalk.jl") + +id = katsura4() +groebnerwalk(id, :degrevlex, :lex ,:standard) +groebnerwalk(id, :degrevlex, :lex ,:pertubed, 2) +groebnerwalk(id, :degrevlex, :lex ,:fractal_combined) +groebnerwalk(id, :degrevlex, :lex ,:generic) +``` + +The following function call performs the standard version of the Gröbner Walk by converting the given Groebner basis w.r.t. the monomial order ``` :degrevlex``` to the Groebner basis w.r.t. the monomial order ```:lex```. + +```julia +using Singular +using Nemo +include("Examples.jl") +include("GroebnerWalk.jl") + +id = katsura4() +R = base_ring(id) +I = Singular.std(id, complete_reduction = true) + +groebnerwalk(I, ordering_as_matrix(:degrevlex, nvars(R)), ordering_as_matrix(:lex, nvars(R)), :standard) + +``` + +### Tests + +The tests can be started with **runtests.jl**. + +### Tracing the Walk + +For informations about the conversions-steps use the parameter ```infoLevel```. diff --git a/src/GroebnerWalk/Testideals.jl b/src/GroebnerWalk/Testideals.jl new file mode 100644 index 000000000..a96ef049b --- /dev/null +++ b/src/GroebnerWalk/Testideals.jl @@ -0,0 +1,171 @@ +include("GroebnerWalk.jl") +include("Examples.jl") +using Test + +@testset "Groebnerwalks" begin + @testset "Testing Groebnerwalks" begin + let id = katsura5() + R = base_ring(id) + dim = nvars(R) + TarOrd = ordering_as_matrix(:lex, dim) + S = change_order(R, TarOrd) + I = Singular.std(id, complete_reduction = true) + + ideals = [] + for i = 2:nvars(S)-2 + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :pertubed, i)) + end + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :standard)) + #push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :tran)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_start_order)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_look_ahead)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_lex)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_combined)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :generic)) + + s = Singular.std( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + complete_reduction = true, + ) + + for id in ideals + @test equalitytest( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + s, + ) + end + end + + let id = katsura4() + R = base_ring(id) + dim = nvars(R) + StartOrd = ordering_as_matrix([1,3,1,3,1],:lex) + TarOrd = ordering_as_matrix([1,0,0,0,2],:lex) + S = change_order(R, TarOrd) + + I = Singular.std(Singular.Ideal(change_order(R,StartOrd), [change_ring(x, change_order(R,StartOrd)) for x in gens(id)]), complete_reduction = true) + + ideals = [] + for i = 1:nvars(S)-1 + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :pertubed, i)) + end + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :standard)) + #push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :tran)) + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :fractal)) + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :fractal_start_order)) + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :fractal_look_ahead)) + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :fractal_lex)) + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :fractal_combined)) + push!(ideals, groebnerwalk(I, ordering_as_matrix([1,3,1,3,1],:lex), ordering_as_matrix([1,0,0,0,2],:lex), :generic)) + + s = Singular.std( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + complete_reduction = true, + ) + + for id in ideals + @test equalitytest( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + s, + ) + end + end + + let id = cyclic4() + R = base_ring(id) + dim = nvars(R) + TarOrd = ordering_as_matrix(:lex, dim) + S = change_order(R, TarOrd) + + I = Singular.std(id, complete_reduction = true) + + ideals = [] + for i = 2:nvars(S)-1 + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :pertubed, i)) + end + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :standard)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :tran)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_start_order)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_look_ahead)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_lex)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_combined)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :generic)) + + s = Singular.std( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + complete_reduction = true, + ) + + + for id in ideals + @test equalitytest( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + s, + ) + end + end + + let id = ex2() + R = base_ring(id) + dim = nvars(R) + TarOrd = ordering_as_matrix(:lex, dim) + S = change_order(R, TarOrd) + I = Singular.std(id, complete_reduction = true) + + ideals = [] + for i = 2:nvars(S) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :pertubed, i)) + end + push!(ideals, groebnerwalk(I, :degrevlex, :lex, :standard)) + #push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :tran)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_look_ahead)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_start_order)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_lex)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_combined)) + push!(ideals, groebnerwalk(I, :degrevlex, :lex, :generic)) + + s = Singular.std( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + complete_reduction = true, + ) + + for id in ideals + @test equalitytest( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + s, + ) + end + end + + + id = redeco7() + R = base_ring(id) + dim = nvars(R) + TarOrd = ordering_as_matrix(:lex, dim) + S = change_order(R, TarOrd) + I = Singular.std(id, complete_reduction = true) + + ideals = [] + for i = 2:nvars(S)-3 + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :pertubed, i)) + end + push!(ideals, groebnerwalk(I, :degrevlex, :lex, :standard)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :fractal_combined)) + push!(ideals, groebnerwalk(I, ordering_as_matrix(:degrevlex, dim), ordering_as_matrix(:lex, dim), :generic)) + + s = Singular.std( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + complete_reduction = true, + ) + + for id in ideals + @test equalitytest( + Singular.Ideal(S, [change_ring(x, S) for x in Singular.gens(id)]), + s, + ) + end + end +end diff --git a/src/GroebnerWalk/TranWalkUtilities.jl b/src/GroebnerWalk/TranWalkUtilities.jl new file mode 100644 index 000000000..e2ebd1600 --- /dev/null +++ b/src/GroebnerWalk/TranWalkUtilities.jl @@ -0,0 +1,43 @@ +# computes the representation of the matrixorder defined by T. +function representation_vector(G::Singular.sideal, T::Matrix{Int}) + n = size(T)[1] + M = 0 + for i = 1:n + for j = 1:n + temp = T[i, j] + if M < temp + M = temp + end + end + end + d0 = 0 + for g in Singular.gens(G) + temp = deg(g, n) + if d0 < temp + d0 = temp + end + end + d = M * (2 * d0^2 + (n + 1) * d0) + w = d^(n - 1) * T[1, :] + for i = 2:n + w = w + d^(n - i) * T[i, :] + end + return w +end + +# checks if Gw is an initialform of an ideal corresponding to a face of the Groebner fan with dimension less than n-1. +function inSeveralCones(Gw::Vector{spoly{L}}) where {L<:Nemo.RingElem} + counter = 0 + for g in Gw + if size(collect(Singular.coefficients(g)))[1] > 2 + return true + end + if size(collect(Singular.coefficients(g)))[1] == 2 + counter = counter + 1 + end + end + if counter > 1 + return true + end + return false +end diff --git a/src/GroebnerWalk/UnitTests.jl b/src/GroebnerWalk/UnitTests.jl new file mode 100644 index 000000000..18bd100af --- /dev/null +++ b/src/GroebnerWalk/UnitTests.jl @@ -0,0 +1,207 @@ +using Test +include("GroebnerWalk.jl") +include("Examples.jl") + +@testset "UnitTests" begin + @testset "Testing GroebnerwalkUtilitys" begin + + R, (x, y, z) = Singular.PolynomialRing( + Singular.QQ, + ["x", "y", "z"], + ordering = Singular.ordering_M(ordering_as_matrix(:degrevlex, 3)), + ) + + f1 = 3 * x^2 + 16 * x^2 * z + 14 * x^2 * y^3 + f2 = y^3 * z + 17 * x^2 * z^2 + 7 * x^2 * y^2 * z^2 + 13 * x^3 * z^2 + I = Singular.Ideal(R, [f1, f2]) + sol = [14 * x^2 * y^3, y^3 * z + 7 * x^2 * y^2 * z^2] + @test initials(R, gens(I), [1, 3, 1]) == sol + + @test difference_lead_tail(I) == + [[0, 3, -1], [0, 3, 0], [-1, 2, 0], [2, -1, 1], [0, 2, 0]] + + F = [ + 13 * x^3 * z^2, + 14 * x^2 * y^3, + 98 * x * y^5 * z^2, + y^7 * z + x^2 * z^3, + 14 * x * y^6 * z, + ] + g = y^7 * z + x^2 * z^3 + 28 * x^2 * y^4 + q = Array{Singular.elem_type(R),1}(undef, 5) + q[1] = R(0) + q[2] = R(2 * y) + q[3] = R(0) + q[4] = R(1) + q[5] = R(0) + @test division_algorithm(g, F, R) == q + + J = Singular.Ideal(R, [f2, f1]) + f1 = 4 * x^2 + 16 * x^2 * z + 14 * x^2 * y^3 + f2 = y^3 * z + 17 * x^2 * z^2 + 7 * x^2 * y^2 * z^2 + 13 * x^3 * z^2 + K = Singular.Ideal(R, [f1, f2]) + @test equalitytest(I, J) == true + @test equalitytest(I, K) == false + + @test deg(f1, 3) == 5 + + id = trinks1() + R = base_ring(id) + dim = nvars(R) + ve = ones(Int, dim) + StartOrd = ordering_as_matrix(:degrevlex, dim) + TarOrd = ordering_as_matrix(:lex, dim) + I = Singular.std(id, complete_reduction = true) + @test pertubed_vector(I, StartOrd, 1) == [1, 1, 1, 1, 1, 1] + @test pertubed_vector(I, StartOrd, 2) == [4, 4, 4, 4, 4, 3] + @test pertubed_vector(I, StartOrd, 3) == [49, 49, 49, 49, 48, 42] + @test pertubed_vector(I, StartOrd, 4) == + [1000, 1000, 1000, 999, 990, 900] + @test pertubed_vector(I, TarOrd, 1) == [1, 0, 0, 0, 0, 0] + @test pertubed_vector(I, TarOrd, 2) == [4, 1, 0, 0, 0, 0] + @test pertubed_vector(I, TarOrd, 3) == [49, 7, 1, 0, 0, 0] + @test pertubed_vector(I, TarOrd, 4) == [1000, 100, 10, 1, 0, 0] + + @test inCone(I, [1000, 1000, 1000, 999, 990, 900]) == true + @test inCone(I, [100, 1000, 1000, 999, 990, 900]) == false + + @test dot([1, 2, 3, 4], [2, 2, 2, 2]) == 20 + + end + + @testset "Testing FraktalWalkUtilitys" begin + R, (x, y, z) = Singular.PolynomialRing( + Singular.QQ, + ["x", "y", "z"], + ordering = Singular.ordering_M(ordering_as_matrix(:degrevlex, 3)), + ) + + f1 = 3 * x^2 + f2 = y^3 * z + 17 * x^2 * z^2 + I = Singular.Ideal(R, [f1, f2]) + f1 = 3 * x^2 + f2 = y^3 * z + J = Singular.Ideal(R, [f1, f2]) + f1 = x^3 * y^2 + z^2 + y^2 + f2 = y^3 + K = Singular.Ideal(R, [f1, f2]) + + @test ismonomial(gens(I)) == false + @test ismonomial(gens(J)) == true + @test isbinomial(gens(I)) == true + @test isbinomial(gens(J)) == true + @test isbinomial(gens(K)) == false + @test ismonomial(gens(K)) == false + end + + @testset "Testing GenericWalkUtilitys" begin + R, (x, y) = Singular.PolynomialRing( + Singular.QQ, + ["x", "y"], + ordering = Singular.ordering_M(ordering_as_matrix(:degrevlex, 2)), + ) + + f1 = x^2 - y^3 + f2 = x^3 - y^2 - x + I = Singular.Ideal(R, [f1, f2]) + G = Singular.std(I, complete_reduction = true) + + @test next_gamma( + gens(G), + [Singular.leading_term(g) for g in gens(G)], + [0], + ordering_as_matrix(:degrevlex, 2), + ordering_as_matrix(:lex, 2), + ) == [-2, 3] + + Rn, (x, y) = Singular.PolynomialRing( + Singular.QQ, + ["x", "y"], + ordering = Singular.ordering_M(ordering_as_matrix(:lex, 2)), + ) + g = [Rn(y^3 - x^2), Rn(x^3)] + @test facet_initials( + [change_ring(x, Rn) for x in gens(G)], + [change_ring(Singular.leading_term(g), Rn) for g in gens(G)], + [-2, 3], + ) == g + + @test difference_lead_tail( + gens(G), + [Singular.leading_term(g) for g in gens(G)], + ) == [[-2, 3], [3, -2], [2, 0]] + + @test isparallel([1, 2], [1, 4]) == false + @test isparallel([1, 2], [2, 4]) == true + @test isparallel([-1, 0], [-2, 1]) == false + @test isparallel([-1, 0], [2, 0]) == true + + @test less_facet( + [-2, 3], + [-1, 4], + ordering_as_matrix(:degrevlex, 2), + ordering_as_matrix(:lex, 2), + ) == true + @test less_facet( + [-1, 7], + [-1, 4], + ordering_as_matrix(:degrevlex, 2), + ordering_as_matrix(:lex, 2), + ) == false + + dim = 5 + ve = [1, 1, 1, 1, 1] + StartOrd = ordering_as_matrix(:degrevlex, dim) + TarOrd = ordering_as_matrix(:lex, dim) + R, (a, b, c, d, e) = Singular.PolynomialRing( + Singular.QQ, + ["a", "b", "c", "d", "e"], + ordering = Singular.ordering_M(StartOrd), + ) + S = change_order(R, TarOrd) + J = Singular.Ideal( + R, + [ + b + 3 * b^3 + 2 * b * c * e + 5 * b * d * e, + 4 + b^2 + 4 * b * c + 5 * b^3 + c * d * e, + d * e + 5 * b^2 * e, + ], + ) + I = Singular.std(J, complete_reduction = true) + f1 = R( + a^3 + + a^2 + + b^5 * a^3 * c^9 + + e^3 + + b^2 * a^2 * c^4 + + d^3 + + e^3 + + b^2 * d^2 * e^7, + ) + f2 = R(a^3 * b^3) + f3 = R(a^2 * b^4 + c^2 + a^3 * 4 + a * e^3) + f4 = R(0) + + @test (reduce(f3, I)) == + reduce_walk(f3, gens(I), [Singular.leading_term(g) for g in gens(I)]) + @test (reduce(f1, I)) == + reduce_walk(f1, gens(I), [Singular.leading_term(g) for g in gens(I)]) + @test (reduce(f2, I)) == + reduce_walk(f2, gens(I), [Singular.leading_term(g) for g in gens(I)]) + @test (reduce(f4, I)) == + reduce_walk(f4, gens(I), [Singular.leading_term(g) for g in gens(I)]) + J = Singular.std(J) + + @test equalitytest( + Singular.std(J, complete_reduction = true), + Singular.Ideal( + R, + interreduce( + collect(gens(J)), + [Singular.leading_term(g) for g in gens(J)], + ), + ), + ) + + end +end diff --git a/src/GroebnerWalk/runtests.jl b/src/GroebnerWalk/runtests.jl new file mode 100644 index 000000000..1ec2efc0d --- /dev/null +++ b/src/GroebnerWalk/runtests.jl @@ -0,0 +1,8 @@ +using Singular +using Nemo +using Test + +@testset "GroebnerWalk" begin + include("Testideals.jl") + include("UnitTests.jl") +end