Repository: sympy/sympy
Pull request: improve identification of imaginary roots for CRootOf instances
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234". See
https://github.com/blog/1506-closing-issues-via-pull-requests .-->

#14293 #14291 

#### Brief description of what is fixed or changed

fixes #14615

Fixes #14291. Because not all imaginary intervals were being identified, some imaginary intervals found their way into code that attempts to refine the interval until both bottom left and upper right were modified. An imaginary interval containing an x-value of 0 at one of the corners will never refine that corner and an infinite loop ensues.

The imaginary identification has now been made robust allowing intervals to be clearly identified as having a real part or being purely imaginary. This allows an `_eval_is_imaginary` routine to be added to CRootOf: so where `CRootOf(x**5+x-1,3).is_imaginary` formerly gave None, it now gives False since it is not real and not imaginary. The cost of doing this is extra refinement at instantiation.

#### Other comments

The sorting of complex roots has been abandoned until a method of determining when two roots share a common real value. (See discussion at #14293.) Instead, we settle for the invariant that the real roots come first followed by complex pairs -- in no guaranteed order -- but with the one having a negative real part coming before the positive real part. So the conjugate of a non-real root will always be an index of +/-1 from a root's index:
```
>>> CRootOf(x**4 + x + 1, 1).conjugate()
CRootOf(x**4 + x + 1, 0)
>>> _.conjugate()
CRootOf(x**4 + x + 1, 1)
```

The docstring of ComplexInterval has been significantly expanded.

An attempt to keep cache entries all PurePoly instances has been made and assertions inserted wherever the cache is accessed.

`_eval_Eq` has been simplified to reflect the fact that refinement has already taken place at instantiation. Also the bounds are included in the check for containment, e.g. `a <= other and other <= b` rather than `a < other and other < b`.

- [x] check that root cache is used when retrieving roots.  The only place that the root cache is being used in rootoftools is when the interval is being retreived.

Problem statement:
poly(((x - 1)**2 + 1)*((x - 1)**2 + 2)*(x - 1)).all_roots() hangs
bisected to 73fafba


wrong simplification of sum of fractions
The result should be zero:
```
In [14]: e1 = (2*b*x+a)/(2*sqrt(c*x**2)*c*x)

In [15]: e2 = (b*x+a)**2/(2*sqrt(c*x**2)*c*x*a)

In [16]: e1.diff(x)-e2.diff(x)
Out[16]: 
                                                               2   
      b           a + 2⋅b⋅x        b⋅(a + b⋅x)        (a + b⋅x)    
───────────── - ────────────── - ─────────────── + ────────────────
       ______           ______            ______             ______
      ╱    2       2   ╱    2            ╱    2         2   ╱    2 
c⋅x⋅╲╱  c⋅x     c⋅x ⋅╲╱  c⋅x     a⋅c⋅x⋅╲╱  c⋅x     a⋅c⋅x ⋅╲╱  c⋅x  

In [17]: _.simplify()
Out[17]: 
   ⎛                    7/2⎞
 2 ⎜   9/2  7     ⎛   2⎞   ⎟
b ⋅⎝- c   ⋅x  + c⋅⎝c⋅x ⎠   ⎠
────────────────────────────
             6  8           
          a⋅c ⋅x            

In [20]: (e1.diff(x)-e2.diff(x)).expand()
Out[20]: 0
```

Hints from issue discussion:
I found that the infinite loop occurs [here](https://github.com/sympy/sympy/blob/master/sympy/polys/rootoftools.py#L323). I am having some difficulty in understanding the code/work of this function. @skirpichev  @jksuom can you show me some examples that what `_refine_complexes` do ?
> what _refine_complexes do

Isn't this the routine that tries to locate each root within a rectangle that is horizontally and vertically distinct from every other root's rectangle?
The problem is that the complex roots of two equations may share the same real value so the roots will never be isolated to satisfy the assumptions of the routine (which wrongly assumed this would not be the case).
@smichr I couldn't locate where in code it is checking the same real value of roots? please help!
And `complexes` in _refine_complexes() is:
```
[((0, 6) x (-6, 0), Poly(x**2 - 2*x + 3, x, domain='ZZ'), 1), ((0, 6) x (0, 6), Poly(x**2 - 2*x + 3, x, domain='ZZ'), 1), ((0, 4) x (-4, 0), Poly(x**2 - 2*x + 2, x, domain='ZZ'), 1), ((0, 4) x (0, 4), Poly(x**2 - 2*x + 2, x, domain='ZZ'), 1)]
```
What are those `(0, 6)`, `(0, -6)`?
> The problem is that the complex roots of two equations may share the same real value

Oh, again wrong (:  Complex roots of one equation may share the same real part (actually, we care about different conjugated pairs of roots, if coefficients of the equation are real).

BTW, @vishalg2235, I wonder what you did to "locate where in code it is"?

@skirpichev I was trying to understand `_refine_complexes()` but couldn't understand what logic is applied here. Can you please get me more familiar with this code ?

Commit message: allow is_disjoint to compare re and cmplx intervals

Patch to review:
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index df22c68955fa..ad7177365ecd 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -17,6 +17,7 @@
 from sympy.core.compatibility import (
     as_int, integer_types, long, string_types, with_metaclass, HAS_GMPY,
     SYMPY_INTS, int_info)
+
 import mpmath
 import mpmath.libmp as mlib
 from mpmath.libmp.backend import MPZ
@@ -1482,35 +1483,38 @@ def __new__(cls, p, q=None, gcd=None):
             if isinstance(p, Rational):
                 return p
 
-            if isinstance(p, string_types):
-                if p.count('/') > 1:
-                    raise TypeError('invalid input: %s' % p)
-                pq = p.rsplit('/', 1)
-                if len(pq) == 2:
-                    p, q = pq
-                    fp = fractions.Fraction(p)
-                    fq = fractions.Fraction(q)
-                    f = fp/fq
-                    return Rational(f.numerator, f.denominator, 1)
-                p = p.replace(' ', '')
-                try:
-                    p = fractions.Fraction(p)
-                except ValueError:
-                    pass  # error will raise below
+            if isinstance(p, SYMPY_INTS):
+                pass
+            else:
+                if isinstance(p, (float, Float)):
+                    return Rational(*_as_integer_ratio(p))
 
-            if not isinstance(p, string_types):
-                try:
-                    if isinstance(p, fractions.Fraction):
+                if not isinstance(p, string_types):
+                    try:
+                        p = sympify(p)
+                    except (SympifyError, SyntaxError):
+                        pass  # error will raise below
+                else:
+                    if p.count('/') > 1:
+                        raise TypeError('invalid input: %s' % p)
+                    p = p.replace(' ', '')
+                    pq = p.rsplit('/', 1)
+                    if len(pq) == 2:
+                        p, q = pq
+                        fp = fractions.Fraction(p)
+                        fq = fractions.Fraction(q)
+                        p = fp/fq
+                    try:
+                        p = fractions.Fraction(p)
+                    except ValueError:
+                        pass  # error will raise below
+                    else:
                         return Rational(p.numerator, p.denominator, 1)
-                except NameError:
-                    pass  # error will raise below
 
-                if isinstance(p, (float, Float)):
-                    return Rational(*_as_integer_ratio(p))
+                if not isinstance(p, Rational):
+                    raise TypeError('invalid input: %s' % p)
 
-            if not isinstance(p, SYMPY_INTS + (Rational,)):
-                raise TypeError('invalid input: %s' % p)
-            q = q or S.One
+            q = 1
             gcd = 1
         else:
             p = Rational(p)
@@ -3817,11 +3821,10 @@ def _mpc_(self):
 
 
 def sympify_fractions(f):
-    return Rational(f.numerator, f.denominator)
+    return Rational(f.numerator, f.denominator, 1)
 
 converter[fractions.Fraction] = sympify_fractions
 
-
 try:
     if HAS_GMPY == 2:
         import gmpy2 as gmpy
@@ -3848,6 +3851,13 @@ def sympify_mpmath(x):
 converter[mpnumeric] = sympify_mpmath
 
 
+def sympify_mpq(x):
+    p, q = x._mpq_
+    return Rational(p, q, 1)
+
+converter[type(mpmath.rational.mpq(1, 2))] = sympify_mpq
+
+
 def sympify_complex(a):
     real, imag = list(map(sympify, (a.real, a.imag)))
     return real + S.ImaginaryUnit*imag
diff --git a/sympy/core/tests/test_numbers.py b/sympy/core/tests/test_numbers.py
index ed2682531234..3ac1b824b2eb 100644
--- a/sympy/core/tests/test_numbers.py
+++ b/sympy/core/tests/test_numbers.py
@@ -10,11 +10,13 @@
 from sympy.core.numbers import (igcd, ilcm, igcdex, seterr, _intcache,
     igcd2, igcd_lehmer, mpf_norm, comp, mod_inverse)
 from sympy.core.mod import Mod
+from sympy.polys.domains.groundtypes import PythonRational
 from sympy.utilities.decorator import conserve_mpmath_dps
 from sympy.utilities.iterables import permutations
 from sympy.utilities.pytest import XFAIL, raises
 
 from mpmath import mpf
+from mpmath.rational import mpq
 import mpmath
 
 
@@ -318,6 +320,8 @@ def test_Rational_new():
     assert Rational(1, 3.0) == Rational(1, 3)
     assert Rational(Float(0.5)) == Rational(1, 2)
     assert Rational('1e2/1e-2') == Rational(10000)
+    assert Rational('1 234') == Rational(1234)
+    assert Rational('1/1 234') == Rational(1, 1234)
     assert Rational(-1, 0) == S.ComplexInfinity
     assert Rational(1, 0) == S.ComplexInfinity
     # Make sure Rational doesn't lose precision on Floats
@@ -332,6 +336,9 @@ def test_Rational_new():
     except ImportError:
         pass
 
+    assert Rational(mpq(2, 6)) == Rational(1, 3)
+    assert Rational(PythonRational(2, 6)) == Rational(1, 3)
+
 
 def test_Number_new():
     """"
diff --git a/sympy/core/tests/test_sympify.py b/sympy/core/tests/test_sympify.py
index 9246d2373d6f..e9d63eee47fe 100644
--- a/sympy/core/tests/test_sympify.py
+++ b/sympy/core/tests/test_sympify.py
@@ -16,6 +16,7 @@
 from sympy.external import import_module
 
 import mpmath
+from mpmath.rational import mpq
 
 
 numpy = import_module('numpy')
@@ -116,6 +117,8 @@ def test_sympify_mpmath():
 
     assert sympify(mpmath.mpc(1.0 + 2.0j)) == Float(1.0) + Float(2.0)*I
 
+    assert sympify(mpq(1, 2)) == S.Half
+
 
 def test_sympify2():
     class A:
diff --git a/sympy/integrals/quadrature.py b/sympy/integrals/quadrature.py
index 2e95e024a858..a096f4901d1d 100644
--- a/sympy/integrals/quadrature.py
+++ b/sympy/integrals/quadrature.py
@@ -53,7 +53,7 @@ def gauss_legendre(n, n_digits):
     >>> x
     [-0.86114, -0.33998, 0.33998, 0.86114]
     >>> w
-    [0.34786, 0.65215, 0.65215, 0.34786]
+    [0.34785, 0.65215, 0.65215, 0.34785]
 
     See Also
     ========
diff --git a/sympy/integrals/tests/test_quadrature.py b/sympy/integrals/tests/test_quadrature.py
index 64487dfc5d31..66453547b135 100644
--- a/sympy/integrals/tests/test_quadrature.py
+++ b/sympy/integrals/tests/test_quadrature.py
@@ -466,7 +466,7 @@ def test_jacobi():
             '0.65486073394528506',
             '0.95949297361449739']
     assert [str(r) for r in w] == [
-            '0.090675770007435371',
+            '0.090675770007435372',
             '0.33391416373675607',
             '0.65248870981926643',
             '0.94525424081394926',
diff --git a/sympy/polys/polyroots.py b/sympy/polys/polyroots.py
index 038415e2a6e0..fa6b0b536504 100644
--- a/sympy/polys/polyroots.py
+++ b/sympy/polys/polyroots.py
@@ -720,6 +720,7 @@ def preprocess_roots(poly):
     """Try to get rid of symbolic coefficients from ``poly``. """
     coeff = S.One
 
+    poly_func = poly.func
     try:
         _, poly = poly.clear_denoms(convert=True)
     except DomainError:
@@ -783,6 +784,8 @@ def func(k, coeff):
             poly = poly.termwise(func)
             coeff *= basis
 
+    if not isinstance(poly, poly_func):
+        poly = poly_func(poly)
     return coeff, poly
 
 
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index d5ac051f636a..3cb472eedfa5 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -3553,6 +3553,7 @@ def nroots(f, n=15, maxsteps=50, cleanup=True):
         [-1.73205080756887729352744634151, 1.73205080756887729352744634151]
 
         """
+        from sympy.functions.elementary.complexes import sign
         if f.is_multivariate:
             raise MultivariatePolynomialError(
                 "can't compute numerical roots of %s" % f)
@@ -3591,7 +3592,7 @@ def nroots(f, n=15, maxsteps=50, cleanup=True):
             # Mpmath puts real roots first, then complex ones (as does all_roots)
             # so we make sure this convention holds here, too.
             roots = list(map(sympify,
-                sorted(roots, key=lambda r: (1 if r.imag else 0, r.real, r.imag))))
+                sorted(roots, key=lambda r: (1 if r.imag else 0, r.real, abs(r.imag), sign(r.imag)))))
         except NoConvergence:
             raise NoConvergence(
                 'convergence to root failed; try n < %s or maxsteps > %s' % (
diff --git a/sympy/polys/rootisolation.py b/sympy/polys/rootisolation.py
index f697dc98e580..733cd43c7747 100644
--- a/sympy/polys/rootisolation.py
+++ b/sympy/polys/rootisolation.py
@@ -565,7 +565,7 @@ def dup_isolate_real_roots_list(polys, K, eps=None, inf=None, sup=None, strict=F
        2. Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. Vigklas: Improving the Performance
        of the Continued Fractions Method Using New Bounds of Positive Roots.
        Nonlinear Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008.
-"""
+    """
     if K.is_QQ:
         K, F, polys = K.get_ring(), K, polys[:]
 
@@ -1724,7 +1724,11 @@ def __repr__(self):
 
     def is_disjoint(self, other):
         """Return ``True`` if two isolation intervals are disjoint. """
-        return (self.b <= other.a or other.b <= self.a)
+        if isinstance(other, RealInterval):
+            return (self.b <= other.a or other.b <= self.a)
+        assert isinstance(other, ComplexInterval)
+        return (self.b <= other.ax or other.bx <= self.a
+            or other.ay*other.by > 0)
 
     def _inner_refine(self):
         """Internal one step real root refinement procedure. """
@@ -1764,14 +1768,146 @@ def refine(self):
         """Perform one step of real root refinement algorithm. """
         return self._inner_refine()
 
+
 class ComplexInterval(object):
     """A fully qualified representation of a complex isolation interval.
-    The printed form is shown as (x1, y1) x (x2, y2): the southwest x northeast
-    coordinates of the interval's rectangle."""
+    The printed form is shown as (ax, bx) x (ay, by) where (ax, ay)
+    and (bx, by) are the coordinates of the southwest and northeast
+    corners of the interval's rectangle, respectively.
+
+    Examples
+    ========
+
+    >>> from sympy import RootOf, Rational, S
+    >>> from sympy.abc import x
+    >>> root = RootOf(x**10 - 2*x + 3, 9)
+    >>> i = root._get_interval(); i
+    (3/64, 3/32) x (9/8, 75/64)
+
+    The real part of the root lies within the range [0, 3/4] while
+    the imaginary part lies within the range [9/8, 3/2]:
+
+    >>> root.n(3)
+    0.0766 + 1.14*I
+
+    The width of the ranges in the x and y directions on the complex
+    plane are:
+
+    >>> i.dx, i.dy
+    (3/64, 3/64)
+
+    The center of the range is
+
+    >>> i.center
+    (9/128, 147/128)
+
+    The northeast coordinate of the rectangle bounding the root in the
+    complex plane is given by attribute b and the x and y components
+    are accessed by bx and by:
+
+    >>> i.b, i.bx, i.by
+    ((3/32, 75/64), 3/32, 75/64)
+
+    The southwest coordinate is similarly given by i.a
+
+    >>> i.a, i.ax, i.ay
+    ((3/64, 9/8), 3/64, 9/8)
+
+    Although the interval prints to show only the real and imaginary
+    range of the root, all the information of the underlying root
+    is contained as properties of the interval.
+
+    For example, an interval with a nonpositive imaginary range is
+    considered to be the conjugate. Since the y values of y are in the
+    range [0, 1/4] it is not the conjugate:
+
+    >>> i.conj
+    False
+
+    The conjugate's interval is
+
+    >>> ic = i.conjugate(); ic
+    (3/64, 3/32) x (-75/64, -9/8)
+
+        NOTE: the values printed still represent the x and y range
+        in which the root -- conjugate, in this case -- is located,
+        but the underlying a and b values of a root and its conjugate
+        are the same:
+
+        >>> assert i.a == ic.a and i.b == ic.b
+
+        What changes are the reported coordinates of the bounding rectangle:
+
+        >>> (i.ax, i.ay), (i.bx, i.by)
+        ((3/64, 9/8), (3/32, 75/64))
+        >>> (ic.ax, ic.ay), (ic.bx, ic.by)
+        ((3/64, -75/64), (3/32, -9/8))
+
+    The interval can be refined once:
+
+    >>> i  # for reference, this is the current interval
+    (3/64, 3/32) x (9/8, 75/64)
+
+    >>> i.refine()
+    (3/64, 3/32) x (9/8, 147/128)
+
+    Several refinement steps can be taken:
+
+    >>> i.refine_step(2)  # 2 steps
+    (9/128, 3/32) x (9/8, 147/128)
+
+    It is also possible to refine to a given tolerance:
+
+    >>> tol = min(i.dx, i.dy)/2
+    >>> i.refine_size(tol)
+    (9/128, 21/256) x (9/8, 291/256)
+
+    A disjoint interval is one whose bounding rectangle does not
+    overlap with another. An interval, necessarily, is not disjoint with
+    itself, but any interval is disjoint with a conjugate since the
+    conjugate rectangle will always be in the lower half of the complex
+    plane and the non-conjugate in the upper half:
+
+    >>> i.is_disjoint(i), i.is_disjoint(i.conjugate())
+    (False, True)
+
+    The following interval j is not disjoint from i:
+
+    >>> close = RootOf(x**10 - 2*x + 300/S(101), 9)
+    >>> j = close._get_interval(); j
+    (75/1616, 75/808) x (225/202, 1875/1616)
+    >>> i.is_disjoint(j)
+    False
+
+    The two can be made disjoint, however:
+
+    >>> newi, newj = i.refine_disjoint(j)
+    >>> newi
+    (39/512, 159/2048) x (2325/2048, 4653/4096)
+    >>> newj
+    (3975/51712, 2025/25856) x (29325/25856, 117375/103424)
+
+    Even though the real ranges overlap, the imaginary do not, so
+    the roots have been resolved as distinct. Intervals are disjoint
+    when either the real or imaginary component of the intervals is
+    distinct. In the case above, the real components have not been
+    resolved (so we don't know, yet, which root has the smaller real
+    part) but the imaginary part of ``close`` is larger than ``root``:
+
+    >>> close.n(3)
+    0.0771 + 1.13*I
+    >>> root.n(3)
+    0.0766 + 1.14*I
+    """
 
     def __init__(self, a, b, I, Q, F1, F2, f1, f2, dom, conj=False):
         """Initialize new complex interval with complete information. """
-        self.a, self.b = a, b  # the southwest and northeast corner: (x1, y1), (x2, y2)
+        # a and b are the SW and NE corner of the bounding interval,
+        # (ax, ay) and (bx, by), respectively, for the NON-CONJUGATE
+        # root (the one with the positive imaginary part); when working
+        # with the conjugate, the a and b value are still non-negative
+        # but the ay, by are reversed and have oppositite sign
+        self.a, self.b = a, b
         self.I, self.Q = I, Q
 
         self.f1, self.F1 = f1, F1
@@ -1780,6 +1916,16 @@ def __init__(self, a, b, I, Q, F1, F2, f1, f2, dom, conj=False):
         self.dom = dom
         self.conj = conj
 
+    @property
+    def args(self):
+        i = self
+        return (i.a, i.b, i.I, i.Q, i.F1, i.F2, i.f1, i.f2, i.dom, i.conj)
+
+    def __eq__(self, other):
+        if type(other) != type(self):
+            return False
+        return self.args == other.args
+
     @property
     def ax(self):
         """Return ``x`` coordinate of south-western corner. """
@@ -1822,7 +1968,8 @@ def center(self):
         return ((self.ax + self.bx)/2, (self.ay + self.by)/2)
 
     def as_tuple(self):
-        """Return tuple representation of complex isolating interval. """
+        """Return tuple representation of the complex isolating
+        interval's SW and NE corners, respectively. """
         return ((self.ax, self.ay), (self.bx, self.by))
 
     def __repr__(self):
@@ -1835,7 +1982,9 @@ def conjugate(self):
 
     def is_disjoint(self, other):
         """Return ``True`` if two isolation intervals are disjoint. """
-        if self.conj != other.conj:
+        if isinstance(other, RealInterval):
+            return other.is_disjoint(self)
+        if self.conj != other.conj:  # above and below real axis
             return True
         re_distinct = (self.bx <= other.ax or other.bx <= self.ax)
         if re_distinct:
@@ -1843,6 +1992,14 @@ def is_disjoint(self, other):
         im_distinct = (self.by <= other.ay or other.by <= self.ay)
         return im_distinct
 
+    def _set_imaginary(self):
+        a1, a2 = self.a
+        b1, b2 = self.b
+        a = a1 - a1, a2
+        b = b1 - b1, b2
+        return ComplexInterval(a, b, self.I, self.Q,
+            self.F1, self.F2, self.f1, self.f2, self.dom, self.conj)
+
     def _inner_refine(self):
         """Internal one step complex root refinement procedure. """
         (u, v), (s, t) = self.a, self.b
diff --git a/sympy/polys/rootoftools.py b/sympy/polys/rootoftools.py
index 9e69f3ff3634..6c4572b05012 100644
--- a/sympy/polys/rootoftools.py
+++ b/sympy/polys/rootoftools.py
@@ -2,8 +2,8 @@
 
 from __future__ import print_function, division
 
-from sympy.core import (S, Expr, Integer, Float, I, Add, Lambda, symbols,
-        sympify, Rational, Dummy)
+from sympy.core import (S, Expr, Integer, Float, I, oo, Add, Lambda,
+    symbols, sympify, Rational, Dummy)
 from sympy.core.cache import cacheit
 from sympy.core.function import AppliedUndef
 from sympy.functions.elementary.miscellaneous import root as _root
@@ -29,22 +29,84 @@
 from sympy.polys.domains import QQ
 
 from mpmath import mpf, mpc, findroot, workprec
-from mpmath.libmp.libmpf import prec_to_dps
+from mpmath.libmp.libmpf import dps_to_prec, prec_to_dps
 
 from sympy.utilities import lambdify, public, sift
 
-from sympy.core.compatibility import range
+from sympy.core.compatibility import range, ordered
 
 from math import log as mathlog
 
 __all__ = ['CRootOf']
 
-def _ispow2(i):
-    v = mathlog(i, 2)
-    return v == int(v)
 
-_reals_cache = {}
-_complexes_cache = {}
+
+class _pure_key_dict(dict):
+    """A minimal dictionary that makes sure that the key is PurePoly.
+
+    Examples
+    ========
+
+    Only the following four actions are guaranteed:
+
+    >>> from sympy.polys.rootoftools import _pure_key_dict
+    >>> from sympy import S, PurePoly
+    >>> from sympy.abc import x, y
+
+    1) creation
+
+    >>> P = _pure_key_dict()
+
+    2) assignment for a PurePoly or univariate polynomial
+
+    >>> P[x] = 1
+    >>> P[PurePoly(x - y, x)] = 2
+
+    3) retrieval based on PurePoly key comparison (use this
+       instead of the get method)
+
+    >>> P[y]
+    1
+
+    4) KeyError when trying to retrieve a nonexisting key
+
+    >>> P[y + 1]
+    Traceback (most recent call last):
+    ...
+    KeyError: PurePoly(y + 1, y, domain='ZZ')
+    """
+    def __getitem__(self, k):
+        if not isinstance(k, PurePoly):
+            assert len(k.free_symbols) == 1
+            k = PurePoly(k, expand=False)
+        return self.__dict__[k]
+
+    def __setitem__(self, k, v):
+        if not isinstance(k, PurePoly):
+            assert len(k.free_symbols) == 1
+            k = PurePoly(k, expand=False)
+        self.__dict__[k] = v
+
+_reals_cache = _pure_key_dict()
+_complexes_cache = _pure_key_dict()
+
+
+def _pure_factors(poly):
+    _, factors = poly.factor_list()
+    return [(PurePoly(f, expand=False), m) for f, m in factors]
+
+
+def _imag_count_of_factor(f):
+    """Return the number of imaginary roots for irreducible
+    univariate polynomial ``f``.
+    """
+    terms = [(i, j) for (i,), j in f.terms()]
+    if any(i % 2 for i, j in terms):
+        return 0
+    # update signs
+    even = [(i, I**i*j) for i, j in terms]
+    even = Poly.from_dict(dict(even), Dummy('x'))
+    return int(even.count_roots(-oo, oo))
 
 
 @public
@@ -55,7 +117,8 @@ def rootof(f, x, index=None, radicals=True, expand=True):
     expression involving radicals.
 
     Parameters
-    ----------
+    ==========
+
     f : Expr
         Univariate polynomial.
     x : Symbol, optional
@@ -91,6 +154,77 @@ class ComplexRootOf(RootOf):
     real or complex intervals and indexed in a fixed order.
     Currently only rational coefficients are allowed.
     Can be imported as ``CRootOf``.
+
+
+    Examples
+    ========
+
+    >>> from sympy import CRootOf, rootof
+    >>> from sympy.abc import x
+
+    CRootOf is a way to reference a particular root of a
+    polynomial. If there is a rational root, it will be returned:
+
+    >>> CRootOf(x**2 - 4, 0)
+    -2
+
+    Whether roots involving radicals are returned or not
+    depends on whether the ``radicals`` flag is true (which is
+    set to True with rootof):
+
+    >>> CRootOf(x**2 - 3, 0)
+    CRootOf(x**2 - 3, 0)
+    >>> CRootOf(x**2 - 3, 0, radicals=True)
+    -sqrt(3)
+    >>> rootof(x**2 - 3, 0)
+    -sqrt(3)
+
+    The following cannot be expressed in terms of radicals:
+
+    >>> r = rootof(4*x**5 + 16*x**3 + 12*x**2 + 7, 0); r
+    CRootOf(4*x**5 + 16*x**3 + 12*x**2 + 7, 0)
+
+    The root bounds can be seen, however, and they are used by the
+    evaluation methods to get numerical approximations for the root.
+
+    >>> interval = r._get_interval(); interval
+    (-1, 0)
+    >>> r.evalf(2)
+    -0.98
+
+    The evalf method refines the width of the root bounds until it
+    guarantees that any decimal approximation within those bounds
+    will satisfy the desired precision. It then stores the refined
+    interval so subsequent requests at or below the requested
+    precision will not have to recompute the root bounds and will
+    return very quickly.
+
+    Before evaluation above, the interval was
+
+    >>> interval
+    (-1, 0)
+
+    After evaluation it is now
+
+    >>. r._get_interval()
+    (-165/169, -206/211)
+
+    The `eval_approx` method will also find the root to a given
+    precision but the interval is not modified unless the search
+    for the root fails to converge within the root bounds.
+
+    To reset all intervals for a given polynomial, the `_reset` method
+    can be called from any CRootOf instance of the polynomial:
+
+    >>> r._reset()
+    >>> r._get_interval()
+    (-1, 0)
+
+    See Also
+    ========
+    eval_approx
+    eval_rational
+    _eval_evalf
     """
 
     __slots__ = ['index']
@@ -181,15 +315,22 @@ def args(self):
 
     @property
     def free_symbols(self):
-        # CRootOf currently only works with univariate expressions and although
-        # the poly attribute is often a PurePoly, sometimes it is a Poly. In
-        # either case no free symbols should be reported.
+        # CRootOf currently only works with univariate expressions
+        # whose poly attribute should be a PurePoly with no free
+        # symbols
         return set()
 
     def _eval_is_real(self):
         """Return ``True`` if the root is real. """
         return self.index < len(_reals_cache[self.poly])
 
+    def _eval_is_imaginary(self):
+        """Return ``True`` if the root is real. """
+        if self.index >= len(_reals_cache[self.poly]):
+            ivl = self._get_interval()
+            return ivl.ax*ivl.bx <= 0  # all others are on one side or the other
+        return False  # XXX is this necessary?
+
     @classmethod
     def real_roots(cls, poly, radicals=True):
         """Get real roots of a polynomial. """
@@ -224,25 +365,41 @@ def _get_complexes_sqf(cls, factor):
         return complex_part
 
     @classmethod
-    def _get_reals(cls, factors):
+    def _get_reals(cls, factors, use_cache=True):
         """Compute real root isolating intervals for a list of factors. """
         reals = []
 
         for factor, k in factors:
-            real_part = cls._get_reals_sqf(factor)
-            reals.extend([(root, factor, k) for root in real_part])
+            try:
+                if not use_cache:
+                    raise KeyError
+                r = _reals_cache[factor]
+                reals.extend([(i, factor, k) for i in r])
+            except KeyError:
+                real_part = cls._get_reals_sqf(factor)
+                new = [(root, factor, k) for root in real_part]
+                reals.extend(new)
 
+        reals = cls._reals_sorted(reals)
         return reals
 
     @classmethod
-    def _get_complexes(cls, factors):
+    def _get_complexes(cls, factors, use_cache=True):
         """Compute complex root isolating intervals for a list of factors. """
         complexes = []
 
-        for factor, k in factors:
-            complex_part = cls._get_complexes_sqf(factor)
-            complexes.extend([(root, factor, k) for root in complex_part])
+        for factor, k in ordered(factors):
+            try:
+                if not use_cache:
+                    raise KeyError
+                c = _complexes_cache[factor]
+                complexes.extend([(i, factor, k) for i in c])
+            except KeyError:
+                complex_part = cls._get_complexes_sqf(factor)
+                new = [(root, factor, k) for root in complex_part]
+                complexes.extend(new)
 
+        complexes = cls._complexes_sorted(complexes)
         return complexes
 
     @classmethod
@@ -271,140 +428,95 @@ def _reals_sorted(cls, reals):
         return reals
 
     @classmethod
-    def _separate_imaginary_from_complex(cls, complexes):
-
-        def is_imag(c):
-            '''
-            return True if all roots are imaginary (ax**2 + b)
-            return False if no roots are imaginary
-            return None if 2 roots are imaginary (ax**N'''
-            u, f, k = c
-            deg = f.degree()
-            if f.length() == 2:
-                if deg == 2:
-                    return True  # both imag
-                elif _ispow2(deg):
-                    if f.LC()*f.TC() < 0:
-                        return None  # 2 are imag
-            return False  # none are imag
-        # separate according to the function
+    def _refine_imaginary(cls, complexes):
         sifted = sift(complexes, lambda c: c[1])
-        del complexes
-        imag = []
         complexes = []
-        for f in sifted:
-            isift = sift(sifted[f], lambda c: is_imag(c))
-            imag.extend(isift.pop(True, []))
-            complexes.extend(isift.pop(False, []))
-            mixed = isift.pop(None, [])
-            assert not isift
-            if not mixed:
-                continue
-            while True:
-                # the non-imaginary ones will be on one side or the other
-                # of the y-axis
-                i = 0
-                while i < len(mixed):
-                    u, f, k = mixed[i]
-                    if u.ax*u.bx > 0:
-                        complexes.append(mixed.pop(i))
-                    else:
-                        i += 1
-                if len(mixed) == 2:
-                    imag.extend(mixed)
-                    break
-                # refine
-                for i, (u, f, k) in enumerate(mixed):
-                    u = u._inner_refine()
-                    mixed[i] = u, f, k
-        return imag, complexes
+        for f in ordered(sifted):
+            nimag = _imag_count_of_factor(f)
+            if nimag == 0:
+                # refine until xbounds are neg or pos
+                for u, f, k in sifted[f]:
+                    while u.ax*u.bx <= 0:
+                        u = u._inner_refine()
+                    complexes.append((u, f, k))
+            else:
+                # refine until all but nimag xbounds are neg or pos
+                potential_imag = list(range(len(sifted[f])))
+                while True:
+                    assert len(potential_imag) > 1
+                    for i in list(potential_imag):
+                        u, f, k = sifted[f][i]
+                        if u.ax*u.bx > 0:
+                            potential_imag.remove(i)
+                        elif u.ax != u.bx:
+                            u = u._inner_refine()
+                            sifted[f][i] = u, f, k
+                    if len(potential_imag) == nimag:
+                        break
+                # set real part to 0 in what are now
+                # known to be the imaginaries
+                for i in potential_imag:
+                    u, f, k = sifted[f][i]
+                    u = u._set_imaginary()
+                    sifted[f][i] = u, f, k
+                complexes.extend(sifted[f])
+        return complexes
 
     @classmethod
     def _refine_complexes(cls, complexes):
         """return complexes such that no bounding rectangles of non-conjugate
-        roots would intersect if slid horizontally or vertically/
+        roots would intersect. In addition, assure that neither ay nor by is
+        0 to guarantee that non-real roots are distinct from real roots in
+        terms of the y-bounds.
         """
-        while complexes:  # break when all are distinct
-            # get the intervals pairwise-disjoint.
-            # If rectangles were drawn around the coordinates of the bounding
-            # rectangles, no rectangles would intersect after this procedure.
-            for i, (u, f, k) in enumerate(complexes):
-                for j, (v, g, m) in enumerate(complexes[i + 1:]):
-                    u, v = u.refine_disjoint(v)
-                    complexes[i + j + 1] = (v, g, m)
-
-                complexes[i] = (u, f, k)
-            # Although there are no intersecting rectangles, a given rectangle
-            # might intersect another when slid horizontally. We have to refine
-            # intervals until this is not true so we can sort the roots
-            # unambiguously. Since complex roots come in conjugate pairs, we
-            # will always have 2 rectangles above each other but we should not
-            # have more than that.
-            N = len(complexes)//2 - 1
-            # check x (real) parts: there must be N + 1 disjoint x ranges, i.e.
-            # the first one must be different from N others
-            uu = set([(u.ax, u.bx) for u, _, _ in complexes])
-            u = uu.pop()
-            if sum([u[1] <= v[0] or v[1] <= u[0] for v in uu]) < N:
-                # refine
-                for i, (u, f, k) in enumerate(complexes):
-                    u = u._inner_refine()
-                    complexes[i] = u, f, k
-            else:
-                # intervals with identical x-values have disjoint y-values or
-                # else they would not be disjoint so there is no need for
-                # further checks
-                break
+        # get the intervals pairwise-disjoint.
+        # If rectangles were drawn around the coordinates of the bounding
+        # rectangles, no rectangles would intersect after this procedure.
+        for i, (u, f, k) in enumerate(complexes):
+            for j, (v, g, m) in enumerate(complexes[i + 1:]):
+                u, v = u.refine_disjoint(v)
+                complexes[i + j + 1] = (v, g, m)
+
+            complexes[i] = (u, f, k)
+
+        # refine until the x-bounds are unambiguously positive or negative
+        # for non-imaginary roots
+        complexes = cls._refine_imaginary(complexes)
+
+        # make sure that all y bounds are off the real axis
+        # and on the same side of the axis
+        for i, (u, f, k) in enumerate(complexes):
+            while u.ay*u.by <= 0:
+                u = u.refine()
+            complexes[i] = u, f, k
         return complexes
 
     @classmethod
     def _complexes_sorted(cls, complexes):
         """Make complex isolating intervals disjoint and sort roots. """
-        if not complexes:
-            return []
-        cache = {}
-
-        # imaginary roots can cause a problem in terms of sorting since
-        # their x-intervals will never refine as distinct from others
-        # so we handle them separately
-        imag, complexes = cls._separate_imaginary_from_complex(complexes)
         complexes = cls._refine_complexes(complexes)
-
-        # sort imaginary roots
-        def key(c):
-            '''return, for ax**n+b, +/-root(abs(b/a), b) according to the
-            apparent sign of the imaginary interval, e.g. if the interval
-            were (0, 3) the positive root would be returned.
-            '''
-            u, f, k = c
-            r = _root(abs(f.TC()/f.LC()), f.degree())
-            if u.ay < 0 or u.by < 0:
-                return -r
-            return r
-        imag = sorted(imag, key=lambda c: key(c))
-
-        # sort complexes and combine with imag
-        if complexes:
-            # key is (x1, y1) e.g. (1, 2)x(3, 4) -> (1,3)
-            complexes = sorted(complexes, key=lambda c: c[0].a)
-            # find insertion point for imaginary
-            for i, c in enumerate(reversed(complexes)):
-                if c[0].bx <= 0:
-                    break
-            i = len(complexes) - i - 1
-            if i:
-                i += 1
-            complexes = complexes[:i] + imag + complexes[i:]
-        else:
-            complexes = imag
+        # XXX don't sort until you are sure that it is compatible
+        # with the indexing method but assert that the desired state
+        # is not broken
+        C, F = 0, 1  # location of ComplexInterval and factor
+        fs = set([i[F] for i in complexes])
+        for i in range(1, len(complexes)):
+            if complexes[i][F] != complexes[i - 1][F]:
+                # if this fails the factors of a root were not
+                # contiguous because a discontinuity should only
+                # happen once
+                fs.remove(complexes[i - 1][F])
+        for i in range(len(complexes)):
+            # negative im part (conj=True) comes before
+            # positive im part (conj=False)
+            assert complexes[i][C].conj is (i % 2 == 0)
 
         # update cache
+        cache = {}
+        # -- collate
         for root, factor, _ in complexes:
-            if factor in cache:
-                cache[factor].append(root)
-            else:
-                cache[factor] = [root]
-
+            cache.setdefault(factor, []).append(root)
+        # -- store
         for factor, roots in cache.items():
             _complexes_cache[factor] = roots
 
@@ -436,8 +548,7 @@ def _complexes_index(cls, complexes, index):
         Map initial complex root index to an index in a factor where
         the root belongs.
         """
-        index, i = index, 0
-
+        i = 0
         for j, (_, factor, k) in enumerate(complexes):
             if index < i + k:
                 poly, index = factor, 0
@@ -460,26 +571,23 @@ def _count_roots(cls, roots):
     @classmethod
     def _indexed_root(cls, poly, index):
         """Get a root of a composite polynomial by index. """
-        (_, factors) = poly.factor_list()
+        factors = _pure_factors(poly)
 
         reals = cls._get_reals(factors)
         reals_count = cls._count_roots(reals)
 
         if index < reals_count:
-            reals = cls._reals_sorted(reals)
             return cls._reals_index(reals, index)
         else:
             complexes = cls._get_complexes(factors)
-            complexes = cls._complexes_sorted(complexes)
             return cls._complexes_index(complexes, index - reals_count)
 
     @classmethod
     def _real_roots(cls, poly):
         """Get real roots of a composite polynomial. """
-        (_, factors) = poly.factor_list()
+        factors = _pure_factors(poly)
 
         reals = cls._get_reals(factors)
-        reals = cls._reals_sorted(reals)
         reals_count = cls._count_roots(reals)
 
         roots = []
@@ -489,13 +597,15 @@ def _real_roots(cls, poly):
 
         return roots
 
+    def _reset(self):
+        self._all_roots(self.poly, use_cache=False)
+
     @classmethod
-    def _all_roots(cls, poly):
+    def _all_roots(cls, poly, use_cache=True):
         """Get real and complex roots of a composite polynomial. """
-        (_, factors) = poly.factor_list()
+        factors = _pure_factors(poly)
 
-        reals = cls._get_reals(factors)
-        reals = cls._reals_sorted(reals)
+        reals = cls._get_reals(factors, use_cache=use_cache)
         reals_count = cls._count_roots(reals)
 
         roots = []
@@ -503,8 +613,7 @@ def _all_roots(cls, poly):
         for index in range(0, reals_count):
             roots.append(cls._reals_index(reals, index))
 
-        complexes = cls._get_complexes(factors)
-        complexes = cls._complexes_sorted(complexes)
+        complexes = cls._get_complexes(factors, use_cache=use_cache)
         complexes_count = cls._count_roots(complexes)
 
         for index in range(0, complexes_count):
@@ -591,26 +700,35 @@ def _eval_subs(self, old, new):
         # don't allow subs to change anything
         return self
 
-    def _eval_evalf(self, prec):
-        """Evaluate this complex root to the given precision. """
+    def _eval_conjugate(self):
+        if self.is_real:
+            return self
+        expr, i = self.args
+        return self.func(expr, i + (1 if self._get_interval().conj else -1))
+
+    def eval_approx(self, n):
+        """Evaluate this complex root to the given precision.
+
+        This uses secant method and root bounds are used to both
+        generate an initial guess and to check that the root
+        returned is valid. If ever the method converges outside the
+        root bounds, the bounds will be made smaller and updated.
+        """
+        prec = dps_to_prec(n)
         with workprec(prec):
             g = self.poly.gen
             if not g.is_Symbol:
                 d = Dummy('x')
+                if self.is_imaginary:
+                    d *= I
                 func = lambdify(d, self.expr.subs(g, d))
             else:
-                func = lambdify(g, self.expr)
+                expr = self.expr
+                if self.is_imaginary:
+                    expr = self.expr.subs(g, I*g)
+                func = lambdify(g, expr)
 
             interval = self._get_interval()
-            if not self.is_real:
-                # For complex intervals, we need to keep refining until the
-                # imaginary interval is disjunct with other roots, that is,
-                # until both ends get refined.
-                ay = interval.ay
-                by = interval.by
-                while interval.ay == ay or interval.by == by:
-                    interval = interval.refine()
-
             while True:
                 if self.is_real:
                     a = mpf(str(interval.a))
@@ -619,31 +737,29 @@ def _eval_evalf(self, prec):
                         root = a
                         break
                     x0 = mpf(str(interval.center))
+                    x1 = x0 + mpf(str(interval.dx))/4
+                elif self.is_imaginary:
+                    a = mpf(str(interval.ay))
+                    b = mpf(str(interval.by))
+                    if a == b:
+                        root = mpc(mpf('0'), a)
+                        break
+                    x0 = mpf(str(interval.center[1]))
+                    x1 = x0 + mpf(str(interval.dy))/4
                 else:
                     ax = mpf(str(interval.ax))
                     bx = mpf(str(interval.bx))
                     ay = mpf(str(interval.ay))
                     by = mpf(str(interval.by))
                     if ax == bx and ay == by:
-                        # the sign of the imaginary part will be assigned
-                        # according to the desired index using the fact that
-                        # roots are sorted with negative imag parts coming
-                        # before positive (and all imag roots coming after real
-                        # roots)
-                        deg = self.poly.degree()
-                        i = self.index  # a positive attribute after creation
-                        if (deg - i) % 2:
-                            if ay < 0:
-                                ay = -ay
-                        else:
-                            if ay > 0:
-                                ay = -ay
                         root = mpc(ax, ay)
                         break
                     x0 = mpc(*map(str, interval.center))
-
+                    x1 = x0 + mpc(*map(str, (interval.dx, interval.dy)))/4
                 try:
-                    root = findroot(func, x0)
+                    # without a tolerance, this will return when (to within
+                    # the given precision) x_i == x_{i-1}
+                    root = findroot(func, (x0, x1))
                     # If the (real or complex) root is not in the 'interval',
                     # then keep refining the interval. This happens if findroot
                     # accidentally finds a different root outside of this
@@ -658,8 +774,11 @@ def _eval_evalf(self, prec):
                     # successful iterations to process (in which case it
                     # will fail to initialize a variable that is tested
                     # after the iterations and raise an UnboundLocalError).
-                    if self.is_real:
-                        if (a <= root <= b):
+                    if self.is_real or self.is_imaginary:
+                        if not bool(root.imag) == self.is_real and (
+                                a <= root <= b):
+                            if self.is_imaginary:
+                                root = mpc(mpf('0'), root.real)
                             break
                     elif (ax <= root.real <= bx and ay <= root.imag <= by):
                         break
@@ -667,40 +786,98 @@ def _eval_evalf(self, prec):
                     pass
                 interval = interval.refine()
 
-        return (Float._new(root.real._mpf_, prec)
-                + I*Float._new(root.imag._mpf_, prec))
-
-    def eval_rational(self, tol):
-        """
-        Return a Rational approximation to ``self`` with the tolerance ``tol``.
+        # update the interval so we at least (for this precision or
+        # less) don't have much work to do to recompute the root
+        self._set_interval(interval)
+        return (Float._new(root.real._mpf_, prec) +
+            I*Float._new(root.imag._mpf_, prec))
 
-        This method uses bisection, which is very robust and it will always
-        converge. The returned Rational instance will be at most 'tol' from the
-        exact root.
+    def _eval_evalf(self, prec, **kwargs):
+        """Evaluate this complex root to the given precision."""
+        # all kwargs are ignored
+        return self.eval_rational(n=prec_to_dps(prec))._evalf(prec)
 
-        The following example first obtains Rational approximation to 1e-7
-        accuracy for all roots of the 4-th order Legendre polynomial, and then
-        evaluates it to 5 decimal digits (so all digits will be correct
-        including rounding):
+    def eval_rational(self, dx=None, dy=None, n=15):
+        """
+        Return a Rational approximation of ``self`` that has real
+        and imaginary component approximations that are within ``dx``
+        and ``dy`` of the true values, respectively. Alternatively,
+        ``n`` digits of precision can be specified.
+
+        The interval is refined with bisection and is sure to
+        converge. The root bounds are updated when the refinement
+        is complete so recalculation at the same or lesser precision
+        will not have to repeat the refinement and should be much
+        faster.
+
+        The following example first obtains Rational approximation to
+        1e-8 accuracy for all roots of the 4-th order Legendre
+        polynomial. Since the roots are all less than 1, this will
+        ensure the decimal representation of the approximation will be
+        correct (including rounding) to 6 digits:
 
         >>> from sympy import S, legendre_poly, Symbol
         >>> x = Symbol("x")
         >>> p = legendre_poly(4, x, polys=True)
-        >>> roots = [r.eval_rational(S(1)/10**7) for r in p.real_roots()]
-        >>> roots = [str(r.n(5)) for r in roots]
-        >>> roots
-        ['-0.86114', '-0.33998', '0.33998', '0.86114']
+        >>> r = p.real_roots()[-1]
+        >>> r.eval_rational(10**-8).n(6)
+        0.861136
 
-        """
+        It is not necessary to a two-step calculation, however: the
+        decimal representation can be computed directly:
 
-        if not self.is_real:
-            raise NotImplementedError(
-                "eval_rational() only works for real polynomials so far")
-        func = lambdify(self.poly.gen, self.expr)
+        >>> r.evalf(17)
+        0.86113631159405258
+
+        """
+        dy = dy or dx
+        if dx:
+            rtol = None
+            dx = dx if isinstance(dx, Rational) else Rational(str(dx))
+            dy = dy if isinstance(dy, Rational) else Rational(str(dy))
+        else:
+            # 5 binary (or 2 decimal) digits are needed to ensure that
+            # a given digit is correctly rounded
+            # prec_to_dps(dps_to_prec(n) + 5) - n <= 2 (tested for
+            # n in range(1000000)
+            rtol = S(10)**-(n + 2)  # +2 for guard digits
         interval = self._get_interval()
-        a = Rational(str(interval.a))
-        b = Rational(str(interval.b))
-        return bisect(func, a, b, tol)
+        while True:
+            if self.is_real:
+                if rtol:
+                    dx = abs(interval.center*rtol)
+                interval = interval.refine_size(dx=dx)
+                c = interval.center
+                real = Rational(c)
+                imag = S.Zero
+                if not rtol or interval.dx < abs(c*rtol):
+                    break
+            elif self.is_imaginary:
+                if rtol:
+                    dy = abs(interval.center[1]*rtol)
+                    dx = 1
+                interval = interval.refine_size(dx=dx, dy=dy)
+                c = interval.center[1]
+                imag = Rational(c)
+                real = S.Zero
+                if not rtol or interval.dy < abs(c*rtol):
+                    break
+            else:
+                if rtol:
+                    dx = abs(interval.center[0]*rtol)
+                    dy = abs(interval.center[1]*rtol)
+                interval = interval.refine_size(dx, dy)
+                c = interval.center
+                real, imag = map(Rational, c)
+                if not rtol or (
+                        interval.dx < abs(c[0]*rtol) and
+                        interval.dy < abs(c[1]*rtol)):
+                    break
+
+        # update the interval so we at least (for this precision or
+        # less) don't have much work to do to recompute the root
+        self._set_interval(interval)
+        return real + I*imag
 
     def _eval_Eq(self, other):
         # CRootOf represents a Root, so if other is that root, it should set
@@ -719,32 +896,22 @@ def _eval_Eq(self, other):
             return S.false
         o = other.is_real, other.is_imaginary
         s = self.is_real, self.is_imaginary
-        if o != s and None not in o and None not in s:
+        assert None not in s  # this is part of initial refinement
+        if o != s and None not in o:
             return S.false
-        i = self._get_interval()
-        was = i.a, i.b
-        need = [True]*2
-        # make sure it would be distinct from others
-        while any(need):
-            i = i.refine()
-            a, b = i.a, i.b
-            if need[0] and a != was[0]:
-                need[0] = False
-            if need[1] and b != was[1]:
-                need[1] = False
         re, im = other.as_real_imag()
-        if not im:
-            if self.is_real:
-                a, b = [Rational(str(i)) for i in (a, b)]
-                return sympify(a < other and other < b)
-            return S.false
         if self.is_real:
-            return S.false
-        z = r1, r2, i1, i2 = [Rational(str(j)) for j in (
+            if im:
+                return S.false
+            i = self._get_interval()
+            a, b = [Rational(str(_)) for _ in (i.a, i.b)]
+            return sympify(a <= other and other <= b)
+        i = self._get_interval()
+        r1, r2, i1, i2 = [Rational(str(j)) for j in (
             i.ax, i.bx, i.ay, i.by)]
         return sympify((
-            r1 < re and re < r2) and (
-            i1 < im and im < i2))
+            r1 <= re and re <= r2) and (
+            i1 <= im and im <= i2))
 
 CRootOf = ComplexRootOf
 
@@ -800,7 +967,7 @@ def __new__(cls, expr, func=None, x=None, auto=True, quadratic=False):
         func = Lambda(var, expr)
 
         rational = cls._is_func_rational(poly, func)
-        (_, factors), terms = poly.factor_list(), []
+        factors, terms = _pure_factors(poly), []
 
         for poly, k in factors:
             if poly.is_linear:
@@ -849,7 +1016,7 @@ def _transform(cls, expr, x):
 
     @classmethod
     def _is_func_rational(cls, poly, func):
-        """Check if a lambda is areational function. """
+        """Check if a lambda is a rational function. """
         var, expr = func.variables[0], func.expr
         return expr.is_rational_function(var)
 
@@ -949,39 +1116,3 @@ def _eval_derivative(self, x):
         var, expr = self.fun.args
         func = Lambda(var, expr.diff(x))
         return self.new(self.poly, func, self.auto)
-
-
-def bisect(f, a, b, tol):
-    """
-    Implements bisection. This function is used in RootOf.eval_rational() and
-    it needs to be robust.
-
-    Examples
-    ========
-
-    >>> from sympy import S
-    >>> from sympy.polys.rootoftools import bisect
-    >>> bisect(lambda x: x**2-1, -10, 0, S(1)/10**2)
-    -1025/1024
-    >>> bisect(lambda x: x**2-1, -10, 0, S(1)/10**4)
-    -131075/131072
-
-    """
-    a = sympify(a)
-    b = sympify(b)
-    fa = f(a)
-    fb = f(b)
-    if fa * fb >= 0:
-        raise ValueError("bisect: f(a) and f(b) must have opposite signs")
-    while (b - a > tol):
-        c = (a + b)/2
-        fc = f(c)
-        if (fc == 0):
-            return c  # We need to make sure f(c) is not zero below
-        if (fa * fc < 0):
-            b = c
-            fb = fc
-        else:
-            a = c
-            fa = fc
-    return (a + b)/2
diff --git a/sympy/polys/specialpolys.py b/sympy/polys/specialpolys.py
index 5dca7a4043dc..3977cf97ad19 100644
--- a/sympy/polys/specialpolys.py
+++ b/sympy/polys/specialpolys.py
@@ -62,7 +62,7 @@ def swinnerton_dyer_poly(n, x=None, polys=False):
         for i in range(2, n + 1):
             p = nextprime(p)
             a.append(sqrt(p))
-        return minimal_polynomial(Add(*a), x, polys=args.get('polys', False))
+        return minimal_polynomial(Add(*a), x, polys=polys)
 
     if n == 1:
         ex = x**2 - 2
diff --git a/sympy/polys/tests/test_polyroots.py b/sympy/polys/tests/test_polyroots.py
index 41865ef98de0..db00d3f1ac16 100644
--- a/sympy/polys/tests/test_polyroots.py
+++ b/sympy/polys/tests/test_polyroots.py
@@ -2,7 +2,7 @@
 
 from sympy import (S, symbols, Symbol, Wild, Rational, sqrt,
     powsimp, sin, cos, pi, I, Interval, re, im, exp, ZZ, Piecewise,
-    acos, root)
+    acos, root, conjugate)
 
 from sympy.polys import Poly, cyclotomic_poly, intervals, nroots, rootof
 
@@ -23,6 +23,20 @@
 a, b, c, d, e, q, t, x, y, z = symbols('a,b,c,d,e,q,t,x,y,z')
 
 
+def _check(roots):
+    # this is the desired invariant for roots returned
+    # by all_roots. It is trivially true for linear
+    # polynomials.
+    nreal = sum([1 if i.is_real else 0 for i in roots])
+    assert list(sorted(roots[:nreal])) == list(roots[:nreal])
+    for ix in range(nreal, len(roots), 2):
+        if not (
+                roots[ix + 1] == roots[ix] or
+                roots[ix + 1] == conjugate(roots[ix])):
+            return False
+    return True
+
+
 def test_roots_linear():
     assert roots_linear(Poly(2*x + 1, x)) == [-Rational(1, 2)]
 
@@ -32,6 +46,7 @@ def test_roots_quadratic():
     assert roots_quadratic(Poly(2*x**2 + 3*x, x)) == [-Rational(3, 2), 0]
     assert roots_quadratic(Poly(2*x**2 + 3, x)) == [-I*sqrt(6)/2, I*sqrt(6)/2]
     assert roots_quadratic(Poly(2*x**2 + 4*x + 3, x)) == [-1 - I*sqrt(2)/2, -1 + I*sqrt(2)/2]
+    _check(Poly(2*x**2 + 4*x + 3, x).all_roots())
 
     f = x**2 + (2*a*e + 2*c*e)/(a - c)*x + (d - b + a*e**2 - c*e**2)/(a - c)
     assert roots_quadratic(Poly(f, x)) == \
@@ -72,29 +87,38 @@ def test_issue_8438():
 
 def test_issue_8285():
     roots = (Poly(4*x**8 - 1, x)*Poly(x**2 + 1)).all_roots()
-    assert roots == _nsort(roots)
+    assert _check(roots)
     f = Poly(x**4 + 5*x**2 + 6, x)
     ro = [rootof(f, i) for i in range(4)]
     roots = Poly(x**4 + 5*x**2 + 6, x).all_roots()
     assert roots == ro
-    assert roots == _nsort(roots)
+    assert _check(roots)
     # more than 2 complex roots from which to identify the
     # imaginary ones
     roots = Poly(2*x**8 - 1).all_roots()
-    assert roots == _nsort(roots)
+    assert _check(roots)
     assert len(Poly(2*x**10 - 1).all_roots()) == 10  # doesn't fail
 
 
 def test_issue_8289():
     roots = (Poly(x**2 + 2)*Poly(x**4 + 2)).all_roots()
-    assert roots == _nsort(roots)
+    assert _check(roots)
     roots = Poly(x**6 + 3*x**3 + 2, x).all_roots()
-    assert roots == _nsort(roots)
+    assert _check(roots)
     roots = Poly(x**6 - x + 1).all_roots()
-    assert roots == _nsort(roots)
-    # all imaginary roots
+    assert _check(roots)
+    # all imaginary roots with multiplicity of 2
     roots = Poly(x**4 + 4*x**2 + 4, x).all_roots()
-    assert roots == _nsort(roots)
+    assert _check(roots)
+
+
+def test_issue_14291():
+    assert Poly(((x - 1)**2 + 1)*((x - 1)**2 + 2)*(x - 1)
+        ).all_roots() == [1, 1 - I, 1 + I, 1 - sqrt(2)*I, 1 + sqrt(2)*I]
+    p = x**4 + 10*x**2 + 1
+    ans = [rootof(p, i) for i in range(4)]
+    assert Poly(p).all_roots() == ans
+    _check(ans)
 
 
 def test_issue_13340():
@@ -241,6 +265,7 @@ def test_roots_binomial():
 
     assert roots_binomial(Poly(5*x**4 + 2, x)) == \
         [-A - A*I, -A + A*I, A - A*I, A + A*I]
+    _check(roots_binomial(Poly(x**8 - 2)))
 
     a1 = Symbol('a1', nonnegative=True)
     b1 = Symbol('b1', nonnegative=True)
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 6c614c87f713..e73ca8d4b6a6 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -2798,6 +2798,13 @@ def test_nroots():
     f = Poly(x**4 - 1)
     assert f.nroots(2) == [w.n(2) for w in f.all_roots()]
 
+    assert str(Poly(x**16 + 32*x**14 + 508*x**12 + 5440*x**10 +
+        39510*x**8 + 204320*x**6 + 755548*x**4 + 1434496*x**2 +
+        877969).nroots(2)) == ('[-1.7 - 1.9*I, -1.7 + 1.9*I, -1.7 '
+        '- 2.5*I, -1.7 + 2.5*I, -1.0*I, 1.0*I, -1.7*I, 1.7*I, -2.8*I, '
+        '2.8*I, -3.4*I, 3.4*I, 1.7 - 1.9*I, 1.7 + 1.9*I, 1.7 - 2.5*I, '
+        '1.7 + 2.5*I]')
+
 
 def test_ground_roots():
     f = x**6 - 4*x**4 + 4*x**3 - x**2
diff --git a/sympy/polys/tests/test_rootoftools.py b/sympy/polys/tests/test_rootoftools.py
index 34bf60007de1..262baf2407e9 100644
--- a/sympy/polys/tests/test_rootoftools.py
+++ b/sympy/polys/tests/test_rootoftools.py
@@ -131,8 +131,9 @@ def test_CRootOf___eq__():
 
 def test_CRootOf___eval_Eq__():
     f = Function('f')
-    r = rootof(x**3 + x + 3, 2)
-    r1 = rootof(x**3 + x + 3, 1)
+    eq = x**3 + x + 3
+    r = rootof(eq, 2)
+    r1 = rootof(eq, 1)
     assert Eq(r, r1) is S.false
     assert Eq(r, r) is S.true
     assert Eq(r, x) is S.false
@@ -141,16 +142,17 @@ def test_CRootOf___eval_Eq__():
     assert Eq(r, I) is S.false
     assert Eq(r, f(0)) is S.false
     assert Eq(r, f(0)) is S.false
-    sol = solve(r.expr)
+    sol = solve(eq)
     for s in sol:
         if s.is_real:
             assert Eq(r, s) is S.false
-    r = rootof(r.expr, 0)
+    r = rootof(eq, 0)
     for s in sol:
         if s.is_real:
             assert Eq(r, s) is S.true
-    eq = (x**3 + x + 1)
-    assert [Eq(rootof(eq,i), j) for i in range(3) for j in solve(eq)] == [
+    eq = x**3 + x + 1
+    sol = solve(eq)
+    assert [Eq(rootof(eq, i), j) for i in range(3) for j in sol] == [
         False, False, True, False, True, False, True, False, False]
     assert Eq(rootof(eq, 0), 1 + S.ImaginaryUnit) == False
 
@@ -191,6 +193,10 @@ def test_CRootOf_evalf():
 
     p = legendre_poly(4, x, polys=True)
     roots = [str(r.n(17)) for r in p.real_roots()]
+    # magnitudes are given by
+    # sqrt(3/S(7) - 2*sqrt(6/S(5))/7)
+    #   and
+    # sqrt(3/S(7) + 2*sqrt(6/S(5))/7)
     assert roots == [
             "-0.86113631159405258",
             "-0.33998104358485626",
@@ -244,7 +250,23 @@ def test_CRootOf_evalf():
 
     # watch out for UnboundLocalError
     c = CRootOf(90720*x**6 - 4032*x**4 + 84*x**2 - 1, 0)
-    assert str(c._eval_evalf(2)) == '-0.e-1'
+    assert c._eval_evalf(2)  # doesn't fail
+
+    # watch out for imaginary parts that don't want to evaluate
+    assert str(RootOf(x**16 + 32*x**14 + 508*x**12 + 5440*x**10 +
+        39510*x**8 + 204320*x**6 + 755548*x**4 + 1434496*x**2 +
+        877969, 10).n(2)) == '-3.4*I'
+
+    # check reset
+    r = [RootOf(x**3 + x + 3, i) for i in range(3)]
+    r[0]._reset()
+    for ri in r:
+        i = ri._get_interval()
+        n = ri.n(2)
+        assert str(i) != str(ri._get_interval())
+        ri._reset()
+        # XXX using strings because intervals don't compare equal
+        assert str(i) == str(ri._get_interval())
 
 
 def test_CRootOf_evalf_caching_bug():
@@ -283,12 +305,9 @@ def test_CRootOf_all_roots():
 
 def test_CRootOf_eval_rational():
     p = legendre_poly(4, x, polys=True)
-    roots = [r.eval_rational(S(1)/10**20) for r in p.real_roots()]
+    roots = [r.eval_rational(n=18) for r in p.real_roots()]
     for r in roots:
         assert isinstance(r, Rational)
-    # All we know is that the Rational instance will be at most 1/10^20 from
-    # the exact root. So if we evaluate to 17 digits, it must be exactly equal
-    # to:
     roots = [str(r.n(17)) for r in roots]
     assert roots == [
             "-0.86113631159405258",
@@ -444,3 +463,43 @@ def test_issue_8316():
     assert len(f.all_roots()) == 8
     f = Poly(7*x**8 - 10)
     assert len(f.all_roots()) == 8
+
+
+def test__imag_count():
+    from sympy.polys.rootoftools import _imag_count_of_factor
+    def imag_count(p):
+        return sum([_imag_count_of_factor(f)*m for f, m in
+        p.factor_list()[1]])
+    assert imag_count(Poly(x**6 + 10*x**2 + 1)) == 2
+    assert imag_count(Poly(x**2)) == 0
+    assert imag_count(Poly([1]*3 + [-1], x)) == 0
+    assert imag_count(Poly(x**3 + 1)) == 0
+    assert imag_count(Poly(x**2 + 1)) == 2
+    assert imag_count(Poly(x**2 - 1)) == 0
+    assert imag_count(Poly(x**4 - 1)) == 2
+    assert imag_count(Poly(x**4 + 1)) == 0
+    assert imag_count(Poly([1, 2, 3], x)) == 0
+    assert imag_count(Poly(x**3 + x + 1)) == 0
+    assert imag_count(Poly(x**4 + x + 1)) == 0
+    def q(r1, r2, p):
+        return Poly(((x - r1)*(x - r2)).subs(x, x**p), x)
+    assert imag_count(q(-1, -2, 2)) == 4
+    assert imag_count(q(-1, 2, 2)) == 2
+    assert imag_count(q(1, 2, 2)) == 0
+    assert imag_count(q(1, 2, 4)) == 4
+    assert imag_count(q(-1, 2, 4)) == 2
+    assert imag_count(q(-1, -2, 4)) == 0
+
+
+def test_set_imaginary():
+    r = RootOf(x**4 + 4*x**2 + 1, 1)
+    i = r._get_interval()
+    assert r.is_imaginary and i.ax == i.bx == 0
+
+
+def test_is_disjoint():
+    eq = x**3 + 5*x + 1
+    ir = rootof(eq, 0)._get_interval()
+    ii = rootof(eq, 1)._get_interval()
+    assert ir.is_disjoint(ii)
+    assert ii.is_disjoint(ir)
diff --git a/sympy/polys/tests/test_specialpolys.py b/sympy/polys/tests/test_specialpolys.py
index ab92607b0e93..943b0bca48ad 100644
--- a/sympy/polys/tests/test_specialpolys.py
+++ b/sympy/polys/tests/test_specialpolys.py
@@ -1,6 +1,7 @@
 """Tests for functions for generating interesting polynomials. """
 
-from sympy import Poly, ZZ, symbols
+from sympy import Poly, ZZ, symbols, sqrt, prime, Add
+from sympy.utilities.iterables import permute_signs
 from sympy.utilities.pytest import raises
 
 from sympy.polys.specialpolys import (
@@ -29,6 +30,12 @@ def test_swinnerton_dyer_poly():
     assert swinnerton_dyer_poly(2, x) == x**4 - 10*x**2 + 1
     assert swinnerton_dyer_poly(
         3, x) == x**8 - 40*x**6 + 352*x**4 - 960*x**2 + 576
+    # we only need to check that the polys arg works but
+    # we may as well test that the roots are correct
+    p = [sqrt(prime(i)) for i in range(1, 5)]
+    assert str([i.n(3) for i in
+        swinnerton_dyer_poly(4, polys=True).all_roots()]
+        ) == str(sorted([Add(*i).n(3) for i in permute_signs(p)]))
 
 
 def test_cyclotomic_poly():
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index cc1fc40e86be..8fa009c72286 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -834,10 +834,10 @@ def solve(f, *symbols, **flags):
         >>> solve(eq, check=False)
         [1/3]
 
-    In the above example there is only a single solution to the equation. Other
-    expressions will yield spurious roots which must be checked manually;
-    roots which give a negative argument to odd-powered radicals will also need
-    special checking:
+    In the above example there is only a single solution to the
+    equation. Other expressions will yield spurious roots which
+    must be checked manually; roots which give a negative argument
+    to odd-powered radicals will also need special checking:
 
         >>> from sympy import real_root, S
         >>> eq = root(x, 3) - root(x, 5) + S(1)/7
@@ -848,15 +848,15 @@ def solve(f, *symbols, **flags):
         >>> [abs(eq.subs(x,i).n(2)) for i in sol]
         [0.48, 0.e-110, 0.e-110, 0.052, 0.052]
 
-        The first solution is negative so real_root must be used to see that
-        it satisfies the expression:
+        The first solution is negative so real_root must be used to see
+        that it satisfies the expression:
 
         >>> abs(real_root(eq.subs(x, sol[0])).n(2))
         0.e-110
 
-    If the roots of the equation are not real then more care will be necessary
-    to find the roots, especially for higher order equations. Consider the
-    following expression:
+    If the roots of the equation are not real then more care will be
+    necessary to find the roots, especially for higher order equations.
+    Consider the following expression:
 
         >>> expr = root(x, 3) - root(x, 5)
 
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index afa5291e1fd4..4ae18490c678 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -537,8 +537,12 @@ def _solve_trig1(f, symbol, domain):
         return ConditionSet(symbol, Eq(f, 0), S.Reals)
 
     solns = solveset_complex(g, y) - solveset_complex(h, y)
+    if isinstance(solns, ConditionSet):
+        raise NotImplementedError
 
     if isinstance(solns, FiniteSet):
+        if any(isinstance(s, RootOf) for s in solns):
+            raise NotImplementedError
         result = Union(*[invert_complex(exp(I*symbol), s, symbol)[1]
                        for s in solns])
         return Intersection(result, domain)
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 14d443d54ecb..eb468dda8003 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -781,18 +781,18 @@ def test_solve_trig():
         imageset(Lambda(n, 2*pi*n + pi/2), S.Integers)
 
     assert solveset_real(cos(x), x) == \
-        Union(imageset(Lambda(n, 2*pi*n - pi/2), S.Integers),
-              imageset(Lambda(n, 2*pi*n + pi/2), S.Integers))
+        Union(imageset(Lambda(n, 2*pi*n + pi/2), S.Integers),
+              imageset(Lambda(n, 2*pi*n + 3*pi/2), S.Integers))
 
     assert solveset_real(sin(x) + cos(x), x) == \
-        Union(imageset(Lambda(n, 2*n*pi - pi/4), S.Integers),
-              imageset(Lambda(n, 2*n*pi + 3*pi/4), S.Integers))
+        Union(imageset(Lambda(n, 2*n*pi + 3*pi/4), S.Integers),
+              imageset(Lambda(n, 2*n*pi + 7*pi/4), S.Integers))
 
     assert solveset_real(sin(x)**2 + cos(x)**2, x) == S.EmptySet
 
     assert solveset_complex(cos(x) - S.Half, x) == \
-        Union(imageset(Lambda(n, 2*n*pi + pi/3), S.Integers),
-              imageset(Lambda(n, 2*n*pi - pi/3), S.Integers))
+        Union(imageset(Lambda(n, 2*n*pi + 5*pi/3), S.Integers),
+              imageset(Lambda(n, 2*n*pi + pi/3), S.Integers))
 
     y, a = symbols('y,a')
     assert solveset(sin(y + a) - sin(y), a, domain=S.Reals) == \


Code review comments:
diff_hunk: @@ -244,7 +250,23 @@ def test_CRootOf_evalf():
 
     # watch out for UnboundLocalError
     c = CRootOf(90720*x**6 - 4032*x**4 + 84*x**2 - 1, 0)
-    assert str(c._eval_evalf(2)) == '-0.e-1'
+    assert c._eval_evalf(2)  # doesn't fail
+
+    # watch out for imaginary parts that don't want to evaluate
+    assert str(RootOf(x**16 + 32*x**14 + 508*x**12 + 5440*x**10 +
+        39510*x**8 + 204320*x**6 + 755548*x**4 + 1434496*x**2 +
+        877969, 10).n(2)) == '-3.4*I'
+
+    # check reset
+    r = [RootOf(x**3 + x + 3, i) for i in range(3)]
+    r[0]._reset()
+    for ri in r:
+        i = ri._get_interval()
+        n = ri.n(2)
+        assert str(i) != str(ri._get_interval())
+        ri._reset()
+        # XXX using strings because intervals don't compare equal
path: sympy/polys/tests/test_rootoftools.py
text: @user1:
Could `__eq__` be redefined to satisfy this?

@author:
I can try. I'm not sure why the equality fails as it is.

@author:
Ah. No `__eq__` for RealInterval. That is now fixed. `func` and `args` are now defined for each, too.

diff_hunk: @@ -181,15 +315,22 @@ def args(self):
 
     @property
     def free_symbols(self):
-        # CRootOf currently only works with univariate expressions and although
-        # the poly attribute is often a PurePoly, sometimes it is a Poly. In
-        # either case no free symbols should be reported.
+        # CRootOf currently only works with univariate expressions
+        # whose poly attribute should be a PurePoly with no free
+        # symbols
         return set()
 
     def _eval_is_real(self):
         """Return ``True`` if the root is real. """
         return self.index < len(_reals_cache[self.poly])
 
+    def _eval_is_imaginary(self):
+        """Return ``True`` if the root is real. """
path: sympy/polys/rootoftools.py
text: @user1:
real -> imaginary

diff_hunk: @@ -29,22 +29,84 @@
 from sympy.polys.domains import QQ
 
 from mpmath import mpf, mpc, findroot, workprec
-from mpmath.libmp.libmpf import prec_to_dps
+from mpmath.libmp.libmpf import dps_to_prec, prec_to_dps
 
 from sympy.utilities import lambdify, public, sift
 
-from sympy.core.compatibility import range
+from sympy.core.compatibility import range, ordered
 
 from math import log as mathlog
 
 __all__ = ['CRootOf']
 
-def _ispow2(i):
-    v = mathlog(i, 2)
-    return v == int(v)
 
-_reals_cache = {}
-_complexes_cache = {}
+
+class _pure_key_dict(dict):
+    """A minimal dictionary that makes sure that the key is PurePoly.
+
+    Examples
+    ========
+
+    Only the following four actions are guaranteed:
+
+    >>> from sympy.polys.rootoftools import _pure_key_dict
+    >>> from sympy import S, PurePoly
+    >>> from sympy.abc import x, y
+
+    1) creation
+
+    >>> P = _pure_key_dict()
+
+    2) assignment for a PurePoly or univariate polynomial
+
+    >>> P[x] = 1
+    >>> P[PurePoly(x - y, x)] = 2
+
+    3) retrieval based on PurePoly key comparison (use this
+       instead of the get method)
+
+    >>> P[y]
+    1
+
+    4) KeyError when trying to retrieve a nonexisting key
+
+    >>> P[y + 1]
+    Traceback (most recent call last):
+    ...
+    KeyError: PurePoly(y + 1, y, domain='ZZ')
+    """
+    def __getitem__(self, k):
+        if not isinstance(k, PurePoly):
+            assert len(k.free_symbols) == 1
+            k = PurePoly(k, expand=False)
+        return self.__dict__[k]
path: sympy/polys/rootoftools.py
text: @user1:
It seems that `__dict__` is not the place where a dict object stores its data.
```
>>> from sympy.polys.rootoftools import _pure_key_dict
>>> from sympy.abc import x
>>> P = _pure_key_dict()
>>> P[x] = 1
>>> x in P
False
>>> P.items()
[]
```
If we want to use other dict methods than `__setitem__` and `__getitem__`, we should probably either define them explicitly or else try to use here something like `dict.__getitem__(self, k)`.

@author:
I am a lightweight when it comes to creating new classes...and this is a lightweight implementation adequate for the needs that I had. I also learned that 'in' doesn't work so I documented what is supported...and if it's not there, laissez-faire :-) Do you think it needs more? See the docstring for examples to see if I have left out anything important or mislead in any way. Maybe a note that `in` is not supported would be good to add?

I will address the other docstring edit at the same time this evening.

@user1:
It is probably not necessary to make lightweight cache a subclass of dict, if its methods are not needed. It seems to me that there could be `self.cache = {}` in `__init__(self)` and `__setitem__(self, k, v)` could then use `self.cache[k] = v`.

@author:
The main need that I have is for there to be a guard that checks that setting and getting is done via PurePoly. Tucking the dict into the `__init__` won't help with that. (At least I don't see how that will solve that problem.)

@author:
LOL. I meant "caveat emptor" not "laissez-faire"!

@user1:
> The main need that I have is for there to be a guard that checks that setting and getting is done via PurePoly.

I think that this could be handled by `__getitem__` and `__setitem__` exactly as done by the code now. There should be no need to change that part. But I would not use `__dict__` for caching the polynomials. `__dict__` is used internally for storing attributes. Python will be confused if its keys are not strings.
```
>>> from sympy.polys.rootoftools import _pure_key_dict
>>> from sympy.abc import x
>>> P = _pure_key_dict()
>>> dir(P)
['__class__', '__cmp__', '__contains__', '__delattr__', '__delitem__', '__dict__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__iter__', '__le__', '__len__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setitem__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'clear', 'copy', 'fromkeys', 'get', 'has_key', 'items', 'iteritems', 'iterkeys', 'itervalues', 'keys', 'pop', 'popitem', 'setdefault', 'update', 'values', 'viewitems', 'viewkeys', 'viewvalues']
>>> P[x] = 1
>>> dir(P)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "sympy/core/expr.py", line 331, in __lt__
    raise TypeError("Invalid comparison %s < %s" % (self, other))
TypeError: Invalid comparison PurePoly(x, x, domain='ZZ') < __setattr__
```
I think that there are two possibilities. Either define `_pure_key_dict` as a subclass of dict (as it is now) but using its (built-in) dictionary instead of `__dict__`, or define it as a class having a dict attribute (that I called `cache` above). In the former case, all methods of dict would be available (with redefined `__getattr__` and `__setattr__`); the latter case would be a really lightweight solution (with only two methods `__getattr__` and `__setattr__` defined).

@author:
OK, now I think I see what you mean. I was thinking that you meant for the data to be stored in RootOf someway...and I was having a hard time imagining how to do that. I have updated using your suggestion about a self.dict attribute and make the class derive from object, not dict.

diff_hunk: @@ -1835,14 +1982,24 @@ def conjugate(self):
 
     def is_disjoint(self, other):
         """Return ``True`` if two isolation intervals are disjoint. """
-        if self.conj != other.conj:
+        if isinstance(other, RealInterval):
+            return other.is_disjoint(self)
+        if self.conj != other.conj:  # above and below real axis
             return True
         re_distinct = (self.bx <= other.ax or other.bx <= self.ax)
         if re_distinct:
             return True
         im_distinct = (self.by <= other.ay or other.by <= self.ay)
         return im_distinct
 
+    def _set_imaginary(self):
+        a1, a2 = self.a
+        b1, b2 = self.b
+        a = a1 - a1, a2
path: sympy/polys/rootisolation.py
text: @user1:
Will this be `0, a2`?

@author:
It will be a `PythonRational(0)`, I believe. I was just avoiding creating it explicitly since PythonRational already knows how to do subtraction.

@user1:
I thought so. It could probably also be GMPYRational. How about adding a comment?

@author:
Changed as follows:

```diff
diff --git a/sympy/polys/rootisolation.py b/sympy/polys/rootisolation.py
index 4d2bde9..8ee07a4 100644
--- a/sympy/polys/rootisolation.py
+++ b/sympy/polys/rootisolation.py
@@ -2013,8 +2013,9 @@ def is_disjoint(self, other):
     def _set_imaginary(self):
         a1, a2 = self.a
         b1, b2 = self.b
-        a = a1 - a1, a2
-        b = b1 - b1, b2
+        zero = a1 - a1  # to preserve type
+        a = zero, a2
+        b = zero, b2
         return ComplexInterval(a, b, self.I, self.Q,
```