Repository: sympy/sympy

Issue:
create Root object
Although some of this is discussed in [here](https://groups.google.com/forum/#!searchin/sympy/Fateman/sympy/40d-NGT_gF4/SJCETmfTfbkJ) I did not want this issue to get buried in the discussion.  When solving a radical expression like the following, SymPy's use of the principle root make solve miss a solution:

```
>>> eq
-x + (x**3 - 3*x**2)**(1/3) + 4/3
>>> solve(_, check=0)
[-8*sqrt(6)/9 + 8/3, 8*sqrt(6)/9 + 8/3]  <--------- both are solutions (as can be seen graphically)
>>> [eq.subs(x,i).n() for i in _]
[1.2659863237109 + 0.730917544784875*I, 0.e-124]  <-- the first one does not appear to be a soln

Here's a case where the root is rational:

>>> eq=eq-S(1)/3; eq
-x + (x**3 - 3*x**2)**(1/3) + 4/3
>>> solve(eq)
[]
>>> solve(eq, check=False)
[1/3]
>>> eq.subs(x,S(1)/3)
2/3 + 2*(-1)**(1/3)/3
```

The 1/3 does not appear to be a solution because the principle value of the root is being used rather than the real value:

```
>>> root(-1,3).n()
0.5 + 0.866025403784439*I
>>> real_root(-1, 3)
-1
```

Note that nsolve cannot find the root either because of this issue:

```
>>> nsolve(eq,.33)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "sympy\solvers\solvers.py", line 2532, in nsolve
    return findroot(f, x0, **kwargs)
  File "c:\Python26\lib\site-packages\mpmath\calculus\optimization.py", line 975
, in findroot
    % (norm(f(*xl))**2, tol))
ValueError: Could not find root within given tolerance. (2.60832 > 2.1684e-19)
Try another starting point or tweak arguments.
>>> nsolve(eq,.33+.01j)  # try an initial guess that is off the real axis
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "sympy\solvers\solvers.py", line 2532, in nsolve
    return findroot(f, x0, **kwargs)
  File "c:\Python26\lib\site-packages\mpmath\calculus\optimization.py", line 975
, in findroot
    % (norm(f(*xl))**2, tol))
ValueError: Could not find root within given tolerance. (2.60832 > 2.1684e-19)
Try another starting point or tweak arguments.
```

How should a formula be written that will allow one to verify that 1/3 is a solution?

Discussion:
asmeurer: Why not use `a**(k/n)`?
smichr: I guess the former can be used for 0<k<n; what I have is correct for arbitrary integer values (so one may index the roots with negative values if desired. If the user specifies an integer for k, I could use your form, otherwise the form that I have. How does sound?
smichr: > Why not use a**(k/n)

I think you mean (for a < 0): `a**(1/n)*(-1)**(k/n)`...or maybe I am misunderstanding something:

```
>>> from sympy.polys.polyutils import _nsort
>>> _nsort([(2*(-1)**(i/S(5))).n(2) for i in range(1,6)])
[-2.0, -1.6 + 1.2*I, -0.62 + 1.9*I, 0.62 + 1.9*I, 1.6 + 1.2*I]
>>> _nsort([(2**(S(1)/5)*(-1)**(i/S(5))).n(2) for i in range(1,6)])
[-1.1, -0.93 + 0.67*I, -0.35 + 1.1*I, 0.35 + 1.1*I, 0.93 + 0.67*I]
>>> _nsort([(2**(S(1)/5)*(-1)**(i/S(5))).n(2) for i in range(1,6)])
[-1.1, -0.93 + 0.67*I, -0.35 + 1.1*I, 0.35 + 1.1*I, 0.93 + 0.67*I]
```
asmeurer: Yes I think I mean that.
smichr: Ahh...kth nth root of x is `root(x, n)*(-1)**(2*k/n)` (note the factor of 2). This has been rewritten in those terms. You can see the doctests for examples.

Patch:
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 1b1c03f781..54bc5031d3 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -157,7 +157,11 @@ def _eval_evalf(self, prec):
         """ Careful! any evalf of polar numbers is flaky """
         from sympy import im, pi, re
         i = im(self.args[0])
-        if i <= -pi or i > pi:
+        try:
+            bad = (i <= -pi or i > pi)
+        except TypeError:
+            bad = True
+        if bad:
             return self  # cannot evalf for this argument
         res = exp(self.args[0])._eval_evalf(prec)
         if i > 0 and im(res) < 0:
diff --git a/sympy/functions/elementary/miscellaneous.py b/sympy/functions/elementary/miscellaneous.py
index 5684ea665f..cbdbe08603 100644
--- a/sympy/functions/elementary/miscellaneous.py
+++ b/sympy/functions/elementary/miscellaneous.py
@@ -162,10 +162,9 @@ def cbrt(arg):
     return C.Pow(arg, C.Rational(1, 3))
 
 
-def root(arg, n):
-    """The n-th root function (a shortcut for ``arg**(1/n)``)
-
-    root(x, n) -> Returns the principal n-th root of x.
+def root(arg, n, k=0):
+    """root(x, n, k) -> Returns the k-th n-th root of x, defaulting to the
+    principle root (k=0).
 
 
     Examples
@@ -186,6 +185,10 @@ def root(arg, n):
     >>> root(x, -Rational(2, 3))
     x**(-3/2)
 
+    To get the k-th n-th root, specify k:
+
+    >>> root(-2, 3, 2)
+    -(-1)**(2/3)*2**(1/3)
 
     To get all n n-th roots you can use the RootOf function.
     The following examples show the roots of unity for n
@@ -193,13 +196,13 @@ def root(arg, n):
 
     >>> from sympy import RootOf, I
 
-    >>> [ RootOf(x**2-1,i) for i in (0,1) ]
+    >>> [ RootOf(x**2 - 1, i) for i in range(2) ]
     [-1, 1]
 
-    >>> [ RootOf(x**3-1,i) for i in (0,1,2) ]
+    >>> [ RootOf(x**3 - 1,i) for i in range(3) ]
     [1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]
 
-    >>> [ RootOf(x**4-1,i) for i in (0,1,2,3) ]
+    >>> [ RootOf(x**4 - 1,i) for i in range(4) ]
     [-1, 1, -I, I]
 
     SymPy, like other symbolic algebra systems, returns the
@@ -211,8 +214,8 @@ def root(arg, n):
     >>> root(-8, 3)
     2*(-1)**(1/3)
 
-    The real_root function can be used to either make such a result
-    real or simply return the real root in the first place:
+    The real_root function can be used to either make the principle
+    result real (or simply to return the real root directly):
 
     >>> from sympy import real_root
     >>> real_root(_)
@@ -220,6 +223,12 @@ def root(arg, n):
     >>> real_root(-32, 5)
     -2
 
+    Alternatively, the n//2-th n-th root of a negative number can be
+    computed with root:
+
+    >>> root(-32, 5, 5//2)
+    -2
+
     See Also
     ========
 
@@ -238,12 +247,16 @@ def root(arg, n):
 
     """
     n = sympify(n)
+    if k:
+        return C.Pow(arg, S.One/n)*S.NegativeOne**(2*k/n)
     return C.Pow(arg, 1/n)
 
 
 def real_root(arg, n=None):
     """Return the real nth-root of arg if possible. If n is omitted then
-    all instances of (-n)**(1/odd) will be changed to -n**(1/odd).
+    all instances of (-n)**(1/odd) will be changed to -n**(1/odd); this
+    will only create a real root of a principle root -- the presence of
+    other factors may cause the result to not be real.
 
     Examples
     ========
@@ -258,6 +271,15 @@ def real_root(arg, n=None):
     >>> real_root(_)
     -2
 
+    If one creates a non-principle root and applies real_root, the
+    result will not be real (so use with caution):
+
+    >>> root(-8, 3, 2)
+    -2*(-1)**(2/3)
+    >>> real_root(_)
+    -2*(-1)**(2/3)
+
+
     See Also
     ========
 
@@ -266,10 +288,20 @@ def real_root(arg, n=None):
     root, sqrt
     """
     if n is not None:
-        n = as_int(n)
-        rv = C.Pow(arg, Rational(1, n))
-        if n % 2 == 0:
-            return rv
+        try:
+            n = as_int(n)
+            arg = sympify(arg)
+            if arg.is_positive or arg.is_negative:
+                rv = root(arg, n)
+            else:
+                raise ValueError
+        except ValueError:
+            return root(arg, n)*C.Piecewise(
+                (S.One, ~C.Equality(C.im(arg), 0)),
+                (C.Pow(S.NegativeOne, S.One/n)**(2*C.floor(n/2)), C.And(
+                    C.Equality(n % 2, 1),
+                    arg < 0)),
+                (S.One, True))
     else:
         rv = sympify(arg)
     n1pow = Transform(lambda x: -(-x.base)**x.exp,
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 306ffa845e..f3c6dc5752 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -741,14 +741,12 @@ def solve(f, *symbols, **flags):
         >>> expr = root(x, 3) - root(x, 5)
 
     We will construct a known value for this expression at x = 3 by selecting
-    a non-principle root for each (i.e. multiplying by a power of the same
-    root of -1):
+    the 1-th root for each radical:
 
-        >>> expr1 = root(x, 3)*root(-1, 3)**2 - root(x, 5)*root(-1, 5)**2
+        >>> expr1 = root(x, 3, 1) - root(x, 5, 1)
         >>> v = expr1.subs(x, -3)
 
-    The solve function is unable to find any exact roots to this equation in
-    either form:
+    The solve function is unable to find any exact roots to this equation:
 
         >>> eq = Eq(expr, v); eq1 = Eq(expr1, v)
         >>> solve(eq, check=False), solve(eq1, check=False)


Test patch:
diff --git a/sympy/functions/elementary/tests/test_miscellaneous.py b/sympy/functions/elementary/tests/test_miscellaneous.py
index 883ad265f8..a2b68fb296 100644
--- a/sympy/functions/elementary/tests/test_miscellaneous.py
+++ b/sympy/functions/elementary/tests/test_miscellaneous.py
@@ -185,6 +185,7 @@ def test_issue_8413():
 def test_root():
     from sympy.abc import x
     n = Symbol('n', integer=True)
+    k = Symbol('k', integer=True)
 
     assert root(2, 2) == sqrt(2)
     assert root(2, 1) == 2
@@ -206,6 +207,8 @@ def test_root():
     assert root(x, n) == x**(1/n)
     assert root(x, -n) == x**(-1/n)
 
+    assert root(x, n, k) == x**(1/n)*(-1)**(2*k/n)
+
 
 def test_real_root():
     assert real_root(-8, 3) == -2
@@ -218,6 +221,16 @@ def test_real_root():
     assert real_root(r1 + r2 + r3) == -1 + r2 + r3
     assert real_root(root(-2, 3)) == -root(2, 3)
     assert real_root(-8., 3) == -2
+    x = Symbol('x')
+    n = Symbol('n')
+    g = real_root(x, n)
+    assert g.subs(dict(x=-8, n=3)) == -2
+    assert g.subs(dict(x=8, n=3)) == 2
+    # give principle root if there is no real root -- if this is not desired
+    # then maybe a Root class is needed to raise an error instead
+    assert g.subs(dict(x=I, n=3)) == cbrt(I)
+    assert g.subs(dict(x=-8, n=2)) == sqrt(-8)
+    assert g.subs(dict(x=I, n=2)) == sqrt(I)
 
 
 def test_rewrite_MaxMin_as_Heaviside():