Opened 5 years ago

Closed 4 years ago

# Improvements to discriminant computation

Reported by: Owned by: gagern major sage-6.2 algebra discriminant pbruin Martin von Gagern Peter Bruin N/A 1892c84 (Commits)

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 `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 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:

```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:1 Changed 5 years ago by gagern

• 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 5 years ago by gagern

• Commit set to 66a8c198b7aa57f4567906baf964b30b88f4048f
• Description modified (diff)
• Status changed from new to needs_review

New commits:

 ​28980e9 `Added test case to expose the segfault of discriminant.` ​0758de1 `Drive-by fix to TeX notation in docstring.` ​d5e0307 `Use simple coefficients if base is multivariate polynomial.` ​66a8c19 `Compute discriminant of power series polynomial using substitution.`

### comment:3 Changed 4 years ago by gagern

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:5 Changed 4 years ago by pbruin

This looks good and passes doctests. However, I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sage-dev 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 an `algorithm` 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 git

• 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 follow-up: ↓ 8 Changed 4 years ago by gagern

I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sage-dev 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 non-public?

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 pbruin

• Reviewers set to Peter Bruin

I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sage-dev 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` 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 non-public?

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 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?

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 time-critical 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 line-by-line 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 pbruin

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 gagern

• Authors set to Martin von Gagern
• Status changed from needs_review to positive_review

### comment:11 Changed 4 years ago by vbraun

• Branch changed from u/gagern/ticket/16014 to 1892c84d6372e882fa0c39667d06ba83cc61f3be
• Resolution set to fixed
• Status changed from positive_review to closed

### comment:12 Changed 2 years ago by embray

• 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`.

Note: See TracTickets for help on using tickets.