sympy/sympy — pull request #7393 by certik (2014-04-09 19:50:11+00:00)
Title: [WIP] Make nroots() robust

Previously:
```py
from sympy import var, legendre_poly
var("x")
n = 64
p = legendre_poly(n, x, polys=True)
print p.nroots(n=3)
```
produced:
```
[-1.77, -1.6 - 0.456*I, -1.6 + 0.456*I, -1.22 - 0.639*I, -1.22 + 0.639*I, -1.20, -0.979 - 0.56*I, -0.979 + 0.56*I, -0.823 + 0.561*I, -0.823 - 0.561*I, -0.701 - 0.466*I, -0.701 + 0.466*I, -0.579 - 0.446*I, -0.579 + 0.446*I, -0.57 + 0.172*I, -0.531 - 0.0702*I, -0.459 - 0.351*I, -0.456 + 0.354*I, -0.44 - 0.109*I, -0.404 + 0.115*I, -0.371 - 0.26*I, -0.371 + 0.26*I, -0.359 - 0.133*I, -0.308 + 0.177*I, -0.272 - 0.112*I, -0.247 + 0.0635*I, -0.209 - 0.0534*I, -0.178 + 0.0231*I, -0.119 - 0.0344*I, -0.0796 - 0.0155*I, -0.0564 + 0.144*I, -0.0366 - 0.00806*I, 0.0233 - 0.00192*I, 0.0736 - 0.000641*I, 0.12 - 0.00184*I, 0.167 + 0.0116*I, 0.18 - 0.0771*I, 0.225 + 0.0617*I, 0.268 + 0.173*I, 0.272 - 0.111*I, 0.296 - 0.0282*I, 0.311 - 0.174*I, 0.311 + 0.174*I, 0.364 - 0.246*I, 0.367 + 0.23*I, 0.401 + 0.0157*I, 0.456 - 0.354*I, 0.457 + 0.354*I, 0.57 - 0.172*I, 0.57 + 0.172*I, 0.579 - 0.446*I, 0.579 + 0.446*I, 0.7 - 0.466*I, 0.703 + 0.488*I, 0.736 + 0.539*I, 0.804 - 0.461*I, 0.979 - 0.56*I, 0.979 + 0.56*I, 1.20, 1.22 - 0.639*I, 1.22 + 0.639*I, 1.6 - 0.456*I, 1.6 + 0.456*I, 1.77]
```
Which is just **wrong**.
Now it not only produces the correct answer:
```
[-0.999, -0.996, -0.991, -0.983, -0.973, -0.961, -0.946, -0.930, -0.911, -0.889, -0.866, -0.841, -0.813, -0.784, -0.753, -0.720, -0.685, -0.649, -0.611, -0.572, -0.531, -0.489, -0.446, -0.402, -0.357, -0.311, -0.265, -0.217, -0.170, -0.121, -0.0730, -0.0243, 0.0243, 0.0730, 0.121, 0.170, 0.217, 0.265, 0.311, 0.357, 0.402, 0.446, 0.489, 0.531, 0.572, 0.611, 0.649, 0.685, 0.720, 0.753, 0.784, 0.813, 0.841, 0.866, 0.889, 0.911, 0.930, 0.946, 0.961, 0.973, 0.983, 0.991, 0.996, 0.999]
```
But it also raises an exception if it can't get convergence to the given accuracy, e.g. with this change:
```py
print p.nroots(n=3, maxsteps=10)
```
it gives:
```
Traceback (most recent call last):
File "a.py", line 7, in
print p.nroots(n=3, maxsteps=10)
File "/auto/nest/nest/u/ondrej/repos/sympy/sympy/polys/polytools.py", line 3397, in nroots
cleanup=cleanup, error=False, extraprec=f.degree()*10)
File "/auto/nest/nest/u/ondrej/repos/sympy/sympy/mpmath/calculus/polynomials.py", line 176, in polyroots
raise Exception("Didn't converge in maxsteps=%d steps." % maxsteps)
Exception: Didn't converge in maxsteps=10 steps.
```

With this better `nroots()` we can now almost use it in `quadrature.py` as follows:
```diff
--- a/sympy/integrals/quadrature.py
+++ b/sympy/integrals/quadrature.py
@@ -69,10 +69,9 @@ def gauss_legendre(n, n_digits):
pd = p.diff(x)
xi = []
w = []
- for r in p.real_roots():
- if isinstance(r, RootOf):
- r = r.eval_rational(S(1)/10**(n_digits+2))
- xi.append(r.n(n_digits))
+ roots = p.nroots(n_digits)
+ for r in roots:
+ xi.append(r)
w.append((2/((1-r**2) * pd.subs(x, r)**2)).n(n_digits))
return xi, w
```
This gets a lot faster. The `xi` points are as accurate as before, but the `w` points are now very inaccurate, due to evaluating the `pd`, which is a polynomial of high order. This is not a small error, i.e. in one of the tests that fail, instead of `0.5555555555555555555555555555555555555556`, we got `0.5555555555574520784163542152657447087753`. One solution is to modify `nroots()` further and allow to also optionally return accurate derivatives, but one has to figure out what the best approach is. I will work on that once this PR is in.

TODO:

* [ ] Write tests that previously fail and now pass