#16014: Improvements to discriminant computation
-------------------------------------+-------------------------------------
Reporter: gagern | Owner:
Type: defect | Status: needs_review
Priority: major | Milestone: sage-6.2
Component: algebra | Resolution:
Keywords: discriminant | Merged in:
Authors: | Reviewers:
Report Upstream: N/A | Work issues:
Branch: | Commit:
u/gagern/ticket/16014 | 66a8c198b7aa57f4567906baf964b30b88f4048f
Dependencies: | Stopgaps:
-------------------------------------+-------------------------------------
Changes (by gagern):
* status: new => needs_review
* commit: => 66a8c198b7aa57f4567906baf964b30b88f4048f
Old description:
> Currently (i.e. on sage 6.1.1), the discriminant computation fails with a
> segmenttation fault when the coefficients are too complicated. For
> example:
>
> {{{
> #!python
> PR.<b,t1,t2,x1,y1,x2,y2> = QQ[]
> PRmu.<mu> = PR[]
> E1 = diagonal_matrix(PR, [1, b^2, -b^2])
> RotScale = matrix(PR, [[1, -t1, 0], [t1, 1, 0], [0, 0, 1]])
> E2 = diagonal_matrix(PR, [1, 1, 1+t1^2])*RotScale.transpose()*E1*RotScale
> Transl = matrix(PR, [[1, 0, -x1], [0, 1, -y1], [0, 0, 1]])
> E3 = Transl.transpose()*E2*Transl
> E4 = E3.subs(t1=t2, x1=x2, y1=y2)
> det(mu*E3 + E4).discriminant()
> }}}
>
> Apparently the Python interpreter reaches some limit on the depth of the
> call stack, since it seems to call `symtable_visit_expr` from
> `Python-2.7.6/Python/symtable.c:1191` over and over and over again. But
> this ticket here is not about fixing Python. (If you want a ticket for
> that, feel free to create one and provide a reference here.)
>
> Instead, one should observe that the above is only a discriminant of a
> third degree polynomial. As such, it has an easy and well-known formula
> for the discriminant. Plugging the complicated coefficients into this
> easy formula leads to the solution rather quickly (considering the
> immense size of said solution). So perhaps we should employ this approach
> in the `discriminant` method.
>
> I've [https://groups.google.com/d/topic/sage-
> support/fABfhHyioa0/discussion discussed this question on sage-support],
> and will probably continue the discussion. Right now, I'm unsure about
> when to choose a new method over the existing one. Should it be based on
> degree, using hardcoded versions of the well-known discriminants for
> lower degrees? Or should it rather be based on base ring, thus employing
> a computation over some simple generators instead of the full
> coefficients if the base ring is a polynomial ring or some other fairly
> complex object. Here is a code snippet to illustrate the latter idea:
>
> {{{
> #!python
> def generic_discriminant(p):
> """
> Compute discriminant of univariate (sparse or dense) 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())
> }}}
>
> Personally I think I prefer a case distinction based on base ring, or a
> combination. If we do the case distinction based on base ring, then I'd
> also value some input as to what rings should make use of this new
> approach, and how to check for them. For example, I just noticed that
> #15061 discusses issues with discriminants of polynomials over power
> series. So I guess those might benefit from this approach as well.
>
> I intend to create a git branch once I know the direction to take, but
> I'd prefer some discussion up front to know what direction to take, how
> to make my case distinctions, and where to place my new code.
New description:
Currently (i.e. on sage 6.1.1), the discriminant computation fails with a
segmenttation fault when the coefficients are too complicated. For
example:
{{{
#!python
PR.<b,t1,t2,x1,y1,x2,y2> = QQ[]
PRmu.<mu> = PR[]
E1 = diagonal_matrix(PR, [1, b^2, -b^2])
RotScale = matrix(PR, [[1, -t1, 0], [t1, 1, 0], [0, 0, 1]])
E2 = diagonal_matrix(PR, [1, 1, 1+t1^2])*RotScale.transpose()*E1*RotScale
Transl = matrix(PR, [[1, 0, -x1], [0, 1, -y1], [0, 0, 1]])
E3 = Transl.transpose()*E2*Transl
E4 = E3.subs(t1=t2, x1=x2, y1=y2)
det(mu*E3 + E4).discriminant()
}}}
Apparently the Python interpreter reaches some limit on the depth of the
call stack, since it seems to call `symtable_visit_expr` from
`Python-2.7.6/Python/symtable.c:1191` over and over and over again. But
this ticket here is not about fixing Python. (If you want a ticket for
that, feel free to create one and provide a reference here.)
Instead, one should observe that the above is only a discriminant of a
third degree polynomial. As such, it has an easy and well-known formula
for the discriminant. Plugging the complicated coefficients into this easy
formula leads to the solution rather quickly (considering the immense size
of said solution). So perhaps we should employ this approach in the
`discriminant` method.
I've [https://groups.google.com/d/topic/sage-
support/fABfhHyioa0/discussion discussed this question on sage-support],
and will probably continue the discussion. Right now, I'm unsure about
when to choose a new method over the existing one. Should it be based on
degree, using hardcoded versions of the well-known discriminants for lower
degrees? Or should it rather be based on base ring, thus employing a
computation over some simple generators instead of the full coefficients
if the base ring is a polynomial ring or some other fairly complex object.
Here is a code snippet to illustrate the latter idea:
{{{
#!python
def generic_discriminant(p):
"""
Compute discriminant of univariate (sparse or dense) 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())
}}}
Personally I think I prefer a case distinction based on base ring, or a
combination. If we do the case distinction based on base ring, then I'd
also value some input as to what rings should make use of this new
approach, and how to check for them. For example, I just noticed that
#15061 discusses issues with discriminants of polynomials over power
series. So I guess those might benefit from this approach as well.
--
Comment:
New commits:
||[http://git.sagemath.org/sage.git/commit/?id=28980e9539c18fe3b9b7751ecd480e11e2089ddb
28980e9]||{{{Added test case to expose the segfault of discriminant.}}}||
||[http://git.sagemath.org/sage.git/commit/?id=0758de125303597e5d207e5c272a003287ebe705
0758de1]||{{{Drive-by fix to TeX notation in docstring.}}}||
||[http://git.sagemath.org/sage.git/commit/?id=d5e03075e3a9a8beac80396027413256d9b39289
d5e0307]||{{{Use simple coefficients if base is multivariate
polynomial.}}}||
||[http://git.sagemath.org/sage.git/commit/?id=66a8c198b7aa57f4567906baf964b30b88f4048f
66a8c19]||{{{Compute discriminant of power series polynomial using
substitution.}}}||
--
Ticket URL: <http://trac.sagemath.org/ticket/16014#comment:2>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-trac.
For more options, visit https://groups.google.com/d/optout.
