Repository: sympy/sympy
Pull request: Clean up flake8 warnings from doctests
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #19374 

#### Brief description of what is fixed or changed
running `flake8 --doctests sympy/` returns a bunch of errors. This PR is to fix them

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* other
  * Improve code quality by cleaning the doctests up
<!-- END RELEASE NOTES -->

Problem statement:
Clean up flake8 warnings from doctests
Currently flake8 runs on the main code but not on doctests. It has a `--doctests` argument to check the doctests as well.

For example:
```console
$ flake8 --doctests sympy/matrices/determinant.py 
sympy/matrices/determinant.py:545:9: F401 'sympy.MatrixSymbol' imported but unused
sympy/matrices/determinant.py:545:9: F401 'sympy.eye' imported but unused
sympy/matrices/determinant.py:545:9: F401 'sympy.det' imported but unused
```
That can be fixed with
```diff
diff --git a/sympy/matrices/determinant.py b/sympy/matrices/determinant.py
index 647541afef..f3bb9a4067 100644
--- a/sympy/matrices/determinant.py
+++ b/sympy/matrices/determinant.py
@@ -542,7 +542,7 @@ def _det(M, method="bareiss", iszerofunc=None):
     Examples
     ========
 
-    >>> from sympy import Matrix, MatrixSymbol, eye, det
+    >>> from sympy import Matrix
     >>> M = Matrix([[1, 2], [3, 4]])
     >>> M.det()
     -2
```

Running flake8 on the whole of sympy with
```
$ flake8 --doctests sympy
```
leads to 698 warnings almost all of which are F401 (unnecessary import) and so easily fixable.

These should all be fixed and then we should add the `--doctest` argument on Travis.

I'm marking this as easy to fix.

Hints from issue discussion:
I'd like to take up on the task of cleaning the imports if I may
Yes, go ahead.
I was working on this cleanup and I found a curious situation:

sympy supports multiple dispatch per argument type, as defined in the docs of `sympy/multipledispatch/core.py`(lines 22-38):
```
Examples
    --------

    >>> from sympy.multipledispatch import dispatch
    >>> @dispatch(int)
    ... def f(x):
    ...     return x + 1

    >>> @dispatch(float)
    ... def f(x):
    ...     return x - 1

    >>> f(3)
    4
    >>> f(3.0)
    2.0

```

`flake8`, on the other hand, reports an error `F811 redefinition of unused 'f' from line 27`
which makes sense from its point of view but it seems that this is actually a wanted behaviour

This could be easily solved by invoking `flake8 --doctests --ignore=F811` in Travis, but I just wanted to make sure that's a correct way, and if it's not I'd like to know a better solution
Also, doctest from `sympy/core/singleton.py` (line 49) throws a warning about comparing Singleton and a literal with `is`, and although technically not an error, I'd like to know is there anything for me to do to remove this warning

This is the output I get:
![Screenshot_20200521_044553](https://user-images.githubusercontent.com/12953598/82517703-fa110e00-9b1d-11ea-9cc5-c7faf7efdb82.png)

I think we can write it like:
```
>>> x = 0
>>> x is S.Zero
```
to avoid the warning.

Multiple dispatch gives lots of flake8 warnings. These have been suppressed with a flake8 flag:
https://github.com/sympy/sympy/blob/0f598f042929f9de97bd8ba25ac06a0e4d36881b/sympy/sets/handlers/intersection.py#L14-L16
I guess we can do the same in the doctests if needed.

Ultimately it would be good to find a better way of using multiple dispatch so we don't repeat the same function name #17795.
I will add the noqa then, since I feel this second solution is out of my scope
I've said this before, but I'm not a fan of when linters force you to write code in a less clear way, and especially when this is done so in documentation. 
Having noqa cluttering the documentation is also not great. Although it isn't so bad if Sphinx automatically removes it like it does for # DOCTEST: +SKIP.
Hey, what would you like me to do with `sympy/utilities/runtests.py`? `flake8 --doctests sympy/utilities/runtests.py` reports an `F403` (which is interesting I guess since it has no doctest code)

edit: same for `benchmarking.py`, `randtest.py`, `tmpfiles.py`, and `quality_unicode.py`.
I think we can avoid using noqa for multiple dispatch by making all the functions called `_` as is done for single dispatch.

Commit message: quality(doc): Test doctests with flake8 on Travis

Mostly removes no longer needed imports.

Co-authored-by: Marijan Smetko <marijansmetko123@gmail.com>

Patch to review:
diff --git a/bin/test_travis.sh b/bin/test_travis.sh
index 51f553cbbf91..104a464bed00 100755
--- a/bin/test_travis.sh
+++ b/bin/test_travis.sh
@@ -7,6 +7,7 @@ set -x
 
 if [[ "${TEST_FLAKE8}" == "true" ]]; then
     flake8 sympy;
+    flake8 --doctests --per-file-ignores='sympy/interactive/session.py:F821' sympy;
 fi
 
 if [[ "${TEST_SETUP}" == "true" ]]; then
diff --git a/sympy/assumptions/assume.py b/sympy/assumptions/assume.py
index 7efa484c0b9c..92e459320d22 100644
--- a/sympy/assumptions/assume.py
+++ b/sympy/assumptions/assume.py
@@ -17,7 +17,7 @@ class to create your own local assumptions contexts. It is basically a thin
     Examples
     ========
 
-    >>> from sympy import AppliedPredicate, Q
+    >>> from sympy import Q
     >>> from sympy.assumptions.assume import global_assumptions
     >>> global_assumptions
     AssumptionsContext()
diff --git a/sympy/assumptions/handlers/calculus.py b/sympy/assumptions/handlers/calculus.py
index 702972eb3091..4bab975a4fbc 100644
--- a/sympy/assumptions/handlers/calculus.py
+++ b/sympy/assumptions/handlers/calculus.py
@@ -16,7 +16,7 @@ class AskFiniteHandler(CommonHandler):
 
     Examples of usage:
 
-    >>> from sympy import Symbol, Q
+    >>> from sympy import Q
     >>> from sympy.assumptions.handlers.calculus import AskFiniteHandler
     >>> from sympy.abc import x
     >>> a = AskFiniteHandler()
@@ -35,7 +35,7 @@ def Symbol(expr, assumptions):
         Examples
         ========
 
-        >>> from sympy import Symbol, Q
+        >>> from sympy import Q
         >>> from sympy.assumptions.handlers.calculus import AskFiniteHandler
         >>> from sympy.abc import x
         >>> a = AskFiniteHandler()
diff --git a/sympy/assumptions/refine.py b/sympy/assumptions/refine.py
index cfe630a63d65..a50d506a8fe0 100644
--- a/sympy/assumptions/refine.py
+++ b/sympy/assumptions/refine.py
@@ -53,7 +53,7 @@ def refine_abs(expr, assumptions):
     Examples
     ========
 
-    >>> from sympy import Symbol, Q, refine, Abs
+    >>> from sympy import Q, Abs
     >>> from sympy.assumptions.refine import refine_abs
     >>> from sympy.abc import x
     >>> refine_abs(Abs(x), Q.real(x))
@@ -89,7 +89,7 @@ def refine_Pow(expr, assumptions):
     """
     Handler for instances of Pow.
 
-    >>> from sympy import Symbol, Q
+    >>> from sympy import Q
     >>> from sympy.assumptions.refine import refine_Pow
     >>> from sympy.abc import x,y,z
     >>> refine_Pow((-1)**x, Q.real(x))
@@ -190,7 +190,7 @@ def refine_atan2(expr, assumptions):
     Examples
     ========
 
-    >>> from sympy import Symbol, Q, refine, atan2
+    >>> from sympy import Q, atan2
     >>> from sympy.assumptions.refine import refine_atan2
     >>> from sympy.abc import x, y
     >>> refine_atan2(atan2(y,x), Q.real(y) & Q.positive(x))
diff --git a/sympy/calculus/finite_diff.py b/sympy/calculus/finite_diff.py
index f09fca706069..c91589c22a06 100644
--- a/sympy/calculus/finite_diff.py
+++ b/sympy/calculus/finite_diff.py
@@ -441,7 +441,7 @@ def differentiate_finite(expr, *symbols,
     Examples
     ========
 
-    >>> from sympy import cos, sin, Function, differentiate_finite
+    >>> from sympy import sin, Function, differentiate_finite
     >>> from sympy.abc import x, y, h
     >>> f, g = Function('f'), Function('g')
     >>> differentiate_finite(f(x)*g(x), x, points=[x-h, x+h])
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index dc893b0531c9..50a9be1c8282 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -720,7 +720,7 @@ def stationary_points(f, symbol, domain=S.Reals):
     Examples
     ========
 
-    >>> from sympy import Symbol, S, sin, log, pi, pprint, stationary_points
+    >>> from sympy import Symbol, S, sin, pi, pprint, stationary_points
     >>> from sympy.sets import Interval
     >>> x = Symbol('x')
 
diff --git a/sympy/categories/baseclasses.py b/sympy/categories/baseclasses.py
index 98fdbab91372..a29fa82e92eb 100644
--- a/sympy/categories/baseclasses.py
+++ b/sympy/categories/baseclasses.py
@@ -563,7 +563,7 @@ class Diagram(Basic):
     ========
 
     >>> from sympy.categories import Object, NamedMorphism, Diagram
-    >>> from sympy import FiniteSet, pprint, default_sort_key
+    >>> from sympy import pprint, default_sort_key
     >>> A = Object("A")
     >>> B = Object("B")
     >>> C = Object("C")
diff --git a/sympy/codegen/array_utils.py b/sympy/codegen/array_utils.py
index 1fca53698284..bae4a9930295 100644
--- a/sympy/codegen/array_utils.py
+++ b/sympy/codegen/array_utils.py
@@ -299,8 +299,8 @@ def _get_contraction_tuples(self):
         Examples
         ========
 
-        >>> from sympy import MatrixSymbol, MatrixExpr, Sum, Symbol
-        >>> from sympy.abc import i, j, k, l, N
+        >>> from sympy import MatrixSymbol
+        >>> from sympy.abc import N
         >>> from sympy.codegen.array_utils import CodegenArrayContraction, CodegenArrayTensorProduct
         >>> A = MatrixSymbol("A", N, N)
         >>> B = MatrixSymbol("B", N, N)
@@ -362,8 +362,8 @@ def sort_args_by_name(self):
         Examples
         ========
 
-        >>> from sympy import MatrixSymbol, MatrixExpr, Sum, Symbol
-        >>> from sympy.abc import i, j, k, l, N
+        >>> from sympy import MatrixSymbol
+        >>> from sympy.abc import N
         >>> from sympy.codegen.array_utils import parse_matrix_expression
         >>> A = MatrixSymbol("A", N, N)
         >>> B = MatrixSymbol("B", N, N)
@@ -402,8 +402,8 @@ def _get_contraction_links(self):
         Examples
         ========
 
-        >>> from sympy import MatrixSymbol, MatrixExpr, Sum, Symbol
-        >>> from sympy.abc import i, j, k, l, N
+        >>> from sympy import MatrixSymbol
+        >>> from sympy.abc import N
         >>> from sympy.codegen.array_utils import parse_matrix_expression
         >>> A = MatrixSymbol("A", N, N)
         >>> B = MatrixSymbol("B", N, N)
@@ -566,7 +566,6 @@ def nest_permutation(self):
 
         >>> from sympy.codegen.array_utils import (CodegenArrayPermuteDims, CodegenArrayTensorProduct, nest_permutation)
         >>> from sympy import MatrixSymbol
-        >>> from sympy.combinatorics import Permutation
 
         >>> M = MatrixSymbol("M", 3, 3)
         >>> N = MatrixSymbol("N", 3, 3)
@@ -1061,7 +1060,6 @@ def parse_indexed_expression(expr, first_indices=None):
 
     >>> from sympy.codegen.array_utils import parse_indexed_expression
     >>> from sympy import MatrixSymbol, Sum, symbols
-    >>> from sympy.combinatorics import Permutation
 
     >>> i, j, k, d = symbols("i j k d")
     >>> M = MatrixSymbol("M", d, d)
@@ -1250,7 +1248,7 @@ def recognize_matrix_expression(expr):
     Examples
     ========
 
-    >>> from sympy import MatrixSymbol, MatrixExpr, Sum, Symbol
+    >>> from sympy import MatrixSymbol, Sum
     >>> from sympy.abc import i, j, k, l, N
     >>> from sympy.codegen.array_utils import CodegenArrayContraction, CodegenArrayTensorProduct
     >>> from sympy.codegen.array_utils import recognize_matrix_expression, parse_indexed_expression, parse_matrix_expression
diff --git a/sympy/codegen/ast.py b/sympy/codegen/ast.py
index b25bf0850bc8..7f2ebc2f4eb0 100644
--- a/sympy/codegen/ast.py
+++ b/sympy/codegen/ast.py
@@ -968,7 +968,6 @@ class Type(Token):
       ...
     ValueError: Casting gives a significantly different value.
     >>> boost_mp50 = Type('boost::multiprecision::cpp_dec_float_50')
-    >>> from sympy import Symbol
     >>> from sympy.printing.cxxcode import cxxcode
     >>> from sympy.codegen.ast import Declaration, Variable
     >>> cxxcode(Declaration(Variable('x', type=boost_mp50)))
@@ -1052,7 +1051,7 @@ def cast_check(self, value, rtol=None, atol=0, limits=None, precision_targets=No
         Examples
         ========
 
-        >>> from sympy.codegen.ast import Type, integer, float32, int8
+        >>> from sympy.codegen.ast import integer, float32, int8
         >>> integer.cast_check(3.0) == 3
         True
         >>> float32.cast_check(1e-40)  # doctest: +ELLIPSIS
@@ -1168,7 +1167,7 @@ class FloatType(FloatBaseType):
     Examples
     ========
 
-    >>> from sympy import S, Float
+    >>> from sympy import S
     >>> from sympy.codegen.ast import FloatType
     >>> half_precision = FloatType('f16', nbits=16, nmant=10, nexp=5)
     >>> half_precision.max
@@ -1564,8 +1563,7 @@ class Declaration(Token):
     Examples
     ========
 
-    >>> from sympy import Symbol
-    >>> from sympy.codegen.ast import Declaration, Type, Variable, NoneToken, integer, untyped
+    >>> from sympy.codegen.ast import Declaration, NoneToken, untyped
     >>> z = Declaration('z')
     >>> z.variable.type == untyped
     True
diff --git a/sympy/codegen/fnodes.py b/sympy/codegen/fnodes.py
index 5f34e3aef46a..18a08ac7ef5a 100644
--- a/sympy/codegen/fnodes.py
+++ b/sympy/codegen/fnodes.py
@@ -454,7 +454,7 @@ def size(array, dim=None, kind=None):
 
     >>> from sympy import Symbol
     >>> from sympy.printing import fcode
-    >>> from sympy.codegen.ast import FunctionDefinition, real, Return, Variable
+    >>> from sympy.codegen.ast import FunctionDefinition, real, Return
     >>> from sympy.codegen.fnodes import array, sum_, size
     >>> a = Symbol('a', real=True)
     >>> body = [Return((sum_(a**2)/size(a))**.5)]
@@ -505,8 +505,8 @@ def bind_C(name=None):
 
     >>> from sympy import Symbol
     >>> from sympy.printing import fcode
-    >>> from sympy.codegen.ast import FunctionDefinition, real, Return, Variable
-    >>> from sympy.codegen.fnodes import array, sum_, size, bind_C
+    >>> from sympy.codegen.ast import FunctionDefinition, real, Return
+    >>> from sympy.codegen.fnodes import array, sum_, bind_C
     >>> a = Symbol('a', real=True)
     >>> s = Symbol('s', integer=True)
     >>> arr = array(a, dim=[s], intent='in')
diff --git a/sympy/codegen/rewriting.py b/sympy/codegen/rewriting.py
index 790a3995c4a2..df2cb513bc31 100644
--- a/sympy/codegen/rewriting.py
+++ b/sympy/codegen/rewriting.py
@@ -73,7 +73,7 @@ class ReplaceOptim(Optimization):
     Examples
     ========
 
-    >>> from sympy import Symbol, Pow
+    >>> from sympy import Symbol
     >>> from sympy.codegen.rewriting import ReplaceOptim
     >>> from sympy.codegen.cfunctions import exp2
     >>> x = Symbol('x')
diff --git a/sympy/combinatorics/coset_table.py b/sympy/combinatorics/coset_table.py
index 9e9b2b0f7ecf..e99d3b6af546 100644
--- a/sympy/combinatorics/coset_table.py
+++ b/sympy/combinatorics/coset_table.py
@@ -1171,7 +1171,7 @@ def modified_coset_enumeration_r(fp_grp, Y, max_cosets=None, draft=None,
     ========
 
     >>> from sympy.combinatorics.free_groups import free_group
-    >>> from sympy.combinatorics.fp_groups import FpGroup, coset_enumeration_r
+    >>> from sympy.combinatorics.fp_groups import FpGroup
     >>> from sympy.combinatorics.coset_table import modified_coset_enumeration_r
     >>> F, x, y = free_group("x, y")
     >>> f = FpGroup(F, [x**3, y**3, x**-1*y**-1*x*y])
diff --git a/sympy/combinatorics/fp_groups.py b/sympy/combinatorics/fp_groups.py
index 1313c3d3ee6d..379a9b8c9bad 100644
--- a/sympy/combinatorics/fp_groups.py
+++ b/sympy/combinatorics/fp_groups.py
@@ -127,7 +127,7 @@ def subgroup(self, gens, C=None, homomorphism=False):
         Examples
         ========
 
-        >>> from sympy.combinatorics.fp_groups import (FpGroup, FpSubgroup)
+        >>> from sympy.combinatorics.fp_groups import FpGroup
         >>> from sympy.combinatorics.free_groups import free_group
         >>> F, x, y = free_group("x, y")
         >>> f = FpGroup(F, [x**3, y**5, (x*y)**2])
diff --git a/sympy/combinatorics/free_groups.py b/sympy/combinatorics/free_groups.py
index 26a3afbef724..6f670a114c74 100644
--- a/sympy/combinatorics/free_groups.py
+++ b/sympy/combinatorics/free_groups.py
@@ -77,7 +77,7 @@ def vfree_group(symbols):
     >>> from sympy.combinatorics.free_groups import vfree_group
     >>> vfree_group("x, y, z")
     <free group on the generators (x, y, z)>
-    >>> x**2*y**-2*z
+    >>> x**2*y**-2*z # noqa: F821
     x**2*y**-2*z
     >>> type(_)
     <class 'sympy.combinatorics.free_groups.FreeGroupElement'>
diff --git a/sympy/combinatorics/generators.py b/sympy/combinatorics/generators.py
index 02d8b5e2c38c..3104c2fc3afc 100644
--- a/sympy/combinatorics/generators.py
+++ b/sympy/combinatorics/generators.py
@@ -11,7 +11,6 @@ def symmetric(n):
     Examples
     ========
 
-    >>> from sympy.combinatorics.permutations import Permutation
     >>> from sympy.combinatorics.generators import symmetric
     >>> list(symmetric(3))
     [(2), (1 2), (2)(0 1), (0 1 2), (0 2 1), (0 2)]
@@ -27,7 +26,6 @@ def cyclic(n):
     Examples
     ========
 
-    >>> from sympy.combinatorics.permutations import Permutation
     >>> from sympy.combinatorics.generators import cyclic
     >>> list(cyclic(5))
     [(4), (0 1 2 3 4), (0 2 4 1 3),
@@ -51,7 +49,6 @@ def alternating(n):
     Examples
     ========
 
-    >>> from sympy.combinatorics.permutations import Permutation
     >>> from sympy.combinatorics.generators import alternating
     >>> list(alternating(3))
     [(2), (0 1 2), (0 2 1)]
@@ -73,7 +70,6 @@ def dihedral(n):
     Examples
     ========
 
-    >>> from sympy.combinatorics.permutations import Permutation
     >>> from sympy.combinatorics.generators import dihedral
     >>> list(dihedral(3))
     [(2), (0 2), (0 1 2), (1 2), (0 2 1), (2)(0 1)]
diff --git a/sympy/combinatorics/homomorphisms.py b/sympy/combinatorics/homomorphisms.py
index 9612121b85da..9e6610eb6b08 100644
--- a/sympy/combinatorics/homomorphisms.py
+++ b/sympy/combinatorics/homomorphisms.py
@@ -447,7 +447,7 @@ def group_isomorphism(G, H, isomorphism=True):
     >>> from sympy.combinatorics.perm_groups import PermutationGroup
     >>> from sympy.combinatorics.free_groups import free_group
     >>> from sympy.combinatorics.fp_groups import FpGroup
-    >>> from sympy.combinatorics.homomorphisms import homomorphism, group_isomorphism
+    >>> from sympy.combinatorics.homomorphisms import group_isomorphism
     >>> from sympy.combinatorics.named_groups import DihedralGroup, AlternatingGroup
 
     >>> D = DihedralGroup(8)
diff --git a/sympy/combinatorics/named_groups.py b/sympy/combinatorics/named_groups.py
index c0c37e40cc0c..dbdb0706d787 100644
--- a/sympy/combinatorics/named_groups.py
+++ b/sympy/combinatorics/named_groups.py
@@ -16,7 +16,6 @@ def AbelianGroup(*cyclic_orders):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.named_groups import AbelianGroup
     >>> AbelianGroup(3, 4)
     PermutationGroup([
diff --git a/sympy/combinatorics/perm_groups.py b/sympy/combinatorics/perm_groups.py
index 2220c4e8a953..de94ddabb4af 100644
--- a/sympy/combinatorics/perm_groups.py
+++ b/sympy/combinatorics/perm_groups.py
@@ -33,7 +33,6 @@ class PermutationGroup(Basic):
     ========
 
     >>> from sympy.combinatorics import Permutation
-    >>> from sympy.combinatorics.permutations import Cycle
     >>> from sympy.combinatorics.polyhedron import Polyhedron
     >>> from sympy.combinatorics.perm_groups import PermutationGroup
 
@@ -254,7 +253,6 @@ def __mul__(self, other):
         Examples
         ========
 
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import CyclicGroup
         >>> G = CyclicGroup(5)
         >>> H = G*G
@@ -912,7 +910,6 @@ def center(self):
         Examples
         ========
 
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import DihedralGroup
         >>> D = DihedralGroup(4)
         >>> G = D.center()
@@ -1464,7 +1461,6 @@ def generate(self, method="coset", af=False):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation
         >>> from sympy.combinatorics import PermutationGroup
         >>> from sympy.combinatorics.polyhedron import tetrahedron
 
@@ -1957,7 +1953,6 @@ def is_alt_sym(self, eps=0.05, _random_prec=None):
         Examples
         ========
 
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import DihedralGroup
         >>> D = DihedralGroup(10)
         >>> D.is_alt_sym()
@@ -2101,7 +2096,6 @@ def is_primitive(self, randomized=True):
         Examples
         ========
 
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import DihedralGroup
         >>> D = DihedralGroup(10)
         >>> D.is_primitive()
@@ -2516,7 +2510,6 @@ def minimal_block(self, points):
         Examples
         ========
 
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import DihedralGroup
         >>> D = DihedralGroup(10)
         >>> D.minimal_block([0, 5])
@@ -2801,8 +2794,6 @@ def orbit_rep(self, alpha, beta, schreier_vector=None):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import AlternatingGroup
         >>> G = AlternatingGroup(5)
         >>> G.orbit_rep(0, 4)
@@ -2843,8 +2834,6 @@ def orbit_transversal(self, alpha, pairs=False):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import DihedralGroup
         >>> G = DihedralGroup(6)
         >>> G.orbit_transversal(0)
@@ -3429,7 +3418,6 @@ def schreier_sims_incremental(self, base=None, gens=None, slp_dict=False):
         ========
 
         >>> from sympy.combinatorics.named_groups import AlternatingGroup
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.testutil import _verify_bsgs
         >>> A = AlternatingGroup(7)
         >>> base = [2, 3]
@@ -3614,7 +3602,6 @@ def schreier_sims_random(self, base=None, gens=None, consec_succ=10,
         Examples
         ========
 
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.testutil import _verify_bsgs
         >>> from sympy.combinatorics.named_groups import SymmetricGroup
         >>> S = SymmetricGroup(5)
@@ -3769,8 +3756,6 @@ def stabilizer(self, alpha):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import DihedralGroup
         >>> G = DihedralGroup(6)
         >>> G.stabilizer(5)
@@ -3872,7 +3857,6 @@ def subgroup_search(self, prop, base=None, strong_gens=None, tests=None,
 
         >>> from sympy.combinatorics.named_groups import (SymmetricGroup,
         ... AlternatingGroup)
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.testutil import _verify_bsgs
         >>> S = SymmetricGroup(7)
         >>> prop_even = lambda x: x.is_even
@@ -4598,7 +4582,6 @@ def strong_presentation(G):
         Examples
         ========
 
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
         >>> from sympy.combinatorics.named_groups import DihedralGroup
         >>> P = DihedralGroup(4)
         >>> G = P.strong_presentation()
@@ -4949,7 +4932,7 @@ def _orbits(degree, generators):
     ========
 
     >>> from sympy.combinatorics.permutations import Permutation
-    >>> from sympy.combinatorics.perm_groups import PermutationGroup, _orbits
+    >>> from sympy.combinatorics.perm_groups import _orbits
     >>> a = Permutation([0, 2, 1])
     >>> b = Permutation([1, 0, 2])
     >>> _orbits(a.size, [a, b])
@@ -4991,7 +4974,6 @@ def _orbit_transversal(degree, generators, alpha, pairs, af=False, slp=False):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.named_groups import DihedralGroup
     >>> from sympy.combinatorics.perm_groups import _orbit_transversal
     >>> G = DihedralGroup(6)
@@ -5043,7 +5025,6 @@ def _stabilizer(degree, generators, alpha):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.perm_groups import _stabilizer
     >>> from sympy.combinatorics.named_groups import DihedralGroup
     >>> G = DihedralGroup(6)
@@ -5119,7 +5100,7 @@ def order(self):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup
+        >>> from sympy.combinatorics import SymmetricPermutationGroup
         >>> G = SymmetricPermutationGroup(4)
         >>> G.order()
         24
@@ -5138,7 +5119,7 @@ def degree(self):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup
+        >>> from sympy.combinatorics import SymmetricPermutationGroup
         >>> G = SymmetricPermutationGroup(4)
         >>> G.degree
         4
@@ -5154,7 +5135,7 @@ def identity(self):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup
+        >>> from sympy.combinatorics import SymmetricPermutationGroup
         >>> G = SymmetricPermutationGroup(4)
         >>> G.identity()
         (3)
diff --git a/sympy/combinatorics/permutations.py b/sympy/combinatorics/permutations.py
index de7bfd918993..da19cbf50601 100644
--- a/sympy/combinatorics/permutations.py
+++ b/sympy/combinatorics/permutations.py
@@ -322,7 +322,6 @@ def __call__(self, *other):
         ========
 
         >>> from sympy.combinatorics.permutations import Cycle as C
-        >>> from sympy.combinatorics.permutations import Permutation as Perm
         >>> C(1, 2)(2, 3)
         (1 3 2)
 
@@ -355,7 +354,6 @@ def list(self, size=None):
         ========
 
         >>> from sympy.combinatorics.permutations import Cycle
-        >>> from sympy.combinatorics.permutations import Permutation
         >>> p = Cycle(2, 3)(4, 5)
         >>> p.list()
         [0, 1, 3, 2, 5, 4]
@@ -1195,7 +1193,7 @@ def rmul(*args):
         Examples
         ========
 
-        >>> from sympy.combinatorics.permutations import _af_rmul, Permutation
+        >>> from sympy.combinatorics.permutations import Permutation
 
         >>> a, b = [1, 0, 2], [0, 2, 1]
         >>> a = Permutation(a); b = Permutation(b)
diff --git a/sympy/combinatorics/polyhedron.py b/sympy/combinatorics/polyhedron.py
index aadb3e8fde20..78176f11fed4 100644
--- a/sympy/combinatorics/polyhedron.py
+++ b/sympy/combinatorics/polyhedron.py
@@ -425,7 +425,6 @@ def array_form(self):
         Examples
         ========
 
-        >>> from sympy.combinatorics import Permutation, Cycle
         >>> from sympy.combinatorics.polyhedron import tetrahedron
         >>> tetrahedron = tetrahedron.copy()
         >>> tetrahedron.array_form
diff --git a/sympy/combinatorics/tensor_can.py b/sympy/combinatorics/tensor_can.py
index 75797b3f20cd..937a25fdf47e 100644
--- a/sympy/combinatorics/tensor_can.py
+++ b/sympy/combinatorics/tensor_can.py
@@ -387,7 +387,6 @@ def double_coset_can_rep(dummies, sym, b_S, sgens, S_transversals, g):
     ========
 
     >>> from sympy.combinatorics.permutations import Permutation
-    >>> from sympy.combinatorics.perm_groups import PermutationGroup
     >>> from sympy.combinatorics.tensor_can import double_coset_can_rep, get_transversals
     >>> gens = [Permutation(x) for x in [[2, 1, 0, 3, 4, 5, 7, 6], [4, 1, 2, 3, 0, 5, 7, 6]]]
     >>> base = [0, 2]
@@ -913,7 +912,6 @@ def bsgs_direct_product(base1, gens1, base2, gens2, signed=True):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.tensor_can import (get_symmetric_group_sgs, bsgs_direct_product)
     >>> base1, gens1 = get_symmetric_group_sgs(1)
     >>> base2, gens2 = get_symmetric_group_sgs(2)
@@ -947,7 +945,6 @@ def get_symmetric_group_sgs(n, antisym=False):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs
     >>> get_symmetric_group_sgs(3)
     ([0, 1], [(4)(0 1), (4)(1 2)])
@@ -1051,7 +1048,6 @@ def tensor_gens(base, gens, list_free_indices, sym=0):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.tensor_can import tensor_gens, get_symmetric_group_sgs
 
     two symmetric tensors with 3 indices without free indices
@@ -1167,7 +1163,6 @@ def gens_products(*v):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs, gens_products
     >>> base, gens = get_symmetric_group_sgs(2)
     >>> gens_products((base, gens, [[], []], 0))
diff --git a/sympy/combinatorics/testutil.py b/sympy/combinatorics/testutil.py
index d9d80db949c5..c34c1665824e 100644
--- a/sympy/combinatorics/testutil.py
+++ b/sympy/combinatorics/testutil.py
@@ -215,7 +215,7 @@ def canonicalize_naive(g, dummies, sym, *v):
 
     >>> from sympy.combinatorics.testutil import canonicalize_naive
     >>> from sympy.combinatorics.tensor_can import get_symmetric_group_sgs
-    >>> from sympy.combinatorics import Permutation, PermutationGroup
+    >>> from sympy.combinatorics import Permutation
     >>> g = Permutation([1, 3, 2, 0, 4, 5])
     >>> base2, gens2 = get_symmetric_group_sgs(2)
     >>> canonicalize_naive(g, [2, 3], 0, (base2, gens2, 2, 0))
diff --git a/sympy/combinatorics/util.py b/sympy/combinatorics/util.py
index 9019f77b3f85..fe7e33c0081a 100644
--- a/sympy/combinatorics/util.py
+++ b/sympy/combinatorics/util.py
@@ -139,7 +139,6 @@ def _distribute_gens_by_base(base, gens):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.named_groups import DihedralGroup
     >>> from sympy.combinatorics.util import _distribute_gens_by_base
     >>> D = DihedralGroup(3)
@@ -206,7 +205,6 @@ def _handle_precomputed_bsgs(base, strong_gens, transversals=None,
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.named_groups import DihedralGroup
     >>> from sympy.combinatorics.util import _handle_precomputed_bsgs
     >>> D = DihedralGroup(3)
@@ -265,16 +263,13 @@ def _orbits_transversals_from_bsgs(base, strong_gens_distr,
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.named_groups import SymmetricGroup
-    >>> from sympy.combinatorics.util import _orbits_transversals_from_bsgs
-    >>> from sympy.combinatorics.util import (_orbits_transversals_from_bsgs,
-    ... _distribute_gens_by_base)
+    >>> from sympy.combinatorics.util import _distribute_gens_by_base
     >>> S = SymmetricGroup(3)
     >>> S.schreier_sims()
     >>> strong_gens_distr = _distribute_gens_by_base(S.base, S.strong_gens)
-    >>> _orbits_transversals_from_bsgs(S.base, strong_gens_distr)
-    ([[0, 1, 2], [1, 2]], [{0: (2), 1: (0 1 2), 2: (0 2 1)}, {1: (2), 2: (1 2)}])
+    >>> (S.base, strong_gens_distr)
+    ([0, 1], [[(0 1 2), (2)(0 1), (1 2)], [(1 2)]])
 
     See Also
     ========
@@ -326,7 +321,6 @@ def _remove_gens(base, strong_gens, basic_orbits=None, strong_gens_distr=None):
     ========
 
     >>> from sympy.combinatorics.named_groups import SymmetricGroup
-    >>> from sympy.combinatorics.perm_groups import PermutationGroup
     >>> from sympy.combinatorics.util import _remove_gens
     >>> from sympy.combinatorics.testutil import _verify_bsgs
     >>> S = SymmetricGroup(15)
@@ -408,7 +402,6 @@ def _strip(g, base, orbits, transversals):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.named_groups import SymmetricGroup
     >>> from sympy.combinatorics.permutations import Permutation
     >>> from sympy.combinatorics.util import _strip
@@ -501,7 +494,6 @@ def _strong_gens_from_distr(strong_gens_distr):
     Examples
     ========
 
-    >>> from sympy.combinatorics import Permutation
     >>> from sympy.combinatorics.named_groups import SymmetricGroup
     >>> from sympy.combinatorics.util import (_strong_gens_from_distr,
     ... _distribute_gens_by_base)
diff --git a/sympy/concrete/delta.py b/sympy/concrete/delta.py
index a53a9fdad606..7c9f28b5035d 100644
--- a/sympy/concrete/delta.py
+++ b/sympy/concrete/delta.py
@@ -267,7 +267,7 @@ def deltasummation(f, limit, no_piecewise=False):
     >>> from sympy.abc import k
     >>> i, j = symbols('i, j', integer=True, finite=True)
     >>> from sympy.concrete.delta import deltasummation
-    >>> from sympy import KroneckerDelta, Piecewise
+    >>> from sympy import KroneckerDelta
     >>> deltasummation(KroneckerDelta(i, k), (k, -oo, oo))
     1
     >>> deltasummation(KroneckerDelta(i, k), (k, 0, oo))
diff --git a/sympy/concrete/gosper.py b/sympy/concrete/gosper.py
index 5b7c3148986e..794a51c5673e 100644
--- a/sympy/concrete/gosper.py
+++ b/sympy/concrete/gosper.py
@@ -173,7 +173,7 @@ def gosper_sum(f, k):
 
     >>> from sympy.concrete.gosper import gosper_sum
     >>> from sympy.functions import factorial
-    >>> from sympy.abc import i, n, k
+    >>> from sympy.abc import n, k
 
     >>> f = (4*k + 1)*factorial(k)/factorial(2*k + 1)
     >>> gosper_sum(f, (k, 0, n))
diff --git a/sympy/concrete/guess.py b/sympy/concrete/guess.py
index b6c2a847c37b..688fe04c4828 100644
--- a/sympy/concrete/guess.py
+++ b/sympy/concrete/guess.py
@@ -248,7 +248,7 @@ def guess_generating_function(v, X=Symbol('x'), types=['all'], maxsqrtn=2):
     >>> ggf([fibonacci(k) for k in range(5, 15)], types=['ogf'])
     {'ogf': (3*x + 5)/(-x**2 - x + 1)}
 
-    >>> from sympy import simplify, factorial
+    >>> from sympy import factorial
     >>> ggf([factorial(k) for k in range(12)], types=['ogf', 'egf', 'lgf'])
     {'egf': 1/(1 - x)}
 
diff --git a/sympy/concrete/products.py b/sympy/concrete/products.py
index b54e7dac0c90..a2c593270c7f 100644
--- a/sympy/concrete/products.py
+++ b/sympy/concrete/products.py
@@ -63,7 +63,7 @@ class Product(ExprWithIntLimits):
     ========
 
     >>> from sympy.abc import a, b, i, k, m, n, x
-    >>> from sympy import Product, factorial, oo
+    >>> from sympy import Product, oo
     >>> Product(k, (k, 1, m))
     Product(k, (k, 1, m))
     >>> Product(k, (k, 1, m)).doit()
@@ -428,7 +428,7 @@ def is_convergent(self):
         Examples
         ========
 
-        >>> from sympy import Interval, S, Product, Symbol, cos, pi, exp, oo
+        >>> from sympy import Product, Symbol, cos, pi, exp, oo
         >>> n = Symbol('n', integer=True)
         >>> Product(n/(n + 1), (n, 1, oo)).is_convergent()
         False
@@ -474,7 +474,7 @@ def reverse_order(expr, *indices):
         Examples
         ========
 
-        >>> from sympy import Product, simplify, RisingFactorial, gamma, Sum
+        >>> from sympy import Product, simplify, Sum
         >>> from sympy.abc import x, y, a, b, c, d
         >>> P = Product(x, (x, a, b))
         >>> Pr = P.reverse_order(x)
diff --git a/sympy/concrete/summations.py b/sympy/concrete/summations.py
index 27633a7c25e2..b4f1d413501c 100644
--- a/sympy/concrete/summations.py
+++ b/sympy/concrete/summations.py
@@ -628,7 +628,7 @@ def is_absolutely_convergent(self):
         Examples
         ========
 
-        >>> from sympy import Sum, Symbol, sin, oo
+        >>> from sympy import Sum, Symbol, oo
         >>> n = Symbol('n', integer=True)
         >>> Sum((-1)**n, (n, 1, oo)).is_absolutely_convergent()
         False
diff --git a/sympy/core/basic.py b/sympy/core/basic.py
index 53beba2396f7..8b024abea205 100644
--- a/sympy/core/basic.py
+++ b/sympy/core/basic.py
@@ -554,7 +554,7 @@ def as_dummy(self):
         ========
 
         >>> from sympy import Integral, Symbol
-        >>> from sympy.abc import x, y
+        >>> from sympy.abc import x
         >>> r = Symbol('r', real=True)
         >>> Integral(r, (r, x)).as_dummy()
         Integral(_0, (_0, x))
@@ -569,7 +569,6 @@ def as_dummy(self):
         Any object that has structurally bound variables should have
         a property, `bound_symbols` that returns those symbols
         appearing in the object.
-
         """
         from sympy.core.symbol import Dummy, Symbol
         def can(x):
diff --git a/sympy/core/cache.py b/sympy/core/cache.py
index 0bad6fbeefcd..7c5045093749 100644
--- a/sympy/core/cache.py
+++ b/sympy/core/cache.py
@@ -77,7 +77,7 @@ def __cacheit(maxsize):
            ...    return a+b
 
            >>> @cacheit
-           ... def f(a, b):
+           ... def f(a, b): # noqa: F811
            ...    return [a, b] # <-- WRONG, returns mutable object
 
            to force cacheit to check returned results mutability and consistency,
@@ -120,7 +120,7 @@ def __cacheit(maxsize):
            ...    return a+b
 
            >>> @cacheit
-           ... def f(a, b):
+           ... def f(a, b): # noqa: F811
            ...    return [a, b] # <-- WRONG, returns mutable object
 
            to force cacheit to check returned results mutability and consistency,
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index beec0ba0b1f6..b5490207ae9c 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -434,7 +434,7 @@ def is_number(self):
         Examples
         ========
 
-        >>> from sympy import log, Integral, cos, sin, pi
+        >>> from sympy import Integral, cos, sin, pi
         >>> from sympy.core.function import Function
         >>> from sympy.abc import x
         >>> f = Function('f')
@@ -2827,7 +2827,7 @@ def series(self, x=None, x0=0, n=6, dir="+", logx=None):
         Examples
         ========
 
-        >>> from sympy import cos, exp, tan, oo, series
+        >>> from sympy import cos, exp, tan
         >>> from sympy.abc import x, y
         >>> cos(x).series()
         1 - x**2/2 + x**4/24 + O(x**6)
@@ -3043,7 +3043,7 @@ def aseries(self, x=None, n=6, bound=0, hir=False):
         ========
 
         >>> from sympy import sin, exp
-        >>> from sympy.abc import x, y
+        >>> from sympy.abc import x
 
         >>> e = sin(1/x + exp(-x)) - sin(1/x)
 
@@ -3688,7 +3688,7 @@ def round(self, n=None):
         Examples
         ========
 
-        >>> from sympy import pi, E, I, S, Add, Mul, Number
+        >>> from sympy import pi, E, I, S, Number
         >>> pi.round()
         3
         >>> pi.round(2)
@@ -3921,7 +3921,7 @@ class UnevaluatedExpr(Expr):
     ========
 
     >>> from sympy import UnevaluatedExpr
-    >>> from sympy.abc import a, b, x, y
+    >>> from sympy.abc import x
     >>> x*(1/x)
     1
     >>> x*UnevaluatedExpr(1/x)
diff --git a/sympy/core/exprtools.py b/sympy/core/exprtools.py
index 820006e447a9..7ff6ef54596d 100644
--- a/sympy/core/exprtools.py
+++ b/sympy/core/exprtools.py
@@ -49,7 +49,7 @@ def _monotonic_sign(self):
     ========
 
     >>> from sympy.core.exprtools import _monotonic_sign as F
-    >>> from sympy import Dummy, S
+    >>> from sympy import Dummy
     >>> nn = Dummy(integer=True, nonnegative=True)
     >>> p = Dummy(integer=True, positive=True)
     >>> p2 = Dummy(integer=True, positive=True)
diff --git a/sympy/core/function.py b/sympy/core/function.py
index d76995116387..c02df2331536 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -216,7 +216,6 @@ def nargs(self):
         ========
 
         >>> from sympy.core.function import Function
-        >>> from sympy.abc import x, y
         >>> f = Function('f')
 
         If the function can take any number of arguments, the set of whole
@@ -1090,7 +1089,7 @@ class Derivative(Expr):
     automatically simplified in a fairly conservative fashion unless the
     keyword ``simplify`` is set to False.
 
-        >>> from sympy import cos, sin, sqrt, diff, Function, symbols
+        >>> from sympy import sqrt, diff, Function, symbols
         >>> from sympy.abc import x, y, z
         >>> f, g = symbols('f,g', cls=Function)
 
@@ -1524,7 +1523,7 @@ def _sort_variable_count(cls, vc):
         Examples
         ========
 
-        >>> from sympy import Derivative, Function, symbols, cos
+        >>> from sympy import Derivative, Function, symbols
         >>> vsort = Derivative._sort_variable_count
         >>> x, y, z = symbols('x y z')
         >>> f, g, h = symbols('f g h', cls=Function)
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index 9a01dd27e1d8..27b671e148c9 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -2583,7 +2583,7 @@ class Zero(IntegerConstant, metaclass=Singleton):
     Examples
     ========
 
-    >>> from sympy import S, Integer, zoo
+    >>> from sympy import S, Integer
     >>> Integer(0) is S.Zero
     True
     >>> 1/S.Zero
@@ -3290,7 +3290,7 @@ class ComplexInfinity(AtomicExpr, metaclass=Singleton):
     Examples
     ========
 
-    >>> from sympy import zoo, oo
+    >>> from sympy import zoo
     >>> zoo + 42
     zoo
     >>> 42/zoo
diff --git a/sympy/core/relational.py b/sympy/core/relational.py
index 33672037e46b..c65577b561b9 100644
--- a/sympy/core/relational.py
+++ b/sympy/core/relational.py
@@ -680,7 +680,7 @@ def as_poly(self, *gens, **kwargs):
         ========
 
         >>> from sympy import Eq
-        >>> from sympy.abc import x, y
+        >>> from sympy.abc import x
         >>> Eq(x**2, 1).as_poly(x)
         Poly(x**2 - 1, x, domain='ZZ')
         '''
diff --git a/sympy/core/singleton.py b/sympy/core/singleton.py
index 4ae3ec3e7131..1f9447021df5 100644
--- a/sympy/core/singleton.py
+++ b/sympy/core/singleton.py
@@ -46,7 +46,8 @@ class SingletonRegistry(Registry):
     When using ``is`` comparison, make sure the argument is sympified. For
     instance,
 
-    >>> 0 is S.Zero
+    >>> x = 0
+    >>> x is S.Zero
     False
 
     This problem is not an issue when using ``==``, which is recommended for
@@ -137,7 +138,6 @@ class is instantiated. Additionally, this instance can be accessed through
 
         >>> from sympy import S, Basic
         >>> from sympy.core.singleton import Singleton
-        >>> from sympy.core.compatibility import with_metaclass
         >>> class MySingleton(Basic, metaclass=Singleton):
         ...     pass
         >>> Basic() is Basic()
diff --git a/sympy/core/symbol.py b/sympy/core/symbol.py
index 1839feccc05e..214bf4ffe4a7 100644
--- a/sympy/core/symbol.py
+++ b/sympy/core/symbol.py
@@ -37,7 +37,7 @@ def _symbol(s, matching_symbol=None, **assumptions):
     Examples
     ========
 
-    >>> from sympy import Symbol, Dummy
+    >>> from sympy import Symbol
     >>> from sympy.core.symbol import _symbol
     >>> _symbol('y')
     y
@@ -116,7 +116,7 @@ def _uniquely_named_symbol(xname, exprs=(), compare=str, modify=None, **assumpti
     Examples
     ========
 
-    >>> from sympy.core.symbol import _uniquely_named_symbol as usym, Dummy
+    >>> from sympy.core.symbol import _uniquely_named_symbol as usym
     >>> from sympy.abc import x
     >>> usym('x', x)
     _x
@@ -720,17 +720,17 @@ def var(names, **args):
 
     >>> var('x')
     x
-    >>> x
+    >>> x # noqa: F821
     x
 
     >>> var('a,ab,abc')
     (a, ab, abc)
-    >>> abc
+    >>> abc # noqa: F821
     abc
 
     >>> var('x,y', real=True)
     (x, y)
-    >>> x.is_real and y.is_real
+    >>> x.is_real and y.is_real # noqa: F821
     True
 
     See :func:`symbols` documentation for more details on what kinds of
diff --git a/sympy/core/sympify.py b/sympy/core/sympify.py
index a75546544696..a3f1ce2851ba 100644
--- a/sympy/core/sympify.py
+++ b/sympy/core/sympify.py
@@ -448,7 +448,7 @@ def kernS(s):
     ========
 
     >>> from sympy.core.sympify import kernS
-    >>> from sympy.abc import x, y, z
+    >>> from sympy.abc import x, y
 
     The 2-arg Mul distributes a number (or minus sign) across the terms
     of an expression, but kernS will prevent that:
diff --git a/sympy/diffgeom/diffgeom.py b/sympy/diffgeom/diffgeom.py
index a4d1ed30b8e3..44edd4596808 100644
--- a/sympy/diffgeom/diffgeom.py
+++ b/sympy/diffgeom/diffgeom.py
@@ -482,8 +482,7 @@ class BaseScalarField(Expr):
     Define boilerplate Manifold, Patch and coordinate systems:
 
     >>> from sympy import symbols, sin, cos, pi, Function
-    >>> from sympy.diffgeom import (
-    ...        Manifold, Patch, CoordSystem, Point, BaseScalarField)
+    >>> from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField
     >>> r0, theta0 = symbols('r0, theta0')
     >>> m = Manifold('M', 2)
     >>> p = Patch('P', m)
@@ -574,7 +573,7 @@ class BaseVectorField(Expr):
 
     Use the predefined R2 manifold, setup some boilerplate.
 
-    >>> from sympy import symbols, pi, Function
+    >>> from sympy import symbols, Function
     >>> from sympy.diffgeom.rn import R2, R2_p, R2_r
     >>> from sympy.diffgeom import BaseVectorField
     >>> from sympy import pprint
@@ -681,7 +680,6 @@ class Commutator(Expr):
 
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.diffgeom import Commutator
-    >>> from sympy import pprint
     >>> from sympy.simplify import simplify
 
     Vector fields:
@@ -1204,7 +1202,7 @@ def intcurve_series(vector_field, param, start_point, n=6, coord_sys=None, coeff
     Use the predefined R2 manifold:
 
     >>> from sympy.abc import t, x, y
-    >>> from sympy.diffgeom.rn import R2, R2_p, R2_r
+    >>> from sympy.diffgeom.rn import R2_p, R2_r
     >>> from sympy.diffgeom import intcurve_series
 
     Specify a starting point and a vector field:
diff --git a/sympy/discrete/convolutions.py b/sympy/discrete/convolutions.py
index 9875c185a492..faf1c18b029a 100644
--- a/sympy/discrete/convolutions.py
+++ b/sympy/discrete/convolutions.py
@@ -296,7 +296,7 @@ def convolution_subset(a, b):
     Examples
     ========
 
-    >>> from sympy import symbols, S, I
+    >>> from sympy import symbols, S
     >>> from sympy.discrete.convolutions import convolution_subset
     >>> u, v, x, y, z = symbols('u v x y z')
 
diff --git a/sympy/functions/combinatorial/factorials.py b/sympy/functions/combinatorial/factorials.py
index 0078ce06f55b..240d4ddf2874 100644
--- a/sympy/functions/combinatorial/factorials.py
+++ b/sympy/functions/combinatorial/factorials.py
@@ -396,7 +396,8 @@ class factorial2(CombinatorialFunction):
     ========
 
     >>> from sympy import factorial2, var
-    >>> var('n')
+    >>> n = var('n')
+    >>> n
     n
     >>> factorial2(n + 1)
     factorial2(n + 1)
@@ -649,7 +650,7 @@ class FallingFactorial(CombinatorialFunction):
     Factorial Factorization and Symbolic Summation", Journal of
     Symbolic Computation, vol. 20, pp. 235-268, 1995.
 
-    >>> from sympy import ff, factorial, rf, gamma, polygamma, binomial, symbols, Poly
+    >>> from sympy import ff, factorial, rf, gamma, binomial, symbols, Poly
     >>> from sympy.abc import x, k
     >>> n, m = symbols('n m', integer=True)
     >>> ff(x, 0)
diff --git a/sympy/functions/combinatorial/numbers.py b/sympy/functions/combinatorial/numbers.py
index bde58ed85a74..2591fa00bf5a 100644
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -1081,8 +1081,8 @@ class catalan(Function):
     Examples
     ========
 
-    >>> from sympy import (Symbol, binomial, gamma, hyper, polygamma,
-    ...             catalan, diff, combsimp, Rational, I)
+    >>> from sympy import (Symbol, binomial, gamma, hyper, catalan,
+    ...                    diff, combsimp, Rational, I)
 
     >>> [catalan(i) for i in range(1,10)]
     [1, 2, 5, 14, 42, 132, 429, 1430, 4862]
diff --git a/sympy/functions/elementary/complexes.py b/sympy/functions/elementary/complexes.py
index 9defa46aa857..7fda45982bc6 100644
--- a/sympy/functions/elementary/complexes.py
+++ b/sympy/functions/elementary/complexes.py
@@ -29,7 +29,7 @@ class re(Function):
     ========
 
     >>> from sympy import re, im, I, E
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
     >>> re(2*E)
     2*E
     >>> re(2*I + 17)
diff --git a/sympy/functions/elementary/integers.py b/sympy/functions/elementary/integers.py
index 908bd4d712ab..960121c8aa73 100644
--- a/sympy/functions/elementary/integers.py
+++ b/sympy/functions/elementary/integers.py
@@ -362,7 +362,7 @@ class frac(Function):
     Examples
     ========
 
-    >>> from sympy import Symbol, frac, Rational, floor, ceiling, I
+    >>> from sympy import Symbol, frac, Rational, floor, I
     >>> frac(Rational(4, 3))
     1/3
     >>> frac(-Rational(4, 3))
diff --git a/sympy/functions/elementary/miscellaneous.py b/sympy/functions/elementary/miscellaneous.py
index a80432889715..3e958c55f8b1 100644
--- a/sympy/functions/elementary/miscellaneous.py
+++ b/sympy/functions/elementary/miscellaneous.py
@@ -250,7 +250,7 @@ def root(arg, n, k=0, evaluate=None):
     The following examples show the roots of unity for n
     equal 2, 3 and 4:
 
-    >>> from sympy import rootof, I
+    >>> from sympy import rootof
 
     >>> [rootof(x**2 - 1, i) for i in range(2)]
     [-1, 1]
@@ -329,8 +329,7 @@ def real_root(arg, n=None, evaluate=None):
     Examples
     ========
 
-    >>> from sympy import root, real_root, Rational
-    >>> from sympy.abc import x, n
+    >>> from sympy import root, real_root
 
     >>> real_root(-8, 3)
     -2
diff --git a/sympy/functions/elementary/piecewise.py b/sympy/functions/elementary/piecewise.py
index 9b5c1862a063..b7ccb6943ace 100644
--- a/sympy/functions/elementary/piecewise.py
+++ b/sympy/functions/elementary/piecewise.py
@@ -75,7 +75,7 @@ class Piecewise(Function):
     Examples
     ========
 
-    >>> from sympy import Piecewise, log, ITE, piecewise_fold
+    >>> from sympy import Piecewise, log, piecewise_fold
     >>> from sympy.abc import x, y
     >>> f = x**2
     >>> g = log(x)
diff --git a/sympy/functions/elementary/trigonometric.py b/sympy/functions/elementary/trigonometric.py
index ac8fddf94a2a..4a8adc4b1aed 100644
--- a/sympy/functions/elementary/trigonometric.py
+++ b/sympy/functions/elementary/trigonometric.py
@@ -137,7 +137,7 @@ def _pi_coeff(arg, cycles=1):
 
     >>> from sympy.functions.elementary.trigonometric import _pi_coeff as coeff
     >>> from sympy import pi, Dummy
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
     >>> coeff(3*x*pi)
     3*x
     >>> coeff(11*pi/7)
@@ -1879,7 +1879,7 @@ class sinc(Function):
     Examples
     ========
 
-    >>> from sympy import sinc, oo, jn, Product, Symbol
+    >>> from sympy import sinc, oo, jn
     >>> from sympy.abc import x
     >>> sinc(x)
     sinc(x)
@@ -2057,11 +2057,15 @@ class asin(InverseTrigonometricFunction):
     Examples
     ========
 
-    >>> from sympy import asin, oo, pi
+    >>> from sympy import asin, oo
     >>> asin(1)
     pi/2
     >>> asin(-1)
     -pi/2
+    >>> asin(-oo)
+    oo*I
+    >>> asin(oo)
+    -oo*I
 
     See Also
     ========
@@ -2227,7 +2231,7 @@ class acos(InverseTrigonometricFunction):
     Examples
     ========
 
-    >>> from sympy import acos, oo, pi
+    >>> from sympy import acos, oo
     >>> acos(1)
     0
     >>> acos(0)
@@ -2395,7 +2399,7 @@ class atan(InverseTrigonometricFunction):
     Examples
     ========
 
-    >>> from sympy import atan, oo, pi
+    >>> from sympy import atan, oo
     >>> atan(0)
     0
     >>> atan(1)
@@ -2758,11 +2762,15 @@ class asec(InverseTrigonometricFunction):
     Examples
     ========
 
-    >>> from sympy import asec, oo, pi
+    >>> from sympy import asec, oo
     >>> asec(1)
     0
     >>> asec(-1)
     pi
+    >>> asec(0)
+    zoo
+    >>> asec(-oo)
+    pi/2
 
     See Also
     ========
@@ -2876,11 +2884,17 @@ class acsc(InverseTrigonometricFunction):
     Examples
     ========
 
-    >>> from sympy import acsc, oo, pi
+    >>> from sympy import acsc, oo
     >>> acsc(1)
     pi/2
     >>> acsc(-1)
     -pi/2
+    >>> acsc(oo)
+    0
+    >>> acsc(-oo) == acsc(oo)
+    True
+    >>> acsc(0)
+    zoo
 
     See Also
     ========
diff --git a/sympy/functions/special/bessel.py b/sympy/functions/special/bessel.py
index db8f55d221f4..30bddffbbddc 100644
--- a/sympy/functions/special/bessel.py
+++ b/sympy/functions/special/bessel.py
@@ -660,7 +660,6 @@ class jn(SphericalBesselBase):
     ========
 
     >>> from sympy import Symbol, jn, sin, cos, expand_func, besselj, bessely
-    >>> from sympy import simplify
     >>> z = Symbol("z")
     >>> nu = Symbol("nu", integer=True)
     >>> print(expand_func(jn(0, z)))
@@ -1742,7 +1741,7 @@ class marcumq(Function):
     ========
 
     >>> from sympy import marcumq
-    >>> from sympy.abc import m, a, b, x
+    >>> from sympy.abc import m, a, b
     >>> marcumq(m, a, b)
     marcumq(m, a, b)
 
diff --git a/sympy/functions/special/delta_functions.py b/sympy/functions/special/delta_functions.py
index 23229a55ca20..be38c348d443 100644
--- a/sympy/functions/special/delta_functions.py
+++ b/sympy/functions/special/delta_functions.py
@@ -58,7 +58,7 @@ class DiracDelta(Function):
     Examples
     ========
 
-    >>> from sympy import DiracDelta, diff, pi, Piecewise
+    >>> from sympy import DiracDelta, diff, pi
     >>> from sympy.abc import x, y
 
     >>> DiracDelta(x)
@@ -160,7 +160,7 @@ class is about to be instantiated and it returns either some simplified
         Examples
         ========
 
-        >>> from sympy import DiracDelta, S, Subs
+        >>> from sympy import DiracDelta, S
         >>> from sympy.abc import x
 
         >>> DiracDelta(x)
@@ -335,7 +335,7 @@ def _eval_rewrite_as_Piecewise(self, *args, **kwargs):
         Examples
         ========
 
-        >>> from sympy import DiracDelta, Piecewise, Symbol, SingularityFunction
+        >>> from sympy import DiracDelta, Piecewise, Symbol
         >>> x = Symbol('x')
 
         >>> DiracDelta(x).rewrite(Piecewise)
@@ -549,7 +549,7 @@ def _eval_rewrite_as_Piecewise(self, arg, H0=None, **kwargs):
         Examples
         ========
 
-        >>> from sympy import Heaviside, Piecewise, Symbol, pprint
+        >>> from sympy import Heaviside, Piecewise, Symbol
         >>> x = Symbol('x')
 
         >>> Heaviside(x).rewrite(Piecewise)
diff --git a/sympy/functions/special/elliptic_integrals.py b/sympy/functions/special/elliptic_integrals.py
index e1a242a9180a..5e6854521925 100644
--- a/sympy/functions/special/elliptic_integrals.py
+++ b/sympy/functions/special/elliptic_integrals.py
@@ -32,7 +32,7 @@ class elliptic_k(Function):
     Examples
     ========
 
-    >>> from sympy import elliptic_k, I, pi
+    >>> from sympy import elliptic_k, I
     >>> from sympy.abc import m
     >>> elliptic_k(0)
     pi/2
@@ -128,7 +128,7 @@ class elliptic_f(Function):
     Examples
     ========
 
-    >>> from sympy import elliptic_f, I, O
+    >>> from sympy import elliptic_f, I
     >>> from sympy.abc import z, m
     >>> elliptic_f(z, m).series(z)
     z + z**5*(3*m**2/40 - m/30) + m*z**3/6 + O(z**6)
@@ -217,7 +217,7 @@ class elliptic_e(Function):
     Examples
     ========
 
-    >>> from sympy import elliptic_e, I, pi, O
+    >>> from sympy import elliptic_e, I
     >>> from sympy.abc import z, m
     >>> elliptic_e(z, m).series(z)
     z + z**5*(-m**2/40 + m/30) - m*z**3/6 + O(z**6)
@@ -338,7 +338,7 @@ class elliptic_pi(Function):
     Examples
     ========
 
-    >>> from sympy import elliptic_pi, I, pi, O, S
+    >>> from sympy import elliptic_pi, I
     >>> from sympy.abc import z, n, m
     >>> elliptic_pi(n, z, m).series(z, n=4)
     z + z**3*(m/6 + n/3) + O(z**4)
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index 5e8a51ad309f..699d39c53e3e 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -608,7 +608,7 @@ class erf2(Function):
     Examples
     ========
 
-    >>> from sympy import I, oo, erf2
+    >>> from sympy import oo, erf2
     >>> from sympy.abc import x, y
 
     Several special values are known:
@@ -748,7 +748,7 @@ class erfinv(Function):
     Examples
     ========
 
-    >>> from sympy import I, oo, erfinv
+    >>> from sympy import erfinv
     >>> from sympy.abc import x
 
     Several special values are known:
@@ -842,7 +842,7 @@ class erfcinv (Function):
     Examples
     ========
 
-    >>> from sympy import I, oo, erfcinv
+    >>> from sympy import erfcinv
     >>> from sympy.abc import x
 
     Several special values are known:
@@ -918,7 +918,7 @@ class erf2inv(Function):
     Examples
     ========
 
-    >>> from sympy import I, oo, erf2inv, erfinv, erfcinv
+    >>> from sympy import erf2inv, oo
     >>> from sympy.abc import x, y
 
     Several special values are known:
@@ -1064,7 +1064,7 @@ class Ei(Function):
     The exponential integral is related to many other special functions.
     For example:
 
-    >>> from sympy import uppergamma, expint, Shi
+    >>> from sympy import expint, Shi
     >>> Ei(x).rewrite(expint)
     -expint(1, x*exp_polar(I*pi)) - I*pi
     >>> Ei(x).rewrite(Shi)
@@ -1536,7 +1536,7 @@ class Li(Function):
     Examples
     ========
 
-    >>> from sympy import I, oo, Li
+    >>> from sympy import Li
     >>> from sympy.abc import z
 
     The following special value is known:
@@ -2213,7 +2213,7 @@ class fresnels(FresnelIntegral):
 
     Defining the Fresnel functions via an integral:
 
-    >>> from sympy import integrate, pi, sin, gamma, expand_func
+    >>> from sympy import integrate, pi, sin, expand_func
     >>> integrate(sin(pi*z**2/2), z)
     3*fresnels(z)*gamma(3/4)/(4*gamma(7/4))
     >>> expand_func(integrate(sin(pi*z**2/2), z))
@@ -2354,7 +2354,7 @@ class fresnelc(FresnelIntegral):
 
     Defining the Fresnel functions via an integral:
 
-    >>> from sympy import integrate, pi, cos, gamma, expand_func
+    >>> from sympy import integrate, pi, cos, expand_func
     >>> integrate(cos(pi*z**2/2), z)
     fresnelc(z)*gamma(1/4)/(4*gamma(5/4))
     >>> expand_func(integrate(cos(pi*z**2/2), z))
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index 4d92e373519e..d9793fd4a971 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -43,7 +43,7 @@ class gamma(Function):
     Examples
     ========
 
-    >>> from sympy import S, I, pi, oo, gamma
+    >>> from sympy import S, I, pi, gamma
     >>> from sympy.abc import x
 
     Several special values are known:
@@ -868,7 +868,7 @@ class loggamma(Function):
 
     For half-integral values:
 
-    >>> from sympy import S, pi
+    >>> from sympy import S
     >>> loggamma(S(5)/2)
     log(3*sqrt(pi)/4)
     >>> loggamma(n/2)
@@ -1241,7 +1241,7 @@ class multigamma(Function):
     Examples
     ========
 
-    >>> from sympy import S, I, pi, oo, gamma, multigamma
+    >>> from sympy import S, multigamma
     >>> from sympy import Symbol
     >>> x = Symbol('x')
     >>> p = Symbol('p', positive=True, integer=True)
diff --git a/sympy/functions/special/polynomials.py b/sympy/functions/special/polynomials.py
index 8d202083c8f4..161d84d5be32 100644
--- a/sympy/functions/special/polynomials.py
+++ b/sympy/functions/special/polynomials.py
@@ -441,7 +441,7 @@ class chebyshevt(OrthogonalPolynomial):
     Examples
     ========
 
-    >>> from sympy import chebyshevt, chebyshevu, diff
+    >>> from sympy import chebyshevt, diff
     >>> from sympy.abc import n,x
     >>> chebyshevt(0, x)
     1
@@ -553,7 +553,7 @@ class chebyshevu(OrthogonalPolynomial):
     Examples
     ========
 
-    >>> from sympy import chebyshevt, chebyshevu, diff
+    >>> from sympy import chebyshevu, diff
     >>> from sympy.abc import n,x
     >>> chebyshevu(0, x)
     1
@@ -1169,7 +1169,7 @@ class assoc_laguerre(OrthogonalPolynomial):
     Examples
     ========
 
-    >>> from sympy import laguerre, assoc_laguerre, diff
+    >>> from sympy import assoc_laguerre, diff
     >>> from sympy.abc import x, n, a
     >>> assoc_laguerre(0, a, x)
     1
diff --git a/sympy/geometry/curve.py b/sympy/geometry/curve.py
index c03a58575b2c..6aab8d79be35 100644
--- a/sympy/geometry/curve.py
+++ b/sympy/geometry/curve.py
@@ -52,7 +52,7 @@ class Curve(GeometrySet):
     Examples
     ========
 
-    >>> from sympy import sin, cos, Symbol, interpolate
+    >>> from sympy import sin, cos, interpolate
     >>> from sympy.abc import t, a
     >>> from sympy.geometry import Curve
     >>> C = Curve((sin(t), cos(t)), (t, 0, 2))
@@ -257,7 +257,6 @@ def length(self):
         ========
 
         >>> from sympy.geometry.curve import Curve
-        >>> from sympy import cos, sin
         >>> from sympy.abc import t
         >>> Curve((t, t), (t, 0, 1)).length
         sqrt(2)
@@ -290,7 +289,7 @@ def plot_interval(self, parameter='t'):
         ========
 
         >>> from sympy import Curve, sin
-        >>> from sympy.abc import x, t, s
+        >>> from sympy.abc import x, s
         >>> Curve((x, sin(x)), (x, 1, 2)).plot_interval()
         [t, 1, 2]
         >>> Curve((x, sin(x)), (x, 1, 2)).plot_interval(s)
@@ -337,7 +336,6 @@ def scale(self, x=1, y=1, pt=None):
         ========
 
         >>> from sympy.geometry.curve import Curve
-        >>> from sympy import pi
         >>> from sympy.abc import x
         >>> Curve((x, x), (x, 0, 1)).scale(2)
         Curve((2*x, x), (x, 0, 1))
@@ -355,7 +353,6 @@ def translate(self, x=0, y=0):
         ========
 
         >>> from sympy.geometry.curve import Curve
-        >>> from sympy import pi
         >>> from sympy.abc import x
         >>> Curve((x, x), (x, 0, 1)).translate(1, 2)
         Curve((x + 1, x + 2), (x, 0, 1))
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index 06ec340fb75a..b8bcdbaf793a 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -635,7 +635,7 @@ def intersection(self, o):
         Examples
         ========
 
-        >>> from sympy import Ellipse, Point, Line, sqrt
+        >>> from sympy import Ellipse, Point, Line
         >>> e = Ellipse(Point(0, 0), 5, 7)
         >>> e.intersection(Point(0, 0))
         []
@@ -879,7 +879,7 @@ def normal_lines(self, p, prec=None):
         Examples
         ========
 
-        >>> from sympy import Line, Point, Ellipse
+        >>> from sympy import Point, Ellipse
         >>> e = Ellipse((0, 0), 2, 3)
         >>> c = e.center
         >>> e.normal_lines(c + Point(1, 0))
@@ -1016,7 +1016,7 @@ def auxiliary_circle(self):
         Examples
         ========
 
-        >>> from sympy import Circle, Ellipse, Point, symbols
+        >>> from sympy import Ellipse, Point, symbols
         >>> c = Point(1, 2)
         >>> Ellipse(c, 8, 7).auxiliary_circle()
         Circle(Point2D(1, 2), 8)
@@ -1040,7 +1040,7 @@ def director_circle(self):
         Examples
         ========
 
-        >>> from sympy import Circle, Ellipse, Point, symbols
+        >>> from sympy import Ellipse, Point, symbols
         >>> c = Point(3,8)
         >>> Ellipse(c, 7, 9).director_circle()
         Circle(Point2D(3, 8), sqrt(130))
@@ -1094,7 +1094,7 @@ def random_point(self, seed=None):
         Examples
         ========
 
-        >>> from sympy import Point, Ellipse, Segment
+        >>> from sympy import Point, Ellipse
         >>> e1 = Ellipse(Point(0, 0), 3, 2)
         >>> e1.random_point() # gives some random point
         Point2D(...)
diff --git a/sympy/geometry/line.py b/sympy/geometry/line.py
index a2773b9a3e15..33805b37e757 100644
--- a/sympy/geometry/line.py
+++ b/sympy/geometry/line.py
@@ -172,7 +172,7 @@ def angle_between(l1, l2):
         Examples
         ========
 
-        >>> from sympy import Point, Line, pi
+        >>> from sympy import Line
         >>> e = Line((0, 0), (1, 0))
         >>> ne = Line((0, 0), (1, 1))
         >>> sw = Line((1, 1), (0, 0))
@@ -224,7 +224,7 @@ def smallest_angle_between(l1, l2):
         Examples
         ========
 
-        >>> from sympy import Point, Line, pi
+        >>> from sympy import Point, Line
         >>> p1, p2, p3 = Point(0, 0), Point(0, 4), Point(2, -2)
         >>> l1, l2 = Line(p1, p2), Line(p1, p3)
         >>> l1.smallest_angle_between(l2)
@@ -319,7 +319,7 @@ def are_concurrent(*lines):
         Examples
         ========
 
-        >>> from sympy import Point, Line, Line3D
+        >>> from sympy import Point, Line
         >>> p1, p2 = Point(0, 0), Point(3, 5)
         >>> p3, p4 = Point(-2, -2), Point(0, 2)
         >>> l1, l2, l3 = Line(p1, p2), Line(p1, p3), Line(p1, p4)
@@ -1969,7 +1969,6 @@ class Line2D(LinearEntity2D, Line):
     ========
 
     >>> from sympy import Point
-    >>> from sympy.abc import L
     >>> from sympy.geometry import Line, Segment
     >>> L = Line(Point(2,3), Point(3,5))
     >>> L
@@ -2496,7 +2495,7 @@ class Line3D(LinearEntity3D, Line):
     ========
 
     >>> from sympy import Point3D
-    >>> from sympy.geometry import Line3D, Segment3D
+    >>> from sympy.geometry import Line3D
     >>> L = Line3D(Point3D(2, 3, 4), Point3D(3, 5, 1))
     >>> L
     Line3D(Point3D(2, 3, 4), Point3D(3, 5, 1))
diff --git a/sympy/geometry/plane.py b/sympy/geometry/plane.py
index 804698250cff..cd4ee63aa3eb 100644
--- a/sympy/geometry/plane.py
+++ b/sympy/geometry/plane.py
@@ -40,7 +40,6 @@ class Plane(GeometryEntity):
     ========
 
     >>> from sympy import Plane, Point3D
-    >>> from sympy.abc import x
     >>> Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2))
     Plane(Point3D(1, 1, 1), (-1, 2, -1))
     >>> Plane((1, 1, 1), (2, 3, 4), (2, 2, 2))
@@ -267,7 +266,7 @@ def distance(self, o):
         Examples
         ========
 
-        >>> from sympy import Point, Point3D, Line, Line3D, Plane
+        >>> from sympy import Point3D, Line3D, Plane
         >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 1, 1))
         >>> b = Point3D(1, 2, 3)
         >>> a.distance(b)
@@ -362,7 +361,7 @@ def intersection(self, o):
         Examples
         ========
 
-        >>> from sympy import Point, Point3D, Line, Line3D, Plane
+        >>> from sympy import Point3D, Line3D, Plane
         >>> a = Plane(Point3D(1, 2, 3), normal_vector=(1, 1, 1))
         >>> b = Point3D(1, 2, 3)
         >>> a.intersection(b)
@@ -438,7 +437,7 @@ def is_coplanar(self, o):
         Examples
         ========
 
-        >>> from sympy import Plane, Point3D
+        >>> from sympy import Plane
         >>> o = (0, 0, 0)
         >>> p = Plane(o, (1, 1, 1))
         >>> p2 = Plane(o, (2, 2, 2))
@@ -621,7 +620,7 @@ def perpendicular_line(self, pt):
         Examples
         ========
 
-        >>> from sympy import Plane, Point3D, Line3D
+        >>> from sympy import Plane, Point3D
         >>> a = Plane(Point3D(1,4,6), normal_vector=(2, 4, 6))
         >>> a.perpendicular_line(Point3D(9, 8, 7))
         Line3D(Point3D(9, 8, 7), Point3D(11, 12, 13))
@@ -656,7 +655,7 @@ def perpendicular_plane(self, *pts):
         Examples
         ========
 
-        >>> from sympy import Plane, Point3D, Line3D
+        >>> from sympy import Plane, Point3D
         >>> a, b = Point3D(0, 0, 0), Point3D(0, 1, 0)
         >>> Z = (0, 0, 1)
         >>> p = Plane(a, normal_vector=Z)
@@ -722,7 +721,7 @@ def projection_line(self, line):
         Examples
         ========
 
-        >>> from sympy import Plane, Line, Line3D, Point, Point3D
+        >>> from sympy import Plane, Line, Line3D, Point3D
         >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 1, 1))
         >>> b = Line(Point3D(1, 1), Point3D(2, 2))
         >>> a.projection_line(b)
@@ -764,7 +763,7 @@ def projection(self, pt):
         Examples
         ========
 
-        >>> from sympy import Plane, Point, Point3D
+        >>> from sympy import Plane, Point3D
         >>> A = Plane(Point3D(1, 1, 2), normal_vector=(1, 1, 1))
 
         The projection is along the normal vector direction, not the z
@@ -828,7 +827,8 @@ def parameter_value(self, other, u, v=None):
         Examples
         ========
 
-        >>> from sympy import Plane, Point, pi
+        >>> from sympy import pi
+        >>> from sympy.geometry import Plane
         >>> from sympy.abc import t, u, v
         >>> p = Plane((2, 0, 0), (0, 0, 1), (0, 1, 0))
 
diff --git a/sympy/geometry/point.py b/sympy/geometry/point.py
index 51a5705123ff..a8e79b277792 100644
--- a/sympy/geometry/point.py
+++ b/sympy/geometry/point.py
@@ -1153,7 +1153,7 @@ def are_collinear(*points):
         Examples
         ========
 
-        >>> from sympy import Point3D, Matrix
+        >>> from sympy import Point3D
         >>> from sympy.abc import x
         >>> p1, p2 = Point3D(0, 0, 0), Point3D(1, 1, 1)
         >>> p3, p4, p5 = Point3D(2, 2, 2), Point3D(x, x, x), Point3D(1, 2, 6)
diff --git a/sympy/geometry/polygon.py b/sympy/geometry/polygon.py
index f3b283a40728..502c0fdbdadb 100644
--- a/sympy/geometry/polygon.py
+++ b/sympy/geometry/polygon.py
@@ -75,7 +75,7 @@ class Polygon(GeometrySet):
     Examples
     ========
 
-    >>> from sympy import Point, Polygon, pi
+    >>> from sympy import Polygon, pi
     >>> p1, p2, p3, p4, p5 = [(0, 0), (1, 0), (5, 1), (0, 1), (3, 0)]
     >>> Polygon(p1, p2, p3, p4)
     Polygon(Point2D(0, 0), Point2D(1, 0), Point2D(5, 1), Point2D(0, 1))
@@ -397,7 +397,7 @@ def second_moment_of_area(self, point=None):
         Examples
         ========
 
-        >>> from sympy import Point, Polygon, symbols
+        >>> from sympy import Polygon, symbols
         >>> a, b = symbols('a, b')
         >>> p1, p2, p3, p4, p5 = [(0, 0), (a, 0), (a, b), (0, b), (a/3, b/3)]
         >>> rectangle = Polygon(p1, p2, p3, p4)
@@ -718,7 +718,6 @@ def encloses_point(self, p):
         ========
 
         >>> from sympy import Polygon, Point
-        >>> from sympy.abc import t
         >>> p = Polygon((0, 0), (4, 0), (4, 4))
         >>> p.encloses_point(Point(2, 1))
         True
@@ -811,7 +810,7 @@ def arbitrary_point(self, parameter='t'):
         Examples
         ========
 
-        >>> from sympy import Polygon, S, Symbol
+        >>> from sympy import Polygon, Symbol
         >>> t = Symbol('t', real=True)
         >>> tri = Polygon((0, 0), (1, 0), (1, 1))
         >>> p = tri.arbitrary_point('t')
@@ -978,7 +977,7 @@ def cut_section(self, line):
         Examples
         ========
 
-        >>> from sympy import Point, Symbol, Polygon, Line
+        >>> from sympy import Polygon, Line
         >>> a, b = 20, 10
         >>> p1, p2, p3, p4 = [(0, b), (0, 0), (a, 0), (a, b)]
         >>> rectangle = Polygon(p1, p2, p3, p4)
@@ -1919,7 +1918,7 @@ def rotate(self, angle, pt=None):
         """Override GeometryEntity.rotate to first rotate the RegularPolygon
         about its center.
 
-        >>> from sympy import Point, RegularPolygon, Polygon, pi
+        >>> from sympy import Point, RegularPolygon, pi
         >>> t = RegularPolygon(Point(1, 0), 1, 3)
         >>> t.vertices[0] # vertex on x-axis
         Point2D(2, 0)
@@ -2625,7 +2624,7 @@ def exradii(self):
         The exradius touches the side of the triangle to which it is keyed, e.g.
         the exradius touching side 2 is:
 
-        >>> from sympy.geometry import Point, Triangle, Segment2D, Point2D
+        >>> from sympy.geometry import Point, Triangle
         >>> p1, p2, p3 = Point(0, 0), Point(6, 0), Point(0, 2)
         >>> t = Triangle(p1, p2, p3)
         >>> t.exradii[t.sides[2]]
diff --git a/sympy/geometry/util.py b/sympy/geometry/util.py
index 809a78dcc68f..92342457912a 100644
--- a/sympy/geometry/util.py
+++ b/sympy/geometry/util.py
@@ -276,7 +276,7 @@ def closest_points(*args):
     Examples
     ========
 
-    >>> from sympy.geometry import closest_points, Point2D, Triangle
+    >>> from sympy.geometry import closest_points, Triangle
     >>> Triangle(sss=(3, 4, 5)).args
     (Point2D(0, 0), Point2D(3, 0), Point2D(3, 4))
     >>> closest_points(*_)
@@ -374,7 +374,7 @@ def convex_hull(*args, polygon=True):
     Examples
     ========
 
-    >>> from sympy.geometry import Point, convex_hull
+    >>> from sympy.geometry import convex_hull
     >>> points = [(1, 1), (1, 2), (3, 1), (-5, 2), (15, 4)]
     >>> convex_hull(*points)
     Polygon(Point2D(-5, 2), Point2D(1, 1), Point2D(3, 1), Point2D(15, 4))
@@ -474,7 +474,7 @@ def farthest_points(*args):
     Examples
     ========
 
-    >>> from sympy.geometry import farthest_points, Point2D, Triangle
+    >>> from sympy.geometry import farthest_points, Triangle
     >>> Triangle(sss=(3, 4, 5)).args
     (Point2D(0, 0), Point2D(3, 0), Point2D(3, 4))
     >>> farthest_points(*_)
diff --git a/sympy/holonomic/holonomic.py b/sympy/holonomic/holonomic.py
index ba804252db94..e00c7064f1ba 100644
--- a/sympy/holonomic/holonomic.py
+++ b/sympy/holonomic/holonomic.py
@@ -150,7 +150,7 @@ class DifferentialOperator:
     ========
 
     >>> from sympy.holonomic.holonomic import DifferentialOperator, DifferentialOperators
-    >>> from sympy.polys.domains import ZZ, QQ
+    >>> from sympy.polys.domains import ZZ
     >>> from sympy import symbols
     >>> x = symbols('x')
     >>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
@@ -397,7 +397,7 @@ class HolonomicFunction:
     ========
 
     >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-    >>> from sympy.polys.domains import ZZ, QQ
+    >>> from sympy.polys.domains import QQ
     >>> from sympy import symbols, S
     >>> x = symbols('x')
     >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
@@ -714,7 +714,7 @@ def integrate(self, limits, initcond=False):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import QQ
         >>> from sympy import symbols
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
@@ -841,7 +841,7 @@ def diff(self, *args, **kwargs):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import ZZ
         >>> from sympy import symbols
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
@@ -1177,7 +1177,7 @@ def composition(self, expr, *args, **kwargs):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import QQ
         >>> from sympy import symbols
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
@@ -1257,7 +1257,7 @@ def to_sequence(self, lb=True):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import QQ
         >>> from sympy import symbols, S
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
@@ -1647,7 +1647,7 @@ def series(self, n=6, coefficient=False, order=True, _recur=None):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import QQ
         >>> from sympy import symbols
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
@@ -1781,7 +1781,7 @@ def evalf(self, points, method='RK4', h=0.05, derivatives=False):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import QQ
         >>> from sympy import symbols
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
@@ -1883,7 +1883,7 @@ def to_hyper(self, as_list=False, _recur=None):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import ZZ
         >>> from sympy import symbols
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
@@ -2056,7 +2056,7 @@ def to_expr(self):
         ========
 
         >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
-        >>> from sympy.polys.domains import ZZ, QQ
+        >>> from sympy.polys.domains import ZZ
         >>> from sympy import symbols, S
         >>> x = symbols('x')
         >>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
@@ -2077,7 +2077,7 @@ def change_ics(self, b, lenics=None):
         ========
 
         >>> from sympy.holonomic import expr_to_holonomic
-        >>> from sympy import symbols, sin, cos, exp
+        >>> from sympy import symbols, sin, exp
         >>> x = symbols('x')
 
         >>> expr_to_holonomic(sin(x)).change_ics(1)
@@ -2148,7 +2148,7 @@ def from_hyper(func, x0=0, evalf=False):
     Examples
     ========
 
-    >>> from sympy.holonomic.holonomic import from_hyper, DifferentialOperators
+    >>> from sympy.holonomic.holonomic import from_hyper
     >>> from sympy import symbols, hyper, S
     >>> x = symbols('x')
     >>> from_hyper(hyper([], [S(3)/2], x**2/4))
@@ -2222,7 +2222,7 @@ def from_meijerg(func, x0=0, evalf=False, initcond=True, domain=QQ):
     Examples
     ========
 
-    >>> from sympy.holonomic.holonomic import from_meijerg, DifferentialOperators
+    >>> from sympy.holonomic.holonomic import from_meijerg
     >>> from sympy import symbols, meijerg, S
     >>> x = symbols('x')
     >>> from_meijerg(meijerg(([], []), ([S(1)/2], [0]), x**2/4))
diff --git a/sympy/holonomic/recurrence.py b/sympy/holonomic/recurrence.py
index 368b5ad1abfb..7202898d98c4 100644
--- a/sympy/holonomic/recurrence.py
+++ b/sympy/holonomic/recurrence.py
@@ -109,7 +109,7 @@ class RecurrenceOperator:
     ========
 
     >>> from sympy.holonomic.recurrence import RecurrenceOperator, RecurrenceOperators
-    >>> from sympy.polys.domains import ZZ, QQ
+    >>> from sympy.polys.domains import ZZ
     >>> from sympy import symbols
     >>> n = symbols('n', integer=True)
     >>> R, Sn = RecurrenceOperators(ZZ.old_poly_ring(n),'Sn')
diff --git a/sympy/integrals/deltafunctions.py b/sympy/integrals/deltafunctions.py
index 1d9b348846e5..66a1d548e307 100644
--- a/sympy/integrals/deltafunctions.py
+++ b/sympy/integrals/deltafunctions.py
@@ -115,7 +115,7 @@ def deltaintegrate(f, x):
 
         >>> from sympy.abc import x, y, z
         >>> from sympy.integrals.deltafunctions import deltaintegrate
-        >>> from sympy import sin, cos, DiracDelta, Heaviside
+        >>> from sympy import sin, cos, DiracDelta
         >>> deltaintegrate(x*sin(x)*cos(x)*DiracDelta(x - 1), x)
         sin(1)*cos(1)*Heaviside(x - 1)
         >>> deltaintegrate(y**2*DiracDelta(x - z)*DiracDelta(y - z), y)
diff --git a/sympy/integrals/heurisch.py b/sympy/integrals/heurisch.py
index 676e4019780c..97b2eb9c6c6b 100644
--- a/sympy/integrals/heurisch.py
+++ b/sympy/integrals/heurisch.py
@@ -47,7 +47,7 @@ def components(f, x):
     minimal, positive exponents.
 
     >>> from sympy import cos, sin
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
     >>> from sympy.integrals.heurisch import components
 
     >>> components(sin(x)*cos(x)**2, x)
diff --git a/sympy/integrals/integrals.py b/sympy/integrals/integrals.py
index b079c5c3fbb3..57ee907a1cb1 100644
--- a/sympy/integrals/integrals.py
+++ b/sympy/integrals/integrals.py
@@ -194,8 +194,8 @@ def transform(self, x, u):
         Examples
         ========
 
-        >>> from sympy.abc import a, b, c, d, x, u, y
-        >>> from sympy import Integral, S, cos, sqrt
+        >>> from sympy.abc import a, x, u
+        >>> from sympy import Integral, cos, sqrt
 
         >>> i = Integral(x*cos(x**2 - 1), (x, 0, 1))
 
@@ -376,7 +376,7 @@ def doit(self, **hints):
         Examples
         ========
 
-        >>> from sympy import Integral, Piecewise, S
+        >>> from sympy import Piecewise, S
         >>> from sympy.abc import x, t
         >>> p = x**2 + Piecewise((0, x/t < 0), (1, True))
         >>> p.integrate((t, S(4)/5, 1), (x, -1, 1))
@@ -1333,10 +1333,9 @@ def principal_value(self, **kwargs):
         Examples
         ========
 
-        >>> from sympy import Dummy, symbols, integrate, limit, oo
+        >>> from sympy import oo
         >>> from sympy.integrals.integrals import Integral
-        >>> from sympy.calculus.singularities import singularities
-        >>> x = symbols('x')
+        >>> from sympy.abc import x
         >>> Integral(x+1, (x, -oo, oo)).principal_value()
         oo
         >>> f = 1 / (x**3)
diff --git a/sympy/integrals/intpoly.py b/sympy/integrals/intpoly.py
index a1ee0957174a..99c05d1259a4 100644
--- a/sympy/integrals/intpoly.py
+++ b/sympy/integrals/intpoly.py
@@ -164,7 +164,6 @@ def main_integrate3d(expr, facets, vertices, hp_params, max_degree=None):
     Examples
     ========
 
-    >>> from sympy.abc import x, y
     >>> from sympy.integrals.intpoly import main_integrate3d, \
     hyperplane_parameters
     >>> cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\
@@ -326,7 +325,6 @@ def polygon_integrate(facet, hp_param, index, facets, vertices, expr, degree):
     Examples
     ========
 
-    >>> from sympy.abc import x, y
     >>> from sympy.integrals.intpoly import polygon_integrate
     >>> cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\
                  (5, 0, 5), (5, 5, 0), (5, 5, 5)],\
@@ -561,8 +559,8 @@ def integration_reduction_dynamic(facets, index, a, b, expr, degree, dims,
     ========
 
     >>> from sympy.abc import x, y
-    >>> from sympy.integrals.intpoly import integration_reduction_dynamic,\
-    hyperplane_parameters, gradient_terms
+    >>> from sympy.integrals.intpoly import (integration_reduction_dynamic, \
+            hyperplane_parameters)
     >>> from sympy.geometry.point import Point
     >>> from sympy.geometry.polygon import Polygon
     >>> triangle = Polygon(Point(0, 3), Point(5, 3), Point(1, 1))
@@ -629,7 +627,6 @@ def left_integral3D(facets, index, expr, vertices, hp_param, degree):
     hp_param : The hyperplane parameters of the face
     degree : Degree of the expr
 
-    >>> from sympy.abc import x, y
     >>> from sympy.integrals.intpoly import left_integral3D
     >>> cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\
                  (5, 0, 5), (5, 5, 0), (5, 5, 5)],\
@@ -664,7 +661,6 @@ def gradient_terms(binomial_power=0, no_of_gens=2):
     Examples
     ========
 
-    >>> from sympy.abc import x, y
     >>> from sympy.integrals.intpoly import gradient_terms
     >>> gradient_terms(2)
     [[1, 0, 0, 0], [y, 0, 1, 0], [y**2, 0, 2, 0], [x, 1, 0, 0],
diff --git a/sympy/integrals/manualintegrate.py b/sympy/integrals/manualintegrate.py
index 6aa5abd999f5..d1d61fb96250 100644
--- a/sympy/integrals/manualintegrate.py
+++ b/sympy/integrals/manualintegrate.py
@@ -1213,7 +1213,7 @@ def integral_steps(integrand, symbol, **options):
     Examples
     ========
 
-    >>> from sympy import exp, sin, cos
+    >>> from sympy import exp, sin
     >>> from sympy.integrals.manualintegrate import integral_steps
     >>> from sympy.abc import x
     >>> print(repr(integral_steps(exp(x) / (1 + exp(2 * x)), x))) \
diff --git a/sympy/integrals/meijerint.py b/sympy/integrals/meijerint.py
index 1f1c1a2a6e08..1871f0743f9c 100644
--- a/sympy/integrals/meijerint.py
+++ b/sympy/integrals/meijerint.py
@@ -381,7 +381,7 @@ def _find_splitting_points(expr, x):
 
     >>> from sympy.integrals.meijerint import _find_splitting_points as fsp
     >>> from sympy import sin
-    >>> from sympy.abc import a, x
+    >>> from sympy.abc import x
     >>> fsp(x, x)
     {0}
     >>> fsp((x-1)**3, x)
@@ -587,7 +587,7 @@ def _condsimp(cond):
     added as need arises rather than following any logical pattern.
 
     >>> from sympy.integrals.meijerint import _condsimp as simp
-    >>> from sympy import Or, Eq, unbranched_argument as arg, And
+    >>> from sympy import Or, Eq, And
     >>> from sympy.abc import x, y, z
     >>> simp(Or(x < y, z, Eq(x, y)))
     z | (x <= y)
diff --git a/sympy/integrals/quadrature.py b/sympy/integrals/quadrature.py
index 27a845bb0a8c..d8866d02ba26 100644
--- a/sympy/integrals/quadrature.py
+++ b/sympy/integrals/quadrature.py
@@ -327,7 +327,6 @@ def gauss_chebyshev_t(n, n_digits):
     Examples
     ========
 
-    >>> from sympy import S
     >>> from sympy.integrals.quadrature import gauss_chebyshev_t
     >>> x, w = gauss_chebyshev_t(3, 5)
     >>> x
@@ -394,7 +393,6 @@ def gauss_chebyshev_u(n, n_digits):
     Examples
     ========
 
-    >>> from sympy import S
     >>> from sympy.integrals.quadrature import gauss_chebyshev_u
     >>> x, w = gauss_chebyshev_u(3, 5)
     >>> x
diff --git a/sympy/integrals/rationaltools.py b/sympy/integrals/rationaltools.py
index 9d7d9b48483f..54bda63ccd3d 100644
--- a/sympy/integrals/rationaltools.py
+++ b/sympy/integrals/rationaltools.py
@@ -320,7 +320,7 @@ def log_to_real(h, q, x, t):
 
         >>> from sympy.integrals.rationaltools import log_to_real
         >>> from sympy.abc import x, y
-        >>> from sympy import Poly, sqrt, S
+        >>> from sympy import Poly, S
         >>> log_to_real(Poly(x + 3*y/2 + S(1)/2, x, domain='QQ[y]'),
         ... Poly(3*y**2 + 1, y, domain='ZZ'), x, y)
         2*sqrt(3)*atan(2*sqrt(3)*x/3 + sqrt(3)/3)/3
diff --git a/sympy/integrals/transforms.py b/sympy/integrals/transforms.py
index 0154c44f1d1b..754fd5808f11 100644
--- a/sympy/integrals/transforms.py
+++ b/sympy/integrals/transforms.py
@@ -1609,7 +1609,7 @@ def inverse_sine_transform(F, k, x, **hints):
     :func:`sympy.integrals.transforms.IntegralTransform.doit`.
     Note that for this transform, by default ``noconds=True``.
 
-    >>> from sympy import inverse_sine_transform, exp, sqrt, gamma, pi
+    >>> from sympy import inverse_sine_transform, exp, sqrt, gamma
     >>> from sympy.abc import x, k, a
     >>> inverse_sine_transform(2**((1-2*a)/2)*k**(a - 1)*
     ...     gamma(-a/2 + 1)/gamma((a+1)/2), k, x)
@@ -1716,7 +1716,7 @@ def inverse_cosine_transform(F, k, x, **hints):
     :func:`sympy.integrals.transforms.IntegralTransform.doit`.
     Note that for this transform, by default ``noconds=True``.
 
-    >>> from sympy import inverse_cosine_transform, exp, sqrt, pi
+    >>> from sympy import inverse_cosine_transform, sqrt, pi
     >>> from sympy.abc import x, k, a
     >>> inverse_cosine_transform(sqrt(2)*a/(sqrt(pi)*(a**2 + k**2)), k, x)
     exp(-a*x)
@@ -1816,7 +1816,7 @@ def hankel_transform(f, r, k, nu, **hints):
     Note that for this transform, by default ``noconds=True``.
 
     >>> from sympy import hankel_transform, inverse_hankel_transform
-    >>> from sympy import gamma, exp, sinh, cosh
+    >>> from sympy import exp
     >>> from sympy.abc import r, k, m, nu, a
 
     >>> ht = hankel_transform(1/r**m, r, k, nu)
@@ -1871,8 +1871,8 @@ def inverse_hankel_transform(F, k, r, nu, **hints):
     :func:`sympy.integrals.transforms.IntegralTransform.doit`.
     Note that for this transform, by default ``noconds=True``.
 
-    >>> from sympy import hankel_transform, inverse_hankel_transform, gamma
-    >>> from sympy import gamma, exp, sinh, cosh
+    >>> from sympy import hankel_transform, inverse_hankel_transform
+    >>> from sympy import exp
     >>> from sympy.abc import r, k, m, nu, a
 
     >>> ht = hankel_transform(1/r**m, r, k, nu)
diff --git a/sympy/integrals/trigonometry.py b/sympy/integrals/trigonometry.py
index adcad981c7ac..009a7c4280b4 100644
--- a/sympy/integrals/trigonometry.py
+++ b/sympy/integrals/trigonometry.py
@@ -28,7 +28,7 @@ def _pat_sincos(x):
 def trigintegrate(f, x, conds='piecewise'):
     """Integrate f = Mul(trig) over x
 
-       >>> from sympy import Symbol, sin, cos, tan, sec, csc, cot
+       >>> from sympy import sin, cos, tan, sec
        >>> from sympy.integrals.trigonometry import trigintegrate
        >>> from sympy.abc import x
 
diff --git a/sympy/interactive/session.py b/sympy/interactive/session.py
index a756b42bebd6..af4edb3e3dbc 100644
--- a/sympy/interactive/session.py
+++ b/sympy/interactive/session.py
@@ -95,8 +95,8 @@ def int_to_Integer(s):
     ========
 
     >>> from __future__ import division
+    >>> from sympy import Integer # noqa: F401
     >>> from sympy.interactive.session import int_to_Integer
-    >>> from sympy import Integer
     >>> s = '1.2 + 1/2 - 0x12 + a1'
     >>> int_to_Integer(s)
     '1.2 +Integer (1 )/Integer (2 )-Integer (0x12 )+a1 '
diff --git a/sympy/logic/algorithms/dpll.py b/sympy/logic/algorithms/dpll.py
index e43441daac43..61f846186ce2 100644
--- a/sympy/logic/algorithms/dpll.py
+++ b/sympy/logic/algorithms/dpll.py
@@ -181,7 +181,6 @@ def unit_propagate(clauses, symbol):
 
     Arguments are expected to be in CNF.
 
-    >>> from sympy import symbols
     >>> from sympy.abc import A, B, D
     >>> from sympy.logic.algorithms.dpll import unit_propagate
     >>> unit_propagate([A | B, D | ~B, B], B)
@@ -223,7 +222,6 @@ def find_pure_symbol(symbols, unknown_clauses):
     Find a symbol and its value if it appears only as a positive literal
     (or only as a negative) in clauses.
 
-    >>> from sympy import symbols
     >>> from sympy.abc import A, B, D
     >>> from sympy.logic.algorithms.dpll import find_pure_symbol
     >>> find_pure_symbol([A, B, D], [A|~B,~B|~D,D|A])
@@ -269,7 +267,6 @@ def find_unit_clause(clauses, model):
     """
     A unit clause has only 1 variable that is not bound in the model.
 
-    >>> from sympy import symbols
     >>> from sympy.abc import A, B, D
     >>> from sympy.logic.algorithms.dpll import find_unit_clause
     >>> find_unit_clause([A | B | D, B | ~D, A | ~B], {A:True})
diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py
index 979d61a68f69..60cee03c90eb 100644
--- a/sympy/logic/boolalg.py
+++ b/sympy/logic/boolalg.py
@@ -27,6 +27,8 @@ def as_Boolean(e):
     >>> from sympy import true, false, nan
     >>> from sympy.logic.boolalg import as_Boolean
     >>> from sympy.abc import x
+    >>> as_Boolean(0) is false
+    True
     >>> as_Boolean(1) is true
     True
     >>> as_Boolean(x)
@@ -35,6 +37,10 @@ def as_Boolean(e):
     Traceback (most recent call last):
     ...
     TypeError: expecting bool or Boolean, not `2`.
+    >>> as_Boolean(nan)
+    Traceback (most recent call last):
+    ...
+    TypeError: expecting bool or Boolean, not `nan`.
 
     """
     from sympy.core.symbol import Symbol
@@ -651,7 +657,6 @@ class And(LatticeOp, BooleanFunction):
     Examples
     ========
 
-    >>> from sympy.core import symbols
     >>> from sympy.abc import x, y
     >>> from sympy.logic.boolalg import And
     >>> x & y
@@ -798,7 +803,6 @@ class Or(LatticeOp, BooleanFunction):
     Examples
     ========
 
-    >>> from sympy.core import symbols
     >>> from sympy.abc import x, y
     >>> from sympy.logic.boolalg import Or
     >>> x | y
@@ -1323,7 +1327,7 @@ class Equivalent(BooleanFunction):
     ========
 
     >>> from sympy.logic.boolalg import Equivalent, And
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
     >>> Equivalent(False, False, False)
     True
     >>> Equivalent(True, False, False)
diff --git a/sympy/logic/inference.py b/sympy/logic/inference.py
index d5aec52eb3ac..e0628f3ec3ab 100644
--- a/sympy/logic/inference.py
+++ b/sympy/logic/inference.py
@@ -144,7 +144,7 @@ def pl_true(expr, model={}, deep=False):
     Examples
     ========
 
-    >>> from sympy.abc import A, B, C
+    >>> from sympy.abc import A, B
     >>> from sympy.logic.inference import pl_true
     >>> pl_true( A & B, {A: True, B: True})
     True
diff --git a/sympy/matrices/common.py b/sympy/matrices/common.py
index e0c7701621ed..a21d4f9d6d9e 100644
--- a/sympy/matrices/common.py
+++ b/sympy/matrices/common.py
@@ -522,7 +522,7 @@ def diagonal(self, k=0):
         Examples
         ========
 
-        >>> from sympy import Matrix, SparseMatrix
+        >>> from sympy import Matrix
         >>> m = Matrix(3, 3, lambda i, j: j - i); m
         Matrix([
         [ 0,  1, 2],
@@ -607,7 +607,7 @@ def todok(self):
         Examples
         ========
 
-        >>> from sympy import Matrix, SparseMatrix
+        >>> from sympy import Matrix
         >>> M = Matrix.eye(3)
         >>> M.todok()
         {(0, 0): 1, (1, 1): 1, (2, 2): 1}
@@ -1246,6 +1246,12 @@ def atoms(self, *types):
         Matrix([[x]])
         >>> _.atoms()
         {x}
+        >>> Matrix([[x, y], [y, x]])
+        Matrix([
+        [x, y],
+        [y, x]])
+        >>> _.atoms()
+        {x, y}
         """
 
         types = tuple(t if isinstance(t, type) else type(t) for t in types)
diff --git a/sympy/matrices/dense.py b/sympy/matrices/dense.py
index e16ccd9379d3..21d836b99430 100644
--- a/sympy/matrices/dense.py
+++ b/sympy/matrices/dense.py
@@ -225,7 +225,6 @@ def equals(self, other, failing_expression=False):
 
         >>> from sympy.matrices import Matrix
         >>> from sympy.abc import x
-        >>> from sympy import cos
         >>> A = Matrix([x*(x - 1), 0])
         >>> B = Matrix([x**2 - x, 0])
         >>> A == B
diff --git a/sympy/matrices/determinant.py b/sympy/matrices/determinant.py
index 647541afef9b..b2cccf0d2ae5 100644
--- a/sympy/matrices/determinant.py
+++ b/sympy/matrices/determinant.py
@@ -542,10 +542,15 @@ def _det(M, method="bareiss", iszerofunc=None):
     Examples
     ========
 
-    >>> from sympy import Matrix, MatrixSymbol, eye, det
+    >>> from sympy import Matrix, eye, det
+    >>> I3 = eye(3)
+    >>> det(I3)
+    1
     >>> M = Matrix([[1, 2], [3, 4]])
-    >>> M.det()
+    >>> det(M)
     -2
+    >>> det(M) == M.det()
+    True
     """
 
     # sanitize `method`
diff --git a/sympy/matrices/expressions/blockmatrix.py b/sympy/matrices/expressions/blockmatrix.py
index 22f0ba693250..4e1bb444e01f 100644
--- a/sympy/matrices/expressions/blockmatrix.py
+++ b/sympy/matrices/expressions/blockmatrix.py
@@ -219,7 +219,7 @@ def transpose(self):
         ========
 
         >>> from sympy import MatrixSymbol, BlockMatrix, ZeroMatrix
-        >>> from sympy.abc import l, m, n
+        >>> from sympy.abc import m, n
         >>> X = MatrixSymbol('X', n, n)
         >>> Y = MatrixSymbol('Y', m ,m)
         >>> Z = MatrixSymbol('Z', n, m)
@@ -287,7 +287,7 @@ class BlockDiagMatrix(BlockMatrix):
     Examples
     ========
 
-    >>> from sympy import MatrixSymbol, BlockDiagMatrix, symbols, Identity
+    >>> from sympy import MatrixSymbol, BlockDiagMatrix, symbols
     >>> n, m, l = symbols('n m l')
     >>> X = MatrixSymbol('X', n, n)
     >>> Y = MatrixSymbol('Y', m ,m)
@@ -408,7 +408,7 @@ def block_collapse(expr):
     """Evaluates a block matrix expression
 
     >>> from sympy import MatrixSymbol, BlockMatrix, symbols, \
-                          Identity, Matrix, ZeroMatrix, block_collapse
+                          Identity, ZeroMatrix, block_collapse
     >>> n,m,l = symbols('n m l')
     >>> X = MatrixSymbol('X', n, n)
     >>> Y = MatrixSymbol('Y', m ,m)
diff --git a/sympy/matrices/expressions/funcmatrix.py b/sympy/matrices/expressions/funcmatrix.py
index 09059d970dbb..5d08fd66cbc6 100644
--- a/sympy/matrices/expressions/funcmatrix.py
+++ b/sympy/matrices/expressions/funcmatrix.py
@@ -31,7 +31,7 @@ class FunctionMatrix(MatrixExpr):
 
     Creating a ``FunctionMatrix`` from ``Lambda``:
 
-    >>> from sympy import FunctionMatrix, symbols, Lambda, MatPow, Matrix
+    >>> from sympy import FunctionMatrix, symbols, Lambda, MatPow
     >>> i, j, n, m = symbols('i,j,n,m')
     >>> FunctionMatrix(n, m, Lambda((i, j), i + j))
     FunctionMatrix(n, m, Lambda((i, j), i + j))
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 4b2bb1d4eb19..31f1226cce69 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -414,7 +414,7 @@ def from_index_summation(expr, first_index=None, last_index=None, dimensions=Non
         Examples
         ========
 
-        >>> from sympy import MatrixSymbol, MatrixExpr, Sum, Symbol
+        >>> from sympy import MatrixSymbol, MatrixExpr, Sum
         >>> from sympy.abc import i, j, k, l, N
         >>> A = MatrixSymbol("A", N, N)
         >>> B = MatrixSymbol("B", N, N)
diff --git a/sympy/matrices/expressions/matmul.py b/sympy/matrices/expressions/matmul.py
index fc19c1fdc643..af9e09695852 100644
--- a/sympy/matrices/expressions/matmul.py
+++ b/sympy/matrices/expressions/matmul.py
@@ -233,7 +233,7 @@ def any_zeros(mul):
 def merge_explicit(matmul):
     """ Merge explicit MatrixBase arguments
 
-    >>> from sympy import MatrixSymbol, eye, Matrix, MatMul, pprint
+    >>> from sympy import MatrixSymbol, Matrix, MatMul, pprint
     >>> from sympy.matrices.expressions.matmul import merge_explicit
     >>> A = MatrixSymbol('A', 2, 2)
     >>> B = Matrix([[1, 1], [1, 1]])
diff --git a/sympy/matrices/expressions/trace.py b/sympy/matrices/expressions/trace.py
index edcfd9ab7ea8..1bfeb44660fa 100644
--- a/sympy/matrices/expressions/trace.py
+++ b/sympy/matrices/expressions/trace.py
@@ -15,6 +15,13 @@ class Trace(Expr):
     >>> A = MatrixSymbol('A', 3, 3)
     >>> Trace(A)
     Trace(A)
+    >>> Trace(eye(3))
+    Trace(Matrix([
+    [1, 0, 0],
+    [0, 1, 0],
+    [0, 0, 1]]))
+    >>> Trace(eye(3)).simplify()
+    3
     """
     is_Trace = True
     is_commutative = True
@@ -117,7 +124,7 @@ def trace(expr):
     Examples
     ========
 
-    >>> from sympy import trace, Symbol, MatrixSymbol, pprint, eye
+    >>> from sympy import trace, Symbol, MatrixSymbol, eye
     >>> n = Symbol('n')
     >>> X = MatrixSymbol('X', n, n)  # A square matrix
     >>> trace(2*X)
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index ae6f9e10d9d7..1e97c753d2d1 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -1517,7 +1517,7 @@ def analytic_func(self, f, x):
         Examples
         ========
 
-        >>> from sympy import Symbol, Matrix, exp, S, log
+        >>> from sympy import Symbol, Matrix, S, log
 
         >>> x = Symbol('x')
         >>> m = Matrix([[S(5)/4, S(3)/4], [S(3)/4, S(5)/4]])
diff --git a/sympy/multipledispatch/core.py b/sympy/multipledispatch/core.py
index f53d4e1516ee..5bebd2a11c01 100644
--- a/sympy/multipledispatch/core.py
+++ b/sympy/multipledispatch/core.py
@@ -28,7 +28,7 @@ def dispatch(*types, **kwargs):
     ...     return x + 1
 
     >>> @dispatch(float)
-    ... def f(x):
+    ... def f(x): # noqa: F811
     ...     return x - 1
 
     >>> f(3)
@@ -50,7 +50,7 @@ def dispatch(*types, **kwargs):
     ...     def __init__(self, data):
     ...         self.data = data
     ...     @dispatch(int)
-    ...     def __init__(self, datum):
+    ...     def __init__(self, datum): # noqa: F811
     ...         self.data = [datum]
     """
     namespace = kwargs.get('namespace', global_namespace)
diff --git a/sympy/multipledispatch/dispatcher.py b/sympy/multipledispatch/dispatcher.py
index 239a25937984..9c5f38651183 100644
--- a/sympy/multipledispatch/dispatcher.py
+++ b/sympy/multipledispatch/dispatcher.py
@@ -57,7 +57,7 @@ class Dispatcher:
     ...     return x + 1
 
     >>> @dispatch(float)
-    ... def f(x):
+    ... def f(x): # noqa: F811
     ...     return x - 1
 
     >>> f(3)
diff --git a/sympy/ntheory/continued_fraction.py b/sympy/ntheory/continued_fraction.py
index ba3ed62a1d7c..8ea00d5eaf0e 100644
--- a/sympy/ntheory/continued_fraction.py
+++ b/sympy/ntheory/continued_fraction.py
@@ -314,7 +314,7 @@ def continued_fraction_convergents(cf):
     Examples
     ========
 
-    >>> from sympy.core import Rational, pi
+    >>> from sympy.core import pi
     >>> from sympy import S
     >>> from sympy.ntheory.continued_fraction import \
             continued_fraction_convergents, continued_fraction_iterator
diff --git a/sympy/ntheory/digits.py b/sympy/ntheory/digits.py
index 57c652dee1a1..1453dcee53ea 100644
--- a/sympy/ntheory/digits.py
+++ b/sympy/ntheory/digits.py
@@ -78,7 +78,7 @@ def count_digits(n, b=10):
     Examples
     ========
 
-    >>> from sympy.ntheory import count_digits, digits
+    >>> from sympy.ntheory import count_digits
 
     >>> count_digits(1111339)
     {1: 4, 3: 2, 9: 1}
diff --git a/sympy/ntheory/elliptic_curve.py b/sympy/ntheory/elliptic_curve.py
index f86055379482..3b51952ae0bb 100644
--- a/sympy/ntheory/elliptic_curve.py
+++ b/sympy/ntheory/elliptic_curve.py
@@ -241,7 +241,6 @@ def order(self):
         Examples
         ========
 
-        >>> from sympy.polys.domains import FF
         >>> from sympy.ntheory.elliptic_curve import EllipticCurve
         >>> e2 = EllipticCurve(1, 0, modulus=19)
         >>> e2.order
diff --git a/sympy/ntheory/factor_.py b/sympy/ntheory/factor_.py
index 9b02fa7bee71..60d4c611f524 100644
--- a/sympy/ntheory/factor_.py
+++ b/sympy/ntheory/factor_.py
@@ -730,7 +730,7 @@ def pollard_pm1(n, B=10, a=2, retries=0, seed=1234):
 
     With the default smoothness bound, this number can't be cracked:
 
-        >>> from sympy.ntheory import pollard_pm1, primefactors
+        >>> from sympy.ntheory import pollard_pm1
         >>> pollard_pm1(21477639576571)
 
     Increasing the smoothness bound helps:
@@ -740,7 +740,6 @@ def pollard_pm1(n, B=10, a=2, retries=0, seed=1234):
 
     Looking at the smoothness of the factors of this number we find:
 
-        >>> from sympy.utilities import flatten
         >>> from sympy.ntheory.factor_ import smoothness_p, factorint
         >>> print(smoothness_p(21477639576571, visual=1))
         p**i=4410317**1 has p-1 B=1787, B-pow=1787
@@ -1048,7 +1047,7 @@ def factorint(n, limit=None, use_trial=True, use_rho=True, use_pm1=True,
     You can easily switch between the two forms by sending them back to
     factorint:
 
-    >>> from sympy import Mul, Pow
+    >>> from sympy import Mul
     >>> regular = factorint(1764); regular
     {2: 2, 3: 2, 7: 2}
     >>> pprint(factorint(regular))
diff --git a/sympy/ntheory/generate.py b/sympy/ntheory/generate.py
index 5ecea6c7e0e1..28a943d857c8 100644
--- a/sympy/ntheory/generate.py
+++ b/sympy/ntheory/generate.py
@@ -766,7 +766,7 @@ def primorial(n, nth=True):
     Examples
     ========
 
-    >>> from sympy.ntheory.generate import primorial, randprime, primerange
+    >>> from sympy.ntheory.generate import primorial, primerange
     >>> from sympy import factorint, Mul, primefactors, sqrt
     >>> primorial(4) # the first 4 primes are 2, 3, 5, 7
     210
diff --git a/sympy/ntheory/modular.py b/sympy/ntheory/modular.py
index c6f28e12f94b..32d3bf13430b 100644
--- a/sympy/ntheory/modular.py
+++ b/sympy/ntheory/modular.py
@@ -42,7 +42,7 @@ def crt(m, v, symmetric=False, check=True):
     As an example consider a set of residues ``U = [49, 76, 65]``
     and a set of moduli ``M = [99, 97, 95]``. Then we have::
 
-       >>> from sympy.ntheory.modular import crt, solve_congruence
+       >>> from sympy.ntheory.modular import crt
 
        >>> crt([99, 97, 95], [49, 76, 65])
        (639985, 912285)
diff --git a/sympy/ntheory/residue_ntheory.py b/sympy/ntheory/residue_ntheory.py
index 69ab7f7964e4..ca944ed4a5f3 100644
--- a/sympy/ntheory/residue_ntheory.py
+++ b/sympy/ntheory/residue_ntheory.py
@@ -975,7 +975,7 @@ def jacobi_symbol(m, n):
     ========
 
     >>> from sympy.ntheory import jacobi_symbol, legendre_symbol
-    >>> from sympy import Mul, S
+    >>> from sympy import S
     >>> jacobi_symbol(45, 77)
     -1
     >>> jacobi_symbol(60, 121)
@@ -1516,7 +1516,6 @@ def polynomial_congruence(expr, m):
     ========
 
     >>> from sympy.ntheory import polynomial_congruence
-    >>> from sympy import Poly
     >>> from sympy.abc import x
     >>> expr = x**6 - 2*x**5 -35
     >>> polynomial_congruence(expr, 6125)
diff --git a/sympy/physics/continuum_mechanics/beam.py b/sympy/physics/continuum_mechanics/beam.py
index d04496717cac..e2312690d7ab 100644
--- a/sympy/physics/continuum_mechanics/beam.py
+++ b/sympy/physics/continuum_mechanics/beam.py
@@ -754,7 +754,7 @@ def solve_for_reaction_loads(self, *reactions):
         being positive.
 
         >>> from sympy.physics.continuum_mechanics.beam import Beam
-        >>> from sympy import symbols, linsolve, limit
+        >>> from sympy import symbols
         >>> E, I = symbols('E, I')
         >>> R1, R2 = symbols('R1, R2')
         >>> b = Beam(30, E, I)
@@ -1513,7 +1513,6 @@ def plot_loading_results(self, subs=None):
 
             >>> from sympy.physics.continuum_mechanics.beam import Beam
             >>> from sympy import symbols
-            >>> from sympy.plotting import PlotGrid
             >>> R1, R2 = symbols('R1, R2')
             >>> b = Beam(8, 200*(10**9), 400*(10**-6))
             >>> b.apply_load(5000, 2, -1)
diff --git a/sympy/physics/hydrogen.py b/sympy/physics/hydrogen.py
index 8ec8e6c2da15..219200c882a9 100644
--- a/sympy/physics/hydrogen.py
+++ b/sympy/physics/hydrogen.py
@@ -21,9 +21,7 @@ def R_nl(n, l, r, Z=1):
     ========
 
     >>> from sympy.physics.hydrogen import R_nl
-    >>> from sympy import var
-    >>> var("r Z")
-    (r, Z)
+    >>> from sympy.abc import r, Z
     >>> R_nl(1, 0, r, Z)
     2*sqrt(Z**3)*exp(-Z*r)
     >>> R_nl(2, 0, r, Z)
@@ -151,10 +149,8 @@ def E_nl(n, Z=1):
     Examples
     ========
 
-    >>> from sympy import var
     >>> from sympy.physics.hydrogen import E_nl
-    >>> var("n Z")
-    (n, Z)
+    >>> from sympy.abc import n, Z
     >>> E_nl(n, Z)
     -Z**2/(2*n**2)
     >>> E_nl(1)
diff --git a/sympy/physics/mechanics/lagrange.py b/sympy/physics/mechanics/lagrange.py
index 7f44352ebf46..148884e7fa61 100644
--- a/sympy/physics/mechanics/lagrange.py
+++ b/sympy/physics/mechanics/lagrange.py
@@ -57,7 +57,7 @@ class LagrangesMethod(object):
 
         >>> from sympy.physics.mechanics import LagrangesMethod, Lagrangian
         >>> from sympy.physics.mechanics import ReferenceFrame, Particle, Point
-        >>> from sympy.physics.mechanics import dynamicsymbols, kinetic_energy
+        >>> from sympy.physics.mechanics import dynamicsymbols
         >>> from sympy import symbols
         >>> q = dynamicsymbols('q')
         >>> qd = dynamicsymbols('q', 1)
diff --git a/sympy/physics/mechanics/rigidbody.py b/sympy/physics/mechanics/rigidbody.py
index 140b64fd9900..b56f9f3b23b2 100644
--- a/sympy/physics/mechanics/rigidbody.py
+++ b/sympy/physics/mechanics/rigidbody.py
@@ -285,7 +285,7 @@ def potential_energy(self, scalar):
         Examples
         ========
 
-        >>> from sympy.physics.mechanics import Particle, Point, outer
+        >>> from sympy.physics.mechanics import Point, outer
         >>> from sympy.physics.mechanics import RigidBody, ReferenceFrame
         >>> from sympy import symbols
         >>> b = ReferenceFrame('b')
diff --git a/sympy/physics/optics/polarization.py b/sympy/physics/optics/polarization.py
index 56ca4d15801d..a0ccd62797ae 100644
--- a/sympy/physics/optics/polarization.py
+++ b/sympy/physics/optics/polarization.py
@@ -307,7 +307,7 @@ def jones_2_stokes(e):
     --------
     The axes on the Poincaré sphere.
 
-    >>> from sympy import pprint, Matrix, pi
+    >>> from sympy import pprint, pi
     >>> from sympy.physics.optics.polarization import jones_vector
     >>> from sympy.physics.optics.polarization import jones_2_stokes
     >>> H = jones_vector(0, 0)
@@ -557,7 +557,7 @@ def mueller_matrix(J):
     --------
     Generic optical components.
 
-    >>> from sympy import pprint, symbols, pi, simplify
+    >>> from sympy import pprint, symbols
     >>> from sympy.physics.optics.polarization import (mueller_matrix,
     ...     linear_polarizer, half_wave_retarder, quarter_wave_retarder)
     >>> theta = symbols("theta", real=True)
diff --git a/sympy/physics/optics/utils.py b/sympy/physics/optics/utils.py
index 02cc3211b5c1..fe644a7f154b 100644
--- a/sympy/physics/optics/utils.py
+++ b/sympy/physics/optics/utils.py
@@ -638,7 +638,6 @@ def hyperfocal_distance(f, N, c):
     Example
     =======
     >>> from sympy.physics.optics import hyperfocal_distance
-    >>> from sympy.abc import f, N, c
     >>> round(hyperfocal_distance(f = 0.5, N = 8, c = 0.0033), 2)
     9.47
     """
diff --git a/sympy/physics/paulialgebra.py b/sympy/physics/paulialgebra.py
index 1b7e68139bd5..a9b19e5c396f 100644
--- a/sympy/physics/paulialgebra.py
+++ b/sympy/physics/paulialgebra.py
@@ -169,7 +169,7 @@ def evaluate_pauli_product(arg):
     >>> evaluate_pauli_product(I*Pauli(1)*Pauli(2))
     -sigma3
 
-    >>> from sympy.abc import x,y
+    >>> from sympy.abc import x
     >>> evaluate_pauli_product(x**2*Pauli(2)*Pauli(1))
     -I*x**2*sigma3
     '''
diff --git a/sympy/physics/pring.py b/sympy/physics/pring.py
index fd5a8a546c7f..797f85568846 100644
--- a/sympy/physics/pring.py
+++ b/sympy/physics/pring.py
@@ -14,7 +14,7 @@ def wavefunction(n, x):
     Examples
     ========
 
-    >>> from sympy.physics.pring import wavefunction, energy
+    >>> from sympy.physics.pring import wavefunction
     >>> from sympy import Symbol, integrate, pi
     >>> x=Symbol("x")
     >>> wavefunction(1, x)
diff --git a/sympy/physics/qho_1d.py b/sympy/physics/qho_1d.py
index b4ab64d6adbd..6bceb345ee30 100644
--- a/sympy/physics/qho_1d.py
+++ b/sympy/physics/qho_1d.py
@@ -23,9 +23,7 @@ def psi_n(n, x, m, omega):
     ========
 
     >>> from sympy.physics.qho_1d import psi_n
-    >>> from sympy import var
-    >>> var("x m omega")
-    (x, m, omega)
+    >>> from sympy.abc import m, x, omega
     >>> psi_n(0, x, m, omega)
     (m*omega)**(1/4)*exp(-m*omega*x**2/(2*hbar))/(hbar**(1/4)*pi**(1/4))
 
@@ -58,9 +56,7 @@ def E_n(n, omega):
     ========
 
     >>> from sympy.physics.qho_1d import E_n
-    >>> from sympy import var
-    >>> var("x omega")
-    (x, omega)
+    >>> from sympy.abc import x, omega
     >>> E_n(x, omega)
     hbar*omega*(x + 1/2)
     """
diff --git a/sympy/physics/quantum/circuitutils.py b/sympy/physics/quantum/circuitutils.py
index 5f05a3ec327b..514c9dcf24c9 100644
--- a/sympy/physics/quantum/circuitutils.py
+++ b/sympy/physics/quantum/circuitutils.py
@@ -339,7 +339,7 @@ def convert_to_real_indices(seq, qubit_map):
 
     >>> from sympy import symbols
     >>> from sympy.physics.quantum.circuitutils import convert_to_real_indices
-    >>> from sympy.physics.quantum.gate import X, Y, Z, H
+    >>> from sympy.physics.quantum.gate import X, Y, H
     >>> i0, i1 = symbols('i:2')
     >>> index_map = {i0 : 0, i1 : 1}
     >>> convert_to_real_indices(X(i0)*Y(i1)*H(i0)*X(i1), index_map)
diff --git a/sympy/physics/quantum/density.py b/sympy/physics/quantum/density.py
index 496d9c3a67b1..d0f2da8b1df8 100644
--- a/sympy/physics/quantum/density.py
+++ b/sympy/physics/quantum/density.py
@@ -228,10 +228,8 @@ def entropy(density):
     ========
 
     >>> from sympy.physics.quantum.density import Density, entropy
-    >>> from sympy.physics.quantum.represent import represent
-    >>> from sympy.physics.quantum.matrixutils import scipy_sparse_matrix
-    >>> from sympy.physics.quantum.spin import JzKet, Jz
-    >>> from sympy import S, log
+    >>> from sympy.physics.quantum.spin import JzKet
+    >>> from sympy import S
     >>> up = JzKet(S(1)/2,S(1)/2)
     >>> down = JzKet(S(1)/2,-S(1)/2)
     >>> d = Density((up,S(1)/2),(down,S(1)/2))
@@ -275,7 +273,7 @@ def fidelity(state1, state2):
     >>> from sympy import S, sqrt
     >>> from sympy.physics.quantum.dagger import Dagger
     >>> from sympy.physics.quantum.spin import JzKet
-    >>> from sympy.physics.quantum.density import Density, fidelity
+    >>> from sympy.physics.quantum.density import fidelity
     >>> from sympy.physics.quantum.represent import represent
     >>>
     >>> up = JzKet(S(1)/2,S(1)/2)
diff --git a/sympy/physics/quantum/hilbert.py b/sympy/physics/quantum/hilbert.py
index 7f57ce90180b..b4bc7e1369fe 100644
--- a/sympy/physics/quantum/hilbert.py
+++ b/sympy/physics/quantum/hilbert.py
@@ -452,7 +452,6 @@ class DirectSumHilbertSpace(HilbertSpace):
     ========
 
     >>> from sympy.physics.quantum.hilbert import ComplexSpace, FockSpace
-    >>> from sympy import symbols
 
     >>> c = ComplexSpace(2)
     >>> f = FockSpace()
diff --git a/sympy/physics/quantum/identitysearch.py b/sympy/physics/quantum/identitysearch.py
index 83f637b9c87c..268ab3430749 100644
--- a/sympy/physics/quantum/identitysearch.py
+++ b/sympy/physics/quantum/identitysearch.py
@@ -725,8 +725,7 @@ def is_reducible(circuit, nqubits, begin, end):
 
     Check if the circuit can be reduced:
 
-    >>> from sympy.physics.quantum.identitysearch import (
-    ...     GateIdentity, is_reducible)
+    >>> from sympy.physics.quantum.identitysearch import is_reducible
     >>> from sympy.physics.quantum.gate import X, Y, Z
     >>> x = X(0); y = Y(0); z = Z(0)
     >>> is_reducible((x, y, z), 1, 0, 3)
@@ -780,7 +779,7 @@ def bfs_identity_search(gate_list, nqubits, max_depth=None,
     Find a list of gate identities:
 
     >>> from sympy.physics.quantum.identitysearch import bfs_identity_search
-    >>> from sympy.physics.quantum.gate import X, Y, Z, H
+    >>> from sympy.physics.quantum.gate import X, Y, Z
     >>> x = X(0); y = Y(0); z = Z(0)
     >>> bfs_identity_search([x], 1, max_depth=2)
     {GateIdentity(X(0), X(0))}
diff --git a/sympy/physics/quantum/innerproduct.py b/sympy/physics/quantum/innerproduct.py
index 8154416bc4e0..f8a726192eb5 100644
--- a/sympy/physics/quantum/innerproduct.py
+++ b/sympy/physics/quantum/innerproduct.py
@@ -35,7 +35,7 @@ class InnerProduct(Expr):
 
     Create an InnerProduct and check its properties:
 
-        >>> from sympy.physics.quantum import Bra, Ket, InnerProduct
+        >>> from sympy.physics.quantum import Bra, Ket
         >>> b = Bra('b')
         >>> k = Ket('k')
         >>> ip = b*k
diff --git a/sympy/physics/quantum/operator.py b/sympy/physics/quantum/operator.py
index 1a8db4b3c30f..8b1e339644c3 100644
--- a/sympy/physics/quantum/operator.py
+++ b/sympy/physics/quantum/operator.py
@@ -51,7 +51,7 @@ class Operator(QExpr):
     Create an operator and examine its attributes::
 
         >>> from sympy.physics.quantum import Operator
-        >>> from sympy import symbols, I
+        >>> from sympy import I
         >>> A = Operator('A')
         >>> A
         A
diff --git a/sympy/physics/quantum/qubit.py b/sympy/physics/quantum/qubit.py
index 536843aeb71c..1c1d4e2ddb86 100644
--- a/sympy/physics/quantum/qubit.py
+++ b/sympy/physics/quantum/qubit.py
@@ -446,7 +446,6 @@ def matrix_to_qubit(matrix):
     Represent a state and then go back to its qubit form:
 
         >>> from sympy.physics.quantum.qubit import matrix_to_qubit, Qubit
-        >>> from sympy.physics.quantum.gate import Z
         >>> from sympy.physics.quantum.represent import represent
         >>> q = Qubit('01')
         >>> matrix_to_qubit(represent(q))
@@ -556,7 +555,7 @@ def measure_all(qubit, format='sympy', normalize=True):
     ========
 
         >>> from sympy.physics.quantum.qubit import Qubit, measure_all
-        >>> from sympy.physics.quantum.gate import H, X, Y, Z
+        >>> from sympy.physics.quantum.gate import H
         >>> from sympy.physics.quantum.qapply import qapply
 
         >>> c = H(0)*H(1)*Qubit('00')
@@ -614,7 +613,7 @@ def measure_partial(qubit, bits, format='sympy', normalize=True):
     ========
 
         >>> from sympy.physics.quantum.qubit import Qubit, measure_partial
-        >>> from sympy.physics.quantum.gate import H, X, Y, Z
+        >>> from sympy.physics.quantum.gate import H
         >>> from sympy.physics.quantum.qapply import qapply
 
         >>> c = H(0)*H(1)*Qubit('00')
diff --git a/sympy/physics/quantum/represent.py b/sympy/physics/quantum/represent.py
index 6ff8be3c867c..3371823c2a6b 100644
--- a/sympy/physics/quantum/represent.py
+++ b/sympy/physics/quantum/represent.py
@@ -306,7 +306,7 @@ def rep_expectation(expr, **options):
     Examples
     ========
 
-    >>> from sympy.physics.quantum.cartesian import XOp, XKet, PxOp, PxKet
+    >>> from sympy.physics.quantum.cartesian import XOp, PxOp, PxKet
     >>> from sympy.physics.quantum.represent import rep_expectation
     >>> rep_expectation(XOp())
     x_1*DiracDelta(x_1 - x_2)
diff --git a/sympy/physics/quantum/state.py b/sympy/physics/quantum/state.py
index 51fc26fce198..0b87cfee5c97 100644
--- a/sympy/physics/quantum/state.py
+++ b/sympy/physics/quantum/state.py
@@ -354,7 +354,7 @@ class for all physical, time-independent Kets in a system. This class
 
     Create a simple Ket and looking at its properties::
 
-        >>> from sympy.physics.quantum import Ket, Bra
+        >>> from sympy.physics.quantum import Ket
         >>> from sympy import symbols, I
         >>> k = Ket('psi')
         >>> k
@@ -418,7 +418,7 @@ class and its subclasses will be the main classes that users will use for
 
     Create a simple Bra and look at its properties::
 
-        >>> from sympy.physics.quantum import Ket, Bra
+        >>> from sympy.physics.quantum import Bra
         >>> from sympy import symbols, I
         >>> b = Bra('psi')
         >>> b
@@ -600,7 +600,6 @@ class TimeDepBra(TimeDepState, BraBase):
     ========
 
         >>> from sympy.physics.quantum import TimeDepBra
-        >>> from sympy import symbols, I
         >>> b = TimeDepBra('psi', 't')
         >>> b
         <psi;t|
@@ -629,7 +628,7 @@ class OrthogonalKet(OrthogonalState, KetBase):
     The inner product of two states with different labels will give zero,
     states with the same label will give one.
 
-        >>> from sympy.physics.quantum import OrthogonalBra, OrthogonalKet, qapply
+        >>> from sympy.physics.quantum import OrthogonalBra, OrthogonalKet
         >>> from sympy.abc import m, n
         >>> (OrthogonalBra(n)*OrthogonalKet(n)).doit()
         1
@@ -968,7 +967,7 @@ def normalize(self):
         ========
 
             >>> from sympy import symbols, pi
-            >>> from sympy.functions import sqrt, sin
+            >>> from sympy.functions import sin
             >>> from sympy.physics.quantum.state import Wavefunction
             >>> x = symbols('x', real=True)
             >>> L = symbols('L', positive=True)
@@ -994,7 +993,7 @@ def prob(self):
         ========
 
             >>> from sympy import symbols, pi
-            >>> from sympy.functions import sqrt, sin
+            >>> from sympy.functions import sin
             >>> from sympy.physics.quantum.state import Wavefunction
             >>> x, L = symbols('x,L', real=True)
             >>> n = symbols('n', integer=True)
diff --git a/sympy/physics/quantum/tensorproduct.py b/sympy/physics/quantum/tensorproduct.py
index 5ebfe9a405cc..6ff0b6a72863 100644
--- a/sympy/physics/quantum/tensorproduct.py
+++ b/sympy/physics/quantum/tensorproduct.py
@@ -71,7 +71,7 @@ class TensorProduct(Expr):
 
     Start with a simple tensor product of sympy matrices::
 
-        >>> from sympy import I, Matrix, symbols
+        >>> from sympy import Matrix
         >>> from sympy.physics.quantum import TensorProduct
 
         >>> m1 = Matrix([[1,2],[3,4]])
diff --git a/sympy/physics/secondquant.py b/sympy/physics/secondquant.py
index 3923c9c68b0a..0e22a282422c 100644
--- a/sympy/physics/secondquant.py
+++ b/sympy/physics/secondquant.py
@@ -2771,8 +2771,8 @@ def wicks(e, **kw_args):
     Examples
     ========
 
-    >>> from sympy import symbols, Function, Dummy
-    >>> from sympy.physics.secondquant import wicks, F, Fd, NO
+    >>> from sympy import symbols, Dummy
+    >>> from sympy.physics.secondquant import wicks, F, Fd
     >>> p, q, r = symbols('p,q,r')
     >>> wicks(Fd(p)*F(q))
     KroneckerDelta(_i, q)*KroneckerDelta(p, q) + NO(CreateFermion(p)*AnnihilateFermion(q))
diff --git a/sympy/physics/sho.py b/sympy/physics/sho.py
index 2f4c2b2fc953..4a38c5cc4a95 100644
--- a/sympy/physics/sho.py
+++ b/sympy/physics/sho.py
@@ -25,9 +25,7 @@ def R_nl(n, l, nu, r):
     ========
 
     >>> from sympy.physics.sho import R_nl
-    >>> from sympy import var
-    >>> var("r nu l")
-    (r, nu, l)
+    >>> from sympy.abc import r, nu, l
     >>> R_nl(0, 0, 1, r)
     2*2**(3/4)*exp(-r**2)/pi**(1/4)
     >>> R_nl(1, 0, 1, r)
diff --git a/sympy/physics/units/prefixes.py b/sympy/physics/units/prefixes.py
index a49e9d13199a..c6995d4d0073 100644
--- a/sympy/physics/units/prefixes.py
+++ b/sympy/physics/units/prefixes.py
@@ -126,7 +126,6 @@ def prefix_unit(unit, prefixes):
 
         >>> from sympy.physics.units.prefixes import (PREFIXES,
         ...                                                 prefix_unit)
-        >>> from sympy.physics.units.systems import MKS
         >>> from sympy.physics.units import m
         >>> pref = {"m": PREFIXES["m"], "c": PREFIXES["c"], "d": PREFIXES["d"]}
         >>> prefix_unit(m, pref)  # doctest: +SKIP
diff --git a/sympy/physics/units/util.py b/sympy/physics/units/util.py
index 3e2ce1a434f5..067893854ad2 100644
--- a/sympy/physics/units/util.py
+++ b/sympy/physics/units/util.py
@@ -76,7 +76,6 @@ def convert_to(expr, target_units, unit_system="SI"):
 
     Conversion to Planck units:
 
-    >>> from sympy.physics.units import gravitational_constant, hbar
     >>> convert_to(atomic_mass_constant, [gravitational_constant, speed_of_light, hbar]).n()
     7.62963085040767e-20*gravitational_constant**(-0.5)*hbar**0.5*speed_of_light**0.5
 
diff --git a/sympy/physics/vector/frame.py b/sympy/physics/vector/frame.py
index bb9d710b8476..57b78b0c873e 100644
--- a/sympy/physics/vector/frame.py
+++ b/sympy/physics/vector/frame.py
@@ -341,7 +341,7 @@ def ang_acc_in(self, otherframe):
         Examples
         ========
 
-        >>> from sympy.physics.vector import ReferenceFrame, Vector
+        >>> from sympy.physics.vector import ReferenceFrame
         >>> N = ReferenceFrame('N')
         >>> A = ReferenceFrame('A')
         >>> V = 10 * N.x
@@ -374,7 +374,7 @@ def ang_vel_in(self, otherframe):
         Examples
         ========
 
-        >>> from sympy.physics.vector import ReferenceFrame, Vector
+        >>> from sympy.physics.vector import ReferenceFrame
         >>> N = ReferenceFrame('N')
         >>> A = ReferenceFrame('A')
         >>> V = 10 * N.x
@@ -942,7 +942,7 @@ def set_ang_acc(self, otherframe, value):
         Examples
         ========
 
-        >>> from sympy.physics.vector import ReferenceFrame, Vector
+        >>> from sympy.physics.vector import ReferenceFrame
         >>> N = ReferenceFrame('N')
         >>> A = ReferenceFrame('A')
         >>> V = 10 * N.x
@@ -978,7 +978,7 @@ def set_ang_vel(self, otherframe, value):
         Examples
         ========
 
-        >>> from sympy.physics.vector import ReferenceFrame, Vector
+        >>> from sympy.physics.vector import ReferenceFrame
         >>> N = ReferenceFrame('N')
         >>> A = ReferenceFrame('A')
         >>> V = 10 * N.x
diff --git a/sympy/physics/vector/point.py b/sympy/physics/vector/point.py
index 2441a6018370..21c86a410b03 100644
--- a/sympy/physics/vector/point.py
+++ b/sympy/physics/vector/point.py
@@ -110,7 +110,7 @@ def a1pt_theory(self, otherpoint, outframe, interframe):
         ========
 
         >>> from sympy.physics.vector import Point, ReferenceFrame
-        >>> from sympy.physics.vector import Vector, dynamicsymbols
+        >>> from sympy.physics.vector import dynamicsymbols
         >>> q = dynamicsymbols('q')
         >>> q2 = dynamicsymbols('q2')
         >>> qd = dynamicsymbols('q', 1)
@@ -387,7 +387,7 @@ def v1pt_theory(self, otherpoint, outframe, interframe):
         ========
 
         >>> from sympy.physics.vector import Point, ReferenceFrame
-        >>> from sympy.physics.vector import Vector, dynamicsymbols
+        >>> from sympy.physics.vector import dynamicsymbols
         >>> q = dynamicsymbols('q')
         >>> q2 = dynamicsymbols('q2')
         >>> qd = dynamicsymbols('q', 1)
diff --git a/sympy/physics/vector/vector.py b/sympy/physics/vector/vector.py
index 888e91a54c6b..89ff52745a2b 100644
--- a/sympy/physics/vector/vector.py
+++ b/sympy/physics/vector/vector.py
@@ -402,7 +402,7 @@ def __xor__(self, other):
         Examples
         ========
 
-        >>> from sympy.physics.vector import ReferenceFrame, Vector
+        >>> from sympy.physics.vector import ReferenceFrame
         >>> from sympy import symbols
         >>> q1 = symbols('q1')
         >>> N = ReferenceFrame('N')
@@ -595,7 +595,7 @@ def express(self, otherframe, variables=False):
         Examples
         ========
 
-        >>> from sympy.physics.vector import ReferenceFrame, Vector, dynamicsymbols
+        >>> from sympy.physics.vector import ReferenceFrame, dynamicsymbols
         >>> q1 = dynamicsymbols('q1')
         >>> N = ReferenceFrame('N')
         >>> A = N.orientnew('A', 'Axis', [q1, N.y])
@@ -625,7 +625,6 @@ def to_matrix(self, reference_frame):
 
         >>> from sympy import symbols
         >>> from sympy.physics.vector import ReferenceFrame
-        >>> from sympy.physics.mechanics.functions import inertia
         >>> a, b, c = symbols('a, b, c')
         >>> N = ReferenceFrame('N')
         >>> vector = a * N.x + b * N.y + c * N.z
diff --git a/sympy/physics/wigner.py b/sympy/physics/wigner.py
index 80a790b91ee1..7e82caaab56c 100644
--- a/sympy/physics/wigner.py
+++ b/sympy/physics/wigner.py
@@ -921,7 +921,6 @@ def wigner_d(J, alpha, beta, gamma):
 
     >>> from sympy.physics.wigner import wigner_d
     >>> from sympy import Integer, symbols, pprint
-    >>> from sympy.physics.wigner import wigner_d_small
     >>> half = 1/Integer(2)
     >>> alpha, beta, gamma = symbols("alpha, beta, gamma", real=True)
     >>> pprint(wigner_d(half, alpha, beta, gamma), use_unicode=True)
diff --git a/sympy/plotting/experimental_lambdify.py b/sympy/plotting/experimental_lambdify.py
index 74c244c9324b..9b3d940ba99b 100644
--- a/sympy/plotting/experimental_lambdify.py
+++ b/sympy/plotting/experimental_lambdify.py
@@ -545,7 +545,7 @@ def tree2str(cls, tree):
         ========
 
         >>> from sympy.abc import x, y, z
-        >>> from sympy import Integral, sin
+        >>> from sympy import sin
         >>> from sympy.plotting.experimental_lambdify import Lambdifier
         >>> str2tree = Lambdifier([x], x).str2tree
         >>> tree2str = Lambdifier([x], x).tree2str
diff --git a/sympy/plotting/plot.py b/sympy/plotting/plot.py
index ae612e144226..f910b8a267d3 100644
--- a/sympy/plotting/plot.py
+++ b/sympy/plotting/plot.py
@@ -2280,7 +2280,7 @@ def check_arguments(args, expr_len, nb_of_free_symbols):
        :format: doctest
        :include-source: True
 
-       >>> from sympy import plot, cos, sin, symbols
+       >>> from sympy import cos, sin, symbols
        >>> from sympy.plotting.plot import check_arguments
        >>> x = symbols('x')
        >>> check_arguments([cos(x), sin(x)], 2, 1)
diff --git a/sympy/plotting/plot_implicit.py b/sympy/plotting/plot_implicit.py
index 2aa257cda8c7..fa75cd2c6f44 100644
--- a/sympy/plotting/plot_implicit.py
+++ b/sympy/plotting/plot_implicit.py
@@ -262,7 +262,7 @@ def plot_implicit(expr, x_var=None, y_var=None, adaptive=True, depth=0,
         :format: doctest
         :include-source: True
 
-        >>> from sympy import plot_implicit, cos, sin, symbols, Eq, And
+        >>> from sympy import plot_implicit, symbols, Eq, And
         >>> x, y = symbols('x y')
 
     Without any ranges for the symbols in the expression:
diff --git a/sympy/plotting/pygletplot/__init__.py b/sympy/plotting/pygletplot/__init__.py
index 7007a80d647e..cd86a505d8c4 100644
--- a/sympy/plotting/pygletplot/__init__.py
+++ b/sympy/plotting/pygletplot/__init__.py
@@ -13,7 +13,6 @@ def PygletPlot(*args, **kwargs):
     See examples/advanced/pyglet_plotting.py for many more examples.
 
     >>> from sympy.plotting.pygletplot import PygletPlot as Plot
-    >>> from sympy import symbols
     >>> from sympy.abc import x, y, z
 
     >>> Plot(x*y**3-y*x**3)
diff --git a/sympy/plotting/pygletplot/plot.py b/sympy/plotting/pygletplot/plot.py
index c8ce706797ef..7fdb1defdf5f 100644
--- a/sympy/plotting/pygletplot/plot.py
+++ b/sympy/plotting/pygletplot/plot.py
@@ -165,7 +165,6 @@ def __init__(self, *fargs, **win_args):
         other words...
 
         >>> from sympy.plotting.pygletplot import PygletPlot as Plot
-        >>> from sympy.core import Symbol
         >>> from sympy.abc import x
         >>> p = Plot(x**2, visible=False)
 
diff --git a/sympy/polys/distributedmodules.py b/sympy/polys/distributedmodules.py
index 5093a2aee41e..6a6cfd78fea1 100644
--- a/sympy/polys/distributedmodules.py
+++ b/sympy/polys/distributedmodules.py
@@ -381,7 +381,7 @@ def sdm_to_vector(f, gens, K, n=None):
 
     >>> from sympy.polys.distributedmodules import sdm_to_vector
     >>> from sympy.abc import x, y, z
-    >>> from sympy.polys import QQ, lex
+    >>> from sympy.polys import QQ
     >>> f = [((1, 0, 0, 1), QQ(2)), ((0, 2, 0, 0), QQ(1)), ((0, 0, 2, 0), QQ(1))]
     >>> sdm_to_vector(f, [x, y, z], QQ)
     [x**2 + y**2, 2*z]
diff --git a/sympy/polys/galoistools.py b/sympy/polys/galoistools.py
index 87c49d109b98..840823930ab5 100644
--- a/sympy/polys/galoistools.py
+++ b/sympy/polys/galoistools.py
@@ -29,7 +29,6 @@ def gf_crt(U, M, K=None):
 
        >>> from sympy.polys.domains import ZZ
        >>> from sympy.polys.galoistools import gf_crt
-       >>> from sympy.ntheory.modular import solve_congruence
 
        >>> gf_crt([49, 76, 65], [99, 97, 95], ZZ)
        639985
@@ -315,7 +314,6 @@ def gf_from_int_poly(f, p):
     Examples
     ========
 
-    >>> from sympy.polys.domains import ZZ
     >>> from sympy.polys.galoistools import gf_from_int_poly
 
     >>> gf_from_int_poly([7, -2, 3], 5)
diff --git a/sympy/polys/monomials.py b/sympy/polys/monomials.py
index 844840d121f5..5df7ed0cfc65 100644
--- a/sympy/polys/monomials.py
+++ b/sympy/polys/monomials.py
@@ -82,8 +82,6 @@ def itermonomials(variables, max_degrees, min_degrees=None):
         >>> from sympy import symbols
         >>> from sympy.polys.monomials import itermonomials
         >>> from sympy.polys.orderings import monomial_key
-        >>> from itertools import product
-        >>> from sympy.core import Mul
         >>> from sympy.abc import x, y
 
         >>> sorted(itermonomials([x, y], [2, 4], [1, 2]), reverse=True, key=monomial_key('lex', [x, y]))
diff --git a/sympy/polys/orderings.py b/sympy/polys/orderings.py
index 59269461d9f3..464eabf144f4 100644
--- a/sympy/polys/orderings.py
+++ b/sympy/polys/orderings.py
@@ -266,7 +266,7 @@ def build_product_order(arg, gens):
 
     For example, build a product of two grlex orders:
 
-    >>> from sympy.polys.orderings import grlex, build_product_order
+    >>> from sympy.polys.orderings import build_product_order
     >>> from sympy.abc import x, y, z, t
 
     >>> O = build_product_order((("grlex", x, y), ("grlex", z, t)), [x, y, z, t])
diff --git a/sympy/polys/partfrac.py b/sympy/polys/partfrac.py
index 2d75731af8eb..ee5ffdc6403c 100644
--- a/sympy/polys/partfrac.py
+++ b/sympy/polys/partfrac.py
@@ -428,7 +428,7 @@ def assemble_partfrac_list(partial_list):
     This example is taken from Bronstein's original paper:
 
     >>> from sympy.polys.partfrac import apart_list, assemble_partfrac_list
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
 
     >>> f = 36 / (x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2)
     >>> pfd = apart_list(f)
diff --git a/sympy/polys/polyroots.py b/sympy/polys/polyroots.py
index d9a9f53249bb..979b27483bee 100644
--- a/sympy/polys/polyroots.py
+++ b/sympy/polys/polyroots.py
@@ -265,7 +265,7 @@ def roots_quartic(f):
     Examples
     ========
 
-        >>> from sympy import Poly, symbols, I
+        >>> from sympy import Poly
         >>> from sympy.polys.polyroots import roots_quartic
 
         >>> r = roots_quartic(Poly('x**4-6*x**3+17*x**2-26*x+20'))
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index f066202e69f3..a0110db39476 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -4484,7 +4484,7 @@ def total_degree(f, *gens):
     Examples
     ========
     >>> from sympy import total_degree, Poly
-    >>> from sympy.abc import x, y, z
+    >>> from sympy.abc import x, y
 
     >>> total_degree(1)
     0
diff --git a/sympy/polys/rings.py b/sympy/polys/rings.py
index d836fd482e02..4e02e9de0f46 100644
--- a/sympy/polys/rings.py
+++ b/sympy/polys/rings.py
@@ -115,7 +115,7 @@ def vring(symbols, domain, order=lex):
 
     >>> vring("x,y,z", ZZ, lex)
     Polynomial ring in x, y, z over ZZ with lex order
-    >>> x + y + z
+    >>> x + y + z # noqa:
     x + y + z
     >>> type(_)
     <class 'sympy.polys.rings.PolyElement'>
@@ -141,8 +141,6 @@ def sring(exprs, *symbols, **options):
 
     >>> from sympy.core import symbols
     >>> from sympy.polys.rings import sring
-    >>> from sympy.polys.domains import ZZ
-    >>> from sympy.polys.orderings import lex
 
     >>> x, y, z = symbols("x,y,z")
     >>> R, f = sring(x + 2*y + 3*z)
diff --git a/sympy/polys/rootisolation.py b/sympy/polys/rootisolation.py
index fccdac6ced45..3108de55304f 100644
--- a/sympy/polys/rootisolation.py
+++ b/sympy/polys/rootisolation.py
@@ -1796,7 +1796,7 @@ class ComplexInterval(object):
     Examples
     ========
 
-    >>> from sympy import CRootOf, Rational, S
+    >>> from sympy import CRootOf, S
     >>> from sympy.abc import x
     >>> CRootOf.clear_cache()  # for doctest reproducibility
     >>> root = CRootOf(x**10 - 2*x + 3, 9)
diff --git a/sympy/polys/rootoftools.py b/sympy/polys/rootoftools.py
index 3a3eaca75632..8d8a199ff826 100644
--- a/sympy/polys/rootoftools.py
+++ b/sympy/polys/rootoftools.py
@@ -43,7 +43,7 @@ class _pure_key_dict(object):
     Only the following actions are guaranteed:
 
     >>> from sympy.polys.rootoftools import _pure_key_dict
-    >>> from sympy import S, PurePoly
+    >>> from sympy import PurePoly
     >>> from sympy.abc import x, y
 
     1) creation
@@ -907,7 +907,7 @@ def eval_rational(self, dx=None, dy=None, n=15):
         ensure the decimal representation of the approximation will be
         correct (including rounding) to 6 digits:
 
-        >>> from sympy import S, legendre_poly, Symbol
+        >>> from sympy import legendre_poly, Symbol
         >>> x = Symbol("x")
         >>> p = legendre_poly(4, x, polys=True)
         >>> r = p.real_roots()[-1]
diff --git a/sympy/printing/dot.py b/sympy/printing/dot.py
index 5ae80e28ed4a..a08d21f817f6 100644
--- a/sympy/printing/dot.py
+++ b/sympy/printing/dot.py
@@ -33,7 +33,8 @@ def purestr(x, with_args=False):
     Examples
     ========
 
-    >>> from sympy import Integer, Float, Symbol, MatrixSymbol
+    >>> from sympy import Float, Symbol, MatrixSymbol
+    >>> from sympy import Integer # noqa: F401
     >>> from sympy.printing.dot import purestr
 
     Applying ``purestr`` for basic symbolic object:
diff --git a/sympy/printing/julia.py b/sympy/printing/julia.py
index ee3dad2fc1d8..0e01c2986fc4 100644
--- a/sympy/printing/julia.py
+++ b/sympy/printing/julia.py
@@ -547,7 +547,7 @@ def julia_code(expr, assign_to=None, **settings):
     >>> julia_code(sin(x).series(x).removeO())
     'x.^5/120 - x.^3/6 + x'
 
-    >>> from sympy import Rational, ceiling, Abs
+    >>> from sympy import Rational, ceiling
     >>> x, y, tau = symbols("x, y, tau")
     >>> julia_code((2*tau)**Rational(7, 2))
     '8*sqrt(2)*tau.^(7/2)'
@@ -628,7 +628,7 @@ def julia_code(expr, assign_to=None, **settings):
     ``contract=False`` will just print the assignment expression that should be
     looped over:
 
-    >>> from sympy import Eq, IndexedBase, Idx, ccode
+    >>> from sympy import Eq, IndexedBase, Idx
     >>> len_y = 5
     >>> y = IndexedBase('y', shape=(len_y,))
     >>> t = IndexedBase('t', shape=(len_y,))
diff --git a/sympy/printing/llvmjitcode.py b/sympy/printing/llvmjitcode.py
index 09facc44be0b..b2623744e457 100644
--- a/sympy/printing/llvmjitcode.py
+++ b/sympy/printing/llvmjitcode.py
@@ -429,7 +429,7 @@ def llvm_callable(args, expr, callback_type=None):
     The 'cubature' callback handles multiple expressions (set `fdim`
     to match in the integration call.)
     >>> import sympy.printing.llvmjitcode as jit
-    >>> from sympy import cse, exp
+    >>> from sympy import cse
     >>> from sympy.abc import x,y
     >>> e1 = x*x + y*y
     >>> e2 = 4*(x*x + y*y) + 8.0
diff --git a/sympy/printing/octave.py b/sympy/printing/octave.py
index 5d3be42ad1ee..66479e4db5ce 100644
--- a/sympy/printing/octave.py
+++ b/sympy/printing/octave.py
@@ -631,7 +631,7 @@ def octave_code(expr, assign_to=None, **settings):
     >>> octave_code(sin(x).series(x).removeO())
     'x.^5/120 - x.^3/6 + x'
 
-    >>> from sympy import Rational, ceiling, Abs
+    >>> from sympy import Rational, ceiling
     >>> x, y, tau = symbols("x, y, tau")
     >>> octave_code((2*tau)**Rational(7, 2))
     '8*sqrt(2)*tau.^(7/2)'
@@ -712,7 +712,7 @@ def octave_code(expr, assign_to=None, **settings):
     ``contract=False`` will just print the assignment expression that should be
     looped over:
 
-    >>> from sympy import Eq, IndexedBase, Idx, ccode
+    >>> from sympy import Eq, IndexedBase, Idx
     >>> len_y = 5
     >>> y = IndexedBase('y', shape=(len_y,))
     >>> t = IndexedBase('t', shape=(len_y,))
diff --git a/sympy/series/formal.py b/sympy/series/formal.py
index acbdab6181be..5f31c26ff82c 100644
--- a/sympy/series/formal.py
+++ b/sympy/series/formal.py
@@ -57,7 +57,7 @@ def rational_algorithm(f, x, k, order=4, full=False):
     Examples
     ========
 
-    >>> from sympy import log, atan, I
+    >>> from sympy import log, atan
     >>> from sympy.series.formal import rational_algorithm as ra
     >>> from sympy.abc import x, k
 
@@ -1172,7 +1172,7 @@ def product(self, other, x=None, n=6):
         Examples
         ========
 
-        >>> from sympy import fps, sin, exp, convolution
+        >>> from sympy import fps, sin, exp
         >>> from sympy.abc import x
         >>> f1 = fps(sin(x))
         >>> f2 = fps(exp(x))
@@ -1264,7 +1264,7 @@ def compose(self, other, x=None, n=6):
         Examples
         ========
 
-        >>> from sympy import fps, sin, exp, bell
+        >>> from sympy import fps, sin, exp
         >>> from sympy.abc import x
         >>> f1 = fps(exp(x))
         >>> f2 = fps(sin(x))
@@ -1337,7 +1337,7 @@ def inverse(self, x=None, n=6):
         Examples
         ========
 
-        >>> from sympy import fps, exp, cos, bell
+        >>> from sympy import fps, exp, cos
         >>> from sympy.abc import x
         >>> f1 = fps(exp(x))
         >>> f2 = fps(cos(x))
@@ -1531,7 +1531,7 @@ def _eval_terms(self, n):
         Examples
         ========
 
-        >>> from sympy import fps, sin, exp, convolution
+        >>> from sympy import fps, sin, exp
         >>> from sympy.abc import x
         >>> f1 = fps(sin(x))
         >>> f2 = fps(exp(x))
@@ -1596,7 +1596,7 @@ def _eval_terms(self, n):
         Examples
         ========
 
-        >>> from sympy import fps, sin, exp, bell
+        >>> from sympy import fps, sin, exp
         >>> from sympy.abc import x
         >>> f1 = fps(exp(x))
         >>> f2 = fps(sin(x))
@@ -1681,7 +1681,7 @@ def _eval_terms(self, n):
         Examples
         ========
 
-        >>> from sympy import fps, exp, cos, bell
+        >>> from sympy import fps, exp, cos
         >>> from sympy.abc import x
         >>> f1 = fps(exp(x))
         >>> f2 = fps(cos(x))
@@ -1751,7 +1751,7 @@ def fps(f, x=None, x0=0, dir=1, hyper=True, order=4, rational=True, full=False):
     Examples
     ========
 
-    >>> from sympy import fps, O, ln, atan, sin
+    >>> from sympy import fps, ln, atan, sin
     >>> from sympy.abc import x, n
 
     Rational Functions
diff --git a/sympy/series/fourier.py b/sympy/series/fourier.py
index ddd5eaf2e286..5030e3ff5c18 100644
--- a/sympy/series/fourier.py
+++ b/sympy/series/fourier.py
@@ -48,7 +48,6 @@ def _process_limits(func, limits):
     Examples
     ========
 
-    >>> from sympy import pi
     >>> from sympy.series.fourier import _process_limits as pari
     >>> from sympy.abc import x
     >>> pari(x**2, (x, -2, 2))
@@ -645,9 +644,8 @@ def fourier_series(f, limits=None, finite=True):
 
     Computing the Fourier series of $f(x) = x^2$:
 
-    >>> from sympy import fourier_series, pi, cos
+    >>> from sympy import fourier_series, pi
     >>> from sympy.abc import x
-
     >>> f = x**2
     >>> s = fourier_series(f, (x, -pi, pi))
     >>> s1 = s.truncate(n=3)
diff --git a/sympy/series/kauers.py b/sympy/series/kauers.py
index bb0d48e32d87..4d8458e8c334 100644
--- a/sympy/series/kauers.py
+++ b/sympy/series/kauers.py
@@ -11,9 +11,8 @@ def finite_diff(expression, variable, increment=1):
     Examples
     ========
 
-    >>> from sympy.abc import x, y, z, k, n
+    >>> from sympy.abc import x, y, z
     >>> from sympy.series.kauers import finite_diff
-    >>> from  sympy import Sum
     >>> finite_diff(x**2, x)
     2*x + 1
     >>> finite_diff(y**3 + 2*y**2 + 3*y + 4, y)
diff --git a/sympy/series/limits.py b/sympy/series/limits.py
index 01b23ea89866..009ae26fdf28 100644
--- a/sympy/series/limits.py
+++ b/sympy/series/limits.py
@@ -38,7 +38,7 @@ def limit(e, z, z0, dir="+"):
     Examples
     ========
 
-    >>> from sympy import limit, sin, Symbol, oo
+    >>> from sympy import limit, sin, oo
     >>> from sympy.abc import x
     >>> limit(sin(x)/x, x, 0)
     1
@@ -138,7 +138,7 @@ class Limit(Expr):
     Examples
     ========
 
-    >>> from sympy import Limit, sin, Symbol
+    >>> from sympy import Limit, sin
     >>> from sympy.abc import x
     >>> Limit(sin(x)/x, x, 0)
     Limit(sin(x)/x, x, 0)
diff --git a/sympy/series/sequences.py b/sympy/series/sequences.py
index eac7276eae5e..1062f0dbd4a7 100644
--- a/sympy/series/sequences.py
+++ b/sympy/series/sequences.py
@@ -188,7 +188,7 @@ def coeff_mul(self, other):
         Examples
         ========
 
-        >>> from sympy import S, oo, SeqFormula
+        >>> from sympy import SeqFormula
         >>> from sympy.abc import n
         >>> SeqFormula(n**2).coeff_mul(2)
         SeqFormula(2*n**2, (n, 0, oo))
@@ -208,7 +208,7 @@ def __add__(self, other):
         Examples
         ========
 
-        >>> from sympy import S, oo, SeqFormula
+        >>> from sympy import SeqFormula
         >>> from sympy.abc import n
         >>> SeqFormula(n**2) + SeqFormula(n**3)
         SeqFormula(n**3 + n**2, (n, 0, oo))
@@ -229,7 +229,7 @@ def __sub__(self, other):
         Examples
         ========
 
-        >>> from sympy import S, oo, SeqFormula
+        >>> from sympy import SeqFormula
         >>> from sympy.abc import n
         >>> SeqFormula(n**2) - (SeqFormula(n))
         SeqFormula(n**2 - n, (n, 0, oo))
@@ -248,7 +248,7 @@ def __neg__(self):
         Examples
         ========
 
-        >>> from sympy import S, oo, SeqFormula
+        >>> from sympy import SeqFormula
         >>> from sympy.abc import n
         >>> -SeqFormula(n**2)
         SeqFormula(-n**2, (n, 0, oo))
@@ -264,7 +264,7 @@ def __mul__(self, other):
         Examples
         ========
 
-        >>> from sympy import S, oo, SeqFormula
+        >>> from sympy import SeqFormula
         >>> from sympy.abc import n
         >>> SeqFormula(n**2) * (SeqFormula(n))
         SeqFormula(n**3, (n, 0, oo))
@@ -379,7 +379,7 @@ class EmptySequence(SeqBase, metaclass=Singleton):
     Examples
     ========
 
-    >>> from sympy import EmptySequence, SeqPer, oo
+    >>> from sympy import EmptySequence, SeqPer
     >>> from sympy.abc import x
     >>> EmptySequence
     EmptySequence
@@ -917,7 +917,7 @@ def sequence(seq, limits=None):
     Examples
     ========
 
-    >>> from sympy import sequence, SeqPer, SeqFormula
+    >>> from sympy import sequence
     >>> from sympy.abc import n
     >>> sequence(n**2, (n, 0, 5))
     SeqFormula(n**2, (n, 0, 5))
diff --git a/sympy/series/series.py b/sympy/series/series.py
index 95ac89b34c79..93bcb61323ff 100644
--- a/sympy/series/series.py
+++ b/sympy/series/series.py
@@ -32,7 +32,7 @@ def series(expr, x=None, x0=0, n=6, dir="+"):
     Examples
     ========
 
-    >>> from sympy import Symbol, series, tan, oo
+    >>> from sympy import series, tan, oo
     >>> from sympy.abc import x
     >>> f = tan(x)
     >>> series(f, x, 2, 6, "+")
diff --git a/sympy/sets/fancysets.py b/sympy/sets/fancysets.py
index 0a94d81db36f..bf013522e14a 100644
--- a/sympy/sets/fancysets.py
+++ b/sympy/sets/fancysets.py
@@ -241,7 +241,7 @@ class Reals(Interval, metaclass=Singleton):
     Examples
     ========
 
-    >>> from sympy import S, Interval, Rational, pi, I
+    >>> from sympy import S, Rational, pi, I
     >>> 5 in S.Reals
     True
     >>> Rational(-1, 2) in S.Reals
diff --git a/sympy/sets/ordinals.py b/sympy/sets/ordinals.py
index 9326163c104f..3b1eadb28880 100644
--- a/sympy/sets/ordinals.py
+++ b/sympy/sets/ordinals.py
@@ -58,7 +58,7 @@ class Ordinal(Basic):
     Internally, this class is just a list of instances of OmegaPower
     Examples
     ========
-    >>> from sympy.sets import Ordinal, ord0, OmegaPower
+    >>> from sympy.sets import Ordinal, OmegaPower
     >>> from sympy.sets.ordinals import omega
     >>> w = omega
     >>> w.is_limit_ordinal
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index 1e3c9d292fca..7788e7a9fb65 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -479,11 +479,10 @@ def powerset(self):
         Examples
         ========
 
-        >>> from sympy import EmptySet, FiniteSet, Interval, PowerSet
+        >>> from sympy import EmptySet, FiniteSet, Interval
 
         A power set of an empty set:
 
-        >>> from sympy import FiniteSet, EmptySet
         >>> A = EmptySet
         >>> A.powerset()
         FiniteSet(EmptySet)
@@ -814,7 +813,7 @@ def is_iterable(self):
         Examples
         ========
 
-        >>> from sympy import FiniteSet, Interval, ProductSet
+        >>> from sympy import FiniteSet, Interval
         >>> I = Interval(0, 1)
         >>> A = FiniteSet(1, 2, 3, 4, 5)
         >>> I.is_iterable
@@ -2200,8 +2199,8 @@ def imageset(*args):
     Examples
     ========
 
-    >>> from sympy import S, Interval, Symbol, imageset, sin, Lambda
-    >>> from sympy.abc import x, y
+    >>> from sympy import S, Interval, imageset, sin, Lambda
+    >>> from sympy.abc import x
 
     >>> imageset(x, 2*x, Interval(0, 2))
     Interval(0, 4)
diff --git a/sympy/simplify/fu.py b/sympy/simplify/fu.py
index d0e85f759719..d564357c0403 100644
--- a/sympy/simplify/fu.py
+++ b/sympy/simplify/fu.py
@@ -451,7 +451,6 @@ def TR4(rv):
     Examples
     ========
 
-    >>> from sympy.simplify.fu import TR4
     >>> from sympy import pi
     >>> from sympy import cos, sin, tan, cot
     >>> for s in (0, pi/6, pi/4, pi/3, pi/2):
@@ -599,7 +598,7 @@ def TR8(rv, first=True):
     Examples
     ========
 
-    >>> from sympy.simplify.fu import TR8, TR7
+    >>> from sympy.simplify.fu import TR8
     >>> from sympy import cos, sin
     >>> TR8(cos(2)*cos(3))
     cos(5)/2 + cos(1)/2
@@ -816,8 +815,8 @@ def TR10i(rv):
     ========
 
     >>> from sympy.simplify.fu import TR10i
-    >>> from sympy import cos, sin, pi, Add, Mul, sqrt, Symbol
-    >>> from sympy.abc import x, y
+    >>> from sympy import cos, sin, sqrt
+    >>> from sympy.abc import x
 
     >>> TR10i(cos(1)*cos(3) + sin(1)*sin(3))
     cos(2)
@@ -1088,7 +1087,6 @@ def TR12(rv, first=True):
     Examples
     ========
 
-    >>> from sympy.simplify.fu import TR12
     >>> from sympy.abc import x, y
     >>> from sympy import tan
     >>> from sympy.simplify.fu import TR12
@@ -1250,7 +1248,7 @@ def TR13(rv):
     ========
 
     >>> from sympy.simplify.fu import TR13
-    >>> from sympy import tan, cot, cos
+    >>> from sympy import tan, cot
     >>> TR13(tan(3)*tan(2))
     -tan(2)/tan(5) - tan(3)/tan(5) + 1
     >>> TR13(cot(3)*cot(2))
@@ -1541,7 +1539,7 @@ def TR15(rv, max=4, pow=False):
 
     >>> from sympy.simplify.fu import TR15
     >>> from sympy.abc import x
-    >>> from sympy import cos, sin
+    >>> from sympy import sin
     >>> TR15(1 - 1/sin(x)**2)
     -cot(x)**2
 
@@ -1570,7 +1568,7 @@ def TR16(rv, max=4, pow=False):
 
     >>> from sympy.simplify.fu import TR16
     >>> from sympy.abc import x
-    >>> from sympy import cos, sin
+    >>> from sympy import cos
     >>> TR16(1 - 1/cos(x)**2)
     -tan(x)**2
 
diff --git a/sympy/simplify/gammasimp.py b/sympy/simplify/gammasimp.py
index f4659c8aa6d4..b65d03e9d9bc 100644
--- a/sympy/simplify/gammasimp.py
+++ b/sympy/simplify/gammasimp.py
@@ -46,7 +46,7 @@ def gammasimp(expr):
     ========
 
     >>> from sympy.simplify import gammasimp
-    >>> from sympy import gamma, factorial, Symbol
+    >>> from sympy import gamma, Symbol
     >>> from sympy.abc import x
     >>> n = Symbol('n', integer = True)
 
diff --git a/sympy/simplify/hyperexpand.py b/sympy/simplify/hyperexpand.py
index bb1efc04ab83..b8b1fdd44181 100644
--- a/sympy/simplify/hyperexpand.py
+++ b/sympy/simplify/hyperexpand.py
@@ -637,7 +637,7 @@ def compute_buckets(self):
 
         >>> from sympy.simplify.hyperexpand import G_Function
         >>> from sympy.abc import y
-        >>> from sympy import S, symbols
+        >>> from sympy import S
 
         >>> a, b = [1, 3, 2, S(3)/2], [1 + y, y, 2, y + 3]
         >>> G_Function(a, b, [2], [y]).compute_buckets()
diff --git a/sympy/simplify/powsimp.py b/sympy/simplify/powsimp.py
index baed0aac8359..280bd9e98fe6 100644
--- a/sympy/simplify/powsimp.py
+++ b/sympy/simplify/powsimp.py
@@ -77,7 +77,7 @@ def powsimp(expr, deep=False, combine='all', force=False, measure=count_ops):
 
     Radicals with Mul bases will be combined if combine='exp'
 
-    >>> from sympy import sqrt, Mul
+    >>> from sympy import sqrt
     >>> x, y = symbols('x y')
 
     Two radicals are automatically joined through Mul:
diff --git a/sympy/simplify/radsimp.py b/sympy/simplify/radsimp.py
index 945792100876..9d96b18e3d0b 100644
--- a/sympy/simplify/radsimp.py
+++ b/sympy/simplify/radsimp.py
@@ -46,7 +46,7 @@ def collect(expr, syms, func=None, evaluate=None, exact=False, distribute_order_
     ========
 
     >>> from sympy import S, collect, expand, factor, Wild
-    >>> from sympy.abc import a, b, c, x, y, z
+    >>> from sympy.abc import a, b, c, x, y
 
     This function can collect symbolic coefficients in polynomials or
     rational expressions. It will manage to find all integer or rational
@@ -634,7 +634,7 @@ def collect_const(expr, *vars, **kwargs):
     ========
 
     >>> from sympy import sqrt
-    >>> from sympy.abc import a, s, x, y, z
+    >>> from sympy.abc import s, x, y, z
     >>> from sympy.simplify.radsimp import collect_const
     >>> collect_const(sqrt(3) + sqrt(3)*(1 + sqrt(2)))
     sqrt(3)*(sqrt(2) + 2)
@@ -757,7 +757,7 @@ def radsimp(expr, symbolic=True, max_terms=4):
     Examples
     ========
 
-    >>> from sympy import radsimp, sqrt, Symbol, denom, pprint, I
+    >>> from sympy import radsimp, sqrt, Symbol, pprint
     >>> from sympy import factor_terms, fraction, signsimp
     >>> from sympy.simplify.radsimp import collect_sqrt
     >>> from sympy.abc import a, b, c
diff --git a/sympy/simplify/simplify.py b/sympy/simplify/simplify.py
index 828db06ac083..b79aba30f508 100644
--- a/sympy/simplify/simplify.py
+++ b/sympy/simplify/simplify.py
@@ -787,7 +787,7 @@ def factor_sum(self, limits=None, radical=False, clear=False, fraction=False, si
     Examples
     ========
 
-    >>> from sympy import Sum, Integral
+    >>> from sympy import Sum
     >>> from sympy.abc import x, y
     >>> from sympy.simplify.simplify import factor_sum
     >>> s = Sum(x*y, (x, 1, 3))
@@ -1169,8 +1169,8 @@ def kroneckersimp(expr):
     The only simplification currently attempted is to identify multiplicative cancellation:
 
     >>> from sympy import KroneckerDelta, kroneckersimp
-    >>> from sympy.abc import i, j
-    >>> kroneckersimp(1 + KroneckerDelta(0, j) * KroneckerDelta(1, j))
+    >>> from sympy.abc import i
+    >>> kroneckersimp(1 + KroneckerDelta(0, i) * KroneckerDelta(1, i))
     1
     """
     def args_cancel(args1, args2):
@@ -1346,7 +1346,7 @@ def nthroot(expr, n, max_len=4, prec=15):
     ========
 
     >>> from sympy.simplify.simplify import nthroot
-    >>> from sympy import Rational, sqrt
+    >>> from sympy import sqrt
     >>> nthroot(90 + 34*sqrt(7), 3)
     sqrt(7) + 3
 
@@ -1414,7 +1414,7 @@ def nsimplify(expr, constants=(), tolerance=None, full=False, rational=None,
     Examples
     ========
 
-    >>> from sympy import nsimplify, sqrt, GoldenRatio, exp, I, exp, pi
+    >>> from sympy import nsimplify, sqrt, GoldenRatio, exp, I, pi
     >>> nsimplify(4/(1+sqrt(5)), [GoldenRatio])
     -2 + 2*GoldenRatio
     >>> nsimplify((1/(exp(3*pi*I/5)+1)))
@@ -1527,7 +1527,6 @@ def _real_to_rational(expr, tolerance=None, rational_conversion='base10'):
     Examples
     ========
 
-    >>> from sympy import Rational
     >>> from sympy.simplify.simplify import _real_to_rational
     >>> from sympy.abc import x
 
diff --git a/sympy/simplify/trigsimp.py b/sympy/simplify/trigsimp.py
index 60f795ca57a2..d1a11d60d03c 100644
--- a/sympy/simplify/trigsimp.py
+++ b/sympy/simplify/trigsimp.py
@@ -444,7 +444,7 @@ def trigsimp(expr, **opts):
     ========
 
     >>> from sympy import trigsimp, sin, cos, log
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
     >>> e = 2*sin(x)**2 + 2*cos(x)**2
     >>> trigsimp(e)
     2
@@ -638,8 +638,8 @@ def trigsimp_old(expr, **opts):
     Examples
     ========
 
-    >>> from sympy import trigsimp, sin, cos, log, cosh, sinh, tan, cot
-    >>> from sympy.abc import x, y
+    >>> from sympy import trigsimp, sin, cos, log, cot
+    >>> from sympy.abc import x
     >>> e = 2*sin(x)**2 + 2*cos(x)**2
     >>> trigsimp(e, old=True)
     2
diff --git a/sympy/solvers/bivariate.py b/sympy/solvers/bivariate.py
index 19994fd43327..4442c9814e8a 100644
--- a/sympy/solvers/bivariate.py
+++ b/sympy/solvers/bivariate.py
@@ -59,7 +59,6 @@ def _mostfunc(lhs, func, X=None):
 
     >>> from sympy.solvers.bivariate import _mostfunc
     >>> from sympy.functions.elementary.exponential import exp
-    >>> from sympy.testing.pytest import raises
     >>> from sympy.abc import x, y
     >>> _mostfunc(exp(x) + exp(exp(x) + 2), exp)
     exp(exp(x) + 2)
diff --git a/sympy/solvers/deutils.py b/sympy/solvers/deutils.py
index 9bdbfad82019..6a45534b9c2a 100644
--- a/sympy/solvers/deutils.py
+++ b/sympy/solvers/deutils.py
@@ -27,7 +27,7 @@ def _preprocess(expr, func=None, hint='_Integral'):
     function to be solved for.
 
     >>> from sympy.solvers.deutils import _preprocess
-    >>> from sympy import Derivative, Function, Integral, sin
+    >>> from sympy import Derivative, Function
     >>> from sympy.abc import x, y, z
     >>> f, g = map(Function, 'fg')
 
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index c67eabfe927f..4e2c7afdbc45 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -867,7 +867,7 @@ def diop_linear(eq, param=symbols("t", integer=True)):
     ========
 
     >>> from sympy.solvers.diophantine.diophantine import diop_linear
-    >>> from sympy.abc import x, y, z, t
+    >>> from sympy.abc import x, y, z
     >>> diop_linear(2*x - 3*y - 5) # solves equation 2*x - 3*y - 5 == 0
     (3*t_0 - 5, 2*t_0 - 5)
 
@@ -2026,7 +2026,6 @@ def transformation_to_DN(eq):
 
     >>> from sympy.abc import x, y
     >>> from sympy.solvers.diophantine.diophantine import transformation_to_DN
-    >>> from sympy.solvers.diophantine import classify_diop
     >>> A, B = transformation_to_DN(x**2 - 3*x*y - y**2 - 2*y + 1)
     >>> A
     Matrix([
@@ -3129,7 +3128,7 @@ def diop_general_sum_of_squares(eq, limit=1):
     ========
 
     >>> from sympy.solvers.diophantine.diophantine import diop_general_sum_of_squares
-    >>> from sympy.abc import a, b, c, d, e, f
+    >>> from sympy.abc import a, b, c, d, e
     >>> diop_general_sum_of_squares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345)
     {(15, 22, 22, 24, 24)}
 
diff --git a/sympy/solvers/inequalities.py b/sympy/solvers/inequalities.py
index 01341ea54b5a..5ce0e102e012 100644
--- a/sympy/solvers/inequalities.py
+++ b/sympy/solvers/inequalities.py
@@ -201,7 +201,7 @@ def reduce_rational_inequalities(exprs, gen, relational=True):
     Examples
     ========
 
-    >>> from sympy import Poly, Symbol
+    >>> from sympy import Symbol
     >>> from sympy.solvers.inequalities import reduce_rational_inequalities
 
     >>> x = Symbol('x', real=True)
@@ -374,7 +374,6 @@ def reduce_abs_inequalities(exprs, gen):
     ========
 
     >>> from sympy import Abs, Symbol
-    >>> from sympy.abc import x
     >>> from sympy.solvers.inequalities import reduce_abs_inequalities
     >>> x = Symbol('x', extended_real=True)
 
@@ -941,7 +940,6 @@ def reduce_inequalities(inequalities, symbols=[]):
     Examples
     ========
 
-    >>> from sympy import sympify as S, Symbol
     >>> from sympy.abc import x, y
     >>> from sympy.solvers.inequalities import reduce_inequalities
 
diff --git a/sympy/solvers/ode/ode.py b/sympy/solvers/ode/ode.py
index c0d3eb1624ed..4a28f813f4a6 100644
--- a/sympy/solvers/ode/ode.py
+++ b/sympy/solvers/ode/ode.py
@@ -1818,9 +1818,8 @@ def classify_sysode(eq, funcs=None, **kwargs):
     Examples
     ========
 
-    >>> from sympy import Function, Eq, symbols, diff, Rational
+    >>> from sympy import Function, Eq, symbols, diff
     >>> from sympy.solvers.ode.ode import classify_sysode
-    >>> from sympy.matrices.dense import Matrix
     >>> from sympy.abc import t
     >>> f, x, y = symbols('f, x, y', cls=Function)
     >>> k, l, m, n = symbols('k, l, m, n', Integer=True)
@@ -2537,7 +2536,7 @@ def ode_sol_simplicity(sol, func, trysolving=True):
     such as ``min(listofsolutions, key=lambda i: ode_sol_simplicity(i,
     f(x)))``.
 
-    >>> from sympy import symbols, Function, Eq, tan, cos, sqrt, Integral
+    >>> from sympy import symbols, Function, Eq, tan, Integral
     >>> from sympy.solvers.ode.ode import ode_sol_simplicity
     >>> x, C1, C2 = symbols('x, C1, C2')
     >>> f = Function('f')
@@ -2838,7 +2837,7 @@ def constant_renumber(expr, variables=None, newconstants=None):
     Examples
     ========
 
-    >>> from sympy import symbols, Eq, pprint
+    >>> from sympy import symbols
     >>> from sympy.solvers.ode.ode import constant_renumber
     >>> x, C1, C2, C3 = symbols('x,C1:4')
     >>> expr = C3 + C2*x + C1*x**2
@@ -3466,7 +3465,7 @@ def ode_2nd_power_series_ordinary(eq, func, order, match):
     ========
 
     >>> from sympy import dsolve, Function, pprint
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
     >>> f = Function("f")
     >>> eq = f(x).diff(x, 2) + f(x)
     >>> pprint(dsolve(eq, hint='2nd_power_series_ordinary'))
@@ -3603,7 +3602,7 @@ def ode_2nd_linear_airy(eq, func, order, match):
     Examples
     ========
 
-    >>> from sympy import dsolve, Function, pprint
+    >>> from sympy import dsolve, Function
     >>> from sympy.abc import x
     >>> f = Function("f")
     >>> eq = f(x).diff(x, 2) - x*f(x)
@@ -3660,7 +3659,7 @@ def ode_2nd_power_series_regular(eq, func, order, match):
     ========
 
     >>> from sympy import dsolve, Function, pprint
-    >>> from sympy.abc import x, y
+    >>> from sympy.abc import x
     >>> f = Function("f")
     >>> eq = x*(f(x).diff(x, 2)) + 2*(f(x).diff(x)) + x*f(x)
     >>> pprint(dsolve(eq, hint='2nd_power_series_regular'))
@@ -3761,11 +3760,11 @@ def ode_2nd_linear_bessel(eq, func, order, match):
     Examples
     ========
 
-    >>> from sympy.abc import x, y, a
+    >>> from sympy.abc import x
     >>> from sympy import Symbol
     >>> v = Symbol('v', positive=True)
-    >>> from sympy.solvers.ode import dsolve, checkodesol
-    >>> from sympy import pprint, Function
+    >>> from sympy.solvers.ode import dsolve
+    >>> from sympy import Function
     >>> f = Function('f')
     >>> y = f(x)
     >>> genform = x**2*y.diff(x, 2) + x*y.diff(x) + (x**2 - v**2)*y
@@ -4104,7 +4103,7 @@ def ode_nth_linear_euler_eq_homogeneous(eq, func, order, match, returns='sol'):
     :py:obj:`~.ComplexRootOf` instance will be returned
     instead.
 
-    >>> from sympy import Function, dsolve, Eq
+    >>> from sympy import Function, dsolve
     >>> from sympy.abc import x
     >>> f = Function('f')
     >>> dsolve(4*x**2*f(x).diff(x, 2) + f(x), f(x),
@@ -4470,8 +4469,8 @@ def ode_linear_coefficients(eq, func, order, match):
     Examples
     ========
 
-    >>> from sympy import Function, Derivative, pprint
-    >>> from sympy.solvers.ode.ode import dsolve, classify_ode
+    >>> from sympy import Function, pprint
+    >>> from sympy.solvers.ode.ode import dsolve
     >>> from sympy.abc import x
     >>> f = Function('f')
     >>> df = f(x).diff(x)
@@ -4508,7 +4507,7 @@ def ode_separable_reduced(eq, func, order, match):
 
     The general solution is:
 
-        >>> from sympy import Function, dsolve, Eq, pprint
+        >>> from sympy import Function, dsolve, pprint
         >>> from sympy.abc import x, n
         >>> f, g = map(Function, ['f', 'g'])
         >>> genform = f(x).diff(x) + (f(x)/x)*g(x**n*f(x))
@@ -4535,8 +4534,8 @@ def ode_separable_reduced(eq, func, order, match):
     Examples
     ========
 
-    >>> from sympy import Function, Derivative, pprint
-    >>> from sympy.solvers.ode.ode import dsolve, classify_ode
+    >>> from sympy import Function, pprint
+    >>> from sympy.solvers.ode.ode import dsolve
     >>> from sympy.abc import x
     >>> f = Function('f')
     >>> d = f(x).diff(x)
@@ -4591,7 +4590,7 @@ def ode_1st_power_series(eq, func, order, match):
     Examples
     ========
 
-    >>> from sympy import Function, Derivative, pprint, exp
+    >>> from sympy import Function, pprint, exp
     >>> from sympy.solvers.ode.ode import dsolve
     >>> from sympy.abc import x
     >>> f = Function('f')
@@ -4672,7 +4671,7 @@ def ode_nth_linear_constant_coeff_homogeneous(eq, func, order, match,
     :py:class:`~sympy.polys.rootoftools.ComplexRootOf` instance will be return
     instead.
 
-    >>> from sympy import Function, dsolve, Eq
+    >>> from sympy import Function, dsolve
     >>> from sympy.abc import x
     >>> f = Function('f')
     >>> dsolve(f(x).diff(x, 5) + 10*f(x).diff(x) - 2*f(x), f(x),
@@ -5508,7 +5507,7 @@ def ode_lie_group(eq, func, order, match):
     Examples
     ========
 
-    >>> from sympy import Function, dsolve, Eq, exp, pprint
+    >>> from sympy import Function, dsolve, exp, pprint
     >>> from sympy.abc import x
     >>> f = Function('f')
     >>> pprint(dsolve(f(x).diff(x) + 2*x*f(x) - x*exp(-x**2), f(x),
@@ -5621,7 +5620,7 @@ def infinitesimals(eq, func=None, order=None, hint='default', match=None):
 
     The infinitesimals can be found by solving the following PDE:
 
-        >>> from sympy import Function, diff, Eq, pprint
+        >>> from sympy import Function, Eq, pprint
         >>> from sympy.abc import x, y
         >>> xi, eta, h = map(Function, ['xi', 'eta', 'h'])
         >>> h = h(x, y)  # dy/dx = h
@@ -5649,7 +5648,7 @@ def infinitesimals(eq, func=None, order=None, hint='default', match=None):
     Examples
     ========
 
-    >>> from sympy import Function, diff
+    >>> from sympy import Function
     >>> from sympy.solvers.ode.ode import infinitesimals
     >>> from sympy.abc import x
     >>> f = Function('f')
diff --git a/sympy/solvers/ode/single.py b/sympy/solvers/ode/single.py
index 3a00e3320cf2..663132dc11ae 100644
--- a/sympy/solvers/ode/single.py
+++ b/sympy/solvers/ode/single.py
@@ -55,7 +55,7 @@ class SingleODEProblem:
     This class is used internally by dsolve. To instantiate an instance
     directly first define an ODE problem:
 
-    >>> from sympy import Eq, Function, Symbol
+    >>> from sympy import Function, Symbol
     >>> x = Symbol('x')
     >>> f = Function('f')
     >>> eq = f(x).diff(x, 2)
@@ -152,7 +152,7 @@ class SingleODESolver:
     You can use a subclass of SingleODEProblem to solve a particular type of
     ODE. We first define a particular ODE problem:
 
-    >>> from sympy import Eq, Function, Symbol
+    >>> from sympy import Function, Symbol
     >>> x = Symbol('x')
     >>> f = Function('f')
     >>> eq = f(x).diff(x, 2)
@@ -466,8 +466,8 @@ class AlmostLinear(SinglePatternODESolver):
     Examples
     ========
 
-    >>> from sympy import Function, Derivative, pprint, sin, cos
-    >>> from sympy.solvers.ode import dsolve, classify_ode
+    >>> from sympy import Function, pprint, sin, cos
+    >>> from sympy.solvers.ode import dsolve
     >>> from sympy.abc import x
     >>> f = Function('f')
     >>> d = f(x).diff(x)
@@ -650,7 +650,7 @@ class Factorable(SingleODESolver):
         Examples
         ========
 
-        >>> from sympy import Function, dsolve, Eq, pprint, Derivative
+        >>> from sympy import Function, dsolve, pprint
         >>> from sympy.abc import x
         >>> f = Function('f')
         >>> eq = (f(x)**2-4)*(f(x).diff(x)+f(x))
@@ -727,7 +727,7 @@ class RiccatiSpecial(SinglePatternODESolver):
     and is valid when neither `a` nor `b` are zero and either `c` or `d` is
     zero.
 
-    >>> from sympy.abc import x, y, a, b, c, d
+    >>> from sympy.abc import x, a, b, c, d
     >>> from sympy.solvers.ode import dsolve, checkodesol
     >>> from sympy import pprint, Function
     >>> f = Function('f')
diff --git a/sympy/solvers/pde.py b/sympy/solvers/pde.py
index b0946d881791..834d1109b885 100644
--- a/sympy/solvers/pde.py
+++ b/sympy/solvers/pde.py
@@ -151,7 +151,7 @@ def pdsolve(eq, func=None, hint='default', dict=False, solvefun=None, **kwargs):
     ========
 
     >>> from sympy.solvers.pde import pdsolve
-    >>> from sympy import Function, diff, Eq
+    >>> from sympy import Function, Eq
     >>> from sympy.abc import x, y
     >>> f = Function('f')
     >>> u = f(x, y)
@@ -256,7 +256,7 @@ def classify_pde(eq, func=None, dict=False, **kwargs):
     ========
 
     >>> from sympy.solvers.pde import classify_pde
-    >>> from sympy import Function, diff, Eq
+    >>> from sympy import Function, Eq
     >>> from sympy.abc import x, y
     >>> f = Function('f')
     >>> u = f(x, y)
@@ -418,7 +418,7 @@ def checkpdesol(pde, sol, func=None, solve_for_func=True):
     Examples
     ========
 
-    >>> from sympy import Function, symbols, diff
+    >>> from sympy import Function, symbols
     >>> from sympy.solvers.pde import checkpdesol, pdsolve
     >>> x, y = symbols('x y')
     >>> f = Function('f')
@@ -525,10 +525,8 @@ def pde_1st_linear_constant_coeff_homogeneous(eq, func, order, match, solvefun):
     Examples
     ========
 
-    >>> from sympy.solvers.pde import (
-    ... pde_1st_linear_constant_coeff_homogeneous)
     >>> from sympy import pdsolve
-    >>> from sympy import Function, diff, pprint
+    >>> from sympy import Function, pprint
     >>> from sympy.abc import x,y
     >>> f = Function('f')
     >>> pdsolve(f(x,y) + f(x,y).diff(x) + f(x,y).diff(y))
@@ -639,7 +637,7 @@ def pde_1st_linear_constant_coeff(eq, func, order, match, solvefun):
     ========
 
     >>> from sympy.solvers.pde import pdsolve
-    >>> from sympy import Function, diff, pprint, exp
+    >>> from sympy import Function, pprint, exp
     >>> from sympy.abc import x,y
     >>> f = Function('f')
     >>> eq = -2*f(x,y).diff(x) + 4*f(x,y).diff(y) + 5*f(x,y) - exp(x + 3*y)
@@ -701,7 +699,6 @@ def pde_1st_linear_variable_coeff(eq, func, order, match, solvefun):
 
     which can be solved using ``dsolve``.
 
-    >>> from sympy.solvers.pde import pdsolve
     >>> from sympy.abc import x, y
     >>> from sympy import Function, pprint
     >>> a, b, c, G, f= [Function(i) for i in ['a', 'b', 'c', 'G', 'f']]
@@ -719,7 +716,7 @@ def pde_1st_linear_variable_coeff(eq, func, order, match, solvefun):
     ========
 
     >>> from sympy.solvers.pde import pdsolve
-    >>> from sympy import Function, diff, pprint, exp
+    >>> from sympy import Function, pprint
     >>> from sympy.abc import x,y
     >>> f = Function('f')
     >>> eq =  x*(u.diff(x)) - y*(u.diff(y)) + y**2*u - y**2
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index b131fb46ae66..914cd8b2c701 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -137,7 +137,6 @@ def denoms(eq, *symbols):
 
     >>> from sympy.solvers.solvers import denoms
     >>> from sympy.abc import x, y, z
-    >>> from sympy import sqrt
 
     >>> denoms(x/y)
     {y}
@@ -2866,7 +2865,6 @@ def nsolve(*args, **kwargs):
     ========
 
     >>> from sympy import Symbol, nsolve
-    >>> import sympy
     >>> import mpmath
     >>> mpmath.mp.dps = 15
     >>> x1 = Symbol('x1')
@@ -2910,7 +2908,6 @@ def nsolve(*args, **kwargs):
     independently verify the solution.
 
     >>> from sympy import cos, cosh
-    >>> from sympy.abc import i
     >>> f = cos(x)*cosh(x) - 1
     >>> nsolve(f, 3.14*100)
     Traceback (most recent call last):
@@ -3283,7 +3280,7 @@ def unrad(eq, *syms, **flags):
 
     >>> from sympy.solvers.solvers import unrad
     >>> from sympy.abc import x
-    >>> from sympy import sqrt, Rational, root, real_roots, solve
+    >>> from sympy import sqrt, Rational, root
 
     >>> unrad(sqrt(x)*x**Rational(1, 3) + 2)
     (x**5 - 64, [])
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index c07a58e2bfef..1e35500be49c 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -134,7 +134,7 @@ def _invert(f_x, y, x, domain=S.Complexes):
 
     >>> from sympy.solvers.solveset import invert_complex, invert_real
     >>> from sympy.abc import x, y
-    >>> from sympy import exp, log
+    >>> from sympy import exp
 
     When does exp(x) == y?
 
@@ -1260,7 +1260,7 @@ def _solve_modular(f, symbol, domain):
     ========
 
     >>> from sympy.solvers.solveset import _solve_modular as solve_modulo
-    >>> from sympy import S, Symbol, sin, Intersection, Range, Interval
+    >>> from sympy import S, Symbol, sin, Intersection, Interval
     >>> from sympy.core.mod import Mod
     >>> x = Symbol('x')
     >>> solve_modulo(Mod(5*x - 8, 7) - 3, x, S.Integers)
@@ -2127,7 +2127,7 @@ def solvify(f, symbol, domain):
     Examples
     ========
 
-    >>> from sympy.solvers.solveset import solvify, solveset
+    >>> from sympy.solvers.solveset import solvify
     >>> from sympy.abc import x
     >>> from sympy import S, tan, sin, exp
     >>> solvify(x**2 - 9, x, S.Reals)
@@ -2443,7 +2443,7 @@ def linsolve(system, *symbols):
     Examples
     ========
 
-    >>> from sympy import Matrix, S, linsolve, symbols
+    >>> from sympy import Matrix, linsolve, symbols
     >>> x, y, z = symbols("x, y, z")
     >>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
     >>> b = Matrix([3, 6, 9])
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index a3e11443b616..ce46b146c8c2 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -252,7 +252,7 @@ def Arcsin(name, a=0, b=1):
     ========
 
     >>> from sympy.stats import Arcsin, density, cdf
-    >>> from sympy import Symbol, simplify
+    >>> from sympy import Symbol
 
     >>> a = Symbol("a", real=True)
     >>> b = Symbol("b", real=True)
@@ -334,7 +334,7 @@ def Benini(name, alpha, beta, sigma):
     ========
 
     >>> from sympy.stats import Benini, density, cdf
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy import Symbol, pprint
 
     >>> alpha = Symbol("alpha", positive=True)
     >>> beta = Symbol("beta", positive=True)
@@ -1403,8 +1403,8 @@ def ExponentialPower(name, mu, alpha, beta):
     Examples
     ========
 
-    >>> from sympy.stats import ExponentialPower, density, E, variance, cdf
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy.stats import ExponentialPower, density, cdf
+    >>> from sympy import Symbol, pprint
     >>> z = Symbol("z")
     >>> mu = Symbol("mu")
     >>> alpha = Symbol("alpha", positive=True)
@@ -1486,7 +1486,7 @@ def FDistribution(name, d1, d2):
     ========
 
     >>> from sympy.stats import FDistribution, density
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy import Symbol, pprint
 
     >>> d1 = Symbol("d1", positive=True)
     >>> d2 = Symbol("d2", positive=True)
@@ -1561,7 +1561,7 @@ def FisherZ(name, d1, d2):
     ========
 
     >>> from sympy.stats import FisherZ, density
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy import Symbol, pprint
 
     >>> d1 = Symbol("d1", positive=True)
     >>> d2 = Symbol("d2", positive=True)
@@ -1644,8 +1644,8 @@ def Frechet(name, a, s=1, m=0):
     Examples
     ========
 
-    >>> from sympy.stats import Frechet, density, E, std, cdf
-    >>> from sympy import Symbol, simplify
+    >>> from sympy.stats import Frechet, density, cdf
+    >>> from sympy import Symbol
 
     >>> a = Symbol("a", positive=True)
     >>> s = Symbol("s", positive=True)
@@ -1828,7 +1828,7 @@ def GammaInverse(name, a, b):
     Examples
     ========
 
-    >>> from sympy.stats import GammaInverse, density, cdf, E, variance
+    >>> from sympy.stats import GammaInverse, density, cdf
     >>> from sympy import Symbol, pprint
 
     >>> a = Symbol("a", positive=True)
@@ -1933,8 +1933,8 @@ def Gumbel(name, beta, mu, minimum=False):
     Examples
     ========
 
-    >>> from sympy.stats import Gumbel, density, E, variance, cdf
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy.stats import Gumbel, density, cdf
+    >>> from sympy import Symbol
     >>> x = Symbol("x")
     >>> mu = Symbol("mu")
     >>> beta = Symbol("beta", positive=True)
@@ -2005,8 +2005,8 @@ def Gompertz(name, b, eta):
     Examples
     ========
 
-    >>> from sympy.stats import Gompertz, density, E, variance
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy.stats import Gompertz, density
+    >>> from sympy import Symbol
 
     >>> b = Symbol("b", positive=True)
     >>> eta = Symbol("eta", positive=True)
@@ -2075,8 +2075,8 @@ def Kumaraswamy(name, a, b):
     Examples
     ========
 
-    >>> from sympy.stats import Kumaraswamy, density, E, variance, cdf
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy.stats import Kumaraswamy, density, cdf
+    >>> from sympy import Symbol, pprint
 
     >>> a = Symbol("a", positive=True)
     >>> b = Symbol("b", positive=True)
@@ -2251,7 +2251,7 @@ def Levy(name, mu, c):
     ========
 
     >>> from sympy.stats import Levy, density, cdf
-    >>> from sympy import Symbol, pprint
+    >>> from sympy import Symbol
 
     >>> mu = Symbol("mu", real=True)
     >>> c = Symbol("c", positive=True)
@@ -2498,7 +2498,7 @@ def LogNormal(name, mean, std):
     ========
 
     >>> from sympy.stats import LogNormal, density
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy import Symbol, pprint
 
     >>> mu = Symbol("mu", real=True)
     >>> sigma = Symbol("sigma", positive=True)
@@ -2897,7 +2897,7 @@ def Normal(name, mean, std):
     ========
 
     >>> from sympy.stats import Normal, density, E, std, cdf, skewness, quantile
-    >>> from sympy import Symbol, simplify, pprint, factor, together, factor_terms
+    >>> from sympy import Symbol, simplify, pprint
 
     >>> mu = Symbol("mu")
     >>> sigma = Symbol("sigma", positive=True)
@@ -3027,7 +3027,7 @@ def GaussianInverse(name, mean, shape):
     Examples
     ========
 
-    >>> from sympy.stats import GaussianInverse, density, cdf, E, std, skewness
+    >>> from sympy.stats import GaussianInverse, density, E, std, skewness
     >>> from sympy import Symbol, pprint
 
     >>> mu = Symbol("mu", positive=True)
@@ -3202,7 +3202,7 @@ def PowerFunction(name, alpha, a, b):
     ========
 
     >>> from sympy.stats import PowerFunction, density, cdf, E, variance
-    >>> from sympy import Symbol, simplify
+    >>> from sympy import Symbol
     >>> alpha = Symbol("alpha", positive=True)
     >>> a = Symbol("a", real=True)
     >>> b = Symbol("b", real=True)
@@ -3295,8 +3295,8 @@ def QuadraticU(name, a, b):
     Examples
     ========
 
-    >>> from sympy.stats import QuadraticU, density, E, variance
-    >>> from sympy import Symbol, simplify, factor, pprint
+    >>> from sympy.stats import QuadraticU, density
+    >>> from sympy import Symbol, pprint
 
     >>> a = Symbol("a", real=True)
     >>> b = Symbol("b", real=True)
@@ -3381,8 +3381,8 @@ def RaisedCosine(name, mu, s):
     Examples
     ========
 
-    >>> from sympy.stats import RaisedCosine, density, E, variance
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy.stats import RaisedCosine, density
+    >>> from sympy import Symbol, pprint
 
     >>> mu = Symbol("mu", real=True)
     >>> s = Symbol("s", positive=True)
@@ -3464,7 +3464,7 @@ def Rayleigh(name, sigma):
     ========
 
     >>> from sympy.stats import Rayleigh, density, E, variance
-    >>> from sympy import Symbol, simplify
+    >>> from sympy import Symbol
 
     >>> sigma = Symbol("sigma", positive=True)
     >>> z = Symbol("z")
@@ -3590,7 +3590,7 @@ def ShiftedGompertz(name, b, eta):
 
     Examples
     ========
-    >>> from sympy.stats import ShiftedGompertz, density, E, variance
+    >>> from sympy.stats import ShiftedGompertz, density
     >>> from sympy import Symbol
 
     >>> b = Symbol("b", positive=True)
@@ -3659,8 +3659,8 @@ def StudentT(name, nu):
     Examples
     ========
 
-    >>> from sympy.stats import StudentT, density, E, variance, cdf
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy.stats import StudentT, density, cdf
+    >>> from sympy import Symbol, pprint
 
     >>> nu = Symbol("nu", positive=True)
     >>> z = Symbol("z")
@@ -3755,7 +3755,7 @@ def Trapezoidal(name, a, b, c, d):
     Examples
     ========
 
-    >>> from sympy.stats import Trapezoidal, density, E
+    >>> from sympy.stats import Trapezoidal, density
     >>> from sympy import Symbol, pprint
 
     >>> a = Symbol("a")
@@ -3854,7 +3854,7 @@ def Triangular(name, a, b, c):
     Examples
     ========
 
-    >>> from sympy.stats import Triangular, density, E
+    >>> from sympy.stats import Triangular, density
     >>> from sympy import Symbol, pprint
 
     >>> a = Symbol("a")
@@ -3966,7 +3966,7 @@ def Uniform(name, left, right):
     Examples
     ========
 
-    >>> from sympy.stats import Uniform, density, cdf, E, variance, skewness
+    >>> from sympy.stats import Uniform, density, cdf, E, variance
     >>> from sympy import Symbol, simplify
 
     >>> a = Symbol("a", negative=True)
@@ -4144,8 +4144,8 @@ def VonMises(name, mu, k):
     Examples
     ========
 
-    >>> from sympy.stats import VonMises, density, E, variance
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy.stats import VonMises, density
+    >>> from sympy import Symbol, pprint
 
     >>> mu = Symbol("mu")
     >>> k = Symbol("k", positive=True)
@@ -4297,7 +4297,7 @@ def WignerSemicircle(name, R):
     ========
 
     >>> from sympy.stats import WignerSemicircle, density, E
-    >>> from sympy import Symbol, simplify
+    >>> from sympy import Symbol
 
     >>> R = Symbol("R", positive=True)
     >>> z = Symbol("z")
diff --git a/sympy/stats/drv_types.py b/sympy/stats/drv_types.py
index 399336b015a3..6a30a2fe12b0 100644
--- a/sympy/stats/drv_types.py
+++ b/sympy/stats/drv_types.py
@@ -592,7 +592,7 @@ def Skellam(name, mu1, mu2):
     ========
 
     >>> from sympy.stats import Skellam, density, E, variance
-    >>> from sympy import Symbol, simplify, pprint
+    >>> from sympy import Symbol, pprint
 
     >>> z = Symbol("z", integer=True)
     >>> mu1 = Symbol("mu1", positive=True)
diff --git a/sympy/stats/joint_rv_types.py b/sympy/stats/joint_rv_types.py
index 756b97f1c490..2438237c7706 100644
--- a/sympy/stats/joint_rv_types.py
+++ b/sympy/stats/joint_rv_types.py
@@ -42,7 +42,7 @@ def JointRV(symbol, pdf, _set=None):
     Examples
     ========
 
-    >>> from sympy import symbols, exp, pi, Indexed, S
+    >>> from sympy import exp, pi, Indexed, S
     >>> from sympy.stats import density
     >>> from sympy.stats.joint_rv_types import JointRV
 
@@ -518,7 +518,6 @@ def GeneralizedMultivariateLogGamma(syms, delta, v, lamda, mu):
     ========
 
     >>> from sympy.stats import density
-    >>> from sympy.stats.joint_rv import marginal_distribution
     >>> from sympy.stats.joint_rv_types import GeneralizedMultivariateLogGamma
     >>> from sympy import symbols, S
     >>> v = 1
@@ -573,7 +572,6 @@ def GeneralizedMultivariateLogGammaOmega(syms, omega, v, lamda, mu):
     ========
 
     >>> from sympy.stats import density
-    >>> from sympy.stats.joint_rv import marginal_distribution
     >>> from sympy.stats.joint_rv_types import GeneralizedMultivariateLogGammaOmega
     >>> from sympy import Matrix, symbols, S
     >>> omega = Matrix([[1, S.Half, S.Half], [S.Half, 1, S.Half], [S.Half, S.Half, 1]])
diff --git a/sympy/stats/random_matrix_models.py b/sympy/stats/random_matrix_models.py
index 02717ee24a0b..a04e2f3bc104 100644
--- a/sympy/stats/random_matrix_models.py
+++ b/sympy/stats/random_matrix_models.py
@@ -234,7 +234,7 @@ class CircularUnitaryEnsemble(CircularEnsemble):
     Examples
     ========
 
-    >>> from sympy.stats import CircularUnitaryEnsemble as CUE, density
+    >>> from sympy.stats import CircularUnitaryEnsemble as CUE
     >>> from sympy.stats import joint_eigen_distribution
     >>> C = CUE('U', 1)
     >>> joint_eigen_distribution(C)
@@ -257,7 +257,7 @@ class CircularOrthogonalEnsemble(CircularEnsemble):
     Examples
     ========
 
-    >>> from sympy.stats import CircularOrthogonalEnsemble as COE, density
+    >>> from sympy.stats import CircularOrthogonalEnsemble as COE
     >>> from sympy.stats import joint_eigen_distribution
     >>> C = COE('O', 1)
     >>> joint_eigen_distribution(C)
@@ -280,7 +280,7 @@ class CircularSymplecticEnsemble(CircularEnsemble):
     Examples
     ========
 
-    >>> from sympy.stats import CircularSymplecticEnsemble as CSE, density
+    >>> from sympy.stats import CircularSymplecticEnsemble as CSE
     >>> from sympy.stats import joint_eigen_distribution
     >>> C = CSE('S', 1)
     >>> joint_eigen_distribution(C)
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 378e5fc0aba8..6537c7bd4df5 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -575,7 +575,6 @@ def pspace(expr):
     ========
 
     >>> from sympy.stats import pspace, Normal
-    >>> from sympy.stats.rv import IndependentProductPSpace
     >>> X = Normal('X', 0, 1)
     >>> pspace(2*X + 1) == X.pspace
     True
@@ -1004,7 +1003,7 @@ def where(condition, given_condition=None, **kwargs):
     ========
 
     >>> from sympy.stats import where, Die, Normal
-    >>> from sympy import symbols, And
+    >>> from sympy import And
 
     >>> D1, D2 = Die('a', 6), Die('b', 6)
     >>> a, b = D1.symbol, D2.symbol
diff --git a/sympy/stats/rv_interface.py b/sympy/stats/rv_interface.py
index 236f774a9441..d725e6f07ddc 100644
--- a/sympy/stats/rv_interface.py
+++ b/sympy/stats/rv_interface.py
@@ -44,7 +44,7 @@ def variance(X, condition=None, **kwargs):
     Examples
     ========
 
-    >>> from sympy.stats import Die, E, Bernoulli, variance
+    >>> from sympy.stats import Die, Bernoulli, variance
     >>> from sympy import simplify, Symbol
 
     >>> X = Die('X', 6)
diff --git a/sympy/stats/stochastic_process_types.py b/sympy/stats/stochastic_process_types.py
index 9380a06b92e0..9cc5ec33f272 100644
--- a/sympy/stats/stochastic_process_types.py
+++ b/sympy/stats/stochastic_process_types.py
@@ -744,7 +744,7 @@ class ContinuousMarkovChain(ContinuousTimeStochasticProcess, MarkovProcess):
     ========
 
     >>> from sympy.stats import ContinuousMarkovChain
-    >>> from sympy import Matrix, S, MatrixSymbol
+    >>> from sympy import Matrix, S
     >>> G = Matrix([[-S(1), S(1)], [S(1), -S(1)]])
     >>> C = ContinuousMarkovChain('C', state_space=[0, 1], gen_mat=G)
     >>> C.limiting_distribution()
@@ -820,7 +820,7 @@ class BernoulliProcess(DiscreteTimeStochasticProcess):
     ========
 
     >>> from sympy.stats import BernoulliProcess, P, E
-    >>> from sympy import Eq, Gt, Lt
+    >>> from sympy import Eq, Gt
     >>> B = BernoulliProcess("B", p=0.7, success=1, failure=0)
     >>> B.state_space
     FiniteSet(0, 1)
diff --git a/sympy/tensor/index_methods.py b/sympy/tensor/index_methods.py
index 08581cfd9230..372263a689ac 100644
--- a/sympy/tensor/index_methods.py
+++ b/sympy/tensor/index_methods.py
@@ -212,7 +212,7 @@ def get_indices(expr):
 
     >>> from sympy.tensor.index_methods import get_indices
     >>> from sympy import symbols
-    >>> from sympy.tensor import IndexedBase, Idx
+    >>> from sympy.tensor import IndexedBase
     >>> x, y, A = map(IndexedBase, ['x', 'y', 'A'])
     >>> i, j, a, z = symbols('i j a z', integer=True)
 
@@ -329,7 +329,7 @@ def get_contraction_structure(expr):
     ========
 
     >>> from sympy.tensor.index_methods import get_contraction_structure
-    >>> from sympy import symbols, default_sort_key
+    >>> from sympy import default_sort_key
     >>> from sympy.tensor import IndexedBase, Idx
     >>> x, y, A = map(IndexedBase, ['x', 'y', 'A'])
     >>> i, j, k, l = map(Idx, ['i', 'j', 'k', 'l'])
diff --git a/sympy/tensor/indexed.py b/sympy/tensor/indexed.py
index a4fde4cb72ea..b8b4da9461e7 100644
--- a/sympy/tensor/indexed.py
+++ b/sympy/tensor/indexed.py
@@ -496,7 +496,7 @@ def shape(self):
         Examples
         ========
 
-        >>> from sympy import IndexedBase, Idx, Symbol
+        >>> from sympy import IndexedBase, Idx
         >>> from sympy.abc import x, y
         >>> IndexedBase('A', shape=(x, y)).shape
         (x, y)
@@ -601,7 +601,7 @@ class Idx(Expr):
     Examples
     ========
 
-    >>> from sympy import IndexedBase, Idx, symbols, oo
+    >>> from sympy import Idx, symbols, oo
     >>> n, i, L, U = symbols('n i L U', integer=True)
 
     If a string is given for the label an integer ``Symbol`` is created and the
diff --git a/sympy/testing/pytest.py b/sympy/testing/pytest.py
index 9843d62ffcbb..174fc590a035 100644
--- a/sympy/testing/pytest.py
+++ b/sympy/testing/pytest.py
@@ -222,8 +222,6 @@ def warns_deprecated_sympy():
 
     >>> from sympy.testing.pytest import warns_deprecated_sympy
     >>> from sympy.utilities.exceptions import SymPyDeprecationWarning
-    >>> import warnings
-
     >>> with warns_deprecated_sympy():
     ...     SymPyDeprecationWarning("Don't use", feature="old thing",
     ...         deprecated_since_version="1.0", issue=123).warn()
diff --git a/sympy/unify/usympy.py b/sympy/unify/usympy.py
index 668bfabd8bec..71a19248a8ae 100644
--- a/sympy/unify/usympy.py
+++ b/sympy/unify/usympy.py
@@ -79,7 +79,7 @@ def unify(x, y, s=None, variables=(), **kwargs):
     ========
 
     >>> from sympy.unify.usympy import unify
-    >>> from sympy import Basic, cos
+    >>> from sympy import Basic
     >>> from sympy.abc import x, y, z, p, q
 
     >>> next(unify(Basic(1, 2), Basic(1, x), variables=[x]))
diff --git a/sympy/utilities/_compilation/compilation.py b/sympy/utilities/_compilation/compilation.py
index 7a72b786c5bb..df640524cd4f 100644
--- a/sympy/utilities/_compilation/compilation.py
+++ b/sympy/utilities/_compilation/compilation.py
@@ -505,13 +505,6 @@ def compile_link_import_py_ext(sources, extname=None, build_dir='.', compile_kwa
     =======
 
     The imported module from of the python extension.
-
-    Examples
-    ========
-
-    >>> mod = compile_link_import_py_ext(['fft.f90', 'conv.cpp', '_fft.pyx'])  # doctest: +SKIP
-    >>> Aprim = mod.fft(A)  # doctest: +SKIP
-
     """
     if extname is None:
         extname = os.path.splitext(os.path.basename(sources[-1]))[0]
diff --git a/sympy/utilities/decorator.py b/sympy/utilities/decorator.py
index d6d290452e1b..97152bd08748 100644
--- a/sympy/utilities/decorator.py
+++ b/sympy/utilities/decorator.py
@@ -192,7 +192,7 @@ def public(obj):
 
     >>> from sympy.utilities.decorator import public
 
-    >>> __all__
+    >>> __all__ # noqa: F821
     Traceback (most recent call last):
     ...
     NameError: name '__all__' is not defined
@@ -201,7 +201,7 @@ def public(obj):
     ... def some_function():
     ...     pass
 
-    >>> __all__
+    >>> __all__ # noqa: F821
     ['some_function']
 
     """
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index 2cc512612919..94c6b700b68f 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -2535,12 +2535,12 @@ def kbins(l, k, ordered=None):
     Examples
     ========
 
-    >>> from sympy.utilities.iterables import kbins
+    >>> from __future__ import print_function
 
     The default is to give the items in the same order, but grouped
     into k partitions without any reordering:
 
-    >>> from __future__ import print_function
+    >>> from sympy.utilities.iterables import kbins
     >>> for p in kbins(list(range(5)), 2):
     ...     print(p)
     ...
@@ -2565,9 +2565,9 @@ def kbins(l, k, ordered=None):
         10 means A == D
         11 means A == A
 
-    >>> for ordered in [None, 0, 1, 10, 11]:
-    ...     print('ordered = %s' % ordered)
-    ...     for p in kbins(list(range(3)), 2, ordered=ordered):
+    >>> for ordered_flag in [None, 0, 1, 10, 11]:
+    ...     print('ordered = %s' % ordered_flag)
+    ...     for p in kbins(list(range(3)), 2, ordered=ordered_flag):
     ...         print('     %s' % p)
     ...
     ordered = None
diff --git a/sympy/utilities/lambdify.py b/sympy/utilities/lambdify.py
index 583864f8d443..4d24aafb14a4 100644
--- a/sympy/utilities/lambdify.py
+++ b/sympy/utilities/lambdify.py
@@ -351,7 +351,6 @@ def lambdify(args, expr, modules=None, printer=None, use_imps=True,
     functions. This may be preferable to using ``evalf`` (which uses mpmath on
     the backend) in some cases.
 
-    >>> import mpmath
     >>> f = lambdify(x, sin(x), 'mpmath')
     >>> f(1)
     0.8414709848078965
@@ -1287,7 +1286,6 @@ def implemented_function(symfunc, implementation):
 
     >>> from sympy.abc import x
     >>> from sympy.utilities.lambdify import lambdify, implemented_function
-    >>> from sympy import Function
     >>> f = implemented_function('f', lambda x: x+1)
     >>> lam_f = lambdify(x, f(x))
     >>> lam_f(4)
diff --git a/sympy/vector/functions.py b/sympy/vector/functions.py
index e3a08c716b7b..8b75445f002d 100644
--- a/sympy/vector/functions.py
+++ b/sympy/vector/functions.py
@@ -354,7 +354,7 @@ def scalar_potential_difference(field, coord_sys, point1, point2):
     Examples
     ========
 
-    >>> from sympy.vector import CoordSys3D, Point
+    >>> from sympy.vector import CoordSys3D
     >>> from sympy.vector import scalar_potential_difference
     >>> R = CoordSys3D('R')
     >>> P = R.origin.locate_new('P', R.x*R.i + R.y*R.j + R.z*R.k)
@@ -483,7 +483,6 @@ def orthogonalize(*vlist, **kwargs):
     ========
 
     >>> from sympy.vector.coordsysrect import CoordSys3D
-    >>> from sympy.vector.vector import Vector, BaseVector
     >>> from sympy.vector.functions import orthogonalize
     >>> C = CoordSys3D('C')
     >>> i, j, k = C.base_vectors()
diff --git a/sympy/vector/point.py b/sympy/vector/point.py
index 2a29bdbdea20..bfcd13e88da8 100644
--- a/sympy/vector/point.py
+++ b/sympy/vector/point.py
@@ -57,7 +57,7 @@ def position_wrt(self, other):
         Examples
         ========
 
-        >>> from sympy.vector import Point, CoordSys3D
+        >>> from sympy.vector import CoordSys3D
         >>> N = CoordSys3D('N')
         >>> p1 = N.origin.locate_new('p1', 10 * N.i)
         >>> N.origin.position_wrt(p1)
@@ -109,7 +109,7 @@ def locate_new(self, name, position):
         Examples
         ========
 
-        >>> from sympy.vector import Point, CoordSys3D
+        >>> from sympy.vector import CoordSys3D
         >>> N = CoordSys3D('N')
         >>> p1 = N.origin.locate_new('p1', 10 * N.i)
         >>> p1.position_wrt(N.origin)
@@ -133,7 +133,7 @@ def express_coordinates(self, coordinate_system):
         Examples
         ========
 
-        >>> from sympy.vector import Point, CoordSys3D
+        >>> from sympy.vector import CoordSys3D
         >>> N = CoordSys3D('N')
         >>> p1 = N.origin.locate_new('p1', 10 * N.i)
         >>> p2 = p1.locate_new('p2', 5 * N.j)
diff --git a/sympy/vector/vector.py b/sympy/vector/vector.py
index 999937b6c3cb..84237b478e2b 100644
--- a/sympy/vector/vector.py
+++ b/sympy/vector/vector.py
@@ -223,7 +223,6 @@ def projection(self, other, scalar=False):
         ========
 
         >>> from sympy.vector.coordsysrect import CoordSys3D
-        >>> from sympy.vector.vector import Vector, BaseVector
         >>> C = CoordSys3D('C')
         >>> i, j, k = C.base_vectors()
         >>> v1 = i + j + k


Code review comments:
diff_hunk: @@ -7,6 +7,7 @@ set -x
 
 if [[ "${TEST_FLAKE8}" == "true" ]]; then
     flake8 sympy;
+    flake8 --doctests --per-file-ignores='sympy/interactive/session.py:F821' sympy;
path: bin/test_travis.sh
text: @user1:
Does this run the normal flake8 checking twice?

@author:
I'm not sure. It seems so, but I can't find anything in the `flake8` documentation that says whether or not the code is also checked with the `--doctests` flag.

edit: I found a reference here: https://bugs.launchpad.net/pyflakes/+bug/1223150 so it seems that the normal checking truly is done twice.