Repository: sympy/sympy

Issue:
`function_range` and `continuous_domain` showing unexpected results
Consider this
```
>>> from sympy import *
>>> from sympy.calculus.util import continuous_domain, function_range

>>> x = Symbol("x", Real = True)
>>> f = abs(x)
>>> function_range(f, x, S.Reals)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-106-a908be9bcc06> in <module>()
----> 1 function_range(f, x, S.Reals)

~/sympy/sympy/calculus/util.py in function_range(f, symbol, domain)
    212                     vals += FiniteSet(f.subs(symbol, limit_point))
    213 
--> 214             solution = solveset(f.diff(symbol), symbol, interval)
    215 
    216             if not iterable(solution):

~/sympy/sympy/solvers/solveset.py in solveset(f, symbol, domain)
   1694             try:
   1695                 r = Dummy('r', **assumptions)
-> 1696                 return solveset(f.xreplace({symbol: r}), r, domain
   1697                     ).xreplace({r: symbol})
   1698             except InconsistentAssumptions:

~/sympy/sympy/solvers/solveset.py in solveset(f, symbol, domain)
   1712     f = piecewise_fold(f)
   1713 
-> 1714     return _solveset(f, symbol, domain, _check=True)
   1715 
   1716 

~/sympy/sympy/solvers/solveset.py in _solveset(f, symbol, domain, _check)
    912                 in_set = in_set.as_set()
    913             solns = solver(expr, symbol, in_set)
--> 914             result += solns
    915     elif isinstance(f, Eq):
    916         result = solver(Add(f.lhs, - f.rhs, evaluate=False), symbol, domain)

~/sympy/sympy/sets/sets.py in __add__(self, other)
    528 
    529     def __add__(self, other):
--> 530         return self.union(other)
    531 
    532     def __or__(self, other):

~/sympy/sympy/sets/sets.py in union(self, other)
     88 
     89         """
---> 90         return Union(self, other)
     91 
     92     def intersect(self, other):

~/sympy/sympy/sets/sets.py in __new__(cls, *args, **kwargs)
   1033         # Reduce sets using known rules
   1034         if evaluate:
-> 1035             return simplify_union(args)
   1036 
   1037         args = list(ordered(args, Set._infimum_key))

~/sympy/sympy/sets/sets.py in simplify_union(args)
   1872             new_args = False
   1873             for t in args - set((s,)):
-> 1874                 new_set = union_sets(s, t)
   1875                 # This returns None if s does not know how to intersect
   1876                 # with t. Returns the newly intersected set otherwise

~/sympy/sympy/multipledispatch/dispatcher.py in __call__(self, *args, **kwargs)
    196             self._cache[types] = func
    197         try:
--> 198             return func(*args, **kwargs)
    199 
    200         except MDNotImplementedError:

~/sympy/sympy/sets/handlers/union.py in union_sets(a, b)
    109 def union_sets(a, b):
    110     # If `b` set contains one of my elements, remove it from `a`
--> 111     if any(b.contains(x) == True for x in a):
    112         return set((
    113             FiniteSet(*[x for x in a if not b.contains(x)]), b))

~/sympy/sympy/sets/handlers/union.py in <genexpr>(.0)
    109 def union_sets(a, b):
    110     # If `b` set contains one of my elements, remove it from `a`
--> 111     if any(b.contains(x) == True for x in a):
    112         return set((
    113             FiniteSet(*[x for x in a if not b.contains(x)]), b))

~/sympy/sympy/sets/conditionset.py in contains(self, other)
    178     def contains(self, other):
    179         return And(Lambda(self.sym, self.condition)(
--> 180             other), self.base_set.contains(other))
    181 
    182     def _eval_subs(self, old, new):

~/sympy/sympy/core/function.py in __call__(self, *args)
   1897                 'plural': 's'*(list(self.nargs)[0] != 1),
   1898                 'given': n})
-> 1899         return self.expr.xreplace(dict(list(zip(self.variables, args))))
   1900 
   1901     def __eq__(self, other):

~/sympy/sympy/core/basic.py in xreplace(self, rule)
   1180 
   1181         """
-> 1182         value, _ = self._xreplace(rule)
   1183         return value
   1184 

~/sympy/sympy/core/basic.py in _xreplace(self, rule)
   1195                 _xreplace = getattr(a, '_xreplace', None)
   1196                 if _xreplace is not None:
-> 1197                     a_xr = _xreplace(rule)
   1198                     args.append(a_xr[0])
   1199                     changed |= a_xr[1]

~/sympy/sympy/core/basic.py in _xreplace(self, rule)
   1202             args = tuple(args)
   1203             if changed:
-> 1204                 return self.func(*args), True
   1205         return self, False
   1206 

~/sympy/sympy/core/function.py in __new__(cls, expr, *variables, **kwargs)
   1278                 raise ValueError(filldedent('''
   1279                     Can't calculate derivative wrt %s.%s''' % (v,
-> 1280                     __)))
   1281 
   1282         # We make a special case for 0th derivative, because there is no

ValueError: 
Can't calculate derivative wrt 0.
```
The ouput should be [0, oo).
Most likely due to `solveset`. Also, expressions like `log(x)*sin(x)`, `sin(x)*cot(x)`, `sin(abs(x))` etc. create such problems. 

__Also See__

#13273 as well
Here, if we pass expression like `sin(x)*exp(x)` to function_range, the problem is coming in solveset(f.diff(x), x, S.Reals) which gives output as `ConditionSet(x, Eq(exp(x)*sin(x) + exp(x)*cos(x), 0), Reals)`. Here, `solveset` is not able to compute `exp(x)*(sin(x) + cos(x)) = 0`. I don't think such simple equations should remain as such.
LambertW has no series expansion at x=0 (nan)
With sympy version 0.7.5:

```
>>> series(LambertW(x),x)
nan
```

This is probably due to the chosen representation of the derivative

```
>>> LambertW(x).diff(x)
LambertW(x)/(x*(LambertW(x) + 1))
```

which requires a 0/0 limiting form with numerator `LambertW(x)`, which requires the derivative at 0 which has a limiting form of 0/0 and so on...

One mathematical solution is to replace `x = LambertW(x)*exp(LambertW(x))` in the derivative `1/(LambertW(x) + 1)/exp(LambertW(x))` which may be safely evaluated at `x=0` where `W(0) = 0`.  Not sure if this type of thing is easy to do in the code though.


Discussion:

Patch:
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 45ee796d63..47a2f73cc2 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -945,6 +945,12 @@ def eval(cls, x, k=None):
                 return S.NegativeOne
             if x == -log(2)/2:
                 return -log(2)
+            if x == 2*log(2):
+                return log(2)
+            if x == -S.Pi/2:
+                return S.ImaginaryUnit*S.Pi/2
+            if x == exp(1 + S.Exp1):
+                return S.Exp1
             if x is S.Infinity:
                 return S.Infinity
 
@@ -1005,3 +1011,21 @@ def _eval_is_algebraic(self):
                 return False
         else:
             return s.is_algebraic
+
+    def _eval_nseries(self, x, n, logx):
+        if len(self.args) == 1:
+            from sympy import Order, ceiling, expand_multinomial
+            arg = self.args[0].nseries(x, n=n, logx=logx)
+            lt = arg.compute_leading_term(x, logx=logx)
+            lte = 1
+            if lt.is_Pow:
+                lte = lt.exp
+            if ceiling(n/lte) >= 1:
+                s = Add(*[(-S.One)**(k - 1)*Integer(k)**(k - 2)/
+                          factorial(k - 1)*arg**k for k in range(1, ceiling(n/lte))])
+                s = expand_multinomial(s)
+            else:
+                s = S.Zero
+
+            return s + Order(x**n, x)
+        return super(LambertW, self)._eval_nseries(x, n, logx)
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index d9bb918e32..40b8873cdf 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -7,7 +7,8 @@
 from sympy.core.power import Pow
 from sympy.core.logic import fuzzy_and, fuzzy_not
 from sympy.functions.special.zeta_functions import zeta
-from sympy.functions.special.error_functions import erf, erfc
+from sympy.functions.special.error_functions import erf, erfc, Ei
+from sympy.functions.elementary.complexes import re
 from sympy.functions.elementary.exponential import exp, log
 from sympy.functions.elementary.integers import ceiling, floor
 from sympy.functions.elementary.miscellaneous import sqrt
@@ -433,16 +434,16 @@ def eval(cls, a, z):
             elif z is S.Infinity:
                 return S.Zero
             elif z is S.Zero:
-                # TODO: Holds only for Re(a) > 0:
-                return gamma(a)
+                if re(a).is_positive:
+                    return gamma(a)
 
         # We extract branching information here. C/f lowergamma.
         nx, n = z.extract_branch_factor()
-        if a.is_integer and (a > 0) == True:
+        if a.is_integer and a.is_positive:
             nx = unpolarify(z)
             if z != nx:
                 return uppergamma(a, nx)
-        elif a.is_integer and (a <= 0) == True:
+        elif a.is_integer and a.is_nonpositive:
             if n != 0:
                 return -2*pi*I*n*(-1)**(-a)/factorial(-a) + uppergamma(a, nx)
         elif n != 0:
@@ -450,7 +451,9 @@ def eval(cls, a, z):
 
         # Special values.
         if a.is_Number:
-            if a is S.One:
+            if a is S.Zero and z.is_positive:
+                return -Ei(-z)
+            elif a is S.One:
                 return exp(-z)
             elif a is S.Half:
                 return sqrt(pi)*erfc(sqrt(z))


Test patch:
diff --git a/sympy/calculus/tests/test_util.py b/sympy/calculus/tests/test_util.py
index bb65362653..d503e9fc52 100644
--- a/sympy/calculus/tests/test_util.py
+++ b/sympy/calculus/tests/test_util.py
@@ -1,13 +1,13 @@
 from sympy import (Symbol, S, exp, log, sqrt, oo, E, zoo, pi, tan, sin, cos,
-                   cot, sec, csc, Abs, symbols, I, re, Lambda, simplify,
-                   ImageSet, expint)
+                   cot, sec, csc, Abs, symbols, I, re, simplify,
+                   expint)
 from sympy.calculus.util import (function_range, continuous_domain, not_empty_in,
                                  periodicity, lcim, AccumBounds, is_convex,
                                  stationary_points, minimum, maximum)
 from sympy.core import Add, Mul, Pow
 from sympy.sets.sets import (Interval, FiniteSet, EmptySet, Complement,
                             Union)
-from sympy.utilities.pytest import raises
+from sympy.utilities.pytest import raises, XFAIL
 from sympy.abc import x
 
 a = Symbol('a', real=True)
@@ -525,3 +525,9 @@ def test_union_AccumBounds():
     assert AccumBounds(0, 3).union(AccumBounds(-1, 2)) == AccumBounds(-1, 3)
     assert AccumBounds(0, 3).union(AccumBounds(-1, 4)) == AccumBounds(-1, 4)
     raises(TypeError, lambda: AccumBounds(0, 3).union(1))
+
+
+def test_issue_16469():
+    x = Symbol("x", real=True)
+    f = abs(x)
+    assert function_range(f, x, S.Reals) == Interval(0, oo, False, True)
diff --git a/sympy/concrete/tests/test_sums_products.py b/sympy/concrete/tests/test_sums_products.py
index eee2beadb5..cd4b1c30f2 100644
--- a/sympy/concrete/tests/test_sums_products.py
+++ b/sympy/concrete/tests/test_sums_products.py
@@ -3,7 +3,8 @@
     factorial, Function, harmonic, I, Integral, KroneckerDelta, log,
     nan, oo, pi, Piecewise, Product, product, Rational, S, simplify, Identity,
     sin, sqrt, Sum, summation, Symbol, symbols, sympify, zeta, gamma, Le,
-    Indexed, Idx, IndexedBase, prod, Dummy, lowergamma, Range, MatrixSymbol)
+    Indexed, Idx, IndexedBase, prod, Dummy, lowergamma, Range, floor,
+    RisingFactorial, MatrixSymbol)
 from sympy.abc import a, b, c, d, k, m, x, y, z
 from sympy.concrete.summations import telescopic, _dummy_with_inherited_properties_concrete
 from sympy.concrete.expr_with_intlimits import ReorderError
@@ -1164,6 +1165,16 @@ def test_exceptions():
     raises(NotImplementedError, lambda: S.is_convergent())
 
 
+@XFAIL
+def test_issue_14313():
+    assert Sum((S(1)/2)**floor(n/2), (n, 1, oo)).is_convergent()
+
+
+@XFAIL
+def test_issue_14871():
+    assert Sum((S(1)/10)**x*RisingFactorial(0, x)/factorial(x), (x, 0, oo)).rewrite(factorial).doit() == 1
+
+
 def test_issue_17165():
     n = symbols("n", integer=True)
     x = symbols('x')
diff --git a/sympy/functions/elementary/tests/test_exponential.py b/sympy/functions/elementary/tests/test_exponential.py
index ebc7d2727f..4d09d3ae45 100644
--- a/sympy/functions/elementary/tests/test_exponential.py
+++ b/sympy/functions/elementary/tests/test_exponential.py
@@ -481,6 +481,9 @@ def test_lambertw():
     assert LambertW(-pi/2, -1) == -I*pi/2
     assert LambertW(-1/E, -1) == -1
     assert LambertW(-2*exp(-2), -1) == -2
+    assert LambertW(2*log(2)) == log(2)
+    assert LambertW(-pi/2) == I*pi/2
+    assert LambertW(exp(1 + E)) == E
 
     assert LambertW(x**2).diff(x) == 2*LambertW(x**2)/x/(1 + LambertW(x**2))
     assert LambertW(x, k).diff(x) == LambertW(x, k)/x/(1 + LambertW(x, k))
diff --git a/sympy/functions/special/tests/test_gamma_functions.py b/sympy/functions/special/tests/test_gamma_functions.py
index 7c1c2f6774..37b02996d2 100644
--- a/sympy/functions/special/tests/test_gamma_functions.py
+++ b/sympy/functions/special/tests/test_gamma_functions.py
@@ -2,7 +2,7 @@
     Symbol, Dummy, gamma, I, oo, nan, zoo, factorial, sqrt, Rational,
     multigamma, log, polygamma, EulerGamma, pi, uppergamma, S, expand_func,
     loggamma, sin, cos, O, lowergamma, exp, erf, erfc, exp_polar, harmonic,
-    zeta, conjugate)
+    zeta, conjugate, Ei)
 
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
@@ -163,6 +163,9 @@ def test_uppergamma():
         uppergamma(x, y)*log(y) + meijerg([], [1, 1], [0, 0, x], [], y)
     assert td(uppergamma(x, randcplx()), x)
 
+    p = Symbol('p', positive=True)
+    assert uppergamma(0, p) == -Ei(-p)
+    assert uppergamma(p, 0) == gamma(p)
     assert uppergamma(S.Half, x) == sqrt(pi)*erfc(sqrt(x))
     assert not uppergamma(S.Half - 3, x).has(uppergamma)
     assert not uppergamma(S.Half + 3, x).has(uppergamma)
@@ -173,6 +176,7 @@ def test_uppergamma():
               uppergamma(S.Half - 3, x), x)
 
     assert unchanged(uppergamma, x, -oo)
+    assert unchanged(uppergamma, x, 0)
 
     assert tn_branch(-3, uppergamma)
     assert tn_branch(-4, uppergamma)
@@ -188,7 +192,6 @@ def test_uppergamma():
     assert uppergamma(-2, x) == expint(3, x)/x**2
 
     assert conjugate(uppergamma(x, y)) == uppergamma(conjugate(x), conjugate(y))
-    assert conjugate(uppergamma(x, 0)) == gamma(conjugate(x))
     assert unchanged(conjugate, uppergamma(x, -oo))
 
     assert uppergamma(x, y).rewrite(expint) == y**x*expint(-x + 1, y)
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 91eea42005..1869a7247f 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -2,10 +2,10 @@
 
 from sympy import (
     limit, exp, oo, log, sqrt, Limit, sin, floor, cos, ceiling,
-    atan, gamma, Symbol, S, pi, Integral, Rational, I, EulerGamma,
+    atan, gamma, Symbol, S, pi, Integral, Rational, I,
     tan, cot, integrate, Sum, sign, Function, subfactorial, symbols,
     binomial, simplify, frac, Float, sec, zoo, fresnelc, fresnels,
-    acos, erfi)
+    acos, erfi, LambertW)
 
 from sympy.calculus.util import AccumBounds
 from sympy.core.add import Add
@@ -13,8 +13,6 @@
 from sympy.series.limits import heuristics
 from sympy.series.order import Order
 from sympy.utilities.pytest import XFAIL, raises
-from sympy.core.numbers import GoldenRatio
-from sympy.functions.combinatorial.numbers import fibonacci
 
 from sympy.abc import x, y, z, k
 n = Symbol('n', integer=True, positive=True)
@@ -558,3 +556,22 @@ def test_issue_17325():
     assert Limit(x**2, x, 0, dir="+-").doit() == 0
     assert Limit(1/x**2, x, 0, dir="+-").doit() == oo
     raises(ValueError, lambda: Limit(1/x, x, 0, dir="+-").doit())
+
+
+def test_issue_10978():
+    assert LambertW(x).limit(x, 0) == 0
+
+
+@XFAIL
+def test_issue_14313_comment():
+    assert limit(floor(n/2), n, oo) == oo
+
+
+@XFAIL
+def test_issue_15323():
+    d = ((1 - 1/x)**x).diff(x)
+    assert limit(d, x, 1, dir='+') == 1
+
+
+def test_issue_12571():
+    assert limit(-LambertW(-log(x))/log(x), x, 1) == 1
diff --git a/sympy/series/tests/test_series.py b/sympy/series/tests/test_series.py
index 8596d4d283..606b521c1c 100644
--- a/sympy/series/tests/test_series.py
+++ b/sympy/series/tests/test_series.py
@@ -1,5 +1,5 @@
 from sympy import sin, cos, exp, E, series, oo, S, Derivative, O, Integral, \
-    Function, log, sqrt, Symbol, Subs, pi, symbols, IndexedBase, atan
+    Function, log, sqrt, Symbol, Subs, pi, symbols, IndexedBase, atan, LambertW
 from sympy.abc import x, y, n, k
 from sympy.utilities.pytest import raises
 from sympy.core.compatibility import range
@@ -213,3 +213,9 @@ def test_issue_15539():
         + O(x**(-6), (x, -oo)))
     assert series(atan(x), x, oo) == (-1/(5*x**5) + 1/(3*x**3) - 1/x + pi/2
         + O(x**(-6), (x, oo)))
+
+
+def test_issue_7259():
+    assert series(LambertW(x), x) == x - x**2 + 3*x**3/2 - 8*x**4/3 + 125*x**5/24 + O(x**6)
+    assert series(LambertW(x**2), x, n=8) == x**2 - x**4 + 3*x**6/2 + O(x**8)
+    assert series(LambertW(sin(x)), x, n=4) == x - x**2 + 4*x**3/3 + O(x**4)