Repository: sympy/sympy

Issue:
collect giving wrong results for derivatives
```py
>>> collect(x*diff(f(x), x, x) - 6*x + f(x).diff(x)*diff(f(x), x, x) + f(x).diff(x), f(x).diff(x, x))
x*Derivative(f(x), x, x) - 6*x + (Derivative(f(x), x, x) + 1)*Derivative(f(x), x)
```

I asked it to collect the second derivative, but it collected the first. I believe this is because it is trying to be smart with derivatives, but it relies on the args order

```py
>>> (diff(f(x), x, x)*diff(f(x), x)).args
(Derivative(f(x), x), Derivative(f(x), x, x))
```

This is a modified example from https://github.com/sympy/sympy/pull/13059#issuecomment-323195614 (in that PR, the args order of `f(x)*diff(f(x), x)` is reversed from master. 

I'm unclear what exactly collect is trying to do with derivatives, but it should respect the symbols given to it and their order (the original example had `collect(x*Derivative(f(x), x) - 6*x + f(x)*Derivative(f(x), x) + f(x), [f(x).diff(x), f(x)])`). Beyond that, I'd say should prefer to collect higher order derivatives. 

Discussion:

Patch:
diff --git a/sympy/combinatorics/coset_table.py b/sympy/combinatorics/coset_table.py
index 7704d76f06..ea1e9486e5 100644
--- a/sympy/combinatorics/coset_table.py
+++ b/sympy/combinatorics/coset_table.py
@@ -117,7 +117,7 @@ def is_complete(self):
         `\alpha \in \Omega` and `x \in A`.
 
         """
-        return not any(None in self.table[coset] for coset in self.omega)
+        return not any([None in self.table[coset] for coset in self.omega])
 
     # Pg. 153 [1]
     def define(self, alpha, x):
diff --git a/sympy/combinatorics/fp_groups.py b/sympy/combinatorics/fp_groups.py
index c0fb07e08b..43b1395ee9 100644
--- a/sympy/combinatorics/fp_groups.py
+++ b/sympy/combinatorics/fp_groups.py
@@ -167,8 +167,7 @@ def standardize_coset_table(self):
         self._coset_table.standardize()
         self._is_standardized = True
 
-    def coset_table(self, H, strategy="relator_based", max_cosets=None,
-                                                 draft=None, incomplete=False):
+    def coset_table(self, H, strategy="relator_based"):
         """
         Return the mathematical coset table of ``self`` in ``H``.
 
@@ -178,14 +177,12 @@ def coset_table(self, H, strategy="relator_based", max_cosets=None,
                 if not self._is_standardized:
                     self.standardize_coset_table()
             else:
-                C = self.coset_enumeration([], strategy, max_cosets=max_cosets,
-                                            draft=draft, incomplete=incomplete)
+                C = self.coset_enumeration([], strategy)
                 self._coset_table = C
                 self.standardize_coset_table()
             return self._coset_table.table
         else:
-            C = self.coset_enumeration(H, strategy, max_cosets=max_cosets,
-                                            draft=draft, incomplete=incomplete)
+            C = self.coset_enumeration(H, strategy)
             C.standardize()
             return C.table
 
@@ -273,25 +270,22 @@ def _finite_index_subgroup(self, s=[]):
         mid = (len(s)+1)//2
         half1 = s[:mid]
         half2 = s[mid:]
-        draft1 = None
-        draft2 = None
         m = 200
         C = None
         while not C and (m/2 < CosetTable.coset_table_max_limit):
             m = min(m, CosetTable.coset_table_max_limit)
-            draft1 = self.coset_enumeration(half1, max_cosets=m,
-                                 draft=draft1, incomplete=True)
-            if draft1.is_complete():
-                C = draft1
+            try:
+                C = self.coset_enumeration(half1, max_cosets=m)
                 half = half1
-            else:
-                draft2 = self.coset_enumeration(half2, max_cosets=m,
-                                 draft=draft2, incomplete=True)
-                if draft2.is_complete():
-                    C = draft2
-                    half = half2
+            except ValueError:
+                pass
             if not C:
-                m *= 2
+                try:
+                    C = self.coset_enumeration(half2, max_cosets=m)
+                    half = half2
+                except ValueError:
+                    m *= 2
+                    continue
         if not C:
             return None, None
         C.compress()
diff --git a/sympy/combinatorics/homomorphisms.py b/sympy/combinatorics/homomorphisms.py
index c94ddce4dc..aa38e6224a 100644
--- a/sympy/combinatorics/homomorphisms.py
+++ b/sympy/combinatorics/homomorphisms.py
@@ -207,39 +207,6 @@ def is_trivial(self):
         '''
         return self.image().order() == 1
 
-    def restrict_to(self, H):
-        '''
-        Return the restriction of the homomorphism to the subgroup `H`
-        of the domain.
-
-        '''
-        if not isinstance(H, PermutationGroup) or not H.is_subgroup(self.domain):
-            raise ValueError("Given H is not a subgroup of the domain")
-        domain = H
-        images = {g: self(g) for g in H.generators}
-        return GroupHomomorphism(domain, self.codomain, images)
-
-    def invert_subgroup(self, H):
-        '''
-        Return the subgroup of the domain that is the inverse image
-        of the subgroup `H` of the homomorphism image
-
-        '''
-        if not H.is_subgroup(self.image()):
-            raise ValueError("Given H is not a subgroup of the image")
-        gens = []
-        P = PermutationGroup(self.image().identity)
-        for h in H.generators:
-            h_i = self.invert(h)
-            if h_i not in P:
-                gens.append(h_i)
-                P = PermutationGroup(gens)
-            for k in self.kernel().generators:
-                if k*h_i not in P:
-                    gens.append(k*h_i)
-                    P = PermutationGroup(gens)
-        return P
-
 def homomorphism(domain, codomain, gens, images=[], check=True):
     '''
     Create (if possible) a group homomorphism from the group `domain`
@@ -336,59 +303,3 @@ def _image(r):
         if not s:
             return False
     return True
-
-def orbit_homomorphism(group, omega):
-    '''
-    Return the homomorphism induced by the action of the permutation
-    group `group` on the set `omega` that is closed under the action.
-
-    '''
-    from sympy.combinatorics import Permutation
-    from sympy.combinatorics.named_groups import SymmetricGroup
-    codomain = SymmetricGroup(len(omega))
-    identity = codomain.identity
-    omega = list(omega)
-    images = {g: identity*Permutation([omega.index(o^g) for o in omega]) for g in group.generators}
-    group._schreier_sims(base=omega)
-    H = GroupHomomorphism(group, codomain, images)
-    if len(group.basic_stabilizers) > len(omega):
-        H._kernel = group.basic_stabilizers[len(omega)]
-    else:
-        H._kernel = PermutationGroup([group.identity])
-    return H
-
-def block_homomorphism(group, blocks):
-    '''
-    Return the homomorphism induced by the action of the permutation
-    group `group` on the block system `blocks`. The latter should be
-    of the same form as returned by the `minimal_block` method for
-    permutation groups, namely a list of length `group.degree` where
-    the i-th entry is a representative of the block i belongs to.
-
-    '''
-    from sympy.combinatorics import Permutation
-    from sympy.combinatorics.named_groups import SymmetricGroup
-
-    n = len(blocks)
-
-    # number the blocks; m is the total number,
-    # b is such that b[i] is the number of the block i belongs to,
-    # p is the list of length m such that p[i] is the representative
-    # of the i-th block
-    m = 0
-    p = []
-    b = [None]*n
-    for i in range(n):
-        if blocks[i] == i:
-            p.append(i)
-            b[i] = m
-            m += 1
-    for i in range(n):
-        b[i] = b[blocks[i]]
-
-    codomain = SymmetricGroup(m)
-    # the list corresponding to the identity permutation in codomain
-    identity = range(m)
-    images = {g: Permutation([b[p[i]^g] for i in identity]) for g in group.generators}
-    H = GroupHomomorphism(group, codomain, images)
-    return H
diff --git a/sympy/combinatorics/perm_groups.py b/sympy/combinatorics/perm_groups.py
index af0e7765e0..a0f29963bb 100644
--- a/sympy/combinatorics/perm_groups.py
+++ b/sympy/combinatorics/perm_groups.py
@@ -1867,95 +1867,12 @@ def is_primitive(self, randomized=True):
         orbits = stab.orbits()
         for orb in orbits:
             x = orb.pop()
-            if x != 0 and any(e != 0 for e in self.minimal_block([0, x])):
+            if x != 0 and self.minimal_block([0, x]) != [0]*n:
                 self._is_primitive = False
                 return False
         self._is_primitive = True
         return True
 
-    def minimal_blocks(self, randomized=True):
-        '''
-        For a transitive group, return the list of all minimal
-        block systems. If a group is intrasitive, return `False`.
-
-        Examples
-        ========
-        >>> from sympy.combinatorics import Permutation
-        >>> from sympy.combinatorics.perm_groups import PermutationGroup
-        >>> from sympy.combinatorics.named_groups import DihedralGroup
-        >>> DihedralGroup(6).minimal_blocks()
-        [[0, 3, 0, 3, 0, 3], [0, 4, 2, 0, 4, 2]]
-        >>> G = PermutationGroup(Permutation(1,2,5))
-        >>> G.minimal_blocks()
-        False
-
-        See Also
-        ========
-
-        minimal_block, is_transitive, is_primitive
-
-        '''
-        def _number_blocks(blocks):
-            # number the blocks of a block system
-            # in order and return the number of
-            # blocks and the tuple with the
-            # reordering
-            n = len(blocks)
-            appeared = {}
-            m = 0
-            b = [None]*n
-            for i in range(n):
-                if blocks[i] not in appeared:
-                    appeared[blocks[i]] = m
-                    b[i] = m
-                    m += 1
-                else:
-                    b[i] = appeared[blocks[i]]
-            return tuple(b), m
-
-        if not self.is_transitive():
-            return False
-        blocks = []
-        num_blocks = []
-        rep_blocks = []
-        if randomized:
-            random_stab_gens = []
-            v = self.schreier_vector(0)
-            for i in range(len(self)):
-                random_stab_gens.append(self.random_stab(0, v))
-            stab = PermutationGroup(random_stab_gens)
-        else:
-            stab = self.stabilizer(0)
-        orbits = stab.orbits()
-        for orb in orbits:
-            x = orb.pop()
-            if x != 0:
-                block = self.minimal_block([0, x])
-                num_block, m = _number_blocks(block)
-                # a representative block (containing 0)
-                rep = set(j for j in range(self.degree) if num_block[j] == 0)
-                # check if the system is minimal with
-                # respect to the already discovere ones
-                minimal = True
-                to_remove = []
-                for i, r in enumerate(rep_blocks):
-                    if len(r) > len(rep) and rep.issubset(r):
-                        # i-th block sytem is not minimal
-                        del num_blocks[i], blocks[i]
-                        to_remove.append(rep_blocks[i])
-                    elif len(r) < len(rep) and r.issubset(rep):
-                        # the system being checked is not minimal
-                        minimal = False
-                        break
-                # remove non-minimal representative blocks
-                rep_blocks = [r for r in rep_blocks if r not in to_remove]
-
-                if minimal and num_block not in num_blocks:
-                    blocks.append(block)
-                    num_blocks.append(num_block)
-                    rep_blocks.append(rep)
-        return blocks
-
     @property
     def is_solvable(self):
         """Test if the group is solvable.
@@ -2759,11 +2676,9 @@ def schreier_sims(self):
         return
 
     def _schreier_sims(self, base=None):
-        schreier = self.schreier_sims_incremental(base=base, slp_dict=True)
-        base, strong_gens = schreier[:2]
+        base, strong_gens = self.schreier_sims_incremental(base=base)
         self._base = base
         self._strong_gens = strong_gens
-        self._strong_gens_slp = schreier[2]
         if not base:
             self._transversals = []
             self._basic_orbits = []
@@ -2783,7 +2698,7 @@ def _schreier_sims(self, base=None):
         self._basic_orbits = [sorted(x) for x in basic_orbits]
         self._transversal_slp = slps
 
-    def schreier_sims_incremental(self, base=None, gens=None, slp_dict=False):
+    def schreier_sims_incremental(self, base=None, gens=None):
         """Extend a sequence of points and generating set to a base and strong
         generating set.
 
@@ -2798,12 +2713,6 @@ def schreier_sims_incremental(self, base=None, gens=None, slp_dict=False):
             relative to the base obtained. Optional parameter with default
             value ``self.generators``.
 
-        slp_dict
-            If `True`, return a dictionary `{g: gens}` for each strong
-            generator `g` where `gens` is a list of strong generators
-            coming before `g` in `strong_gens`, such that the product
-            of the elements of `gens` is equal to `g`.
-
         Returns
         =======
 
@@ -2852,8 +2761,7 @@ def schreier_sims_incremental(self, base=None, gens=None, slp_dict=False):
         id_af = list(range(degree))
         # handle the trivial group
         if len(gens) == 1 and gens[0].is_Identity:
-            if slp_dict:
-                return base, gens, {gens[0]: [gens[0]]}
+            self._strong_gens_slp = {gens[0]: [gens[0]]}
             return base, gens
         # prevent side effects
         _base, _gens = base[:], gens[:]
@@ -2946,27 +2854,22 @@ def schreier_sims_incremental(self, base=None, gens=None, slp_dict=False):
                 continue
             i -= 1
 
+        # create the list of the strong generators strong_gens and rewrite
+        # the indices of strong_gens_slp in terms of the elements of strong_gens
         strong_gens = _gens[:]
-
-        if slp_dict:
-            # create the list of the strong generators strong_gens and
-            # rewrite the indices of strong_gens_slp in terms of the
-            # elements of strong_gens
-            for k, slp in strong_gens_slp:
-                strong_gens.append(k)
-                for i in range(len(slp)):
-                    s = slp[i]
-                    if isinstance(s[1], tuple):
-                        slp[i] = strong_gens_distr[s[0]][s[1][0]]**-1
-                    else:
-                        slp[i] = strong_gens_distr[s[0]][s[1]]
-            strong_gens_slp = dict(strong_gens_slp)
-            # add the original generators
-            for g in _gens:
-                strong_gens_slp[g] = [g]
-            return (_base, strong_gens, strong_gens_slp)
-
-        strong_gens.extend([k for k, _ in strong_gens_slp])
+        for k, slp in strong_gens_slp:
+            strong_gens.append(k)
+            for i in range(len(slp)):
+                s = slp[i]
+                if isinstance(s[1], tuple):
+                    slp[i] = strong_gens_distr[s[0]][s[1][0]]**-1
+                else:
+                    slp[i] = strong_gens_distr[s[0]][s[1]]
+        strong_gens_slp = dict(strong_gens_slp)
+        # add the original generators
+        for g in _gens:
+            strong_gens_slp[g] = [g]
+        self._strong_gens_slp = strong_gens_slp
         return _base, strong_gens
 
     def schreier_sims_random(self, base=None, gens=None, consec_succ=10,
@@ -3457,7 +3360,7 @@ def update_nu(l):
                 else:
                     nu[l] = sorted_orbits[l][temp_index]
             # line 29: set the next element from the current branch and update
-            # accordingly
+            # accorndingly
             c[l] += 1
             if l == 0:
                 gamma  = sorted_orbits[l][c[l]]
@@ -3518,177 +3421,6 @@ def transitivity_degree(self):
         else:
             return self._transitivity_degree
 
-    def _p_elements_group(G, p):
-        '''
-        For an abelian p-group G return the subgroup consisting of
-        all elements of order p (and the identity)
-
-        '''
-        gens = G.generators[:]
-        gens = sorted(gens, key=lambda x: x.order(), reverse=True)
-        gens_p = [g**(g.order()/p) for g in gens]
-        gens_r = []
-        for i in range(len(gens)):
-            x = gens[i]
-            x_order = x.order()
-            # x_p has order p
-            x_p = x**(x_order/p)
-            if i > 0:
-                P = PermutationGroup(gens_p[:i])
-            else:
-                P = PermutationGroup(G.identity)
-            if x**(x_order/p) not in P:
-                gens_r.append(x**(x_order/p))
-            else:
-                # replace x by an element of order (x.order()/p)
-                # so that gens still generates G
-                g = P.generator_product(x_p, original=True)
-                for s in g:
-                    x = x*s**-1
-                x_order = x_order/p
-                # insert x to gens so that the sorting is preserved
-                del gens[i]
-                del gens_p[i]
-                j = i - 1
-                while j < len(gens) and gens[j].order() >= x_order:
-                    j += 1
-                gens = gens[:j] + [x] + gens[j:]
-                gens_p = gens_p[:j] + [x] + gens_p[j:]
-        return PermutationGroup(gens_r)
-
-    def sylow_subgroup(self, p):
-        '''
-        Return a p-Sylow subgroup of the group.
-
-        The algorithm is described in [1], Chapter 4, Section 7
-
-        Examples
-        ========
-        >>> from sympy.combinatorics.named_groups import DihedralGroup
-        >>> from sympy.combinatorics.named_groups import SymmetricGroup
-        >>> from sympy.combinatorics.named_groups import AlternatingGroup
-
-        >>> D = DihedralGroup(6)
-        >>> S = D.sylow_subgroup(2)
-        >>> S.order()
-        4
-        >>> G = SymmetricGroup(6)
-        >>> S = G.sylow_subgroup(5)
-        >>> S.order()
-        5
-
-        >>> G1 = AlternatingGroup(3)
-        >>> G2 = AlternatingGroup(5)
-        >>> G3 = AlternatingGroup(9)
-
-        >>> S1 = G1.sylow_subgroup(3)
-        >>> S2 = G2.sylow_subgroup(3)
-        >>> S3 = G3.sylow_subgroup(3)
-
-        >>> len1 = len(S1.lower_central_series())
-        >>> len2 = len(S2.lower_central_series())
-        >>> len3 = len(S3.lower_central_series())
-
-        >>> len1 == len2
-        True
-        >>> len1 < len3
-        True
-
-        '''
-        from sympy.combinatorics.homomorphisms import (homomorphism,
-            orbit_homomorphism, block_homomorphism)
-        from sympy.ntheory.primetest import isprime
-
-        if not isprime(p):
-            raise ValueError("p must be a prime")
-
-        def is_p_group(G):
-            # check if the order of G is a power of p
-            # and return the power
-            m = G.order()
-            n = 0
-            while m % p == 0:
-                m = m/p
-                n += 1
-                if m == 1:
-                    return True, n
-            return False, n
-
-        def _sylow_reduce(mu, nu):
-            # reduction based on two homomorphisms
-            # mu and nu with trivially intersecting
-            # kernels
-            Q = mu.image().sylow_subgroup(p)
-            Q = mu.invert_subgroup(Q)
-            nu = nu.restrict_to(Q)
-            R = nu.image().sylow_subgroup(p)
-            return nu.invert_subgroup(R)
-
-        order = self.order()
-        if order % p != 0:
-            return PermutationGroup([self.identity])
-        p_group, n = is_p_group(self)
-        if p_group:
-            return self
-
-        # if there is a non-trivial orbit with size not divisible
-        # by p, the sylow subgroup is contained in its stabilizer
-        # (by orbit-stabilizer theorem)
-        orbits = self.orbits()
-        non_p_orbits = [o for o in orbits if len(o) % p != 0 and len(o) != 1]
-        if non_p_orbits:
-            G = self.stabilizer(list(non_p_orbits[0]).pop())
-            return G.sylow_subgroup(p)
-
-        if not self.is_transitive():
-            # apply _sylow_reduce to orbit actions
-            orbits = sorted(orbits, key = lambda x: len(x))
-            omega1 = orbits.pop()
-            omega2 = orbits[0].union(*orbits)
-            mu = orbit_homomorphism(self, omega1)
-            nu = orbit_homomorphism(self, omega2)
-            return _sylow_reduce(mu, nu)
-
-        blocks = self.minimal_blocks()
-        if len(blocks) > 1:
-            # apply _sylow_reduce to block system actions
-            mu = block_homomorphism(self, blocks[0])
-            nu = block_homomorphism(self, blocks[1])
-            return _sylow_reduce(mu, nu)
-        elif len(blocks) == 1:
-            block = list(blocks)[0]
-            if any(e != 0 for e in block):
-                # self is imprimitive
-                mu = block_homomorphism(self, block)
-                if not is_p_group(mu.image())[0]:
-                    S = mu.image().sylow_subgroup(p)
-                    return mu.invert_subgroup(S).sylow_subgroup(p)
-
-        # find an element of order p
-        g = self.random()
-        g_order = g.order()
-        while g_order % p != 0 or g_order == 0:
-            g = self.random()
-            g_order = g.order()
-        g = g**(g_order // p)
-        if order % p**2 != 0:
-            return PermutationGroup(g)
-
-        C = self.centralizer(g)
-        while C.order() % p**n != 0:
-            S = C.sylow_subgroup(p)
-            s_order = S.order()
-            Z = S.center()
-            P = Z._p_elements_group(p)
-            h = P.random()
-            C_h = self.centralizer(h)
-            while C_h.order() % p*s_order != 0:
-                h = P.random()
-                C_h = self.centralizer(h)
-            C = C_h
-
-        return C.sylow_subgroup(p)
-
     def presentation(G):
         '''
         Return an `FpGroup` presentation of the group.
@@ -3803,7 +3535,6 @@ def _rewrite(group, w):
 
         return simplify_presentation(G_p)
 
-
 def _orbit(degree, generators, alpha, action='tuples'):
     r"""Compute the orbit of alpha `\{g(\alpha) | g \in G\}` as a set.
 
diff --git a/sympy/geometry/entity.py b/sympy/geometry/entity.py
index ef9f3d0fcd..6fece70a4d 100644
--- a/sympy/geometry/entity.py
+++ b/sympy/geometry/entity.py
@@ -100,7 +100,6 @@ def __contains__(self, other):
         raise NotImplementedError()
 
     def __getnewargs__(self):
-        """Returns a tuple that will be passed to __new__ on unpickling."""
         return tuple(self.args)
 
     def __ne__(self, o):
@@ -120,11 +119,9 @@ def is_seq_and_not_point(a):
         return Basic.__new__(cls, *args)
 
     def __radd__(self, a):
-        """Implementation of reverse add method."""
         return a.__add__(self)
 
     def __rdiv__(self, a):
-        """Implementation of reverse division method."""
         return a.__div__(self)
 
     def __repr__(self):
@@ -133,11 +130,9 @@ def __repr__(self):
         return type(self).__name__ + repr(self.args)
 
     def __rmul__(self, a):
-        """Implementation of reverse multiplication method."""
         return a.__mul__(self)
 
     def __rsub__(self, a):
-        """Implementation of reverse substraction method."""
         return a.__sub__(self)
 
     def __str__(self):
@@ -349,32 +344,6 @@ def is_similar(self, other):
         raise NotImplementedError()
 
     def reflect(self, line):
-        """
-        Reflects an object across a line.
-
-        Parameters
-        ==========
-
-        line: Line
-
-        Examples
-        ========
-
-        >>> from sympy import pi, sqrt, Line, RegularPolygon
-        >>> l = Line((0, pi), slope=sqrt(2))
-        >>> pent = RegularPolygon((1, 2), 1, 5)
-        >>> rpent = pent.reflect(l)
-        >>> rpent
-        RegularPolygon(Point2D(-2*sqrt(2)*pi/3 - 1/3 + 4*sqrt(2)/3, 2/3 + 2*sqrt(2)/3 + 2*pi/3), -1, 5, -pi/5 + acos(1/3))
-
-        >>> from sympy import pi, Line, Circle, Point
-        >>> l = Line((0, pi), slope=1)
-        >>> circ = Circle(Point(0, 0), 5)
-        >>> rcirc = circ.reflect(l)
-        >>> rcirc
-        Circle(Point2D(-pi, pi), -5)
-
-        """
         from sympy import atan, Point, Dummy, oo
 
         g = self
diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py
index 7b5d943048..93535adb9b 100644
--- a/sympy/logic/boolalg.py
+++ b/sympy/logic/boolalg.py
@@ -1516,7 +1516,7 @@ def _check_pair(minterm1, minterm2):
 
 def _convert_to_varsSOP(minterm, variables):
     """
-    Converts a term in the expansion of a function from binary to its
+    Converts a term in the expansion of a function from binary to it's
     variable form (for SOP).
     """
     temp = []
@@ -1532,7 +1532,7 @@ def _convert_to_varsSOP(minterm, variables):
 
 def _convert_to_varsPOS(maxterm, variables):
     """
-    Converts a term in the expansion of a function from binary to its
+    Converts a term in the expansion of a function from binary to it's
     variable form (for POS).
     """
     temp = []
diff --git a/sympy/physics/units/quantities.py b/sympy/physics/units/quantities.py
index c452df0840..7ab29c2687 100644
--- a/sympy/physics/units/quantities.py
+++ b/sympy/physics/units/quantities.py
@@ -114,12 +114,10 @@ def get_dimensional_expr(expr):
             return dim
         elif isinstance(expr, Function):
             args = [Quantity.get_dimensional_expr(arg) for arg in expr.args]
-            if all(i == 1 for i in args):
-                return S.One
             return expr.func(*args)
         elif isinstance(expr, Quantity):
             return expr.dimension.name
-        return S.One
+        return 1
 
     @staticmethod
     def _collect_factor_and_dimension(expr):
diff --git a/sympy/physics/units/systems/si.py b/sympy/physics/units/systems/SI.py
similarity index 95%
rename from sympy/physics/units/systems/si.py
rename to sympy/physics/units/systems/SI.py
index 5b2f564f4f..b5de2c345e 100644
--- a/sympy/physics/units/systems/si.py
+++ b/sympy/physics/units/systems/SI.py
@@ -7,7 +7,7 @@
 
 Example:
 
-    >>> from sympy.physics.units.systems.si import SI
+    >>> from sympy.physics.units.systems.SI import SI
     >>> from sympy.physics.units import avogadro, boltzmann, lux
     >>> SI.print_unit_base(avogadro)
     6.022140857e+23/mole
diff --git a/sympy/physics/units/systems/__init__.py b/sympy/physics/units/systems/__init__.py
index 7a25846f52..eb421e400f 100644
--- a/sympy/physics/units/systems/__init__.py
+++ b/sympy/physics/units/systems/__init__.py
@@ -3,4 +3,4 @@
 from sympy.physics.units.systems.mks import MKS, _mks_dim
 from sympy.physics.units.systems.mksa import MKSA, _mksa_dim
 from sympy.physics.units.systems.natural import _natural_dim, natural
-from sympy.physics.units.systems.si import SI, _si_dim
+from sympy.physics.units.systems.SI import SI, _si_dim
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index ff6fedc17b..f88f76083c 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -4383,7 +4383,7 @@ def _update_args(args, key, value):
 
 
 @public
-def degree(f, gen=0):
+def degree(f, *gens, **args):
     """
     Return the degree of ``f`` in the given variable.
 
@@ -4402,32 +4402,15 @@ def degree(f, gen=0):
     >>> degree(0, x)
     -oo
 
-    See also
-    ========
-    total_degree
-    degree_list
     """
+    options.allowed_flags(args, ['gen', 'polys'])
 
-    f = sympify(f, strict=True)
-    if f.is_Poly:
-        p = f
-        isNum = p.as_expr().is_Number
-    else:
-        isNum = f.is_Number
-        if not isNum:
-            p, _ = poly_from_expr(f)
-
-    if isNum:
-        return S.Zero if f else S.NegativeInfinity
-
-    if not sympify(gen, strict=True).is_Number:
-        if f.is_Poly and gen not in p.gens:
-            # try recast without explicit gens
-            p, _ = poly_from_expr(f.as_expr())
-        if gen not in p.gens:
-            return S.Zero
+    try:
+        F, opt = poly_from_expr(f, *gens, **args)
+    except PolificationFailed as exc:
+        raise ComputationFailed('degree', 1, exc)
 
-    return Integer(p.degree(gen))
+    return sympify(F.degree(opt.gen))
 
 
 @public
diff --git a/sympy/printing/str.py b/sympy/printing/str.py
index bef6d4271f..0f0edcaeae 100644
--- a/sympy/printing/str.py
+++ b/sympy/printing/str.py
@@ -20,7 +20,6 @@ class StrPrinter(Printer):
     _default_settings = {
         "order": None,
         "full_prec": "auto",
-        "sympy_integers": False,
     }
 
     _relationals = dict()
@@ -550,8 +549,6 @@ def _print_Rational(self, expr):
         if expr.q == 1:
             return str(expr.p)
         else:
-            if self._settings.get("sympy_integers", False):
-                return "S(%s)/%s" % (expr.p, expr.q)
             return "%s/%s" % (expr.p, expr.q)
 
     def _print_PythonRational(self, expr):
diff --git a/sympy/simplify/radsimp.py b/sympy/simplify/radsimp.py
index 91c29dd4cb..a65693d270 100644
--- a/sympy/simplify/radsimp.py
+++ b/sympy/simplify/radsimp.py
@@ -157,6 +157,9 @@ def collect(expr, syms, func=None, evaluate=None, exact=False, distribute_order_
     ========
     collect_const, collect_sqrt, rcollect
     """
+    expr = sympify(expr)
+    syms = list(syms) if iterable(syms) else [syms]
+
     if evaluate is None:
         evaluate = global_evaluate[0]
 
@@ -331,12 +334,8 @@ def parse_expression(terms, pattern):
                 expr.base, syms, func, True, exact, distribute_order_term)
             return Pow(b, expr.exp)
 
-    if iterable(syms):
-        syms = [expand_power_base(i, deep=False) for i in syms]
-    else:
-        syms = [expand_power_base(syms, deep=False)]
+    syms = [expand_power_base(i, deep=False) for i in syms]
 
-    expr = sympify(expr)
     order_term = None
 
     if distribute_order_term:
@@ -352,7 +351,10 @@ def parse_expression(terms, pattern):
 
     collected, disliked = defaultdict(list), S.Zero
     for product in summa:
-        terms = [parse_term(i) for i in Mul.make_args(product)]
+        c, nc = product.args_cnc(split_1=False)
+        args = list(ordered(c)) + nc
+        terms = [parse_term(i) for i in args]
+        small_first = True
 
         for symbol in syms:
             if SYMPY_DEBUG:
@@ -360,6 +362,9 @@ def parse_expression(terms, pattern):
                     str(terms), str(symbol))
                 )
 
+            if isinstance(symbol, Derivative) and small_first:
+                terms = list(reversed(terms))
+                small_first = not small_first
             result = parse_expression(terms, symbol)
 
             if SYMPY_DEBUG:
@@ -371,12 +376,14 @@ def parse_expression(terms, pattern):
                 # when there was derivative in current pattern we
                 # will need to rebuild its expression from scratch
                 if not has_deriv:
-                    index = 1
+                    margs = []
                     for elem in elems:
-                        e = elem[1]
-                        if elem[2] is not None:
-                            e *= elem[2]
-                        index *= Pow(elem[0], e)
+                        if elem[2] is None:
+                            e = elem[1]
+                        else:
+                            e = elem[1]*elem[2]
+                        margs.append(Pow(elem[0], e))
+                    index = Mul(*margs)
                 else:
                     index = make_expression(elems)
                 terms = expand_power_base(make_expression(terms), deep=False)
diff --git a/sympy/vector/functions.py b/sympy/vector/functions.py
index aa2cf7064c..ddc60ab7e9 100644
--- a/sympy/vector/functions.py
+++ b/sympy/vector/functions.py
@@ -157,7 +157,7 @@ def directional_derivative(field, direction_vector):
     coord_sys = _get_coord_sys_from_expr(field)
     if len(coord_sys) > 0:
         # TODO: This gets a random coordinate system in case of multiple ones:
-        coord_sys = next(iter(coord_sys))
+        coord_sys = coord_sys.pop()
         field = express(field, coord_sys, variables=True)
         i, j, k = coord_sys.base_vectors()
         x, y, z = coord_sys.base_scalars()
diff --git a/sympy/vector/operators.py b/sympy/vector/operators.py
index d7d4357c0a..abd6eaeddc 100644
--- a/sympy/vector/operators.py
+++ b/sympy/vector/operators.py
@@ -1,24 +1,14 @@
 import collections
 from sympy.core.expr import Expr
-from sympy.core import sympify, S, preorder_traversal
+from sympy.core import sympify, S
 from sympy.vector.coordsysrect import CoordSys3D
-from sympy.vector.vector import Vector, VectorMul, VectorAdd, Cross, Dot, dot
+from sympy.vector.vector import Vector, VectorMul, VectorAdd
 from sympy.vector.scalar import BaseScalar
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.core.function import Derivative
 from sympy import Add, Mul
 
 
-def _get_coord_systems(expr):
-    g = preorder_traversal(expr)
-    ret = set([])
-    for i in g:
-        if isinstance(i, CoordSys3D):
-            ret.add(i)
-            g.skip()
-    return frozenset(ret)
-
-
 def _get_coord_sys_from_expr(expr, coord_sys=None):
     """
     expr : expression
@@ -32,13 +22,16 @@ def _get_coord_sys_from_expr(expr, coord_sys=None):
             issue=12884,
         ).warn()
 
-    return _get_coord_systems(expr)
+    try:
+        return expr.atoms(CoordSys3D)
+    except:
+        return set([])
 
 
 def _split_mul_args_wrt_coordsys(expr):
     d = collections.defaultdict(lambda: S.One)
     for i in expr.args:
-        d[_get_coord_systems(i)] *= i
+        d[frozenset(i.atoms(CoordSys3D))] *= i
     return list(d.values())
 
 
@@ -156,7 +149,7 @@ def curl(vect, coord_sys=None, doit=True):
     if len(coord_sys) == 0:
         return Vector.zero
     elif len(coord_sys) == 1:
-        coord_sys = next(iter(coord_sys))
+        coord_sys = coord_sys.pop()
         i, j, k = coord_sys.base_vectors()
         x, y, z = coord_sys.base_scalars()
         h1, h2, h3 = coord_sys.lame_coefficients()
@@ -176,7 +169,7 @@ def curl(vect, coord_sys=None, doit=True):
         return outvec
     else:
         # TODO: use some of the vector calculus properties for this:
-        coord_sys = next(iter(coord_sys))  # get one random coord_sys
+        coord_sys = coord_sys.pop()  # get one random coord_sys
         i, j, k = coord_sys.base_vectors()
         x, y, z = coord_sys.base_scalars()
         h1, h2, h3 = coord_sys.lame_coefficients()
@@ -235,15 +228,13 @@ def divergence(vect, coord_sys=None, doit=True):
     2*R.z
 
     """
-    # TODO: Remove this line when warning from issue #12884 will be removed
     coord_sys = _get_coord_sys_from_expr(vect, coord_sys)
+
     if len(coord_sys) == 0:
         return S.Zero
-    elif len(coord_sys) == 1:
-        if isinstance(vect, (Cross, Curl, Gradient)):
-            return Divergence(vect)
+    else:
         # TODO: is case of many coord systems, this gets a random one:
-        coord_sys = next(iter(coord_sys))
+        coord_sys = coord_sys.pop()
         i, j, k = coord_sys.base_vectors()
         x, y, z = coord_sys.base_scalars()
         h1, h2, h3 = coord_sys.lame_coefficients()
@@ -253,27 +244,9 @@ def divergence(vect, coord_sys=None, doit=True):
              / (h1 * h2 * h3)
         vz = _diff_conditional(vect.dot(k), z, h1, h2) \
              / (h1 * h2 * h3)
-        res = vx + vy + vz
         if doit:
-            return res.doit()
-        return res
-    else:
-        if isinstance(vect, (Add, VectorAdd)):
-            return Add.fromiter(divergence(i, doit=doit) for i in vect.args)
-        elif isinstance(vect, (Mul, VectorMul)):
-            vector = [i for i in vect.args if isinstance(i, (Vector, Cross, Gradient))][0]
-            scalar = Mul.fromiter(i for i in vect.args if not isinstance(i, (Vector, Cross, Gradient)))
-            res = Dot(vector, gradient(scalar)) + scalar*divergence(vector, doit=doit)
-            if doit:
-                try:
-                    return res.doit()
-                except:
-                    return res
-            return res
-        elif isinstance(vect, (Cross, Curl, Gradient)):
-            return Divergence(vect)
-        else:
-            raise Divergence(vect)
+            return (vx + vy + vz).doit()
+        return vx + vy + vz
 
 
 def gradient(scalar_field, coord_sys=None, doit=True):
@@ -314,7 +287,7 @@ def gradient(scalar_field, coord_sys=None, doit=True):
     if len(coord_sys) == 0:
         return Vector.zero
     elif len(coord_sys) == 1:
-        coord_sys = next(iter(coord_sys))
+        coord_sys = coord_sys.pop()
         h1, h2, h3 = coord_sys.lame_coefficients()
         i, j, k = coord_sys.base_vectors()
         x, y, z = coord_sys.base_scalars()
diff --git a/sympy/vector/vector.py b/sympy/vector/vector.py
index ad96cae192..09ded3ca44 100755
--- a/sympy/vector/vector.py
+++ b/sympy/vector/vector.py
@@ -1,5 +1,5 @@
 from sympy.core.assumptions import StdFactKB
-from sympy.core import S, Pow, sympify
+from sympy.core import S, Pow, Symbol, sympify
 from sympy.core.expr import AtomicExpr, Expr
 from sympy.core.compatibility import range
 from sympy.matrices import det
@@ -488,6 +488,8 @@ class Dot(Expr):
     """
 
     def __new__(cls, expr1, expr2):
+        if not (isinstance(expr1, Vector) and isinstance(expr2, Vector)):
+            raise ValueError('Arguments must be the vectors.')
         expr1 = sympify(expr1)
         expr2 = sympify(expr2)
         obj = Expr.__new__(cls, expr1, expr2)


Test patch:
diff --git a/sympy/combinatorics/tests/test_perm_groups.py b/sympy/combinatorics/tests/test_perm_groups.py
index 87da3838a9..56d4e461c9 100644
--- a/sympy/combinatorics/tests/test_perm_groups.py
+++ b/sympy/combinatorics/tests/test_perm_groups.py
@@ -460,15 +460,6 @@ def test_minimal_block():
     assert P1.minimal_block([0, 2]) == [0, 3, 0, 3, 0, 3]
     assert P2.minimal_block([0, 2]) == [0, 3, 0, 3, 0, 3]
 
-def test_minimal_blocks():
-    P = PermutationGroup(Permutation(1, 5)(2, 4), Permutation(0, 1, 2, 3, 4, 5))
-    assert P.minimal_blocks() == [[0, 3, 0, 3, 0, 3], [0, 4, 5, 0, 4, 5]]
-
-    P = SymmetricGroup(5)
-    assert P.minimal_blocks() == [[0]*5]
-
-    P = PermutationGroup(Permutation(0, 3))
-    assert P.minimal_blocks() == False
 
 def test_max_div():
     S = SymmetricGroup(10)
@@ -787,47 +778,6 @@ def test_generator_product():
         w = g*w
     assert w == p
 
-def test_sylow_subgroup():
-    P = PermutationGroup(Permutation(1, 5)(2, 4), Permutation(0, 1, 2, 3, 4, 5))
-    S = P.sylow_subgroup(2)
-    assert S.order() == 4
-
-    P = DihedralGroup(12)
-    S = P.sylow_subgroup(3)
-    assert S.order() == 3
-
-    P = PermutationGroup(Permutation(1, 5)(2, 4), Permutation(0, 1, 2, 3, 4, 5), Permutation(0, 2))
-    S = P.sylow_subgroup(3)
-    assert S.order() == 9
-    S = P.sylow_subgroup(2)
-    assert S.order() == 8
-
-    P = SymmetricGroup(10)
-    S = P.sylow_subgroup(2)
-    assert S.order() == 256
-    S = P.sylow_subgroup(3)
-    assert S.order() == 81
-    S = P.sylow_subgroup(5)
-    assert S.order() == 25
-
-    # the length of the lower central series
-    # of a p-Sylow subgroup of Sym(n) grows with
-    # the highest exponent exp of p such
-    # that n >= p**exp
-    exp = 1
-    length = 0
-    for i in range(2, 9):
-        P = SymmetricGroup(i)
-        S = P.sylow_subgroup(2)
-        ls = S.lower_central_series()
-        if i // 2**exp > 0:
-            # length increases with exponent
-            assert len(ls) > length
-            length = len(ls)
-            exp += 1
-        else:
-            assert len(ls) == length
-
 def test_presentation():
     def _test(P):
         G = P.presentation()
diff --git a/sympy/physics/units/tests/test_quantities.py b/sympy/physics/units/tests/test_quantities.py
index 35fc0818da..82bb8c190e 100644
--- a/sympy/physics/units/tests/test_quantities.py
+++ b/sympy/physics/units/tests/test_quantities.py
@@ -341,12 +341,6 @@ def test_get_dimensional_expr_with_function():
         sin(Quantity.get_dimensional_expr(v_w1))
 
 
-def test_get_dimensional_expr_with_function_1():
-    v_w1 = Quantity('v_w1', length / time, meter / second)
-    v_w2 = Quantity('v_w2', length / time, meter / second)
-    assert Quantity.get_dimensional_expr(sin(v_w1/v_w2)) == 1
-
-
 def test_binary_information():
     assert convert_to(kibibyte, byte) == 1024*byte
     assert convert_to(mebibyte, byte) == 1024**2*byte
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 5e2d651ffa..32d1780880 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -1158,26 +1158,17 @@ def test_Poly_degree():
     assert Poly(2*y, x, y).degree(gen=y) == 1
     assert Poly(x*y, x, y).degree(gen=y) == 1
 
-    assert degree(0, x) == -oo
     assert degree(1, x) == 0
     assert degree(x, x) == 1
 
-    assert degree(x*y**2, x) == 1
-    assert degree(x*y**2, y) == 2
-    assert degree(x*y**2, z) == 0
+    assert degree(x*y**2, gen=x) == 1
+    assert degree(x*y**2, gen=y) == 2
 
-    assert degree(y**2+x**3) == 3
-    assert degree(y**2+x**3, 1) == 2
-    assert degree(pi) == 1
+    assert degree(x*y**2, x, y) == 1
+    assert degree(x*y**2, y, x) == 2
 
-    raises(PolynomialError, lambda: degree(x, 1.1))
+    raises(ComputationFailed, lambda: degree(1))
 
-    assert degree(Poly(0,x),z) == -oo
-    assert degree(Poly(1,x),z) == 0
-    assert degree(Poly(x**2+y**3,y)) == 3
-    assert degree(Poly(y**2 + x**3, y, x), 1) == 3
-    assert degree(Poly(y**2 + x**3, x), z) == 0
-    assert degree(Poly(y**2 + x**3 + z**4, x), z) == 4
 
 def test_Poly_degree_list():
     assert Poly(0, x).degree_list() == (-oo,)
diff --git a/sympy/printing/tests/test_str.py b/sympy/printing/tests/test_str.py
index c894039104..d297790204 100644
--- a/sympy/printing/tests/test_str.py
+++ b/sympy/printing/tests/test_str.py
@@ -489,9 +489,6 @@ def test_Rational():
     assert str(sqrt(-4)) == str(2*I)
     assert str(2**Rational(1, 10**10)) == "2**(1/10000000000)"
 
-    assert sstr(Rational(2, 3), sympy_integers=True) == "S(2)/3"
-    assert sstr(Symbol("x")**Rational(2, 3), sympy_integers=True) == "x**(S(2)/3)"
-
 
 def test_Float():
     # NOTE dps is the whole number of decimal digits
diff --git a/sympy/simplify/tests/test_radsimp.py b/sympy/simplify/tests/test_radsimp.py
index 601ea29649..4b0563d8d8 100644
--- a/sympy/simplify/tests/test_radsimp.py
+++ b/sympy/simplify/tests/test_radsimp.py
@@ -158,6 +158,7 @@ def test_radsimp_issue_3214():
 def test_collect_1():
     """Collect with respect to a Symbol"""
     x, y, z, n = symbols('x,y,z,n')
+    assert collect(1, x) == 1
     assert collect( x + y*x, x ) == x * (1 + y)
     assert collect( x + x**2, x ) == x + x**2
     assert collect( x**2 + y*x**2, x ) == (x**2)*(1 + y)
@@ -254,6 +255,10 @@ def test_collect_func():
         x**3: 1
     }
 
+    assert collect(f, x, factor, evaluate=False) == {
+        S.One: (a + 1)**3, x: 3*(a + 1)**2,
+        x**2: umul(S(3), a + 1), x**3: 1}
+
 
 def test_collect_order():
     a, b, x, t = symbols('a,b,x,t')
@@ -283,13 +288,6 @@ def test_rcollect():
     assert rcollect(sqrt(-((x + 1)*(y + 1))), z) == sqrt(-((x + 1)*(y + 1)))
 
 
-@XFAIL
-def test_collect_func_xfail():
-    # XXX: this test will pass when automatic constant distribution is removed (issue 4596)
-    assert collect(f, x, factor, evaluate=False) == {S.One: (a + 1)**3,
-                   x: 3*(a + 1)**2, x**2: 3*(a + 1), x**3: 1}
-
-
 @XFAIL
 def test_collect_issues():
     D = Derivative
@@ -353,6 +351,19 @@ def test_collect_const():
         2*sqrt(sqrt(2) + 3)*(sqrt(5)*x + y) + 2
 
 
+def test_issue_13143():
+    fx = f(x).diff(x)
+    e = f(x) + fx + f(x)*fx
+    # collect function before derivative
+    assert collect(e, Wild('w')) == f(x)*(fx + 1) + fx
+    e = f(x) + f(x)*fx + x*fx*f(x)
+    assert collect(e, fx) == (x*f(x) + f(x))*fx + f(x)
+    assert collect(e, f(x)) == (x*fx + fx + 1)*f(x)
+    e = f(x) + fx + f(x)*fx
+    assert collect(e, [f(x), fx]) == f(x)*(1 + fx) + fx
+    assert collect(e, [fx, f(x)]) == fx*(1 + f(x)) + f(x)
+
+
 def test_issue_6097():
     assert collect(a*y**(2.0*x) + b*y**(2.0*x), y**x) == y**(2.0*x)*(a + b)
     assert collect(a*2**(2.0*x) + b*2**(2.0*x), 2**x) == 2**(2.0*x)*(a + b)
diff --git a/sympy/vector/tests/test_field_functions.py b/sympy/vector/tests/test_field_functions.py
index 46d25d80c0..6adabc7b2c 100644
--- a/sympy/vector/tests/test_field_functions.py
+++ b/sympy/vector/tests/test_field_functions.py
@@ -5,8 +5,7 @@
 from sympy.core.symbol import symbols
 from sympy.core import S
 from sympy import sin, cos
-from sympy.vector.vector import Dot
-from sympy.vector.operators import curl, divergence, gradient, Gradient, Divergence, Cross
+from sympy.vector.operators import curl, divergence, gradient, Gradient
 from sympy.vector.deloperator import Del
 from sympy.vector.functions import (is_conservative, is_solenoidal,
                                     scalar_potential, directional_derivative,
@@ -277,16 +276,3 @@ def test_mixed_coordinates():
            (3*a.y*a.z*b.x*cos(q) - b.z*sin(a.x + a.y))*a.i + \
            (3*a.x*a.z*b.x*cos(q) - b.z*sin(a.x + a.y))*a.j + (3*a.x*a.y*b.x*cos(q))*a.k + \
            (3*a.x*a.y*a.z*cos(q))*b.i + b.z*b.j + (b.y + cos(a.x + a.y))*b.k
-    # divergence
-    assert divergence(a.i*a.x+a.j*a.y+a.z*a.k + b.i*b.x+b.j*b.y+b.z*b.k + c.i*c.x+c.j*c.y+c.z*c.k) == S(9)
-    # assert divergence(3*a.i*a.x*cos(a.x+b.z) + a.j*b.x*c.z) == None
-    assert divergence(3*a.i*a.x*a.z + b.j*b.x*c.z + 3*a.j*a.z*a.y) == \
-            6*a.z + Dot(b.j, c.z*b.i + b.x*c.k)
-    assert divergence(3*cos(q)*a.x*b.x*b.i*c.x) == \
-            Dot(b.i, (3*b.x*c.x*cos(q))*a.i + (3*a.x*c.x*cos(q))*b.i + (3*a.x*b.x*cos(q))*c.i)
-    assert divergence(a.x*b.x*c.x*Cross(a.x*a.i, a.y*b.j)) ==\
-           a.x*b.x*c.x*Divergence(Cross(a.x*a.i, a.y*b.j)) + \
-           Dot(Cross(a.x*a.i, a.y*b.j), b.x*c.x*a.i + a.x*c.x*b.i + a.x*b.x*c.i)
-    assert divergence(a.x*b.x*c.x*(a.x*a.i + b.x*b.i)) ==\
-           Dot(a.i, 2*a.x*b.x*c.x*a.i + a.x**2*c.x*b.i + a.x**2*b.x*c.i) + \
-           Dot(b.i, b.x**2*c.x*a.i + 2*a.x*b.x*c.x*b.i + a.x*b.x**2*c.i)