#16014 closed defect (fixed)
Improvements to discriminant computation
Reported by:  gagern  Owned by:  

Priority:  major  Milestone:  sage6.2 
Component:  algebra  Keywords:  discriminant 
Cc:  pbruin  Merged in:  
Authors:  Martin von Gagern  Reviewers:  Peter Bruin 
Report Upstream:  N/A  Work issues:  
Branch:  1892c84 (Commits)  Commit:  
Dependencies:  Stopgaps: 
Description (last modified by )
Currently (i.e. on sage 6.1.1), the discriminant computation fails with a segmenttation fault when the coefficients are too complicated. For example:
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 Python2.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 wellknown 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 discussed this question on sagesupport, 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 wellknown 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:
def generic_discriminant(p): """ Compute discriminant of univariate (sparse or dense) polynomial. The discriminant is computed in two steps. First all nonzero 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.
Change History (12)
comment:1 Changed 4 years ago by
 Branch set to u/gagern/ticket/16014
 Created changed from 03/26/14 14:11:15 to 03/26/14 14:11:15
 Modified changed from 03/26/14 14:11:15 to 03/26/14 14:11:15
comment:2 Changed 4 years ago by
 Commit set to 66a8c198b7aa57f4567906baf964b30b88f4048f
 Description modified (diff)
 Status changed from new to needs_review
comment:3 Changed 4 years ago by
What I implemented so far is a case distinction based on the base ring, with the new implementation being used for multivariate polynomials and power series.
I implemented the new code in the existing method for now, in part because I wasn't sure whether I should add new methods (e.g. discriminant_by_resultant
and discriminant_by_substitution
) or rather an algorithm
argument for the discriminant
method. In the latter case, the code would stay in that method but be reorganized because the substitution version could call the resultant version instead of executing that computation directly.
I'm not sure about this whole ZeroDivisionError
stuff from #11782 resp. a224fe3. When doing the substitution approach, there should be no zero divisors, so replacing self
by poly
in that exception handler is only for the sake of consistency, but not really required. I guess if I were to write code for the substitution approach exclusively, I might consider using // poly[n]
instead of * an
to avoid the detour to the fraction field.
comment:4 Changed 4 years ago by
 Cc pbruin added
comment:5 Changed 4 years ago by
This looks good and passes doctests. However, I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sagedev to have a separate function generic_univariate_discriminant()
(or similar), which caches its output, and hardcodes it for small degrees. This would be useful for two reasons:
 computing the universal discriminant takes a substantial amount of time (several seconds for degree 7, several minutes for degree 8, probably extremely long for higher degrees);
 it would make the
discriminant()
method cleaner and therefore more attractive to add analgorithm
flag and/or a more sophisticated algorithm selection heuristic.
I would think sage.rings.polynomial.polynomial_element
would be a logical place for such a function, although sage.rings.invariant_theory
would also be defensible.
Here is a simple function to compute the universal discriminant:
def universal_discriminant(n): r"""Return the discriminant of the 'universal' univariate polynomial `a_n x^n + \cdots + a_1 x + a_0` in `\ZZ[a_0, \ldots, a_n][x]`.""" R = PolynomialRing(ZZ, n + 1, 'a') S = PolynomialRing(R, 'x') return S(list(R.gens())).discriminant()
comment:6 Changed 4 years ago by
 Commit changed from 66a8c198b7aa57f4567906baf964b30b88f4048f to 1892c84d6372e882fa0c39667d06ba83cc61f3be
Branch pushed to git repo; I updated commit sha1. New commits:
1892c84  Introduce a caching function universal_discriminant.

comment:7 followup: ↓ 8 Changed 4 years ago by
Replying to pbruin:
I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sagedev to have a separate function
generic_univariate_discriminant()
(or similar), which caches its output, and hardcodes it for small degrees.
Sure I can be convinced, particularly if two people request something. I've implemented the caching, but not the hardcoding so far. I wonder how small “small” actually should be. Up to 7, things compute fairly quickly so there seems to be little benefit from hardcoding, particularly since we do caching. For higher degrees, computing them for hardcoding takes considerable time, and writing down the hard code will take up quite a lot of space in the source file. Should this really be implemented in Cython code, or should we make some effort to obtain or compute discriminants up to say 10 or so and store them in some kind of persistent data representation instead of source code?
I also wonder how public this function should be. I prefer polynomial_element
over invariant_theory
simply because it keeps all my changes in one place. I've written a full docstring for the code, and inside that docstring I had to import the function. Is this a reasonable approach? Should I add the function to sage.rings.polynomial.all
? Should I prepend the function name with an underscore to mark it nonpublic?
I also have a question regarding my late import. I noticed the cdef void late_import
near the beginning of the module. I could add is_PowerSeriesRing
to that. But on the other hand, it seems as if that function only gets called by Polynomial.roots
. Which suggests that I might not be able to rely on is_PowertSeriesRing
if I were to follow that approach. Do you agree with keeping the import where I have it so far? Should the function late_import
be renamed to something like late_import_for_root
to make it clearer that this does no general late importing but specific to that method? Or should I add my imports to that function and make sure to call late_import
in discriminant
as well?
Should I add myself as an author of that module, at the beginning of the file?
comment:8 in reply to: ↑ 7 Changed 4 years ago by
 Reviewers set to Peter Bruin
Replying to gagern:
Replying to pbruin:
I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sagedev to have a separate function
generic_univariate_discriminant()
(or similar), which caches its output, and hardcodes it for small degrees.Sure I can be convinced, particularly if two people request something. I've implemented the caching, but not the hardcoding so far. I wonder how small “small” actually should be. Up to 7, things compute fairly quickly so there seems to be little benefit from hardcoding, particularly since we do caching.
Thanks, this looks very good. You're probably right about not hardcoding the polynomials.
For higher degrees, computing them for hardcoding takes considerable time, and writing down the hard code will take up quite a lot of space in the source file. Should this really be implemented in Cython code, or should we make some effort to obtain or compute discriminants up to say 10 or so and store them in some kind of persistent data representation instead of source code?
Maybe, but I doubt it is very important. If at all, it should be done in a later ticket.
I also wonder how public this function should be. I prefer
polynomial_element
overinvariant_theory
simply because it keeps all my changes in one place. I've written a full docstring for the code, and inside that docstring I had to import the function. Is this a reasonable approach? Should I add the function tosage.rings.polynomial.all
? Should I prepend the function name with an underscore to mark it nonpublic?
I actually agree completely with the way you did all these things. I wouldn't add anything to sage.rings.polynomial.all
; this would add it to the global namespace, which is already very crowded.
I also have a question regarding my late import. I noticed the
cdef void late_import
near the beginning of the module. I could addis_PowerSeriesRing
to that. But on the other hand, it seems as if that function only gets called byPolynomial.roots
. Which suggests that I might not be able to rely onis_PowertSeriesRing
if I were to follow that approach. Do you agree with keeping the import where I have it so far? Should the functionlate_import
be renamed to something likelate_import_for_root
to make it clearer that this does no general late importing but specific to that method? Or should I add my imports to that function and make sure to calllate_import
indiscriminant
as well?
Here, too, I think the way you did it is absolutely fine. In fact one can argue that modules should generally be imported locally in the function that needs it. On the other hand, global imports (including the late_import
approach) are useful if the imported module is used in many places, or if the function is timecritical enough that a local import would add too much overhead.
Should I add myself as an author of that module, at the beginning of the file?
If you want; as you will have seen, authors for some other additions/changes are listed, but it is certainly not done in a very systematic way througout Sage. In any case, your changes and authorship will be remembered by the revision control system, and using git annotate
gives a linebyline overview of who wrote/changed what.
Could you also fill in your name (i.e. real name, not Trac username) in the "Author" field of this ticket?
I will test this and then set it to positive review.
comment:9 Changed 4 years ago by
All tests pass and I am happy with the patch as it is. Please add your name and set the ticket to positive review.
comment:10 Changed 4 years ago by
 Status changed from needs_review to positive_review
comment:11 Changed 4 years ago by
 Branch changed from u/gagern/ticket/16014 to 1892c84d6372e882fa0c39667d06ba83cc61f3be
 Resolution set to fixed
 Status changed from positive_review to closed
comment:12 Changed 17 months ago by
 Commit 1892c84d6372e882fa0c39667d06ba83cc61f3be deleted
FWIW the segfault in Python was fixed in newer Python versions: http://bugs.python.org/issue5765
For some reason (it seems just developer resources) a version of the patch was never backported to 2.7. I might take a stab at that (and open a separate ticket if so). Note the patch doesn't magically allow unlimited stack depth and one can still blow the stack. But with the default settings it will raise a RuntimeError
instead, and give one the opportunity to increase the limit with sys.setrecursionlimit
.
New commits:
Added test case to expose the segfault of discriminant.
Driveby fix to TeX notation in docstring.
Use simple coefficients if base is multivariate polynomial.
Compute discriminant of power series polynomial using substitution.