segmentation fault primary decomposition over number field

[wdecker]

julia> T, t = PolynomialRing(QQ, "t")
(Univariate Polynomial Ring in t over Rational Field, t)

julia> K, a = NumberField(t^2 + 1, "a")
(Number field over Rational Field with defining polynomial t^2 + 1, a)

julia> R, (y, z) = PolynomialRing(K, ["y", "z"])
(Multivariate Polynomial Ring in y, z over K, AbstractAlgebra.Generic.MPoly{nf_elem}[y, z])

julia> p = z^2+1
z^2 + 1

julia> q = z^3+2
z^3 + 2

julia> I = ideal(R, [p*q^2, y-z^2])
ideal generated by: z^8 + z^6 + 4*z^5 + 4*z^3 + 4*z^2 + 4, y - z^2

julia> L = primary_decomposition(I)

signal (11): Segmentation fault: 11
in expression starting at REPL[9]:1
unknown function (ip: 0x0)
Allocations: 73532342 (Pool: 73482814; Big: 49528); GC: 50
zsh: segmentation fault  julia

[thofma]

Seems to be a Singular issue. Here is how to reproduce it using only Singular.jl:

julia> using Singular

julia> T, t = Singular.Nemo.PolynomialRing(Singular.Nemo.QQ, "t")
(Univariate Polynomial Ring in t over Rational Field, t)

julia> K, a = Singular.Nemo.NumberField(t^2 + 1, "a")
(Number field over Rational Field with defining polynomial t^2 + 1, a)

julia> R, (y, z) = Singular.PolynomialRing(K, ["y", "z"])
(Singular Polynomial Ring (Coeffs(17)),(y,z),(dp(2),C), spoly{Singular.n_unknown{Nemo.nf_elem}}[y, z])

julia> p = z^2+1
z^2 + 1

julia> q = z^3+2
z^3 + 2

julia> I = Singular.Ideal(R, p*q^2, y-z^2)
Singular Ideal over Singular Polynomial Ring (Coeffs(17)),(y,z),(dp(2),C) with generators (z^8 + z^6 + 4*z^5 + 4*z^3 + 4*z^2 + 4, -z^2 + y)

julia> Singular.LibPrimdec.primdecGTZ(R, I)

signal (11): Speicherzugriffsfehler
in expression starting at REPL[16]:1
unknown function (ip: (nil))
Allocations: 25972479 (Pool: 25965025; Big: 7454); GC: 28
[1]    30623 segmentation fault (core dumped)  julia

@tthsqe12 Any idea?

On a related note, do you know why we don't create a Singular number field when doing computations with number fields @tthsqe12? Why is singular_ring(::AnticNumberField) not creating a Singular number field?

[wdecker]

In Singular itself you get:

> LIB "primdec.lib";
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/primdec.lib (4.0.1.2,Jan_2015)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/ring.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/absfact.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/triang.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/matrix.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/nctools.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/inout.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/random.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/poly.lib (4.0.0.0,Jun_2013)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/elim.lib (4.0.0.1,Jan_2014)
// ** loaded /Applications/Singular.app/Contents/bin/../share/singular/LIB/general.lib (4.0.0.1,Jan_2014)
> ring  r = (0,a),(x,y,z),lp;
> minpoly = a2+1;
> poly  p = z2+1;
> poly  q = z3+2;
> ideal I = p*q^2,y-z2;
> primdecGTZ(I);
[1]:
   [1]:
      _[1]=-y+z2
      _[2]=y3+4*yz+4
   [2]:
      _[1]=-y+z2
      _[2]=yz+2
      _[3]=y2+2*z
[2]:
   [1]:
      _[1]=z+(a)
      _[2]=y+1
   [2]:
      _[1]=z+(a)
      _[2]=y+1
[3]:
   [1]:
      _[1]=z+(-a)
      _[2]=y+1
   [2]:
      _[1]=z+(-a)
      _[2]=y+1

[tthsqe12]

Well, of course it crashes in singular. Singular is going to get the n_unknown{nf_elem} and you can imagine what happened from there.

[thofma]

Isn't there the generic fields functionality?

[tthsqe12]

Singular's library procedures are not coded with generic field functionality in mind.

[thofma]

Makes sense. There is no primary decomposition for generic fields without more primitives.

[wdecker]

S o what can we do?

[thofma]

I discussed the fix with @tthsqe12. He will take care.

[thofma]

Do you already have code for this @joschmitt? Going back and forth between Singular and Nemo number fields?

[joschmitt]

Not yet. The plan is to start with this on Wednesday (assuming nobody else is faster or insists on doing it themselves).

[hannes14]

The coefficients in Singular have an array of function pointers to define their behaviour.
primdecGTZ needs factorization, so it calls factory for that.
First action is converting Singular polynomials to factory polynomials, it uses convSingNFactoryN
for conversion of the coeffcients. If this field is not set (i.e. is ndConvSingNFactoryN), an error should occur ("not implemented").

[tthsqe12]

@joschmitt I insist on doing this myself in singular.jl tmr (Tuesday). It is a low-level operation that should try to be "efficient".

[joschmitt]

Sure.

[tthsqe12]

So the fixed version probably outputs an output isomorphic to the singular output:

julia> T, t = PolynomialRing(QQ, "t");

julia> K, a = NumberField(t^2 + 1, "a");

julia> R, (y, z) = PolynomialRing(K, ["y", "z"]);

julia> p = z^2+1;

julia> q = z^3+2;

julia> I = ideal(R, [p*q^2, y-z^2]);

julia> L = primary_decomposition(I)
3-element Vector{Tuple{MPolyIdeal{AbstractAlgebra.Generic.MPoly{nf_elem}}, MPolyIdeal{AbstractAlgebra.Generic.MPoly{nf_elem}}}}:
 (ideal generated by: -y + z^2, y^3 + 4*y*z + 4, ideal generated by: -y + z^2, y*z + 2, y^2 + 2*z)
 (ideal generated by: z - a, y + 1, ideal generated by: z - a, y + 1)
 (ideal generated by: z + a, y + 1, ideal generated by: z + a, y + 1)