Solving solvable quintics: First implementation #1746
Conversation
| self.p, self.q, self.r, self.s = poly.all_coeffs()[2:] | ||
| self.F = self.getF() | ||
|
|
||
| def getf20(self): |
There was a problem hiding this comment.
It would be more pythonic to drop the "get" from all these function names and decorate them with @Property. For example,
@property
def f20(self):
...And the access it as PolyQuintic().f20.
There was a problem hiding this comment.
For the ones that take arguments, keep them as functions, but still drop the "get".
|
Add a reference to the the paper and Mathematica file in a comment at the top of the quintics file. |
|
Just so you know, this does need tests. I'm not sure how best to do it. We need to make sure the formulas are correctly input, and that all the lines are covered. You can do the latter by running the coverage tool |
|
If you substitute those solutions back into the quintic and simplify, does it go to zero? |
| Res[3] = solve_five(sol**5 - R3, sol) | ||
| Res[4] = solve_five(sol**5 - R4, sol) | ||
|
|
||
| for i in range(1, 5, 1): |
There was a problem hiding this comment.
The third argument is not necessary.
|
I'll make the changes and commit soon. |
|
As far as checking the solutions is concerned, here's what I did: As you can see, the roots are same. |
|
|
||
| r1 = Res[1][0] | ||
| for root in Res[4]: | ||
| if Abs(im(r1*root)) < S(0.00001): |
There was a problem hiding this comment.
perhaps pure_complex should be used here. That 1e-5 test looks like it would be easy to break.
There was a problem hiding this comment.
I don't quite see the point of using pure_complex. It would just give me (re, im) of r1*root while I only want the imaginary part for which I am using im().
As an example:
In [42]: z = 23.234542525 + I*3653.34653636
In [43]: pure_complex(z)
Out[43]: (23.2345425250000, 3653.34653636000)
In [44]: re
Out[44]: sympy.functions.elementary.complexes.re
In [45]: re(z), im(z)
Out[45]: (23.2345425250000, 3653.34653636000)
There was a problem hiding this comment.
You are just testing that the imaginary part is insignificant, right? Then perhaps check that im(foo.n(chop=True)) == 0
>>> im(S(3))
0
>>> im(S(3)+I)
1
>>> im(S(3)+I*1e-22)
1.00000000000000e-22
>>> im((S(3)+I*1e-22).n(chop=True))
0
>>> im((S(3)+I*1e-12).n(chop=True))
1.00000000000000e-12
But again, you get into the "how small is small" issue so perhaps the magnitude of the im part should be compared to the real part ith the comp function again
>>> im((S(3)*1e-22+I*1e-22).n(chop=True)) # this shouldn't have chopped, IMO
0
There was a problem hiding this comment.
Well, that is an issue. In this case, I think that 10**-6 is small enough though I suppose we can reduce it further to about 10**-10. The thing is that we don't want to have too small a value to compare against since then the difference might not small enough.
While debugging, I observed that if the difference isn't exactly zero, then it was of the order 10**-15 to 10**-17. So, I think 10**-10 should be small enough but not too small that it undermines the purpose if comparing.
|
I think there's some problem somewhere in the PolyQuintic class. In case of some other solvable equations, for which f20 doesn't have a zero root, we are not getting correct results. I'll look into it. |
|
I've implemented the You can match the solutions from here: Also, I tried to check the numerical solution with SymPy: As you can clearly see, I've stumbled upon another bug. The |
|
|
||
| for r2temp in Res[2]: | ||
| for r3temp in Res[3]: | ||
| # Again storig away exact number and using |
|
Now, Clearly, if the quintic is not solvable, or, if upon dividing all the coefficients by the coeff of the leading term we do not get all integers, then it falls back to ordinary solve. |
|
Also, the method seems much faster on a personal machine than on a shared machine. On shared machines, it took ~ 1 minute to solve |
|
Doesn't the algorithm work with rational coefficients? |
|
Let me check. |
|
Yes. It should work with rational coefficients too. But the quintic The solutions match. http://www.wolframalpha.com/input/?i=solve+x**5+-+110*x**3+-+55*x**2+%2B+2310*x+%2B+979&dataset= |
|
It seems we can now add tests. |
|
But that's a really long time. Why does it take so long? |
|
Well, as I said before, when the quintics are like: x5 + a*x + b =0, then we can get away with 15 - 20 secs of solution time. Let me call this equation (1)
In the end, I suppose we can't have too much speed for (2) at least. If this were a complied language, I suppose the solution time could have been reduced to about 60 secs. |
|
If you comment out those simplify calls, how long does it take? And can you give me an example of a before and after with one of them, so I can see what kind of simplification it is doing? |
|
I guess this is coming from the design of solve that is already there (so maybe @smichr can answer this better), but why not always use this algorithm when it can be applied (in other words, why do we have the quintics flag)? Ditto for cubic and quartic. |
| @@ -501,6 +641,7 @@ def roots(f, *gens, **flags): | |||
| auto = flags.pop('auto', True) | |||
| cubics = flags.pop('cubics', True) | |||
| quartics = flags.pop('quartics', True) | |||
| quintics = flags.pop('quintics', False) | |||
There was a problem hiding this comment.
Oh, so I guess cubics and quartics default to True. So why did you make this default to False? It's not very useful if it's never enabled.
There was a problem hiding this comment.
This is intentional. As you have seen, while there are quintics which can be evaluated quickly, there are also those which can take ~ 7 minutes to solve. Surely, we don't want the user to wait for 7 minutes to get a solution? This way, we can document this and if the user wants, then he can go for exact solutions knowing full well that the solution may take a long time to compute.
Also, another reason is this: Suppose that the user has a script. In that script somewhere, there's a solvable quintic. So, to the user, the script may seem to be hanged because of the quintic. I'm sure that won't be a good thing as far as the user satisfaction is concerned.
There was a problem hiding this comment.
Well, the solution here is obviously to make things fast, so that it never takes a long time. I am convinced that it can be done. We just need to figure out what is slow, and how to make it faster.
And anyway, my philosophy is that if it is off by default, then it might as well not exist, because 99% of users will never know it was there. I also believe that it's better to have a slow function that returns a solution than a fast function that doesn't. So I think even if it is slow that this should be True by default.
There was a problem hiding this comment.
Though we can always start with it being off, merge it, and then after it gets more tested, flip the default to True.
There was a problem hiding this comment.
Except if it's off, no one will use it, and it will never be tested.
|
In short, we need to reduce the time of this by a lot. At this rate, it will be impossible to even test the algorithm. |
| if coeff_5 != 1: | ||
| l = [p/coeff_5, q/coeff_5, r/coeff_5, s/coeff_5] | ||
| for coeff in l: | ||
| if not coeff.is_rational: |
There was a problem hiding this comment.
You should use is_Rational (with a capital R). Otherwise, this will work with something like Symbol('a', rational=True), which we don't want.
There was a problem hiding this comment.
You can simplify statements like this by using all or any and a list comprehension. This can be written as
if not all(coeff.is_Rational for coeff in l):
return result
|
You might want to try using line_profilier/kernprof to get an idea of exactly what parts of the algorithm are slow. |
|
Here are the comparisons for f(x) = x*_5 + 15_x + 12: For the same function: I subtracted these outputs and evaluated them. I got |
| d = discriminant(f) | ||
| delta = sqrt(d) | ||
| # zeta = a fifth root of unity | ||
| zeta = cos(2*pi/5) + I*sin(2*pi/5) |
There was a problem hiding this comment.
OK, here is one thing. If zeta = exp(2*I/5), then zeta**2 = exp(4*I/5) (and so on). So I think we can more efficiently compute the powers of zeta in R1, R2, R3, and R4.
And even if we couldn't, we could at least expand them only once and then use that in each, removing the duplicated computations in lines 395-398.
There was a problem hiding this comment.
And I'm curious why you use exp(2*pi/5) as your primitive root instead of just exp(pi/5). Does the paper call for this?
There was a problem hiding this comment.
exp(I*pi/5) isn't a fifth root of unity; exp(I*pi/5)**5 = -1.
There was a problem hiding this comment.
Oh right, I forgot it has to be exp(2*pi*I*k/n) (I forgot the 2).
|
I don't understand why it is happening, but the Travis error is is legitimately something that you broke with this branch. |
|
SymPy Bot Summary: ❗ There were merge conflicts (could not merge prasoon2211/quintic (6efff32) into master (49eea7e)); could not test the branch. |
|
Okay then. I'll fix the merge conflict and fix the tests. Thanks. |
|
Sometimes when I get a messy rebase I prefer to just re-add my changes by hand. You have to make your changes in context of the updated solver. |
|
This doesn't look messy. It's just one conflict. It's just that it might not be clear how to correctly resolve it if you don't know what the code is doing. |
|
By the way @certik, I just realized that your method for bisecting cache bugs is no longer relevant, since the cache is always cleared for each test file (and it has been that way for some time). See https://github.com/sympy/sympy/blob/master/sympy/utilities/runtests.py#L845 and 92ef41f. |
So go to master; copy solvers.py; go to your quintics branch; open solvers.py and paste what was in master and save it. You should now be in a green merge state. Now make your changes to the solvers.py file. Look for the word 'quartic' and make your changes there. |
|
btw, there are two places in solvers.py where quartics are being allowed...I'm not sure if they should both be modified or not, but you should look at both of them and raise any questions you may have. |
|
I have just merged from master and the conflicts have been resolved. And yes, quintics will go at both places where quartics is. I am still trying to find a way for the tests to work. As I mentioned somewhere above, if we check for two different quintics, then the tests hang. But if we test just one of them, tests run fine. So, I'm trying to make a workaround using this. Thanks for you help. |
|
The test pass ok, here. I guess having it XFAILed is ok, but won't this hang in practice for the same reason that it is hanging during tests? (So does that mean that someone will only be able to solve a quintic once with SymPy and the second time, in a given session, it will hang?) Since you are just finishing this, I wonder if you would be a good candidate to look into issue 1890, especially my comment 1: I suspect RootOf instances should be returned rather than symbolic forms that may be wrong. Alternatively Piecewise could be used to returned a solution where different pieces correspond to the different forms of the roots for different conditions on the coefficients. |
|
As to your first question, no, it won't be a problem. A user can solve any number of quintics with the interactive session, as well as with standalone scripts. The problem here seems to be explicitly with the way the testing framework works. |
|
SymPy Bot Summary: ✳️ Passed after merging prasoon2211/quintic (4f45c5a) into master (07933cd). |
| @@ -38,7 +38,6 @@ def clashing(): | |||
|
|
|||
| ns = {} | |||
| exec 'from sympy import *' in ns | |||
|
|
|||
There was a problem hiding this comment.
I would leave that blank line in place.
And the constants file line lengths are egregiously long. If you write them in parentheses they can be broken at will, typically after a + or -.
y = (
x +
y +
z)
There was a problem hiding this comment.
I'll add the line (it got removed when I was trying to fix the clashing bug earlier). And yes, the question about length of lines in the constants file did bother me initially. In the end I decided to keep one constant on one line. The reason for this being - the constant definitions don't contain any programming logic. They are just definitions. For someone reading the code, I think, it improves readability to keep one long constant definition on a separate line. This way, the user doesn't get bogged down with 100s of lines of constant definitions that don't have anything to with the actual implementation of the solver logic.
That being said, I have nothing against the 80 char limit in particular (I am obviously talking in reference of the constants file). If indeed wrapping lines to 80 chars is required, please let me know.
|
How were the constants obtained and how do we know there is no typo in them? |
|
The constants we obtained from a mathematica notebook file that is mentioned in the docstring. I ran a script to convert them to SymPy specific format.
Those 2 points in conjunction prove that there exist no typos in the code :) |
Symbolically or numerically checked? And are all constants used for every calculation? |
|
Checked numerically. Mathematica doesn't provide exact solutions to solvable quintics. And yes, every single constant is used for a solution. |
|
OK, I changed |
|
The tests are passing because f20 isn't used directly in the solution. f20 is just used to check whether quintics is solvable or not. So, when you did the change, it just changed the f(x) defined in test to became non-solvable. Then, the polynomial is just solved numerically within SymPy and if you look closely at the tests, the test just see whether the solution make f(x) zero at the required values, whether these solutions be numerical or exact. |
|
But the tests should be structured so that changes like that do cause them to fail. |
|
Yes. I'll do that. On 3/28/13, Aaron Meurer notifications@github.com wrote:
|
|
|
||
| f = x**5 - 15*x**3 - 5*x**2 + 10*x + 20 | ||
| s = solve(f) | ||
| for root in s: |
There was a problem hiding this comment.
I'm not able to get the exact solutions of f, what I get is unevaluated RootOf object. I doubt this is the expected behaviour.
In [31]: solve(f, quintics=True)
Out[31]:
[RootOf(x**5 - 15*x**3 - 5*x**2 + 10*x + 20, 0),
RootOf(x**5 - 15*x**3 - 5*x**2 + 10*x + 20, 1),
RootOf(x**5 - 15*x**3 - 5*x**2 + 10*x + 20, 2),
RootOf(x**5 - 15*x**3 - 5*x**2 + 10*x + 20, 3),
RootOf(x**5 - 15*x**3 - 5*x**2 + 10*x + 20, 4)]
There was a problem hiding this comment.
If you see line 216, you shall find that this is indeed the expected behaviour. I'm sorry that I can't be of more help regarding this issue - I'm terribly bogged down by other matters. But, I should like to refer you to the docstring at the top of sympy/polys/polyquinticconst.py. The links provided there will be of much help, I think.
|
Changes Unknown when pulling 4f45c5a on prasoon2211:quintic into * on sympy:master*. |
Please see this:
http://code.google.com/p/sympy/issues/detail?id=3548
Also, see this:
https://groups.google.com/forum/?fromgroups=#!topic/sympy/xXyiOUWH0SU
According to the discussions on above links, this PR introduces solving quintics exactly (not numerically) when they are solvable.
There are a lot of constants involved. So, I created a new file and made a new class for quintics containing required constants.
Also note that solving quintics exactly takes( for x5 + 15*x +12 =0 ) ~ 40 secs. The long time is attributed to solving five equations of the form x5 - R =0 somewhere in the middle of the code.
Also, I have used simplify a few times. It seems that the code became almost 3 times as fast with using simplify, than without it.
There is an extra flag; quintics. When set to true, the new roots_quintic method is called and if the quintic is solvable, the solution is returned. I am attaching a screenshot :
Here you can see the solution of x*_5 + 15_x +12 = 0.
I have not added any tests.
@asmeurer Please look at the code. Let me know for changes required.