segmentation fault computing discriminants

[martin....@gmx.net]

The following sage code dies on me, after computing for quite some time:

┌────────────────────────────────────────────────────────────────────┐
│ Sage Version 6.1.1, Release Date: 2014-02-04                       │
│ Type "notebook()" for the browser-based notebook interface.        │
│ Type "help()" for help.                                            │
└────────────────────────────────────────────────────────────────────┘
sage: PR.<b,t1,t2,x1,y1,x2,y2> = QQ[]
sage: PR2.<mu> = PR[]
sage: E1 = diagonal_matrix(PR, [1, b^2, -b^2])
sage: RotScale = matrix(PR, [[1, -t1, 0], [t1, 1, 0], [0, 0, 1]])
sage: E2 = diagonal_matrix(PR, [1, 1, 1+t1^2])*RotScale.transpose()*E1*RotScale
sage: Transl = matrix(PR, [[1, 0, -x1], [0, 1, -y1], [0, 0, 1]])
sage: E3 = Transl.transpose()*E2*Transl
sage: E4 = E3.subs(t1=t2, x1=x2, y1=y2)
sage: Transl = matrix(PR, [[1, 0, 0], [0, 1, 1+b], [0, 0, 1]])
sage: E2 = Transl.transpose()*diagonal_matrix(PR, [b^2, 1, -b^2])*Transl
sage: eqs = [det(mu*Ei + Ej).discriminant() for Ei, Ej in
....:        [(E1, E3), (E1, E4), (E2, E3), (E3, E4)]]
Segmentation fault

Attaching a debugger, I got thousands of stack frames from within the python interpreter:
symtable_visit_expr from Python-2.7.6/Python/symtable.c:1191 seems to be calling itself till the call stack overflows.

The above is from Sage on Gentoo, but the Sage 6.1.1 official release for Mac segfaults as well.

I'm a bit surprised both at the segfault and at the fact that this takes so long. After all, I'm just computing the discriminant of a bunch of cubic polynomials. Wikipedia has a formula for this, and using that formula, things complete in no time at all. Should discriminants for polynomials of small degree be special-cased to avoid such issues? Or perhaps the discriminant should be computed first for a generic polynomial (i.e. one variable for every non-zero coefficient of the input), and the actual coefficients of the input substituted into that in a second step.

I guess I'll be using that hardcoded formula for cubic discriminants for now, but it would be nice if such things could work out of the box.

[John Cremona]

I certainly agree that this should work. I found that if I set
M=mu*E1+E3 I could compute det(M).discriminant() immediately, and
similarly for 2 of the other 3.
So there is something funny happening, as the basic operation of
computing the discriminant of a cubic with coefficients in your
7-variable polynomial ring is ok; it's the (E3,E4) case which is
hanging. When I substituted the coefficients of det(mu*E3+E4) into
the generic discriminant forumla the reult was an expression 3759965
characters long. The first three are only [14360, 14360, 50720] so
the last one is more complicated,

Looking at the code used, it uses the resultant formula which in turn
evaluates a determinant. I agree with you that for small degrees it
would be better (almost certainly in a lot of cases) be better to
substitute into the generic formula.


John Cremona

> --
> You received this message because you are subscribed to the Google Groups
> "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-support...@googlegroups.com.
> To post to this group, send email to sage-s...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-support.
> For more options, visit https://groups.google.com/d/optout.

[luisfe]

On Wednesday, March 26, 2014 10:34:35 AM UTC+1, John Cremona wrote:

Looking at the code used, it uses the resultant formula which in turn
evaluates a determinant.  I agree with you that for small degrees it
would be better (almost certainly in a lot of cases) be better to
substitute into the generic formula.

Even if it is not the best code around, it should not segfault. I will try the latest beta...

[martin....@gmx.net]

On Wednesday, March 26, 2014 10:34:35 AM UTC+1, John Cremona wrote:

Looking at the code used, it uses the resultant formula which in turn
evaluates a determinant.  I agree with you that for small degrees it
would be better (almost certainly in a lot of cases) be better to
substitute into the generic formula.

Here is some code which does that:

def generic_discriminant(p):
    """
    Compute discriminant of univariate polynomial.
   
    The discriminant is computed in two steps. First all non-zero
    coefficients of the polynomial are replaced with variables, and
    the discriminant is computed using these variables. Then the
    original coefficients are substituted into the result of this
    computation. This should be faster in many cases, particularly
    with big coefficients like complicated polynomials. In those
    cases, the matrix computation on simple generators is a lot
    faster, and the final substitution is cheaper than carrying all
    that information along at all times.
    """
    pr = p.parent()
    if not pr.ngens() == 1:
        raise ValueError("Must be univariate")
    ex = p.exponents()
    pr1 = PolynomialRing(ZZ, "x", len(ex))
    pr2.<y> = pr1[]
    py = sum(v*y^e for v, e in zip(pr1.gens(), ex))
    d = py.discriminant()
    return d(p.coefficients())

Should I file a trac ticket for this? If so, should I try to determine the cases when to use it, or should I simply get the function in there and ask someone else to take care of the required connections to existing code?

What is the correct ring to choose here? I thought about using pr.base_ring() instead of ZZ, but that would mean that in the case of stacked polynomial rings, I'd be doing the computation in a more difficult setup than needed. It seems that sage will automatically coerce things correctly: if the input is a polynomial over a polynomial ring over QQ or GF(p), then so is the output. But I'm not sure whether that will always be the case, or whether there might be cases where conversion from ZZ to something else fails to do the right thing.

[John Cremona]

You should definitely open a trac ticket. Don't expect that if you
post some code as text on a ticket that someone else will magically do
the right thing though! If you have not before contributed to Sage
you can make this your chance to learn how that works, and people will
help.

I would separate out two things here: one is a generic computation of
discriminants of univariate polynomials of degree d, the result being
a polynomial in d+1 variables (the generic coefficients) over ZZ (that
is correct -- there's a unique map from ZZ to any other ring). This
should be part of the current invariant_theory module, though when I
just looked at that it only appeared to implement invariants for
binary forms of degrees 2 and 4 (not 3) as well as some ternary stuff.
As a stop-gap, a global function generic_univariate_discriminant(deg)
whose output is a polynomial in ZZ[a0,a1,...,ad] and which caches its
results would do.

Then your function should be a method for the class which is already
computing your discriminants, namely
sage.rings.polynomial.polynomial_element.Polynomial_generic_dense
but as there is already a discriminant function there you don't need a
new function, just special case code to catch small degrees (up to 3
or 4 at least) which gets the generic expression and substitutes
coefficients. This would be as simple as

if f.degree()<5:
return generic_univariate_discriminant(f.list())

As for the generic function, something like

def generic_univariate_discriminant(d):
R = PolynomialRing(ZZ,d,names='a')
a = R.gens()
if d==1: return 1
if d==2: return a[1]^2 - 4*a[0]*a[2]

if d==3: return a[1]**2*a[2]**2 - 4*a[0]*a[2]**3 - 4*a[1]**3*a[3] +
18*a[0]*a[1]*a[2]*a[3] - 27*a[0]**2*a[3]**2

(and so on) would do.

Go for it!

John

[martin....@gmx.net]

On Wednesday, March 26, 2014 1:32:13 PM UTC+1, John Cremona wrote:

You should definitely open a trac ticket.

 
I just opened #16014 for this. And would make this a link to http://trac.sagemath.org/ticket/16014 if the Google Groups editor weren't too broken to let me do so…

Don't expect that if you
post some code as text on a ticket that someone else will magically do
the right thing though!  If you have not before contributed to Sage
you can make this your chance to learn how that works, and people will
help.

I've submitted a couple of patches, e.g. for #13608 and #14239, and the one for #14403 which even got committed.

But I'll take that opportunity to come to grips with the new git-based workflow. After pulling recent changes into my existing git checkout, I'm currently running “./sage --upgrade”, and hope that after that and a “./sage -b” I have an up-to-date snapshot of the current sage master which I can then use as the basis for this contribution. If I should rather start from a released version, and avoid updating prior to contribution, please tell me so and also let me know why that's preferred.
 

I would separate out two things here: one is a generic computation of
discriminants of univariate polynomials of degree d, the result being
a polynomial in d+1 variables (the generic coefficients) over ZZ (that
is correct -- there's a unique map from ZZ to any other ring).  This
should be part of the current invariant_theory module, though when I
just looked at that it only appeared to implement invariants for
binary forms of degrees 2 and 4 (not 3) as well as some ternary stuff.

Uh… I had never seen that module before. But the way I see it, one would have to encode the degree of the polynomial into the name of some identifier, which would be possible if I did some case distinction as you indicated, but which feels weird nevertheless.
 

 As a stop-gap, a global function generic_univariate_discriminant(deg)
whose output is a polynomial in ZZ[a0,a1,...,ad] and which caches its
results would do.

Sounds more reasonable to me. Would you suggest that global to be in invariant_theory nevertheless, or rather in polynomial_element.pyx where it is used?
 

Then your function should be a method for the class which is already
computing your discriminants, namely
sage.rings.polynomial.polynomial_element.Polynomial_generic_dense
but as there is already a discriminant function there you don't need a
new function, just special case code to catch small degrees (up to 3
or 4 at least) which gets the generic expression and substitutes
coefficients.

 
I wonder whether doing the case distinction on degree is really the right thing to do. On the one hand, if the computation already takes ages for degree 3, it will be even worse for higher degrees. On the other hand, this computation is only bad due to the fact that the determinant computation on big polynomials is so expensive. If the base ring were just some finite field or some such, then I doubt we'd have to worry. So I wonder whether it might make more sense to do the case distinction based on the underlying ring. And in that case, I'm far from sure which cases should be handled which way. I guess I'd start by only detecting polynomial rings, and wait for other use cases to report problems. I don't know the most appropriate way to check whether a ring is a polynomial ring. Is some form of isinstance check the right way to do this, or should I combine checks is_PolynomialRing and is_MPolynomialRing?

Perhaps I could also add an argument “algorithm” to the discriminant method. It would default to None (meaning auto-selection) but could be used to force either direct computation or substitution. So possible values might be “direct” and “substitution”, perhaps others later on as well. That way, users could try alternatives to work out when a given method is faster than a different one, and could file bug reports requesting a change of default for certain cases.

This would be as simple as

if f.degree()<5:
    return generic_univariate_discriminant(f.list())

As for the generic function, something like

def generic_univariate_discriminant(d):
   R = PolynomialRing(ZZ,d,names='a')
   a = R.gens()
   if d==1: return 1
  if d==2: return a[1]^2 - 4*a[0]*a[2]

   if d==3: return a[1]**2*a[2]**2 - 4*a[0]*a[2]**3 - 4*a[1]**3*a[3] +
18*a[0]*a[1]*a[2]*a[3] - 27*a[0]**2*a[3]**2

(and so on) would do.

Hmmm… If we do the case distinction not on degree but on base ring, this doesn't look as attractive. Plus my own code would have been particularly efficient for sparse polynomials of high degree, since it would only do the generic computation for the non-zero terms. Hard-coding a few simple cases might still be faster than doing the code I posted for every case, but it also might be slower than using cached values. So perhaps we want a case distinction on degree to choose between cached polynomials for dense low degree cases and on-the fly computation for possibly sparse computation of high degree. The former using code like you posted, the latter code like I did, and both used for subsequent substitution.

Of course, all of this could as well be discussed on the trac ticket.

[John Cremona]

Let's continue the discussion on the trac ticket (your link did work),
ideally there should be alink back from there to this thread too.

Make a new branch from the trac develop branch (currently at
6.2.beta5), not the master (still at 6.1.1). Sorry for not
recognising your name as a seasoned contributor!

John

On 26 March 2014 14:13, <martin....@gmx.net> wrote:
> On Wednesday, March 26, 2014 1:32:13 PM UTC+1, John Cremona wrote:
>>
>> You should definitely open a trac ticket.
>
>
> I just opened #16014 for this. And would make this a link to
> http://trac.sagemath.org/ticket/16014 if the Google Groups editor weren't

> too broken to let me do so...

>
>> Don't expect that if you
>> post some code as text on a ticket that someone else will magically do
>> the right thing though! If you have not before contributed to Sage
>> you can make this your chance to learn how that works, and people will
>> help.
>
>
> I've submitted a couple of patches, e.g. for #13608 and #14239, and the one
> for #14403 which even got committed.
>
> But I'll take that opportunity to come to grips with the new git-based
> workflow. After pulling recent changes into my existing git checkout, I'm
> currently running "./sage --upgrade", and hope that after that and a "./sage
> -b" I have an up-to-date snapshot of the current sage master which I can
> then use as the basis for this contribution. If I should rather start from a
> released version, and avoid updating prior to contribution, please tell me
> so and also let me know why that's preferred.
>
>>
>> I would separate out two things here: one is a generic computation of
>> discriminants of univariate polynomials of degree d, the result being
>> a polynomial in d+1 variables (the generic coefficients) over ZZ (that
>> is correct -- there's a unique map from ZZ to any other ring). This
>> should be part of the current invariant_theory module, though when I
>> just looked at that it only appeared to implement invariants for
>> binary forms of degrees 2 and 4 (not 3) as well as some ternary stuff.
>
>

> Uh... I had never seen that module before. But the way I see it, one would

> Hmmm... If we do the case distinction not on degree but on base ring, this

> doesn't look as attractive. Plus my own code would have been particularly
> efficient for sparse polynomials of high degree, since it would only do the
> generic computation for the non-zero terms. Hard-coding a few simple cases
> might still be faster than doing the code I posted for every case, but it
> also might be slower than using cached values. So perhaps we want a case
> distinction on degree to choose between cached polynomials for dense low
> degree cases and on-the fly computation for possibly sparse computation of
> high degree. The former using code like you posted, the latter code like I
> did, and both used for subsequent substitution.
>
> Of course, all of this could as well be discussed on the trac ticket.
>