diff --git a/.travis.yml b/.travis.yml
index 470b48827b8..ad13994a242 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,7 @@ python:
   - 3.7
 env:
   global:
-    - PYTEST_ADDOPTS="-n auto --cov-config setup.cfg --cov-append --cov-report="
+    - PYTEST_ADDOPTS="-n auto --cov-append --cov-report="
     - COVOPTS="--cov diofant -m \"not slow and not xfail\""
     - PYTHONHASHSEED=0
     - secure: "VNgpd/TFAkLGjWZzRBWxujDHktdvp2Fg9D8PbGmr/pbmLcJS3a5pf6po1T/\
@@ -40,6 +40,15 @@ stages:
   - deploy
 matrix:
   include:
+    - stage: test
+      python: 3.8-dev
+      env: SPLIT=1/3
+    - stage: test
+      python: 3.8-dev
+      env: SPLIT=2/3
+    - stage: test
+      python: 3.8-dev
+      env: SPLIT=3/3
     - stage: quality
       python: 3.7
       env: COVERAGE='on'
@@ -52,9 +61,11 @@ matrix:
             python setup.py build_sphinx -W -b linkcheck
           fi
         - pip uninstall -y gmpy2
+        - pip uninstall -y numpy
         - PYTEST_ADDOPTS="${PYTEST_ADDOPTS} ${COVOPTS}"
         - DIOFANT_GROUND_TYPES='gmpy' py.test diofant/tests/domains
         - travis_wait py.test diofant/tests/polys
+        - travis_wait py.test diofant/tests/external
     - stage: deploy
       env:
         -
@@ -107,10 +118,15 @@ matrix:
             tags: true
             condition: "$TRAVIS_TAG =~ ^v[0-9]+\\.[0-9]+\\.[0-9]+.*$"
 before_install:
-  - pip uninstall -y pytest
+  - pip uninstall -y pytest attrs
   - pip install --pre -U gmpy2
 install:
-  - travis_wait pip install .[exports,plot,interactive,develop,gmpy,docs]
+  - |
+    INSTALL_DEPS="interactive,develop,gmpy"
+    if [ ${TRAVIS_PYTHON_VERSION} == '3.7' ]; then
+      INSTALL_DEPS="${INSTALL_DEPS},exports,plot,docs"
+    fi
+    travis_wait pip install -U .[${INSTALL_DEPS}]
 script:
   - |
     PYTEST_ADDOPTS="${PYTEST_ADDOPTS} --split=${SPLIT}"
@@ -120,6 +136,6 @@ script:
       PYTEST_ADDOPTS="${PYTEST_ADDOPTS} --doctest-modules"
     fi
     python setup.py test
-after_success: test -n "${COVERAGE}" && codecov
+after_success: test -n "${COVERAGE}" && codecov --required
 notifications:
   email: false
diff --git a/conftest.py b/conftest.py
index 500afb0d7d9..a94e6884012 100644
--- a/conftest.py
+++ b/conftest.py
@@ -9,6 +9,13 @@
 
 
 collect_ignore = ["setup.py"]
+try:
+    import matplotlib
+    matplotlib.rc('figure', max_open_warning=0)
+    del matplotlib
+except ImportError:
+    collect_ignore.extend(["diofant/plotting/plot.py",
+                           "diofant/plotting/plot_implicit.py"])
 
 sp = re.compile(r'([0-9]+)/([1-9][0-9]*)')
 
diff --git a/diofant/combinatorics/permutations.py b/diofant/combinatorics/permutations.py
index e95ab868e52..91b55f6f0d4 100644
--- a/diofant/combinatorics/permutations.py
+++ b/diofant/combinatorics/permutations.py
@@ -1980,7 +1980,7 @@ def inversions(self):
         References
         ==========
 
-        [1] https://www.cp.eng.chula.ac.th/~piak/teaching/algo/algo2008/count-inv.htm
+        * https://www.cp.eng.chula.ac.th/~prabhas//teaching/algo/algo2008/count-inv.htm
 
         Examples
         ========
diff --git a/diofant/concrete/gosper.py b/diofant/concrete/gosper.py
index b293f7b05b6..d5cbdfbe5ee 100644
--- a/diofant/concrete/gosper.py
+++ b/diofant/concrete/gosper.py
@@ -2,7 +2,7 @@
 
 from ..core import Dummy, Integer, nan, symbols
 from ..core.compatibility import is_sequence
-from ..polys import Poly, factor, parallel_poly_from_expr
+from ..polys import factor, parallel_poly_from_expr
 from ..simplify import hypersimp
 from ..solvers import solve
 
@@ -34,8 +34,9 @@ def gosper_normal(f, g, n, polys=True):
     Examples
     ========
 
-    >>> gosper_normal(4*n + 5, 2*(4*n + 1)*(2*n + 3), n, polys=False)
-    (1/4, n + 3/2, n + 1/4)
+    >>> gosper_normal(4*n + 5, 2*(4*n + 1)*(2*n + 3), n)
+    (Poly(1/4, n, domain='QQ'), Poly(n + 3/2, n, domain='QQ'),
+     Poly(n + 1/4, n, domain='QQ'))
 
     """
     (p, q), opt = parallel_poly_from_expr((f, g), n, field=True)
@@ -56,14 +57,7 @@ def gosper_normal(f, g, n, polys=True):
         for j in range(1, i + 1):
             C *= d.shift(-j)
 
-    A *= Z
-
-    if not polys:
-        A = A.as_expr()
-        B = B.as_expr()
-        C = C.as_expr()
-
-    return A, B, C
+    return A*Z, B, C
 
 
 def gosper_term(f, n):
@@ -111,29 +105,24 @@ def gosper_term(f, n):
             D.remove(d)
 
     if not D:
-        return  # 'f(n)' is *not* Gosper-summable
+        return
 
     d = max(D)
 
     coeffs = symbols('c:%s' % (d + 1), cls=Dummy)
     domain = A.domain.inject(*coeffs)
 
-    x = Poly(coeffs, n, domain=domain)
+    x = sum(c*n**i for i, c in enumerate(reversed(coeffs)))
+    x = x.as_poly(n, domain=domain)
     H = A*x.shift(1) - B*x - C
 
     solution = solve(H.coeffs(), coeffs)
     if solution:
         solution = solution[0]
 
-    x = x.as_expr().subs(solution)
+    x = x.as_expr().subs(solution).subs({_: 0 for _ in coeffs})
 
-    for coeff in coeffs:
-        if coeff not in solution:
-            x = x.subs({coeff: 0})
-
-    if x == 0:
-        return  # 'f(n)' is *not* Gosper-summable
-    else:
+    if x != 0:
         return B.as_expr()*x/C.as_expr()
 
 
@@ -152,17 +141,8 @@ def gosper_sum(f, k):
     Examples
     ========
 
-    >>> from diofant.abc import i
-
-    >>> f = (4*k + 1)*factorial(k)/factorial(2*k + 1)
-    >>> gosper_sum(f, (k, 0, n))
+    >>> gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n))
     (-factorial(n) + 2*factorial(2*n + 1))/factorial(2*n + 1)
-    >>> _.subs({n: 2}) == sum(f.subs({k: i}) for i in [0, 1, 2])
-    True
-    >>> gosper_sum(f, (k, 3, n))
-    (-60*factorial(n) + factorial(2*n + 1))/(60*factorial(2*n + 1))
-    >>> _.subs({n: 5}) == sum(f.subs({k: i}) for i in [3, 4, 5])
-    True
 
     References
     ==========
@@ -188,9 +168,6 @@ def gosper_sum(f, k):
         result = (f*(g + 1)).subs({k: b}) - (f*g).subs({k: a})
 
         if result is nan:
-            try:
-                result = (f*(g + 1)).limit(k, b) - (f*g).limit(k, a)
-            except NotImplementedError:  # pragma: no cover
-                result = None
+            result = (f*(g + 1)).limit(k, b) - (f*g).limit(k, a)
 
     return factor(result)
diff --git a/diofant/core/add.py b/diofant/core/add.py
index d3ecb77b2b2..d59e6bb8147 100644
--- a/diofant/core/add.py
+++ b/diofant/core/add.py
@@ -413,16 +413,9 @@ def _eval_is_extended_real(self):
     def _eval_is_complex(self):
         return _fuzzy_group((a.is_complex for a in self.args), quick_exit=True)
 
-    def _eval_is_antihermitian(self):
-        return _fuzzy_group((a.is_antihermitian for a in self.args),
-                            quick_exit=True)
-
     def _eval_is_finite(self):
         return _fuzzy_group((a.is_finite for a in self.args), quick_exit=True)
 
-    def _eval_is_hermitian(self):
-        return _fuzzy_group((a.is_hermitian for a in self.args), quick_exit=True)
-
     def _eval_is_integer(self):
         return _fuzzy_group((a.is_integer for a in self.args), quick_exit=True)
 
diff --git a/diofant/core/assumptions.py b/diofant/core/assumptions.py
index 221d34ec50d..956fe01398c 100644
--- a/diofant/core/assumptions.py
+++ b/diofant/core/assumptions.py
@@ -56,11 +56,11 @@
     'real           ==  extended_real & finite',
     'rational       ->  algebraic',
     'algebraic      ->  complex',
-    'real           ->  complex & hermitian',
-    'imaginary      ->  complex & antihermitian',
+    'real           ->  complex',
+    'imaginary      ->  complex',
     'complex        ->  finite & commutative',
     'extended_real  ->  commutative',
-    'extended_real  -> real | infinite',
+    'extended_real  ->  real | infinite',
 
     'odd            ==  integer & ~even',
     'even           ==  integer & ~odd',
@@ -321,10 +321,6 @@
     """Test if self can have only positive values.""",
     'nonpositive':
     """Test if self can have only nonpositive values.""",
-    'hermitian':
-    """Test if self belongs to the field of hermitian operators.""",
-    'antihermitian':
-    """Test if self belongs to the field of antihermitian operators.""",
 }
 
 
diff --git a/diofant/core/basic.py b/diofant/core/basic.py
index ce23fe51296..4fec222b50c 100644
--- a/diofant/core/basic.py
+++ b/diofant/core/basic.py
@@ -1,6 +1,7 @@
 """Base class for all the objects in Diofant."""
 
-from collections import Mapping, defaultdict
+from collections import defaultdict
+from collections.abc import Mapping
 from itertools import zip_longest
 
 from .cache import cacheit
@@ -161,7 +162,7 @@ def _repr_pretty_(self, p, cycle):
 
     def _repr_latex_(self):
         from ..printing import latex
-        return latex(self, mode='equation')
+        return latex(self, mode='equation*')
 
     def atoms(self, *types):
         """Returns the atoms that form the current object.
diff --git a/diofant/core/containers.py b/diofant/core/containers.py
index da0ef8eefd3..d0e3b1d9386 100644
--- a/diofant/core/containers.py
+++ b/diofant/core/containers.py
@@ -6,7 +6,7 @@
 They are supposed to work seamlessly within the Diofant framework.
 """
 
-from collections import Mapping
+from collections.abc import Mapping
 
 from .basic import Basic
 from .compatibility import as_int, iterable
diff --git a/diofant/core/expr.py b/diofant/core/expr.py
index 35469854c37..e07c2cd1c6c 100644
--- a/diofant/core/expr.py
+++ b/diofant/core/expr.py
@@ -704,13 +704,8 @@ def conjugate(self):
         return c(self)
 
     def _eval_transpose(self):
-        from ..functions.elementary.complexes import conjugate
         if self.is_complex or self.is_extended_real:
             return self
-        elif self.is_hermitian:
-            return conjugate(self)
-        elif self.is_antihermitian:
-            return -conjugate(self)
 
     def transpose(self):
         """Transpose self.
@@ -726,10 +721,6 @@ def transpose(self):
 
     def _eval_adjoint(self):
         from ..functions.elementary.complexes import conjugate, transpose
-        if self.is_hermitian:
-            return self
-        elif self.is_antihermitian:
-            return -self
         obj = self._eval_conjugate()
         if obj is not None:
             return transpose(obj)
diff --git a/diofant/core/mul.py b/diofant/core/mul.py
index 9394886e7dd..fa34f897bf4 100644
--- a/diofant/core/mul.py
+++ b/diofant/core/mul.py
@@ -1034,38 +1034,6 @@ def _eval_is_imaginary(self):
         else:
             return obj.is_real
 
-    def _eval_is_hermitian(self):
-        hermitian = True
-        one_nc = zero = False
-
-        for t in self.args:
-            if not t.is_commutative:
-                if one_nc:
-                    return
-                one_nc = True
-
-            if t.is_antihermitian or t.is_hermitian:
-                if t.is_antihermitian:
-                    hermitian = not hermitian
-                z = t.is_zero
-                if not z and zero is False:
-                    zero = z
-                elif z:
-                    if self.is_finite:
-                        return True
-                    return
-            else:
-                return
-
-        if zero is False or hermitian:
-            return hermitian
-
-    def _eval_is_antihermitian(self):
-        if self.is_zero:
-            return False
-        elif self.is_nonzero:
-            return (I*self).is_hermitian
-
     def _eval_is_irrational(self):
         for t in self.args:
             a = t.is_irrational
diff --git a/diofant/core/numbers.py b/diofant/core/numbers.py
index 02b353e3e4d..3ea1c6ddc3e 100644
--- a/diofant/core/numbers.py
+++ b/diofant/core/numbers.py
@@ -1074,6 +1074,9 @@ def _eval_power(self, expt):
     def _as_mpf_val(self, prec):
         return mlib.from_rational(self.numerator, self.denominator, prec, rnd)
 
+    def _mpmath_(self, prec, rnd):
+        return mpmath.make_mpf(mlib.from_rational(self.numerator, self.denominator, prec, rnd))
+
     def __abs__(self):
         return Rational(abs(self.numerator), self.denominator)
 
@@ -2563,6 +2566,16 @@ def as_base_exp(self):
 
 converter[float] = converter[decimal.Decimal] = Float
 
+try:
+    import numpy as np
+
+    def _sympify_numpy(x):
+        return Float(mpmath.mpmathify(x))
+
+    converter[np.floating] = _sympify_numpy
+except ImportError:
+    pass
+
 converter[int] = Integer
 converter[fractions.Fraction] = Rational
 
diff --git a/diofant/functions/elementary/complexes.py b/diofant/functions/elementary/complexes.py
index b95628fe980..598c9464edf 100644
--- a/diofant/functions/elementary/complexes.py
+++ b/diofant/functions/elementary/complexes.py
@@ -97,6 +97,10 @@ def _eval_rewrite_as_im(self, arg):
     def _eval_is_algebraic(self):
         return self.args[0].is_algebraic
 
+    def _eval_is_real(self):
+        if self.args[0].is_complex:
+            return True
+
 
 class im(Function):
     """Returns imaginary part of expression.
@@ -189,6 +193,10 @@ def _eval_rewrite_as_re(self, arg):
     def _eval_is_algebraic(self):
         return self.args[0].is_algebraic
 
+    def _eval_is_real(self):
+        if self.args[0].is_complex:
+            return True
+
 
 ###############################################################################
 # ############# SIGN, ABSOLUTE VALUE, ARGUMENT and CONJUGATION ############## #
diff --git a/diofant/functions/elementary/hyperbolic.py b/diofant/functions/elementary/hyperbolic.py
index 1a44889276d..3d501f2c3db 100644
--- a/diofant/functions/elementary/hyperbolic.py
+++ b/diofant/functions/elementary/hyperbolic.py
@@ -919,6 +919,16 @@ def inverse(self, argindex=1):
         """Returns the inverse of this function."""
         return cosh
 
+    def _eval_rewrite_as_log(self, x):
+        return log(x + sqrt(x - 1)*sqrt(x + 1))
+
+    def _eval_nseries(self, x, n, logx):
+        x0 = self.args[0].limit(x, 0)
+        if x0 == 1:
+            return self._eval_rewrite_as_log(self.args[0])._eval_nseries(x, n, logx)
+        else:
+            return super()._eval_nseries(x, n, logx)
+
 
 class atanh(Function):
     """
diff --git a/diofant/functions/special/delta_functions.py b/diofant/functions/special/delta_functions.py
index e9429d6e15c..7d5b3951986 100644
--- a/diofant/functions/special/delta_functions.py
+++ b/diofant/functions/special/delta_functions.py
@@ -190,6 +190,8 @@ class Heaviside(Function):
 
     """
 
+    is_real = True
+
     def fdiff(self, argindex=1):
         if argindex == 1:
             return DiracDelta(self.args[0])
@@ -215,7 +217,3 @@ def _eval_rewrite_as_Piecewise(self, arg):
     def _eval_rewrite_as_sign(self, arg):
         if arg.is_extended_real:
             return (sign(arg)+1)/2
-
-    def _eval_is_real(self):
-        if self.args[0].is_extended_real:
-            return True
diff --git a/diofant/integrals/integrals.py b/diofant/integrals/integrals.py
index fd170e260b5..da6856f631e 100644
--- a/diofant/integrals/integrals.py
+++ b/diofant/integrals/integrals.py
@@ -483,7 +483,7 @@ def try_meijerg(function, xab):
                 if len(xab) >= 2:
                     if (any(b.is_extended_real for b in xab[1:]) and
                             not any(b.is_extended_real is False for b in xab[1:])):
-                        r = Dummy('r', extended_real=True)
+                        r = Dummy('r', real=True)
                         function = function.subs({xab[0]: r})
                         function = function.rewrite(Piecewise)
                         function = function.subs({r: xab[0]})
diff --git a/diofant/plotting/plot.py b/diofant/plotting/plot.py
index cd0c47a749a..948d46fe9a5 100644
--- a/diofant/plotting/plot.py
+++ b/diofant/plotting/plot.py
@@ -23,7 +23,7 @@
 """
 
 import warnings
-from collections import Callable
+from collections.abc import Callable
 from inspect import getfullargspec
 
 from ..core import Dummy, Expr, Symbol, Tuple, sympify
diff --git a/diofant/polys/compatibility.py b/diofant/polys/compatibility.py
index c518dcf425f..57f8c797bc8 100644
--- a/diofant/polys/compatibility.py
+++ b/diofant/polys/compatibility.py
@@ -32,7 +32,7 @@
                             dup_isolate_all_roots, dup_isolate_all_roots_sqf,
                             dup_isolate_complex_roots_sqf,
                             dup_isolate_real_roots,
-                            dup_isolate_real_roots_list,
+                            dup_isolate_real_roots_pair,
                             dup_isolate_real_roots_sqf, dup_refine_real_root,
                             dup_root_upper_bound, dup_sign_variations,
                             dup_sturm)
@@ -424,8 +424,8 @@ def dup_isolate_real_roots_sqf(self, f, eps=None, inf=None, sup=None, blackbox=F
     def dup_isolate_real_roots(self, f, eps=None, inf=None, sup=None):
         return dup_isolate_real_roots(self.to_dense(f), self.domain, eps=eps, inf=inf, sup=sup)
 
-    def dup_isolate_real_roots_list(self, polys, eps=None, inf=None, sup=None, strict=False, basis=False):
-        return dup_isolate_real_roots_list(list(map(self.to_dense, polys)), self.domain, eps=eps, inf=inf, sup=sup, strict=strict, basis=basis)
+    def dup_isolate_real_roots_pair(self, f, g, eps=None, inf=None, sup=None, strict=False, basis=False):
+        return dup_isolate_real_roots_pair(*(map(self.to_dense, [f, g])), self.domain, eps=eps, inf=inf, sup=sup, strict=strict, basis=basis)
 
     def dup_count_real_roots(self, f, inf=None, sup=None):
         return dup_count_real_roots(self.to_dense(f), self.domain, inf=inf, sup=sup)
diff --git a/diofant/polys/polytools.py b/diofant/polys/polytools.py
index 1e17e8a5db3..9b3a792becb 100644
--- a/diofant/polys/polytools.py
+++ b/diofant/polys/polytools.py
@@ -4488,43 +4488,12 @@ def intervals(F, all=False, eps=None, inf=None, sup=None, strict=False, sqf=Fals
     [((-26/15, -19/11), 1), ((19/11, 26/15), 1)]
 
     """
-    if not hasattr(F, '__iter__'):
-        try:
-            F = Poly(F)
-        except GeneratorsNeeded:
-            return []
-
-        return F.intervals(all=all, eps=eps, inf=inf, sup=sup, sqf=sqf)
-    else:
-        polys, opt = parallel_poly_from_expr(F, domain='QQ')
-
-        if len(opt.gens) > 1:
-            raise MultivariatePolynomialError
-
-        R = polys[0].rep.ring
-        polys = [p.rep for p in polys]
-
-        if eps is not None:
-            eps = opt.domain.convert(eps)
-
-            if eps <= 0:
-                raise ValueError("'eps' must be a positive rational")
-
-        if inf is not None:
-            inf = opt.domain.convert(inf)
-        if sup is not None:
-            sup = opt.domain.convert(sup)
-
-        intervals = R.dup_isolate_real_roots_list(polys, eps=eps, inf=inf, sup=sup,
-                                                  strict=strict)
-
-        result = []
-
-        for (s, t), indices in intervals:
-            s, t = opt.domain.to_expr(s), opt.domain.to_expr(t)
-            result.append(((s, t), indices))
+    try:
+        F = Poly(F)
+    except GeneratorsNeeded:
+        return []
 
-        return result
+    return F.intervals(all=all, eps=eps, inf=inf, sup=sup, sqf=sqf)
 
 
 def refine_root(f, s, t, eps=None, steps=None, check_sqf=False):
diff --git a/diofant/polys/rootisolation.py b/diofant/polys/rootisolation.py
index 83a6d09fa99..c91e9027ed5 100644
--- a/diofant/polys/rootisolation.py
+++ b/diofant/polys/rootisolation.py
@@ -458,14 +458,8 @@ def dup_isolate_real_roots_sqf(f, K, eps=None, inf=None, sup=None, blackbox=Fals
     f = dmp_clear_denoms(f, 0, K)[1]
 
     if K.is_ComplexAlgebraicField and not K.is_RealAlgebraicField:
-        A, K = K, K.domain
-        polys = [dmp_eval_in(_, K.zero, 1, 1, K) for _ in dup_real_imag(f, A)]
-        if not polys[1]:
-            f = polys[0]
-        else:
-            roots = dup_isolate_real_roots_list(polys, K, eps=eps, inf=inf, sup=sup, strict=True)
-            roots = [_[0] for _ in roots if _[1].keys() == {0, 1}]
-            return [RealInterval((a, b), f, K) for (a, b) in roots] if blackbox else roots
+        roots = [r for r, _ in dup_isolate_real_roots(f, K, eps=eps, inf=inf, sup=sup)]
+        return [RealInterval((a, b), f, K) for (a, b) in roots] if blackbox else roots
 
     if dmp_degree_in(f, 0, 0) <= 0:
         return []
@@ -496,6 +490,12 @@ def dup_isolate_real_roots(f, K, eps=None, inf=None, sup=None):
     if not (K.is_ComplexAlgebraicField or K.is_RationalField):
         raise DomainError("isolation of real roots not supported over %s" % K)
 
+    if K.is_ComplexAlgebraicField and not K.is_RealAlgebraicField:
+        A, K = K, K.domain
+        polys = [dmp_eval_in(_, K.zero, 1, 1, K) for _ in dup_real_imag(f, A)]
+        roots = dup_isolate_real_roots_pair(*polys, K, eps=eps, inf=inf, sup=sup, strict=True)
+        return [(_[0], _[1][0]) for _ in roots if _[1].keys() == {0, 1}]
+
     _, factors = dmp_sqf_list(f, 0, K)
     factors = [(dmp_clear_denoms(f, 0, K)[1], k) for f, k in factors]
 
@@ -508,7 +508,7 @@ def dup_isolate_real_roots(f, K, eps=None, inf=None, sup=None):
         return sorted(I_neg + I_zero + I_pos)
 
 
-def dup_isolate_real_roots_list(polys, K, eps=None, inf=None, sup=None, strict=False, basis=False):
+def dup_isolate_real_roots_pair(f, g, K, eps=None, inf=None, sup=None, strict=False, basis=False):
     """Isolate real roots of a list of polynomials."""
     R, K = K, K.field
 
@@ -520,6 +520,7 @@ def dup_isolate_real_roots_list(polys, K, eps=None, inf=None, sup=None, strict=F
     else:
         zeros = False
 
+    polys = [f, g]
     gcd = []
 
     for i, p in enumerate(polys):
@@ -534,16 +535,18 @@ def dup_isolate_real_roots_list(polys, K, eps=None, inf=None, sup=None, strict=F
 
     polys = [dmp_quo(p, gcd, 0, K) for p in polys]
 
-    if len(polys) != 2:
-        raise NotImplementedError
-
     factors = collections.defaultdict(dict)
 
     for i, p in enumerate(polys):
+        ni = (i + 1) % 2
+
+        if not p and zeros and ni in zero_indices:
+            zero_indices[i] = zero_indices[ni]
+
         for f, _ in dmp_sqf_list(dmp_gcd(p, gcd, 0, K), 0, K)[1]:
             k1 = dmp_trial_division(gcd, [f], 0, K)[0][1]
             k2 = dmp_trial_division(p, [f], 0, K)[0][1]
-            factors[tuple(f)] = {i: k1 + k2, (i + 1) % 2: k1}
+            factors[tuple(f)] = {i: k1 + k2, ni: k1}
 
             gcd = dmp_quo(gcd, dmp_pow(f, k1, 0, K), 0, K)
             p = dmp_quo(p, dmp_pow(f, k2, 0, K), 0, K)
@@ -1173,15 +1176,10 @@ def _get_rectangle(f1, f2, F, inf, sup, exclude=None):
     f1L4 = dmp_eval_in(f1, u, 0, 1, F)
     f2L4 = dmp_eval_in(f2, u, 0, 1, F)
 
-    S_L1 = [f1L1, f2L1]
-    S_L2 = [f1L2, f2L2]
-    S_L3 = [f1L3, f2L3]
-    S_L4 = [f1L4, f2L4]
-
-    I_L1 = dup_isolate_real_roots_list(S_L1, F, inf=u, sup=s, basis=True, strict=True)
-    I_L2 = dup_isolate_real_roots_list(S_L2, F, inf=v, sup=t, basis=True, strict=True)
-    I_L3 = dup_isolate_real_roots_list(S_L3, F, inf=u, sup=s, basis=True, strict=True)
-    I_L4 = dup_isolate_real_roots_list(S_L4, F, inf=v, sup=t, basis=True, strict=True)
+    I_L1 = dup_isolate_real_roots_pair(f1L1, f2L1, F, inf=u, sup=s, basis=True, strict=True)
+    I_L2 = dup_isolate_real_roots_pair(f1L2, f2L2, F, inf=v, sup=t, basis=True, strict=True)
+    I_L3 = dup_isolate_real_roots_pair(f1L3, f2L3, F, inf=u, sup=s, basis=True, strict=True)
+    I_L4 = dup_isolate_real_roots_pair(f1L4, f2L4, F, inf=v, sup=t, basis=True, strict=True)
 
     I_L3 = _reverse_intervals(I_L3)
     I_L4 = _reverse_intervals(I_L4)
@@ -1253,7 +1251,7 @@ def _vertical_bisection(N, a, b, I, Q, F1, F2, f1, f2, F):
     f1V = dmp_eval_in(f1, x, 0, 1, F)
     f2V = dmp_eval_in(f2, x, 0, 1, F)
 
-    I_V = dup_isolate_real_roots_list([f1V, f2V], F, inf=v, sup=t, strict=True, basis=True)
+    I_V = dup_isolate_real_roots_pair(f1V, f2V, F, inf=v, sup=t, strict=True, basis=True)
 
     I_L1_L, I_L1_R = [], []
     I_L2_L, I_L2_R = I_V, I_L2
@@ -1360,7 +1358,7 @@ def _horizontal_bisection(N, a, b, I, Q, F1, F2, f1, f2, F):
     f1H = dmp_eval_in(f1, y, 1, 1, F)
     f2H = dmp_eval_in(f2, y, 1, 1, F)
 
-    I_H = dup_isolate_real_roots_list([f1H, f2H], F, inf=u, sup=s, strict=True, basis=True)
+    I_H = dup_isolate_real_roots_pair(f1H, f2H, F, inf=u, sup=s, strict=True, basis=True)
 
     I_L1_B, I_L1_U = I_L1, I_H
     I_L2_B, I_L2_U = [], []
diff --git a/diofant/printing/__init__.py b/diofant/printing/__init__.py
index e670192021b..aee84783e5e 100644
--- a/diofant/printing/__init__.py
+++ b/diofant/printing/__init__.py
@@ -1,6 +1,6 @@
 """Printing subsystem"""
 
-from .ccode import ccode, print_ccode
+from .ccode import ccode
 from .dot import dotprint
 from .fcode import fcode
 from .latex import latex
@@ -8,7 +8,7 @@
 from .mathml import mathml
 from .octave import octave_code
 from .pretty import pprint, pprint_use_unicode, pretty, pretty_print
-from .python import print_python, python
+from .python import python
 from .repr import srepr
 from .str import StrPrinter, sstr, sstrrepr
 
diff --git a/diofant/printing/ccode.py b/diofant/printing/ccode.py
index 070bff21402..9f1f9ba5ddc 100644
--- a/diofant/printing/ccode.py
+++ b/diofant/printing/ccode.py
@@ -378,8 +378,3 @@ def ccode(expr, assign_to=None, **settings):
     """
 
     return CCodePrinter(settings).doprint(expr, assign_to)
-
-
-def print_ccode(expr, **settings):
-    """Prints C representation of the given expression."""
-    print(ccode(expr, **settings))
diff --git a/diofant/printing/mathematica.py b/diofant/printing/mathematica.py
index 63c8f774062..876a3bc2582 100644
--- a/diofant/printing/mathematica.py
+++ b/diofant/printing/mathematica.py
@@ -19,10 +19,15 @@
     "asin": [(lambda x: True, "ArcSin")],
     "acos": [(lambda x: True, "ArcCos")],
     "atan": [(lambda x: True, "ArcTan")],
+    "acot": [(lambda x: True, "ArcCot")],
     "sinh": [(lambda x: True, "Sinh")],
     "cosh": [(lambda x: True, "Cosh")],
     "tanh": [(lambda x: True, "Tanh")],
     "coth": [(lambda x: True, "Coth")],
+    "asinh": [(lambda x: True, "ArcSinh")],
+    "acosh": [(lambda x: True, "ArcCosh")],
+    "atanh": [(lambda x: True, "ArcTanh")],
+    "acoth": [(lambda x: True, "ArcCoth")],
     "sech": [(lambda x: True, "Sech")],
     "csch": [(lambda x: True, "Csch")],
     "sign": [(lambda x: True, "Sign")],
@@ -40,12 +45,10 @@
     "factorial2": [(lambda *x: True, "Factorial2")],
     "RisingFactorial": [(lambda x, k: True, "Pochhammer")],
     "gamma": [(lambda x: True, "Gamma")],
-    "asinh": [(lambda x: True, "ArcSinh")],
     "zeta": [(lambda x: True, "Zeta")],
     "Heaviside": [(lambda x: True, "UnitStep")],
     "fibonacci": [(lambda x: True, "Fibonacci")],
     "polylog": [(lambda x, y: True, "PolyLog")],
-    "atanh": [(lambda x: True, "ArcTanh")],
 }
 
 
diff --git a/diofant/printing/pretty/pretty.py b/diofant/printing/pretty/pretty.py
index 4068af5e9db..8abbf722eda 100644
--- a/diofant/printing/pretty/pretty.py
+++ b/diofant/printing/pretty/pretty.py
@@ -1292,12 +1292,12 @@ def _print_nth_root(self, base, expt):
         _zZ = xobj('/', 1)
         rootsign = xobj('\\', 1) + _zZ
         # Make exponent number to put above it
-        if isinstance(expt, Rational):
+        if expt.is_Rational:
             exp = str(expt.denominator)
             if exp == '2':
                 exp = ''
         else:
-            exp = str(expt.args[0])
+            exp = str(self._print(expt.args[0]))
         exp = exp.ljust(2)
         if len(exp) > 2:
             rootsign = ' '*(len(exp) - 2) + rootsign
diff --git a/diofant/printing/pretty/stringpict.py b/diofant/printing/pretty/stringpict.py
index 9daee5b9d31..243c0b6006b 100644
--- a/diofant/printing/pretty/stringpict.py
+++ b/diofant/printing/pretty/stringpict.py
@@ -1,4 +1,5 @@
 """Prettyprinter by Jurjen Bos.
+
 (I hate spammers: mail me at pietjepuk314 at the reverse of ku.oc.oohay).
 All objects have a method that create a "stringPict",
 that can be used in the str method for pretty printing.
@@ -10,6 +11,7 @@
 TODO:
     - Allow left/center/right alignment options for above/below and
       top/center/bottom alignment options for left/right
+
 """
 
 from .pretty_symbology import hobj, pretty_use_unicode, vobj, xsym
diff --git a/diofant/printing/python.py b/diofant/printing/python.py
index 7aa708d001a..a4883cd4026 100644
--- a/diofant/printing/python.py
+++ b/diofant/printing/python.py
@@ -32,7 +32,7 @@ def _print_Function(self, expr):
             self.functions.append(func)
         return StrPrinter._print_Function(self, expr)
 
-    # procedure (!) for defining symbols which have be defined in print_python()
+    # procedure (!) for defining symbols which have be defined in python()
     def _print_Symbol(self, expr):
         symbol = self._str(expr)
         if symbol not in self.symbols:
@@ -80,8 +80,3 @@ def python(expr, **settings):
         exprp = expr.subs(renamings)
     result += 'e = ' + printer._str(exprp)
     return result
-
-
-def print_python(expr, **settings):
-    """Print output of python() function"""
-    print(python(expr, **settings))
diff --git a/diofant/series/limits.py b/diofant/series/limits.py
index 61f756199e0..9c7294ce9b8 100644
--- a/diofant/series/limits.py
+++ b/diofant/series/limits.py
@@ -140,10 +140,12 @@ def doit(self, **hints):
 
         use_heuristics = hints.get('heuristics', True)
 
-        has_Floats = e.has(Float)
+        has_Floats = e.has(Float) or z0.has(Float)
         if has_Floats:
             e = e.subs({k: Rational(k) for k in e.atoms(Float)},
                        simultaneous=True)
+            z0 = z0.subs({k: Rational(k) for k in z0.atoms(Float)},
+                         simultaneous=True)
 
         if z0.has(z):
             newz = z.as_dummy()
diff --git a/diofant/solvers/diophantine.py b/diofant/solvers/diophantine.py
index 079cab407a1..07a9f5ade37 100644
--- a/diofant/solvers/diophantine.py
+++ b/diofant/solvers/diophantine.py
@@ -734,7 +734,7 @@ def diop_quadratic(eq, param=symbols("t", integer=True)):
     ==========
 
     * Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online],
-      Available: https://www.alpertron.com.ar/METHODS.HTM
+      Available: https://web.archive.org/web/20181231080858/https://www.alpertron.com.ar/METHODS.HTM
     * Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
       Available: https://web.archive.org/web/20180831180321/http://www.jpr2718.org/ax2p.pdf
 
@@ -769,7 +769,7 @@ def _diop_quadratic(var, coeff, t):
     # (1) Simple-Hyperbolic case: A = C = 0, B != 0
     # In this case equation can be converted to (Bx + E)(By + D) = DE - BF
     # We consider two cases; DE - BF = 0 and DE - BF != 0
-    # More details, https://www.alpertron.com.ar/METHODS.HTM#SHyperb
+    # More details, https://web.archive.org/web/20181231080858/https://www.alpertron.com.ar/METHODS.HTM#SHyperb
 
     sol = set()
     discr = B**2 - 4*A*C
@@ -797,7 +797,7 @@ def _diop_quadratic(var, coeff, t):
     # (2) Parabolic case: B**2 - 4*A*C = 0
     # There are two subcases to be considered in this case.
     # sqrt(c)D - sqrt(a)E = 0 and sqrt(c)D - sqrt(a)E != 0
-    # More Details, https://www.alpertron.com.ar/METHODS.HTM#Parabol
+    # More Details, https://web.archive.org/web/20181231080858/https://www.alpertron.com.ar/METHODS.HTM#Parabol
 
     elif discr == 0:
 
diff --git a/diofant/solvers/solvers.py b/diofant/solvers/solvers.py
index 1511fd3359a..d5d5fc08efd 100644
--- a/diofant/solvers/solvers.py
+++ b/diofant/solvers/solvers.py
@@ -310,7 +310,7 @@ def _sympified_list(w):
     ###########################################################################
     for i, fi in enumerate(f):
         if isinstance(fi, Equality):
-            if 'ImmutableMatrix' in [type(a).__name__ for a in fi.args]:
+            if 'ImmutableMatrix' in (type(a).__name__ for a in fi.args):
                 f[i] = fi.lhs - fi.rhs
             else:
                 f[i] = Add(fi.lhs, -fi.rhs, evaluate=False)
diff --git a/diofant/tests/concrete/test_gosper.py b/diofant/tests/concrete/test_gosper.py
index e6d93bcb572..655cc1643e5 100644
--- a/diofant/tests/concrete/test_gosper.py
+++ b/diofant/tests/concrete/test_gosper.py
@@ -12,8 +12,6 @@
 def test_gosper_normal():
     assert gosper_normal(4*n + 5, 2*(4*n + 1)*(2*n + 3), n) == \
         (Poly(Rational(1, 4), n), Poly(n + Rational(3, 2)), Poly(n + Rational(1, 4)))
-    assert gosper_normal(4*n + 5, 2*(4*n + 1)*(2*n + 3), n, polys=False) == \
-        (Rational(1, 4), n + Rational(3, 2), n + Rational(1, 4))
 
 
 def test_gosper_term():
@@ -43,8 +41,11 @@ def test_gosper_sum():
     assert gosper_sum((
         -1)**k*binomial(n, k), (k, 0, m)) == -(-1)**m*(m - n)*binomial(n, m)/n
 
-    assert gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n)) == \
-        (2*factorial(2*n + 1) - factorial(n))/factorial(2*n + 1)
+    f = (4*k + 1)*factorial(k)/factorial(2*k + 1)
+    assert gosper_sum(f, (k, 0, n)) == (2*factorial(2*n + 1) -
+                                        factorial(n))/factorial(2*n + 1)
+    assert gosper_sum(f, (k, 3, n)) == (-60*factorial(n) +
+                                        factorial(2*n + 1))/(60*factorial(2*n + 1))
 
     # issue sympy/sympy#6033:
     assert gosper_sum(
diff --git a/diofant/tests/core/test_arit.py b/diofant/tests/core/test_arit.py
index c2e53fe5bfe..966b59749c7 100644
--- a/diofant/tests/core/test_arit.py
+++ b/diofant/tests/core/test_arit.py
@@ -1260,35 +1260,6 @@ def test_Mul_is_imaginary_real():
     assert (b*nr).is_extended_real is None
 
 
-def test_Mul_hermitian_antihermitian():
-    a = Symbol('a', hermitian=True, zero=False)
-    b = Symbol('b', hermitian=True)
-    c = Symbol('c', hermitian=False)
-    d = Symbol('d', antihermitian=True)
-    e1 = Mul(a, b, c, evaluate=False)
-    e2 = Mul(b, a, c, evaluate=False)
-    e3 = Mul(a, b, c, d, evaluate=False)
-    e4 = Mul(b, a, c, d, evaluate=False)
-    e5 = Mul(a, c, evaluate=False)
-    e6 = Mul(a, c, d, evaluate=False)
-    assert e1.is_hermitian is None
-    assert e2.is_hermitian is None
-    assert e1.is_antihermitian is None
-    assert e2.is_antihermitian is None
-    assert e3.is_antihermitian is None
-    assert e4.is_antihermitian is None
-    assert e5.is_antihermitian is None
-    assert e6.is_antihermitian is None
-
-    z = Symbol('z', zero=True)
-    e = Symbol('e', antihermitian=True, finite=True)
-    assert (z*e).is_antihermitian is False
-    assert (z*e).is_hermitian is True
-    A = Symbol('A', hermitian=True, commutative=False)
-    B = Symbol('B', hermitian=True, commutative=False)
-    assert (A*B).is_hermitian is None
-
-
 def test_Add_is_comparable():
     assert (x + y).is_comparable is False
     assert (x + 1).is_comparable is False
diff --git a/diofant/tests/core/test_assumptions.py b/diofant/tests/core/test_assumptions.py
index 852ea60d32d..78bad821dc3 100644
--- a/diofant/tests/core/test_assumptions.py
+++ b/diofant/tests/core/test_assumptions.py
@@ -907,3 +907,15 @@ def test_sympyissue_16530():
     e = 1/abs(x)
     assert e.is_real is None
     assert e.is_extended_real is None
+
+
+def test_sympyissue_17555():
+    x = Symbol('x', infinite=True, extended_real=True)
+    assert x.is_positive is None
+    assert (-x).is_positive is None
+
+
+def test_sympyissue_17556():
+    z = I*oo
+    assert z.is_imaginary is False
+    assert z.is_finite is False
diff --git a/diofant/tests/core/test_noncommutative.py b/diofant/tests/core/test_noncommutative.py
index 83ad70922f6..095c42258a7 100644
--- a/diofant/tests/core/test_noncommutative.py
+++ b/diofant/tests/core/test_noncommutative.py
@@ -11,8 +11,6 @@
 __all__ = ()
 
 A, B, C = symbols("A B C", commutative=False)
-X = symbols("X", commutative=False, hermitian=True)
-Y = symbols("Y", commutative=False, antihermitian=True)
 
 
 def test_adjoint():
@@ -23,16 +21,6 @@ def test_adjoint():
     assert adjoint(A*B - B*A) == adjoint(B)*adjoint(A) - adjoint(A)*adjoint(B)
     assert adjoint(A + I*B) == adjoint(A) - I*adjoint(B)
 
-    assert adjoint(X) == X
-    assert adjoint(-I*X) == I*X
-    assert adjoint(Y) == -Y
-    assert adjoint(-I*Y) == -I*Y
-
-    assert adjoint(X) == conjugate(transpose(X))
-    assert adjoint(Y) == conjugate(transpose(Y))
-    assert adjoint(X) == transpose(conjugate(X))
-    assert adjoint(Y) == transpose(conjugate(Y))
-
     assert adjoint(2**x) == 2**adjoint(x)
     assert adjoint(x**pi) == adjoint(x**pi, evaluate=False)
 
@@ -124,13 +112,6 @@ def test_transpose():
         transpose(B)*transpose(A) - transpose(A)*transpose(B)
     assert transpose(A + I*B) == transpose(A) + I*transpose(B)
 
-    assert transpose(X) == conjugate(X)
-    assert transpose(-I*X) == -I*conjugate(X)
-    assert transpose(Y) == -conjugate(Y)
-    assert transpose(-I*Y) == I*conjugate(Y)
-
-    assert transpose(X**pi) == transpose(X**pi, evaluate=False)
-
 
 def test_trigsimp():
     assert trigsimp(A*sin(x)**2 + A*cos(x)**2) == A
diff --git a/diofant/tests/core/test_numbers.py b/diofant/tests/core/test_numbers.py
index add87dd7162..cb6560d0665 100644
--- a/diofant/tests/core/test_numbers.py
+++ b/diofant/tests/core/test_numbers.py
@@ -1235,6 +1235,8 @@ def test_conversion_to_mpmath():
     assert mpmath.mpmathify(Rational(1, 2)) == mpmath.mpf(0.5)
     assert mpmath.mpmathify(Float('1.23', 15)) == mpmath.mpf('1.23')
 
+    assert mpmath.mpf(Rational(1, 3)) == mpmath.mpf('0.33333333333333331')
+
 
 def test_relational():
     # real
diff --git a/diofant/tests/external/test_autowrap.py b/diofant/tests/external/test_autowrap.py
index d338b5f3ba6..ae4a39b5046 100644
--- a/diofant/tests/external/test_autowrap.py
+++ b/diofant/tests/external/test_autowrap.py
@@ -172,21 +172,25 @@ def test_wrap_twice_c_cython():
 
 
 @with_cython
+@with_numpy
 def test_autowrap_trace_C_Cython():
     runtest_autowrap_trace('C', 'cython')
 
 
 @with_cython
+@with_numpy
 def test_autowrap_matrix_vector_C_cython():
     runtest_autowrap_matrix_vector('C', 'cython')
 
 
 @with_cython
+@with_numpy
 def test_autowrap_matrix_matrix_C_cython():
     runtest_autowrap_matrix_matrix('C', 'cython')
 
 
 @with_cython
+@with_numpy
 def test_ufuncify_C_Cython():
     runtest_ufuncify('C', 'cython')
 
@@ -199,6 +203,7 @@ def test_sympyissue_10274_C_cython():
 # Numpy
 
 @with_cython
+@with_numpy
 def test_ufuncify_numpy():
     # This test doesn't use Cython, but if Cython works, then there is a valid
     # C compiler, which is needed.
diff --git a/diofant/tests/external/test_numpy.py b/diofant/tests/external/test_numpy.py
index cb581356a46..741045b16f7 100644
--- a/diofant/tests/external/test_numpy.py
+++ b/diofant/tests/external/test_numpy.py
@@ -13,7 +13,7 @@
 import diofant
 from diofant import (DeferredVector, Float, Integer, Matrix, MatrixSymbol,
                      Rational, Symbol, lambdify, list2numpy, matrix2numpy, sin,
-                     symarray, symbols)
+                     symarray, symbols, sympify)
 from diofant.abc import x, y, z
 from diofant.matrices.expressions.matexpr import MatrixElement
 from diofant.utilities.decorator import conserve_mpmath_dps
@@ -256,6 +256,7 @@ def test_lambdify_matrix_vec_input():
     assert numpy.allclose(actual, expected)
 
 
+@pytest.mark.skipif(sys.version_info >= (3, 8), reason="Broken on 3.8")
 def test_lambdify_transl():
     for sym, mat in NUMPY_TRANSLATIONS.items():
         assert sym in dir(diofant)
@@ -312,3 +313,7 @@ def test_from_ndarray():
     assert Matrix(numpy.array([x, y, z])) == Matrix([x, y, z])
     pytest.raises(NotImplementedError, lambda: Matrix(numpy.array([[
         [1, 2], [3, 4]], [[5, 6], [7, 8]]])))
+
+
+def test_sympify():
+    assert sympify(numpy.float128(1.1)) == Float(1.1)
diff --git a/diofant/tests/functions/elementary/test_complexes.py b/diofant/tests/functions/elementary/test_complexes.py
index 87952ade745..1b2c6d0d50f 100644
--- a/diofant/tests/functions/elementary/test_complexes.py
+++ b/diofant/tests/functions/elementary/test_complexes.py
@@ -526,15 +526,6 @@ def test_arg_rewrite():
 
 
 def test_adjoint():
-    a = Symbol('a', antihermitian=True)
-    b = Symbol('b', hermitian=True)
-    assert adjoint(a) == -a
-    assert adjoint(I*a) == I*a
-    assert adjoint(b) == b
-    assert adjoint(I*b) == -I*b
-    assert adjoint(a*b) == -b*a
-    assert adjoint(I*a*b) == I*b*a
-
     x, y = symbols('x y')
     assert adjoint(adjoint(x)) == x
     assert adjoint(x + y) == adjoint(x) + adjoint(y)
diff --git a/diofant/tests/functions/special/test_delta_functions.py b/diofant/tests/functions/special/test_delta_functions.py
index daa885adc85..3658ef775ac 100644
--- a/diofant/tests/functions/special/test_delta_functions.py
+++ b/diofant/tests/functions/special/test_delta_functions.py
@@ -61,7 +61,6 @@ def test_heaviside():
     assert Heaviside(nan) == nan
 
     assert Heaviside(x).is_real
-    assert Heaviside(z).is_real is None
 
     assert adjoint(Heaviside(x)) == Heaviside(x)
     assert adjoint(Heaviside(x - y)) == Heaviside(x - y)
diff --git a/diofant/tests/geometry/test_geometry.py b/diofant/tests/geometry/test_geometry.py
index f3a6bd2a0a0..7c9905f9ca9 100644
--- a/diofant/tests/geometry/test_geometry.py
+++ b/diofant/tests/geometry/test_geometry.py
@@ -850,5 +850,5 @@ def test_idiff():
     assert idiff(circ, y, x, 3).simplify() == ans
     explicit = 12*x/sqrt(-x**2 + 4)**5
     assert ans.subs(solve(circ, y)[0]).equals(explicit)
-    assert True in [sol[y].diff(x, 3).equals(explicit) for sol in solve(circ, y)]
+    assert True in (sol[y].diff(x, 3).equals(explicit) for sol in solve(circ, y))
     assert idiff(x + t + y, [y, t], x) == -Derivative(t, x) - 1
diff --git a/diofant/tests/integrals/test_integrals.py b/diofant/tests/integrals/test_integrals.py
index 274b6a3dc29..d3568d476f2 100644
--- a/diofant/tests/integrals/test_integrals.py
+++ b/diofant/tests/integrals/test_integrals.py
@@ -381,6 +381,14 @@ def test_evalf_integrals():
 
     assert Integral(pi, (x, y, z)).evalf() == Integral(pi, (x, y, z))
     assert Integral(pi, (x, y, y + z)).evalf() == Integral(pi, (x, y, y + z))
+    #
+    # Endpoints causing trouble (rounding error in integration points -> complex log)
+    assert NS(
+        2 + Integral(log(2*cos(x/2)), (x, -pi, pi)), 17, chop=True) == NS(2, 17)
+    assert NS(
+        2 + Integral(log(2*cos(x/2)), (x, -pi, pi)), 20, chop=True) == NS(2, 20)
+    assert NS(
+        2 + Integral(log(2*cos(x/2)), (x, -pi, pi)), 22, chop=True) == NS(2, 22)
 
 
 @pytest.mark.slow
@@ -411,14 +419,6 @@ def test_evalf_integrals_slow():
     # http://mathworld.wolfram.com/VardisIntegral.html
     assert NS(Integral(log(log(sin(x)/cos(x))), (x, pi/4, pi/2)), 15, chop=True) == \
         NS('pi/4*log(4*pi**3/gamma(1/4)**4)', 15)
-    #
-    # Endpoints causing trouble (rounding error in integration points -> complex log)
-    assert NS(
-        2 + Integral(log(2*cos(x/2)), (x, -pi, pi)), 17, chop=True) == NS(2, 17)
-    assert NS(
-        2 + Integral(log(2*cos(x/2)), (x, -pi, pi)), 20, chop=True) == NS(2, 20)
-    assert NS(
-        2 + Integral(log(2*cos(x/2)), (x, -pi, pi)), 22, chop=True) == NS(2, 22)
     # Needs zero handling
     assert NS(pi - 4*Integral(sqrt(1 - x**2), (x, 0, 1)),
               15, maxn=30, chop=True, strict=False) in ('0.0', '0')
diff --git a/diofant/tests/integrals/test_transforms.py b/diofant/tests/integrals/test_transforms.py
index d21ef08034c..3592c75c9d9 100644
--- a/diofant/tests/integrals/test_transforms.py
+++ b/diofant/tests/integrals/test_transforms.py
@@ -143,6 +143,7 @@ def test_mellin_transform():
 
 @pytest.mark.slow
 def test_mellin_transform_bessel():
+    a, b = symbols('a, b', complex=True)
     MT = mellin_transform
 
     # 8.4.19
diff --git a/diofant/tests/interactive/test_printing.py b/diofant/tests/interactive/test_printing.py
new file mode 100644
index 00000000000..1a09797d5bf
--- /dev/null
+++ b/diofant/tests/interactive/test_printing.py
@@ -0,0 +1,34 @@
+import sys
+
+from diofant import init_printing, sqrt
+from diofant.abc import theta, x, y
+
+
+def test_init_printing(capsys):
+    init_printing()
+    sys.displayhook(sqrt(5))
+    assert capsys.readouterr().out == 'sqrt(5)\n'
+    sys.displayhook("xyz")
+    assert capsys.readouterr().out == "'xyz'\n"
+    sys.displayhook(None)
+    assert capsys.readouterr().out == ''
+
+    init_printing(pretty_print=True, no_global=True)
+    sys.displayhook(sqrt(5))
+    assert capsys.readouterr().out == \
+        """\
+  ___\n\
+╲╱ 5 \n"""
+    sys.displayhook(theta)
+    assert capsys.readouterr().out == 'θ\n'
+
+    init_printing(pretty_print=True, use_unicode=False, no_global=True)
+    sys.displayhook(theta)
+    assert capsys.readouterr().out == 'theta\n'
+
+    init_printing(pretty_print=True, order='grevlex', no_global=True)
+    sys.displayhook(y + x + y**2 + x**2)
+    assert capsys.readouterr().out == \
+        """\
+ 2    2        \n\
+x  + y  + x + y\n"""
diff --git a/diofant/tests/polys/test_polytools.py b/diofant/tests/polys/test_polytools.py
index b2fede5d1fc..058a7b8914d 100644
--- a/diofant/tests/polys/test_polytools.py
+++ b/diofant/tests/polys/test_polytools.py
@@ -2402,13 +2402,9 @@ def test_intervals():
     assert intervals(x) == [((0, 0), 1)]
 
     assert intervals(x**128) == [((0, 0), 128)]
-    assert intervals([x**2, x**4]) == [((0, 0), {0: 2, 1: 4})]
-    assert intervals([x**2, x**4], eps=0.1) == [((0, 0), {0: 2, 1: 4})]
 
-    pytest.raises(MultivariatePolynomialError,
-                  lambda: intervals([x**2 - y, x*y + 1]))
+    pytest.raises(MultivariatePolynomialError, lambda: intervals(x**2 - y))
     pytest.raises(MultivariatePolynomialError, lambda: Poly(x**2 - y).intervals())
-    pytest.raises(ValueError, lambda: intervals([x**2, x**4], eps=-0.1))
 
     f = Poly((2*x/5 - Rational(17, 3))*(4*x + Rational(1, 257)))
 
@@ -2448,17 +2444,7 @@ def test_intervals():
          ((-1, -1), 1), ((-1, 0), 3),
          ((1, Rational(3, 2)), 1), ((Rational(3, 2), 2), 7)]
 
-    assert intervals([x**5 - 200, x**5 - 201]) == \
-        [((Rational(75, 26), Rational(101, 35)), {0: 1}), ((Rational(309, 107), Rational(26, 9)), {1: 1})]
-
-    assert intervals([x**5 - 200, x**5 - 201]) == \
-        [((Rational(75, 26), Rational(101, 35)), {0: 1}), ((Rational(309, 107), Rational(26, 9)), {1: 1})]
-
-    assert intervals([x**2 - 200, x**2 - 201]) == \
-        [((-Rational(71, 5), -Rational(85, 6)), {1: 1}), ((-Rational(85, 6), -14), {0: 1}),
-         ((14, Rational(85, 6)), {0: 1}), ((Rational(85, 6), Rational(71, 5)), {1: 1})]
-
-    f, g, h = x**2 - 2, x**4 - 4*x**2 + 4, x - 1
+    f, g = x**2 - 2, x**4 - 4*x**2 + 4
 
     assert intervals(f, inf=Rational(7, 4), sqf=True) == []
     assert intervals(f, inf=Rational(7, 5), sqf=True) == [(Rational(7, 5), Rational(3, 2))]
@@ -2470,17 +2456,6 @@ def test_intervals():
     assert intervals(g, sup=Rational(7, 4)) == [((-2, -1), 2), ((1, Rational(3, 2)), 2)]
     assert intervals(g, sup=Rational(7, 5)) == [((-2, -1), 2)]
 
-    assert intervals([g, h], inf=Rational(7, 4)) == []
-    assert intervals([g, h], inf=Rational(7, 5)) == [((Rational(7, 5), Rational(3, 2)), {0: 2})]
-    assert intervals([g, h], sup=Rational(7, 4)) == [((-2, -1), {0: 2}), ((1, 1), {1: 1}), ((1, Rational(3, 2)), {0: 2})]
-    assert intervals(
-        [g, h], sup=Rational(7, 5)) == [((-2, -1), {0: 2}), ((1, 1), {1: 1})]
-
-    assert intervals([x + 2, x**2 - 2]) == \
-        [((-2, -2), {0: 1}), ((-2, -1), {1: 1}), ((1, 2), {1: 1})]
-    assert intervals([x + 2, x**2 - 2], strict=True) == \
-        [((-2, -2), {0: 1}), ((-Rational(3, 2), -1), {1: 1}), ((1, 2), {1: 1})]
-
     f = 7*z**4 - 19*z**3 + 20*z**2 + 17*z + 20
 
     assert intervals(f) == []
@@ -2502,15 +2477,6 @@ def test_intervals():
 
     pytest.raises(ValueError, lambda: intervals(x**2 - 2, eps=10**-100000))
     pytest.raises(ValueError, lambda: Poly(x**2 - 2).intervals(eps=10**-100000))
-    pytest.raises(
-        ValueError, lambda: intervals([x**2 - 2, x**2 - 3], eps=10**-100000))
-
-
-@pytest.mark.xfail
-def test_intervals_xfail():
-    assert intervals([x + 1, x + 2, x - 1, x + 1, 1, x - 1, x - 1, (x - 2)**2]) == \
-        [((-2, -2), {1: 1}), ((-1, -1), {0: 1, 3: 1}),
-         ((1, 1), {2: 1, 5: 1, 6: 1}), ((2, 2), {7: 2})]
 
 
 def test_refine_root():
diff --git a/diofant/tests/polys/test_rootisolation.py b/diofant/tests/polys/test_rootisolation.py
index 80d9872e1d5..516e1e88b4e 100644
--- a/diofant/tests/polys/test_rootisolation.py
+++ b/diofant/tests/polys/test_rootisolation.py
@@ -429,90 +429,99 @@ def test_dup_isolate_real_roots():
     pytest.raises(DomainError, lambda: R.dup_isolate_real_roots(x + 3))
     pytest.raises(DomainError, lambda: R.dup_isolate_real_roots((x + 2)*(x + 3)**2))
 
+    R, x = ring("x", QQ.algebraic_field(I))
+
+    f = (x**2 - I)**2*(x - 2*I)**3
+
+    assert R.dup_isolate_real_roots(f) == []  # issue diofant/diofant#789
+    assert R.dup_isolate_real_roots(f*(x - 1)**3) == [((1, 1), 3)]
+
+    f = x**4*(x - 1)**3*(x**2 - 2)**2
+
+    assert R.dup_isolate_real_roots(f) == \
+        [((-2, -1), 2), ((0, 0), 4), ((1, 1), 3), ((QQ(4, 3), QQ(3, 2)), 2)]
+
 
-def test_dup_isolate_real_roots_list():
+def test_dup_isolate_real_roots_pair():
     R, x = ring("x", ZZ)
 
-    assert R.dup_isolate_real_roots_list([x*(x + 1), x]) == \
+    assert R.dup_isolate_real_roots_pair(x*(x + 1), x) == \
         [((-1, -1), {0: 1}), ((0, 0), {0: 1, 1: 1})]
-    assert R.dup_isolate_real_roots_list([x*(x - 1), x]) == \
+    assert R.dup_isolate_real_roots_pair(x*(x - 1), x) == \
         [((0, 0), {0: 1, 1: 1}), ((1, 1), {0: 1})]
 
     f, g = (x**2 - 2)**2, x - 1
 
-    assert R.dup_isolate_real_roots_list([f, g], inf=QQ(7, 4)) == []
-    assert R.dup_isolate_real_roots_list([f, g], inf=QQ(7, 5)) == \
+    assert R.dup_isolate_real_roots_pair(f, g, inf=QQ(7, 4)) == []
+    assert R.dup_isolate_real_roots_pair(f, g, inf=QQ(7, 5)) == \
         [((QQ(7, 5), QQ(3, 2)), {0: 2})]
-    assert R.dup_isolate_real_roots_list([f, g], sup=QQ(7, 5)) == \
+    assert R.dup_isolate_real_roots_pair(f, g, sup=QQ(7, 5)) == \
         [((-2, -1), {0: 2}), ((1, 1), {1: 1})]
-    assert R.dup_isolate_real_roots_list([f, g], sup=QQ(7, 4)) == \
+    assert R.dup_isolate_real_roots_pair(f, g, sup=QQ(7, 4)) == \
         [((-2, -1), {0: 2}), ((1, 1), {1: 1}), ((1, QQ(3, 2)), {0: 2})]
-    assert R.dup_isolate_real_roots_list([f, g], sup=-QQ(7, 4)) == []
-    assert R.dup_isolate_real_roots_list([f, g], sup=-QQ(7, 5)) == \
+    assert R.dup_isolate_real_roots_pair(f, g, sup=-QQ(7, 4)) == []
+    assert R.dup_isolate_real_roots_pair(f, g, sup=-QQ(7, 5)) == \
         [((-QQ(3, 2), -QQ(7, 5)), {0: 2})]
-    assert R.dup_isolate_real_roots_list([f, g], inf=-QQ(7, 5)) == \
+    assert R.dup_isolate_real_roots_pair(f, g, inf=-QQ(7, 5)) == \
         [((1, 1), {1: 1}), ((1, 2), {0: 2})]
-    assert R.dup_isolate_real_roots_list([f, g], inf=-QQ(7, 4)) == \
+    assert R.dup_isolate_real_roots_pair(f, g, inf=-QQ(7, 4)) == \
         [((-QQ(3, 2), -1), {0: 2}), ((1, 1), {1: 1}), ((1, 2), {0: 2})]
 
     f, g = 2*x**2 - 1, x**2 - 2
 
-    assert R.dup_isolate_real_roots_list([f, g]) == \
+    assert R.dup_isolate_real_roots_pair(f, g) == \
         [((-2, -1), {1: 1}), ((-1, 0), {0: 1}),
          ((0, 1), {0: 1}), ((1, 2), {1: 1})]
-    assert R.dup_isolate_real_roots_list([f, g], strict=True) == \
+    assert R.dup_isolate_real_roots_pair(f, g, strict=True) == \
         [((-QQ(3, 2), -QQ(4, 3)), {1: 1}), ((-1, -QQ(2, 3)), {0: 1}),
          ((QQ(2, 3), 1), {0: 1}), ((QQ(4, 3), QQ(3, 2)), {1: 1})]
 
     f, g = x**2 - 2, (x - 1)*(x**2 - 2)
 
-    assert R.dup_isolate_real_roots_list([f, g]) == \
+    assert R.dup_isolate_real_roots_pair(f, g) == \
         [((-2, -1), {1: 1, 0: 1}), ((1, 1), {1: 1}), ((1, 2), {1: 1, 0: 1})]
 
     f, g = x*(x**2 - 2), x**2*(x - 1)*(x**2 - 2)
 
-    assert R.dup_isolate_real_roots_list([f, g]) == \
+    assert R.dup_isolate_real_roots_pair(f, g) == \
         [((-2, -1), {1: 1, 0: 1}), ((0, 0), {0: 1, 1: 2}),
          ((1, 1), {1: 1}), ((1, 2), {1: 1, 0: 1})]
 
     f, g = x**2*(x - 1)**3*(x**2 - 2)**2, x*(x - 1)**2*(x**2 + 2)
 
-    assert R.dup_isolate_real_roots_list([f, g]) == \
+    assert R.dup_isolate_real_roots_pair(f, g) == \
         [((-2, -1), {0: 2}), ((0, 0), {0: 2, 1: 1}),
          ((1, 1), {0: 3, 1: 2}), ((1, 2), {0: 2})]
-    assert R.dup_isolate_real_roots_list([f, g], basis=True) == \
+    assert R.dup_isolate_real_roots_pair(f, g, basis=True) == \
         [((-2, -1), {0: 2}, [1, 0, -2]), ((0, 0), {0: 2, 1: 1}, [1, 0]),
          ((1, 1), {0: 3, 1: 2}, [1, -1]), ((1, 2), {0: 2}, [1, 0, -2])]
 
-    R, x = ring("x", EX)
-
-    pytest.raises(DomainError, lambda: R.dup_isolate_real_roots_list([x + 3]))
+    f, g = x, R.zero
 
+    assert R.dup_isolate_real_roots_pair(f, g) == \
+        R.dup_isolate_real_roots_pair(g, f) == [((0, 0), {0: 1, 1: 1})]
 
-@pytest.mark.xfail
-def test_dup_isolate_real_roots_list_xfail():
-    R, x = ring("x", ZZ)
+    f *= x**2
 
-    assert R.dup_isolate_real_roots_list([x + 1, x + 2, x - 1, x + 1, x - 1, x - 1]) == \
-        [((-2, -2), {1: 1}), ((-1, -1), {0: 1, 3: 1}), ((1, 1), {2: 1, 4: 1, 5: 1})]
+    assert R.dup_isolate_real_roots_pair(f, g) == \
+        R.dup_isolate_real_roots_pair(g, f) == [((0, 0), {0: 3, 1: 3})]
 
-    assert R.dup_isolate_real_roots_list([x + 1, x + 2, x - 1, x + 1, x - 1, x + 2]) == \
-        [((-2, -2), {1: 1, 5: 1}), ((-1, -1), {0: 1, 3: 1}), ((1, 1), {2: 1, 4: 1})]
+    R, x = ring("x", EX)
 
+    pytest.raises(DomainError, lambda: R.dup_isolate_real_roots_pair(x, x + 3))
 
-def test_dup_isolate_real_roots_list_QQ():
     R, x = ring("x", ZZ)
 
     f, g = x**5 - 200, x**5 - 201
 
-    assert R.dup_isolate_real_roots_list([f, g]) == \
+    assert R.dup_isolate_real_roots_pair(f, g) == \
         [((QQ(75, 26), QQ(101, 35)), {0: 1}), ((QQ(309, 107), QQ(26, 9)), {1: 1})]
 
     R, x = ring("x", QQ)
 
     f, g = -x**5/200 + 1, -x**5/201 + 1
 
-    assert R.dup_isolate_real_roots_list([f, g]) == \
+    assert R.dup_isolate_real_roots_pair(f, g) == \
         [((QQ(75, 26), QQ(101, 35)), {0: 1}), ((QQ(309, 107), QQ(26, 9)), {1: 1})]
 
 
diff --git a/diofant/tests/printing/test_mathematica.py b/diofant/tests/printing/test_mathematica.py
index 943382361dc..9e51d5dfc5b 100644
--- a/diofant/tests/printing/test_mathematica.py
+++ b/diofant/tests/printing/test_mathematica.py
@@ -7,12 +7,12 @@
 from diofant.core import (Catalan, Derivative, Dummy, E, Eq, EulerGamma,
                           Function, Gt, Integer, Lambda, Le, Ne, Rational,
                           Tuple, oo, pi, symbols)
-from diofant.functions import (Heaviside, Max, Min, Piecewise, acos, asin,
-                               asinh, atan, atanh, binomial, conjugate, cos,
-                               cosh, cot, coth, csch, erfc, exp, factorial,
-                               factorial2, fibonacci, gamma, hyper, im, log,
-                               meijerg, polygamma, polylog, re, rf, sech, sign,
-                               sin, sinh, tan, tanh, zeta)
+from diofant.functions import (Heaviside, Max, Min, Piecewise, acos, acosh,
+                               acot, acoth, asin, asinh, atan, atanh, binomial,
+                               conjugate, cos, cosh, cot, coth, csch, erfc,
+                               exp, factorial, factorial2, fibonacci, gamma,
+                               hyper, im, log, meijerg, polygamma, polylog, re,
+                               rf, sech, sign, sin, sinh, tan, tanh, zeta)
 from diofant.integrals import Integral
 from diofant.logic import Or, false, true
 from diofant.matrices import Matrix, SparseMatrix
@@ -69,10 +69,15 @@ def test_Function():
     assert mcode(asin(x)) == "ArcSin[x]"
     assert mcode(acos(x)) == "ArcCos[x]"
     assert mcode(atan(x)) == "ArcTan[x]"
+    assert mcode(acot(x)) == "ArcCot[x]"
     assert mcode(sinh(x)) == "Sinh[x]"
     assert mcode(cosh(x)) == "Cosh[x]"
     assert mcode(tanh(x)) == "Tanh[x]"
     assert mcode(coth(x)) == "Coth[x]"
+    assert mcode(asinh(x)) == "ArcSinh[x]"
+    assert mcode(acosh(x)) == "ArcCosh[x]"
+    assert mcode(atanh(x)) == "ArcTanh[x]"
+    assert mcode(acoth(x)) == "ArcCoth[x]"
     assert mcode(sech(x)) == "Sech[x]"
     assert mcode(csch(x)) == "Csch[x]"
     assert mcode(erfc(x)) == "Erfc[x]"
@@ -85,11 +90,9 @@ def test_Function():
     assert mcode(rf(x, y)) == "Pochhammer[x, y]"
     assert mcode(gamma(x)) == "Gamma[x]"
     assert mcode(zeta(x)) == "Zeta[x]"
-    assert mcode(asinh(x)) == "ArcSinh[x]"
     assert mcode(Heaviside(x)) == "UnitStep[x]"
     assert mcode(fibonacci(x)) == "Fibonacci[x]"
     assert mcode(polylog(x, y)) == "PolyLog[x, y]"
-    assert mcode(atanh(x)) == "ArcTanh[x]"
 
     class myfunc1(Function):
         @classmethod
diff --git a/diofant/tests/printing/test_pretty.py b/diofant/tests/printing/test_pretty.py
index be4b10180a3..5ff58b2344f 100644
--- a/diofant/tests/printing/test_pretty.py
+++ b/diofant/tests/printing/test_pretty.py
@@ -5,14 +5,15 @@
 import pytest
 
 from diofant import (FF, QQ, RR, ZZ, Add, And, Basic, Complement, Contains,
-                     Derivative, Dict, Eq, Equivalent, EulerGamma, FiniteSet,
-                     Float, Function, Ge, GoldenRatio, Gt, I, Implies, Integer,
-                     Integral, Intersection, Interval, Lambda, Le, Limit, Lt,
-                     Matrix, MatrixSymbol, Mod, Mul, Nand, Ne, Nor, Not, O, Or,
-                     Pow, Product, Range, Rational, Ray, RealField, RootOf,
-                     RootSum, S, Segment, Subs, Sum, Symbol,
-                     SymmetricDifference, Trace, Tuple, Union, Xor, cbrt,
-                     conjugate, grlex, groebner, ilex, oo, pi, root, symbols)
+                     Derivative, Dict, E, Eq, Equivalent, EulerGamma,
+                     FiniteSet, Float, Function, Ge, GoldenRatio, Gt, I,
+                     Implies, Integer, Integral, Intersection, Interval,
+                     Lambda, Le, Limit, Lt, Matrix, MatrixSymbol, Mod, Mul,
+                     Nand, Ne, Nor, Not, O, Or, Pow, Product, Range, Rational,
+                     Ray, RealField, RootOf, RootSum, S, Segment, Subs, Sum,
+                     Symbol, SymmetricDifference, Trace, Tuple, Union, Xor,
+                     cbrt, conjugate, grlex, groebner, ilex, oo, pi, root,
+                     symbols)
 from diofant.abc import a, b, c, d, e, f, k, l, lamda, m, n, t, w, x, y, z
 from diofant.core.trace import Tr
 from diofant.diffgeom import BaseVectorField
@@ -2535,6 +2536,21 @@ def test_pretty_integrals():
     assert pretty(expr) == ascii_str
     assert upretty(expr) == ucode_str
 
+    expr = Integral((x**4 + x**2*exp(x) - x**2 - 2*x*exp(x) - 2*x -
+                     exp(x))*exp(x)/((x - 1)**2*(x + 1)**2*(exp(x) + 1)))
+
+    ucode_str = \
+        """\
+⌠                                            \n\
+⎮  x ⎛ x  2      x      x    4    2      ⎞   \n\
+⎮ ℯ ⋅⎝ℯ ⋅x  - 2⋅ℯ ⋅x - ℯ  + x  - x  - 2⋅x⎠   \n\
+⎮ ──────────────────────────────────────── dx\n\
+⎮        ⎛ x    ⎞        2        2          \n\
+⎮        ⎝ℯ  + 1⎠⋅(x - 1) ⋅(x + 1)           \n\
+⌡                                            \
+"""
+    assert upretty(expr) == ucode_str
+
 
 def test_pretty_matrix():
     # Empty Matrix
@@ -5305,3 +5321,70 @@ def test_pretty_Mod():
     assert upretty(Mod(x, 7) + 1) == ucode_str4
     assert pretty(2 * Mod(x, 7)) == ascii_str5
     assert upretty(2 * Mod(x, 7)) == ucode_str5
+
+
+def test_Pow_roots():
+    expr = root(pi, E)
+
+    assert pretty(expr) == \
+        """\
+E ____\n\
+\\/ pi \
+"""
+    assert upretty(expr) == \
+        """\
+ℯ ___\n\
+╲╱ π \
+"""
+
+    expr = root(pi, pi)
+
+    assert pretty(expr) == \
+        """\
+pi____\n\
+\\/ pi \
+"""
+    assert upretty(expr) == \
+        """\
+π ___\n\
+╲╱ π \
+"""
+
+    expr = root(7, pi)  # issue diofant/diofant#888
+
+    assert pretty(expr) == \
+        """\
+pi___\n\
+\\/ 7 \
+"""
+    assert upretty(expr) == \
+        """\
+π ___\n\
+╲╱ 7 \
+"""
+
+    expr = root(pi, EulerGamma)
+
+    assert upretty(expr) == \
+        """\
+γ ___\n\
+╲╱ π \
+"""
+    assert pretty(expr) == \
+        """\
+EulerGamma____\n\
+        \\/ pi \
+"""
+
+    expr = root(7, Symbol("x_17"))
+
+    assert upretty(expr) == \
+        """\
+x₁₇___\n\
+ ╲╱ 7 \
+"""
+    assert pretty(expr) == \
+        """\
+x_17___\n\
+  \\/ 7 \
+"""
diff --git a/diofant/tests/printing/test_str.py b/diofant/tests/printing/test_str.py
index e65bb1a9e0b..bac9d59ef5f 100644
--- a/diofant/tests/printing/test_str.py
+++ b/diofant/tests/printing/test_str.py
@@ -722,6 +722,10 @@ def test_MatrixSlice():
     assert str(MatrixSymbol('X', 10, 10)[5, :5:2]) == 'X[5, :5:2]'
 
 
+def test_MatrixInverse():
+    assert str(MatrixSymbol('X', 3, 3).inverse()) == 'X^-1'
+
+
 def test_true_false():
     assert str(true) == repr(true) == sstr(true) == "true"
     assert str(false) == repr(false) == sstr(false) == "false"
diff --git a/diofant/tests/printing/test_stringpict.py b/diofant/tests/printing/test_stringpict.py
index 76b455e64a9..c0a623a216e 100644
--- a/diofant/tests/printing/test_stringpict.py
+++ b/diofant/tests/printing/test_stringpict.py
@@ -36,3 +36,8 @@ def test_stringpict():
     assert p1 != 0
 
     assert len(p1) == 1
+
+
+def test_dumb_term(capsys):
+    print(pretty(sin(x)))
+    assert capsys.readouterr().out == 'sin(x)\n'
diff --git a/diofant/tests/series/test_gruntz.py b/diofant/tests/series/test_gruntz.py
index 1bd3c55e26c..59c34a03b68 100644
--- a/diofant/tests/series/test_gruntz.py
+++ b/diofant/tests/series/test_gruntz.py
@@ -383,6 +383,7 @@ def test_intractable():
     assert gruntz(gamma(Rational(1, 7) + 1/x), x) == gamma(Rational(1, 7))
     assert gruntz(log(x**x)/log(gamma(x)), x) == 1
     assert gruntz(log(gamma(gamma(x)))/exp(x), x) == oo
+    assert gruntz(acosh(1 + 1/x)*sqrt(x), x) == sqrt(2)
 
     # issue sympy/sympy#10804
     assert gruntz(2*airyai(x)*root(x, 4) *
diff --git a/diofant/tests/series/test_limits.py b/diofant/tests/series/test_limits.py
index b3439428c0e..198d63da098 100644
--- a/diofant/tests/series/test_limits.py
+++ b/diofant/tests/series/test_limits.py
@@ -3,8 +3,8 @@
 import pytest
 
 from diofant import (E, Float, Function, I, Integral, Limit, Matrix, Piecewise,
-                     PoleError, Rational, Sum, Symbol, acos, atan, cbrt,
-                     ceiling, cos, cot, diff, digamma, erf, erfi, exp,
+                     PoleError, Rational, Sum, Symbol, acos, atan, binomial,
+                     cbrt, ceiling, cos, cot, diff, digamma, erf, erfi, exp,
                      factorial, floor, gamma, integrate, limit, log, nan, oo,
                      pi, polygamma, root, sign, simplify, sin, sinh, sqrt,
                      subfactorial, symbols, tan)
@@ -655,3 +655,46 @@ def test_sympyissue_16222():
 def test_sympyissue_16714():
     e = ((n**(n + 1) + (n + 1)**n)/n**(n + 1))**n
     assert limit(e, n, oo) == E**E
+
+
+@pytest.mark.timeout(20)
+def test_sympyissue_15282():
+    assert limit((x**2000 - (x + 1)**2000)/x**1999, x, oo) == -2000
+
+
+def test_sympyissue_16722():
+    n, z = symbols('n z')
+    assert isinstance(limit(binomial(n + z, n)*n**-z, n, oo), Limit)
+
+    z = symbols('z', positive=True)
+    assert limit(binomial(n + z, n)*n**-z, n, oo) == 1/gamma(z + 1)
+
+    n = symbols('n', positive=True, integer=True)
+    z = symbols('z', positive=True)
+    assert isinstance(limit(binomial(n + z, n)*n**-z, n, oo), Limit)
+
+    n, z = symbols('n z', positive=True, integer=True)
+    assert isinstance(limit(binomial(n + z, n)*n**-z, n, oo), Limit)
+
+
+def test_sympyissue_15673():
+    p = symbols('p')
+    alpha = symbols('α', positive=True)
+
+    e = Limit(4*pi*p**(-alpha)*(p**3 - p**alpha)/(alpha - 3), p, 0)
+    assert isinstance(e.doit(), Limit)  # but see diofant/diofant#425
+
+
+def test_sympyissue_17380():
+    assert limit(x*(((x + 1)**2 + 1)/(x**2 + 1) - 1), x, oo) == 2
+
+
+def test_sympyissue_17431():
+    assert limit(((n + 1) + 1)/(((n + 1) + 2)*factorial(n + 1)) *
+                 (n + 2)*factorial(n)/(n + 1), n, oo) == 0
+    assert limit((n + 2)**2*factorial(n)/((n + 1)*(n + 3)*factorial(n + 1)),
+                 n, oo) == 0
+
+    # test from sympy/sympy#17434 (see also diofant/diofant#425):
+    y = symbols('y', integer=True, positive=True)
+    assert isinstance(limit(x*factorial(x)/factorial(x + y), x, oo), Limit)
diff --git a/diofant/tests/simplify/test_simplify.py b/diofant/tests/simplify/test_simplify.py
index 73e99267587..504fb1d356d 100644
--- a/diofant/tests/simplify/test_simplify.py
+++ b/diofant/tests/simplify/test_simplify.py
@@ -352,6 +352,8 @@ def test_nsimplify():
         assert nsimplify(i) == ans
         assert nsimplify(i + x) == x + ans
 
+    assert nsimplify(Sum(1/n**2, (n, 1, oo)), [pi]) == pi**2/6
+
 
 def test_sympyissue_9448():
     expr = (1/(1 - (-1)**Rational(2, 3) - cbrt(-1)) +
diff --git a/diofant/tests/solvers/test_diophantine.py b/diofant/tests/solvers/test_diophantine.py
index 911bcdc4057..89d4245d04a 100644
--- a/diofant/tests/solvers/test_diophantine.py
+++ b/diofant/tests/solvers/test_diophantine.py
@@ -742,6 +742,7 @@ def test_diopcoverage():
     assert cornacchia(1, 1, 5) == {(2, 1)}
     assert cornacchia(1, 2, 17) == {(3, 2)}
     assert cornacchia(2, 3, 31) == set()
+    assert cornacchia(1, 4, 52) == {(4, 3)}
 
     pytest.raises(ValueError, lambda: reconstruct(4, 20, 1))
 
diff --git a/diofant/tests/vector/test_printing.py b/diofant/tests/vector/test_printing.py
index 4c61348c2c4..dd7613af8bc 100644
--- a/diofant/tests/vector/test_printing.py
+++ b/diofant/tests/vector/test_printing.py
@@ -109,6 +109,7 @@ def test_str_printing():
     assert str(d[5]) == 'a*(N.i|N.k) + (-b)*(N.j|N.k)'
     assert str(d[8]) == ('(N.j|N.k) + (C.x**2 - ' +
                          'Integral(f(b), b))*(N.k|N.k)')
+    assert str(N.origin) == 'N.origin'
 
 
 def test_pretty_print_unicode():
diff --git a/diofant/tests/vector/test_vector.py b/diofant/tests/vector/test_vector.py
index 393ac421e86..64cbebb17aa 100644
--- a/diofant/tests/vector/test_vector.py
+++ b/diofant/tests/vector/test_vector.py
@@ -100,6 +100,10 @@ def test_vector():
     assert VectorAdd(v1, Vector.zero) == v1
     assert VectorMul(0, Vector.zero) == Vector.zero
 
+    v = (a*b + a*c + b**2 + b*c)*i + j
+
+    assert v.factor() == (a + b)*(b + c)*i + j
+
 
 def test_vector_magnitude_normalize():
     assert Vector.zero.magnitude() == 0
diff --git a/docs/aboutus.rst b/docs/aboutus.rst
index 1d1acdc2e23..d595e2c1616 100644
--- a/docs/aboutus.rst
+++ b/docs/aboutus.rst
@@ -2,8 +2,8 @@ About
 =====
 
 The `Diofant`_ project is a fork of the `SymPy`_, started by Sergey B
-Kirpichev, last regular SymPy's commit is :commit:`cbdd072`, (22 Feb
-2015).  The git history goes back to 2007, when development was in svn.
+Kirpichev, last regular `SymPy`_'s commit is :commit:`cbdd072`, (22
+Feb 2015).  The git history goes back to 2007.
 
 The project was named after the Diophantus of Alexandria.  His
 "Arithmetica" is one of the earliest known texts that use symbols in
@@ -12,55 +12,36 @@ equations.  "Diofant" is a transliteration of Диофант, from Russian.
 License
 -------
 
-Unless stated otherwise, all files in the Diofant project are licensed
-using the new BSD license:
+Unless stated otherwise, all files in the `Diofant`_ project are
+licensed using the new BSD license:
 
 .. include:: ../LICENSE.rst
 
-SymPy Development Team
-----------------------
+SymPy
+-----
 
-SymPy is a team project and it was developed by a lot of people.
+`SymPy`_ is a team project and it was developed by a lot of people.
 Please refer to the `AUTHORS
 <https://github.com/sympy/sympy/blob/master/AUTHORS>`_ file in the
 `SymPy repository <https://github.com/sympy/sympy/>`_ for up-to-date
 list of contributors.
 
-Please note that that list is incomplete intentionally, because some
-former SymPy developers don't want to be mentioned in the context of
-the SymPy or derived projects.
-
-Brief History
--------------
-
-`SymPy`_ was started by Ondřej Čertík in 2005, he wrote some code during the
-summer, then he wrote some more code during the summer 2006. In February 2007,
-Fabian Pedregosa joined the project and helped fixed many things, contributed
-documentation and made it alive again. 5 students (Mateusz Paprocki, Brian
-Jorgensen, Jason Gedge, Robert Schwarz and Chris Wu) improved SymPy incredibly
-during the summer 2007 as part of the Google Summer of Code (GSoC). Pearu Peterson
-joined the development during the summer 2007 and he has made SymPy much more
-competitive and fast (from 10x to 100x) by rewriting the core from scratch.
-Jurjen N.E. Bos has contributed pretty printing and other patches.
-Fredrik Johansson has wrote mpmath and contributed a lot of patches.
-
-SymPy has participated in every GSoC since 2007.  Moderate amount
-of SymPy's development has come from GSoC students.
-
-In 2011, Ondřej Čertík stepped down as lead developer, with Aaron Meurer, who
-also started as a GSoC student, taking his place.
-
-Ondřej Čertík is still active in the community, but is too busy with work
-and family to play a lead development role.  Unfortunately, his remaining
-activity neither constructive nor productive anymore and SymPy just
-slowly dies now: most former contributors are inactive now or
-explicitly leaving this "friendly and welcoming" project.
-
-This unfortunate situation was major reason to fork the
-SymPy as the `Diofant`_ project.  Development in the new project will
-be open and public, without hidden double standards, centered about
-good, proved code and not project popularity counters.  Here we believe
-that mathematical correctness is more important than political one.
+Please refer to the `SymPy`_'s `README.rst
+<https://github.com/sympy/sympy/blob/master/README.rst>`_ for a brief
+history of the `SymPy`_ project.  It was started by Ondřej Čertík in
+2005, he is still active in the community, but is too busy with work
+and family to play a lead development role.  Unfortunately, his
+remaining activity neither constructive nor productive anymore and
+`SymPy`_ just slowly dies now: most former contributors are inactive
+or explicitly leaving this "friendly and welcoming" project, some
+developers even don't want to be mentioned in the context of this
+project and its derivatives.  Current maintainer continue bad practice
+of blocking unpleasant opinions and even delete bugreports.
+
+This unfortunate situation was major reason to fork the `SymPy`_ as
+the `Diofant`_ project.  Development in the new project will be open
+and public, without hidden double standards, centered about good,
+proved code and not the project popularity counters.
 
 .. _SymPy : https://www.sympy.org/
 .. _Diofant : https://diofant.readthedocs.io/en/latest/
diff --git a/docs/conf.py b/docs/conf.py
index 387df934c13..a2ef39aa980 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -9,6 +9,9 @@
 # removed automatically).
 #
 
+import inspect
+import os
+import sys
 import warnings
 
 import diofant
@@ -18,7 +21,7 @@
 warnings.simplefilter('error', UserWarning)
 
 # Add any Sphinx extension module names here, as strings.
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.viewcode', 'sphinx.ext.mathjax',
+extensions = ['sphinx.ext.autodoc', 'sphinx.ext.linkcode', 'sphinx.ext.mathjax',
               'sphinx.ext.graphviz', 'sphinx.ext.intersphinx',
               'sphinx.ext.extlinks', 'sphinx.ext.napoleon',
               'sphinxcontrib.bibtex']
@@ -57,7 +60,7 @@
 latex_use_xindy = False
 
 # This value determines how to group the document tree into LaTeX source
-# files. It must be a list of tuples (startdocname, targetname, title,
+# files.  It must be a list of tuples (startdocname, targetname, title,
 # author, documentclass, toctree_only),
 latex_documents = [('index', 'diofant.tex', 'Diofant Documentation',
                     'Diofant Development Team', 'manual', True)]
@@ -111,11 +114,23 @@
 # set, is inherited form parents.
 autodoc_inherit_docstrings = False
 
-# A list of paths that contain custom static files. Relative paths are taken as
-# relative to the configuration directory. They are copied to the output’s
+# A list of paths that contain custom static files.  Relative paths are taken as
+# relative to the configuration directory.  They are copied to the output’s
 # _static directory.
 html_static_path = ['_static']
 
+# These paths are either relative to html_static_path
+# or fully qualified paths (eg. https://...)
+html_css_files = [
+    'custom.css',
+]
+
+# A list of paths that contain extra files not directly related to the
+# documentation, such as robots.txt or .htaccess.  Relative paths are taken
+# as relative to the configuration directory.  They are copied to the
+# output directory. They will overwrite any existing file of the same name.
+html_extra_path = ['robots.txt']
+
 # Should we show "Created using Sphinx" in the HTML footer?
 html_show_sphinx = False
 
@@ -128,8 +143,11 @@
 html_theme_options = {
     'logo_only': True,
     'display_version': False,
+    'analytics_id': 'UA-147167689-1',
 }
 
+# The inline configuration options for mathjax.  The value is used as
+# a parameter of MathJax.Hub.Config().
 mathjax_config = {
     'CommonHTML': {'linebreaks': {'automatic': True}},
     'HTML-CSS': {'linebreaks': {'automatic': True}},
@@ -137,6 +155,55 @@
 }
 
 
-# https://docs.readthedocs.io/en/latest/guides/adding-custom-css.html
-def setup(app):
-    app.add_stylesheet('custom.css')
+def linkcode_resolve(domain, info):
+    """Determine the URL corresponding to Python object. """
+    if domain != 'py':
+        return
+
+    modname = info['module']
+    fullname = info['fullname']
+
+    submod = sys.modules.get(modname)
+    if submod is None:
+        return
+
+    obj = submod
+    for part in fullname.split('.'):
+        try:
+            obj = getattr(obj, part)
+        except Exception:
+            return
+
+    # strip decorators, which would resolve to the source of the decorator
+    # possibly an upstream bug in getsourcefile, bpo-1764286
+    try:
+        unwrap = inspect.unwrap
+    except AttributeError:
+        pass
+    else:
+        obj = unwrap(obj)
+
+    try:
+        fn = inspect.getsourcefile(obj)
+    except Exception:
+        fn = None
+    if not fn:
+        return
+
+    try:
+        source, lineno = inspect.getsourcelines(obj)
+    except Exception:
+        lineno = None
+
+    if lineno:
+        linespec = "#L%d-L%d" % (lineno, lineno + len(source) - 1)
+    else:
+        linespec = ""
+
+    fn = os.path.relpath(fn, start=os.path.dirname(diofant.__file__))
+
+    blobpath = "https://github.com/diofant/diofant/blob/"
+    if 'dev' in version:
+        return blobpath + "master/diofant/%s%s" % (fn, linespec)
+    else:
+        return blobpath + "v%s/diofant/%s%s" % (version, fn, linespec)
diff --git a/docs/guide.rst b/docs/guide.rst
index da398f82378..9a081ffb163 100644
--- a/docs/guide.rst
+++ b/docs/guide.rst
@@ -37,6 +37,8 @@ Contributing Code
 All work should be submitted via `Pull Requests (PR)`_.
 
 1. PR can be submitted as soon as there is code worth discussing.
+   Please make a draft PR, if one is not intended to be merged
+   in its present shape even if all checks pass.
 
 2. Please put your work on the branch of your fork, not in the
    master branch.  PR should generally be made against master.
@@ -52,7 +54,10 @@ All work should be submitted via `Pull Requests (PR)`_.
 
 5. PR should include tests:
 
-   1. Bugfixes should include regression tests.
+   1. Bugfixes should include regression tests.  Please format
+      them accordingly, to include references for fixed
+      issues (e.g. by naming test like ``test_diofantissue_123``
+      or adding comment with issue number).
    2. All new functionality should be tested, every new line
       should be covered by tests.
    3. Optionally, provide doctests to illustrate usage.  But keep in
@@ -135,4 +140,4 @@ and publish this release tag::
 .. _PEP 257: https://www.python.org/dev/peps/pep-0257/
 .. _flake8: http://flake8.rtfd.io/
 .. _No Code of Conduct: https://github.com/domgetter/NCoC
-.. _mention closed issues: https://help.github.com/articles/closing-issues-via-commit-messages
+.. _mention closed issues: https://help.github.com/en/articles/closing-issues-using-keywords
diff --git a/docs/index.rst b/docs/index.rst
index b5758a381c0..039a1d31418 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,7 +1,10 @@
+.. meta::
+    :google-site-verification: MJQX0mr79AK_SDl4JzriJBKrebj_maVDMKFTpsPhGiU
+
 Diofant's documentation
 =======================
 
-.. image:: https://mybinder.org/badge_logo.svg
+.. image:: https://static.mybinder.org/badge_logo.svg
     :target: https://mybinder.org/v2/gh/diofant/diofant/master
 
 `Diofant <https://diofant.readthedocs.io/en/latest/>`_ is a Python
diff --git a/docs/internals/hyperexpand.rst b/docs/internals/hyperexpand.rst
index 52e7a86a326..b920a65e23b 100644
--- a/docs/internals/hyperexpand.rst
+++ b/docs/internals/hyperexpand.rst
@@ -325,6 +325,7 @@ and
 Hence
 
 .. math ::
+    \begin{aligned}
     res_{s = b + l_u + t} =& -z^{b + l_u}
        \frac{(-1)^{l_u}}{l_u!}
        \prod_{i=1}^{u} \frac{(-1)^{\delta_i}}{(l_u - k_i + 1)_{\delta_i}}
@@ -340,6 +341,7 @@ Hence
              \prod_{j=n+1}^p (-1)^t (l_u + b + 1 - a_j)_t^*}
             {\prod_{j=1}^m (-1)^t (l_u + b + 1 - b_j)_t^*
              \prod_{j=m+1}^q (1 - b_j + l_u + b)_t},
+    \end{aligned}
 
 where the `*` means to omit the terms we treated specially.
 
@@ -426,11 +428,13 @@ We get these functions mainly by guessing and testing the result. Hence
 we proceed by computing `f(g(w))` (and simplifying naively)
 
 .. math ::
+   \begin{aligned}
    f(g(w)) &= -\frac{\pi^2 g(w)^4}{16} \\
            &= -\frac{\pi^2 g\left(\frac{2}{\sqrt{\pi}} \exp\left(\frac{i \pi}{4}\right) w^{\frac{1}{4}}\right)^4}{16} \\
            &= -\frac{\pi^2 \frac{2^4}{\sqrt{\pi}^4} \exp\left(\frac{i \pi}{4}\right)^4 {w^{\frac{1}{4}}}^4}{16} \\
            &= -\exp\left(i \pi\right) w \\
            &= w
+   \end{aligned}
 
 and indeed get back `w`. (In case of branched functions we have to be
 aware of branch cuts. In that case we take `w` to be a positive real
diff --git a/docs/modules/printing.rst b/docs/modules/printing.rst
index b83e2c45abb..365284e3849 100644
--- a/docs/modules/printing.rst
+++ b/docs/modules/printing.rst
@@ -54,11 +54,11 @@ to strings of C code).
 
 Usage::
 
-    >>> print_ccode(sin(x)**2 + cos(x)**2)
+    >>> print(ccode(sin(x)**2 + cos(x)**2))
     pow(sin(x), 2) + pow(cos(x), 2)
-    >>> print_ccode(2*x + cos(x), assign_to="result")
+    >>> print(ccode(2*x + cos(x), assign_to="result"))
     result = 2*x + cos(x);
-    >>> print_ccode(Abs(x**2))
+    >>> print(ccode(Abs(x**2)))
     fabs(pow(x, 2))
 
 .. autoclass:: diofant.printing.ccode.CCodePrinter
@@ -69,8 +69,6 @@ Usage::
 
 .. autofunction:: diofant.printing.ccode.ccode
 
-.. autofunction:: diofant.printing.ccode.print_ccode
-
 Fortran Printing
 ----------------
 
@@ -262,7 +260,7 @@ PythonPrinter
 
 This class implements Python printing. Usage::
 
-    >>> print_python(5*x**3 + sin(x))
+    >>> print(python(5*x**3 + sin(x)))
     x = Symbol('x')
     e = 5*x**3 + sin(x)
 
diff --git a/docs/modules/solvers/diophantine.rst b/docs/modules/solvers/diophantine.rst
index c22711d19ce..1bb55a86965 100644
--- a/docs/modules/solvers/diophantine.rst
+++ b/docs/modules/solvers/diophantine.rst
@@ -316,7 +316,7 @@ References
 * Diophantine Equation, Wolfram Mathworld, [online]. Available:
   http://mathworld.wolfram.com/DiophantineEquation.html
 * Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0,[online],
-  Available: https://www.alpertron.com.ar/METHODS.HTM
+  Available: https://web.archive.org/web/20181231080858/https://www.alpertron.com.ar/METHODS.HTM
 * Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
   Available: https://web.archive.org/web/20180831180321/http://www.jpr2718.org/ax2p.pdf
 
diff --git a/docs/release/notes-0.11.rst b/docs/release/notes-0.11.rst
index 87c42240667..4c418c262e0 100644
--- a/docs/release/notes-0.11.rst
+++ b/docs/release/notes-0.11.rst
@@ -29,6 +29,9 @@ Compatibility breaks
 * Change order of keyword arguments for :meth:`~diofant.polys.rings.PolyElement.integrate`, see :pull:`834`.
 * Removed support for ``dps=''`` in :class:`~diofant.core.numbers.Float`.  Significant digits automatically counted for :class:`int` and :class:`str` inputs, see :pull:`797`.
 * Removed ``numer/denom`` properties of :class:`~diofant.polys.fields.FracElement`, see :pull:`851`.
+* Removed ``is_hermitian/is_antihermitian`` core properties, see :pull:`873`.
+* :meth:`~diofant.polys.polytools.Poly.intervals` support only scalar arguments, see :pull:`874`.
+* Removed ``print_python()`` and `print_ccode()` functions, see :pull:`891`.
 
 Minor changes
 =============
@@ -99,3 +102,11 @@ These Sympy issues also were addressed:
 * :sympyissue:`17034` isqrt gives incorrect results
 * :sympyissue:`17044` is_square gives incorrect answers
 * :sympyissue:`10996` Bug in polynomial GCD computation
+* :sympyissue:`15282` Works too long on some limits with big powers
+* :sympyissue:`16722` limit(binomial(n + z, n)*n**-z, n, oo) gives different answers based on assumptions of n and z
+* :sympyissue:`15673` Wrong results. (Limit, Integral, sphere(Space polar coordinates))
+* :sympyissue:`17380` Incorrect results given by some limit expressions
+* :sympyissue:`17431` Wrong results. (Limit, factorial, Power)
+* :sympyissue:`17492` Add link to GitHub in the Sphinx documentation
+* :sympyissue:`17555` (-x).is_extended_positive fails for extended_real and infinite
+* :sympyissue:`17556` Mul.is_imaginary fails for infinite values
diff --git a/docs/release/notes-0.7.3.rst b/docs/release/notes-0.7.3.rst
index 44bb9e4bcdb..172555141dd 100644
--- a/docs/release/notes-0.7.3.rst
+++ b/docs/release/notes-0.7.3.rst
@@ -37,7 +37,7 @@ Major changes
 
 * Context manager for New Assumptions
 
-  - Added the ``with assuming(*facts)`` context manager for new assumptions.  See `blogpost <https://matthewrocklin.com/blog/work/2013/02/05/Assuming>`_.
+  - Added the ``with assuming(*facts)`` context manager for new assumptions.  See `blogpost <https://web.archive.org/web/20181111092915/https://matthewrocklin.com/blog/work/2013/02/05/Assuming>`_.
 
 Compatibility breaks
 ====================
diff --git a/docs/robots.txt b/docs/robots.txt
new file mode 100644
index 00000000000..59f353c04e9
--- /dev/null
+++ b/docs/robots.txt
@@ -0,0 +1,4 @@
+User-agent: *
+Disallow: /en/latest/
+Disallow: /en/v0.9.0/
+Disallow: /en/v0.8.0/
diff --git a/docs/sources.bib b/docs/sources.bib
index b5d13de6ab0..4fb968c05b5 100644
--- a/docs/sources.bib
+++ b/docs/sources.bib
@@ -1511,7 +1511,7 @@ @book{Petkovsek1997AeqB
     year          = {1997},
     publisher     = pub:peters,
     address       = adr:wellesley,
-    url           = {http://www.math.rutgers.edu/~zeilberg/AeqB.pdf},
+    url           = {http://sites.math.rutgers.edu/~zeilberg/AeqB.pdf},
     keywords      = {},
 }
 
@@ -2582,7 +2582,7 @@ @inproceedings{Roach1996
     location      = adr:zurich,
     pages         = {301--308},
     numpages      = {8},
-    url           = {http://doi.acm.org/10.1145/236869.237088},
+    url           = {https://dl.acm.org/citation.cfm?doid=236869.237088},
     doi           = {10.1145/236869.237088},
     acmid         = {237088},
     publisher     = pub:acm,
@@ -2599,7 +2599,7 @@ @inproceedings{Roach1997
     location      = adr:kihei,
     pages         = {205--211},
     numpages      = {7},
-    url           = {http://doi.acm.org/10.1145/258726.258784},
+    url           = {https://dl.acm.org/citation.cfm?doid=258726.258784},
     doi           = {10.1145/258726.258784},
     acmid         = {258784},
     publisher     = pub:acm,
@@ -2627,7 +2627,7 @@ @article{Collins1967subr
     issn          = {0004-5411},
     pages         = {128--142},
     numpages      = {15},
-    url           = {http://doi.acm.org/10.1145/321371.321381},
+    url           = {https://dl.acm.org/citation.cfm?doid=321371.321381},
     doi           = {10.1145/321371.321381},
     acmid         = {321381},
     publisher     = pub:acm,
@@ -2674,7 +2674,7 @@ @inproceedings{Man1994disp
     location      = adr:oxford,
     pages         = {175--180},
     numpages      = {6},
-    url           = {http://doi.acm.org/10.1145/190347.190413},
+    url           = {https://dl.acm.org/citation.cfm?doid=190347.190413},
     doi           = {10.1145/190347.190413},
     acmid         = {190413},
     publisher     = pub:acm,
@@ -2693,7 +2693,7 @@ @article{Man1993indefsum
     issn          = {0747-7171},
     pages         = {355--376},
     numpages      = {22},
-    url           = {http://dx.doi.org/10.1006/jsco.1993.1053},
+    url           = {https://www.sciencedirect.com/science/article/pii/S0747717183710539},
     doi           = {10.1006/jsco.1993.1053},
     acmid         = {181787},
     publisher     = pub:academic,
@@ -2724,7 +2724,7 @@ @article{Kredel1988indep
     year          = {1988},
     issn          = {0747-7171},
     doi           = {10.1016/S0747-7171(88)80045-2},
-    url           = {http://www.sciencedirect.com/science/article/pii/S0747717188800452},
+    url           = {https://www.sciencedirect.com/science/article/pii/S0747717188800452},
     author        = {Heinz Kredel and Volker Weispfenning},
 }
 
@@ -2738,7 +2738,7 @@ @inproceedings{Monagan2000Brown
     location      = {St. Andrews, Scotland},
     pages         = {225--233},
     numpages      = {9},
-    url           = {http://doi.acm.org/10.1145/345542.345639},
+    url           = {https://dl.acm.org/citation.cfm?doid=345542.345639},
     doi           = {10.1145/345542.345639},
     acmid         = {345639},
     publisher     = pub:acm,
@@ -2757,7 +2757,7 @@ @article{Faugere1993groebner
     issn          = {0747-7171},
     pages         = {329--344},
     numpages      = {16},
-    url           = {http://dx.doi.org/10.1006/jsco.1993.1051},
+    url           = {https://www.sciencedirect.com/science/article/pii/S0747717183710515},
     doi           = {10.1006/jsco.1993.1051},
     acmid         = {181785},
     publisher     = pub:academic,
@@ -2770,7 +2770,7 @@ @article{SunWang2010f5
     journal       = {CoRR},
     volume        = {abs/1004.0084},
     year          = {2010},
-    url           = {http://arxiv.org/abs/1004.0084},
+    url           = {https://arxiv.org/abs/1004.0084},
     archivePrefix = {arXiv},
     eprint        = {1004.0084},
     timestamp     = {Mon, 13 Aug 2018 16:47:47 +0200},
@@ -2787,7 +2787,7 @@ @article{Yokoyama1989primitive
     year          = {1989},
     issn          = {0747-7171},
     doi           = {10.1016/S0747-7171(89)80061-6},
-    url           = {http://www.sciencedirect.com/science/article/pii/S0747717189800616},
+    url           = {https://www.sciencedirect.com/science/article/pii/S0747717189800616},
     author        = {Kazuhiro Yokoyama and Masayuki Noro and Taku Takeshima},
 }
 
@@ -2801,7 +2801,7 @@ @inproceedings{Shoup1991ffactor
     location      = adr:bonn,
     pages         = {14--21},
     numpages      = {8},
-    url           = {http://doi.acm.org/10.1145/120694.120697},
+    url           = {https://dl.acm.org/citation.cfm?doid=120694.120697},
     doi           = {10.1145/120694.120697},
     acmid         = {120697},
     publisher     = pub:acm,
@@ -2820,7 +2820,7 @@ @article{Kaltofen1998subquadratic
     issn          = {0025-5718},
     pages         = {1179--1197},
     numpages      = {19},
-    url           = {http://dx.doi.org/10.1090/S0025-5718-98-00944-2},
+    url           = {http://www.ams.org/journals/mcom/1998-67-223/S0025-5718-98-00944-2/home.html},
     doi           = {10.1090/S0025-5718-98-00944-2},
     acmid         = {294748},
     publisher     = pub:ams,
@@ -2836,7 +2836,7 @@ @article{Abramov71rat
     year          = {1971},
     issn          = {0041-5553},
     doi           = {10.1016/0041-5553(71)90028-0},
-    url           = {http://www.sciencedirect.com/science/article/pii/0041555371900280},
+    url           = {https://www.sciencedirect.com/science/article/abs/pii/0041555371900280},
     author        = {S.A. Abramov}
 }
 
@@ -2849,7 +2849,7 @@ @article{Arno1996alg
     year          = {1996},
     issn          = {0022-314X},
     doi           = {10.1006/jnth.1996.0049},
-    url           = {http://www.sciencedirect.com/science/article/pii/S0022314X96900499},
+    url           = {https://www.sciencedirect.com/science/article/pii/S0022314X96900499},
     author        = {Steven Arno and M.L. Robinson and Ferell S. Wheeler},
 }
 
@@ -2872,7 +2872,7 @@ @article{Brown1978prs
     issn          = {0098-3500},
     pages         = {237--249},
     numpages      = {13},
-    url           = {http://doi.acm.org/10.1145/355791.355795},
+    url           = {https://dl.acm.org/citation.cfm?doid=355791.355795},
     doi           = {10.1145/355791.355795},
     acmid         = {355795},
     publisher     = pub:acm,
@@ -2891,7 +2891,7 @@ @article{Hoeven02
     issn          = {0747-7171},
     pages         = {479--542},
     numpages      = {64},
-    url           = {http://dx.doi.org/10.1006/jsco.2002.0562},
+    url           = {https://www.sciencedirect.com/science/article/pii/S0747717102905626},
     doi           = {10.1006/jsco.2002.0562},
     acmid         = {766453},
     publisher     = pub:academic,
diff --git a/setup.cfg b/setup.cfg
index 6de4321a0da..192ca6075cc 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -28,7 +28,7 @@ classifiers = Development Status :: 4 - Beta
 [options]
 zip_safe = True
 packages = find:
-python_requires = >=3.5
+python_requires = >=3.7
 setup_requires = setuptools>=36.7.0
                  setuptools_scm
                  pip>=9.0.1
@@ -39,7 +39,7 @@ install_requires = mpmath>=0.19
 tests_require = diofant[develop]
 [options.extras_require]
 exports =
-          numpy>=1.12.1
+          numpy>=1.12.1,<=1.16.5
           scipy
           cython
 gmpy =
@@ -55,8 +55,7 @@ tests = pytest>=3.10.1
         pytest-timeout
 develop = %(tests)s
           flake8>=2.5.5,!=3.1.0
-          flake8-docstrings>=1.2.0
-          pydocstyle!=4.0.0
+          flake8-docstrings>=1.3.1
           pep8-naming
           flake8-comprehensions
           flake8-isort
@@ -86,8 +85,9 @@ doctest_optionflags = ELLIPSIS NORMALIZE_WHITESPACE IGNORE_EXCEPTION_DETAIL
 addopts = --durations=20
           -r X
           --doctest-glob='*.rst'
+          --cov-config=setup.cfg
 norecursedirs = build .eggs .git
-timeout = 600
+timeout = 700
 xfail_strict = true
 filterwarnings = ignore::UserWarning
                  error::DeprecationWarning