Repository: sympy/sympy

Issue:
solve_linear_system contains duplicate rref algorithm
```
It seems that we have two different rref (reduced row echelon form) algorithms, one in 
Matrix.rref and one in solve_linear_system in solvers.py.  I don't know too much about these 
algorithms, so I don't know if they really are the same algorithm, but in either case, I think 
solve should just be calling Matrix.rref.

This might speed things up too (examples taken from Lay's Linear Algebra, Third Edition):

In [15]: c = Matrix([[1, -2, 1, 0],[0, 2, -8, 8],[-4, 5, 9, -9]])

In [16]: c
Out[16]: 
⎡1   -2  1   0 ⎤
⎢              ⎥
⎢0   2   -8  8 ⎥
⎢              ⎥
⎣-4  5   9   -9⎦

In [17]: c.rref()
Out[17]: 
⎛⎡1  0  0  29⎤, [0, 1, 2]⎞
⎜⎢           ⎥           ⎟
⎜⎢0  1  0  16⎥           ⎟
⎜⎢           ⎥           ⎟
⎝⎣0  0  1  3 ⎦           ⎠

In [19]: solve_linear_system(c, x, y, z)
Out[19]: {x: 29, y: 16, z: 3}

In [25]: %timeit solve_linear_system(c, x, y, z)
100 loops, best of 3: 11.2 ms per loop

In [26]: %timeit c.rref()
1000 loops, best of 3: 3.59 ms per loop

In [30]: e = Matrix([[4, -1, -1, 0, 0, 0, 0, 0, 5],[-1, 4, 0, -1, 0, 0, 0, 0, 15],[-1, 0, 4, -1, -1, 0, 0, 
0, 0],[0, -1, -1, 4, 0, -1, 0, 0, 10],[0, 0, -1, 0, 4, -1, -1, 0, 0],[0, 0, 0, -1, -1, 4, 0, -1, 10],[0, 0, 
0, 0, -1, 0, 4, -1, 20],[0, 0, 0, 0, 0, -1, -1, 4, 30]])

In [32]: e.rref()
Out[32]: 
<output omitted>

In [33]: var('x1 x2 x3 x4 x5 x6 x6 x7 x8')
Out[33]: (x₁, x₂, x₃, x₄, x₅, x₆, x₆, x₇, x₈)

In [34]: solve_linear_system(e, x1, x2, x3, x4, x5, x6, x7, x8)
Out[34]: 
⎧    827      1377      886      1546      1171      1831      1967      2517⎫
⎨x₁: ───, x₂: ────, x₃: ───, x₄: ────, x₅: ────, x₆: ────, x₇: ────, x₈: 
────⎬
⎩    209      209       209      209       209       209       209       209 ⎭

In [35]: %timeit e.rref()
10 loops, best of 3: 65 ms per loop

In [36]: %timeit solve_linear_system(e, x1, x2, x3, x4, x5, x6, x7, x8)
10 loops, best of 3: 84.5 ms per loop
```

Original issue for #4949: http://code.google.com/p/sympy/issues/detail?id=1850
Original author: https://code.google.com/u/asmeurer@gmail.com/
Original owner: https://code.google.com/u/asmeurer@gmail.com/


Discussion:

Patch:
diff --git a/sympy/simplify/ratsimp.py b/sympy/simplify/ratsimp.py
index 63c52a1822..080e634a69 100644
--- a/sympy/simplify/ratsimp.py
+++ b/sympy/simplify/ratsimp.py
@@ -40,7 +40,7 @@ def ratsimpmodprime(expr, G, *gens, **args):
     >>> from sympy.abc import x, y
     >>> eq = (x + y**5 + y)/(x - y)
     >>> ratsimpmodprime(eq, [x*y**5 - x - y], x, y, order='lex')
-    (x**2 + x*y + x + y)/(x**2 - x*y)
+    (-x**2 - x*y - x - y)/(-x**2 + x*y)
 
     If ``polynomial`` is False, the algorithm computes a rational
     simplification which minimizes the sum of the total degrees of
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 54023e4682..cfb982e12f 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -2301,162 +2301,22 @@ def solve_linear_system(system, *symbols, **flags):
     {}
 
     """
-    do_simplify = flags.get('simplify', True)
+    from sympy.solvers.solveset import linsolve
+    from sympy.sets import FiniteSet
 
-    if system.rows == system.cols - 1 == len(symbols):
-        try:
-            # well behaved n-equations and n-unknowns
-            inv = inv_quick(system[:, :-1])
-            rv = dict(zip(symbols, inv*system[:, -1]))
-            if do_simplify:
-                for k, v in rv.items():
-                    rv[k] = simplify(v)
-            if not all(i.is_zero for i in rv.values()):
-                # non-trivial solution
-                return rv
-        except ValueError:
-            pass
-
-    matrix = system[:, :]
-    syms = list(symbols)
-
-    i, m = 0, matrix.cols - 1  # don't count augmentation
-
-    while i < matrix.rows:
-        if i == m:
-            # an overdetermined system
-            if any(matrix[i:, m]):
-                return None   # no solutions
-            else:
-                # remove trailing rows
-                matrix = matrix[:i, :]
-                break
-
-        if not matrix[i, i]:
-            # there is no pivot in current column
-            # so try to find one in other columns
-            for k in range(i + 1, m):
-                if matrix[i, k]:
-                    break
-            else:
-                if matrix[i, m]:
-                    # We need to know if this is always zero or not. We
-                    # assume that if there are free symbols that it is not
-                    # identically zero (or that there is more than one way
-                    # to make this zero). Otherwise, if there are none, this
-                    # is a constant and we assume that it does not simplify
-                    # to zero XXX are there better (fast) ways to test this?
-                    # The .equals(0) method could be used but that can be
-                    # slow; numerical testing is prone to errors of scaling.
-                    if not matrix[i, m].free_symbols:
-                        return None  # no solution
-
-                    # A row of zeros with a non-zero rhs can only be accepted
-                    # if there is another equivalent row. Any such rows will
-                    # be deleted.
-                    nrows = matrix.rows
-                    rowi = matrix.row(i)
-                    ip = None
-                    j = i + 1
-                    while j < matrix.rows:
-                        # do we need to see if the rhs of j
-                        # is a constant multiple of i's rhs?
-                        rowj = matrix.row(j)
-                        if rowj == rowi:
-                            matrix.row_del(j)
-                        elif rowj[:-1] == rowi[:-1]:
-                            if ip is None:
-                                _, ip = rowi[-1].as_content_primitive()
-                            _, jp = rowj[-1].as_content_primitive()
-                            if not (simplify(jp - ip) or simplify(jp + ip)):
-                                matrix.row_del(j)
-
-                        j += 1
-
-                    if nrows == matrix.rows:
-                        # no solution
-                        return None
-                # zero row or was a linear combination of
-                # other rows or was a row with a symbolic
-                # expression that matched other rows, e.g. [0, 0, x - y]
-                # so now we can safely skip it
-                matrix.row_del(i)
-                if not matrix:
-                    # every choice of variable values is a solution
-                    # so we return an empty dict instead of None
-                    return dict()
-                continue
-
-            # we want to change the order of columns so
-            # the order of variables must also change
-            syms[i], syms[k] = syms[k], syms[i]
-            matrix.col_swap(i, k)
-
-        pivot_inv = S.One/matrix[i, i]
-
-        # divide all elements in the current row by the pivot
-        matrix.row_op(i, lambda x, _: x * pivot_inv)
-
-        for k in range(i + 1, matrix.rows):
-            if matrix[k, i]:
-                coeff = matrix[k, i]
-
-                # subtract from the current row the row containing
-                # pivot and multiplied by extracted coefficient
-                matrix.row_op(k, lambda x, j: simplify(x - matrix[i, j]*coeff))
+    assert system.shape[1] == len(symbols) + 1
 
-        i += 1
+    # This is just a wrapper for linsolve:
+    sol = linsolve(system, *symbols)
 
-    # if there weren't any problems, augmented matrix is now
-    # in row-echelon form so we can check how many solutions
-    # there are and extract them using back substitution
-
-    if len(syms) == matrix.rows:
-        # this system is Cramer equivalent so there is
-        # exactly one solution to this system of equations
-        k, solutions = i - 1, {}
-
-        while k >= 0:
-            content = matrix[k, m]
-
-            # run back-substitution for variables
-            for j in range(k + 1, m):
-                content -= matrix[k, j]*solutions[syms[j]]
-
-            if do_simplify:
-                solutions[syms[k]] = simplify(content)
-            else:
-                solutions[syms[k]] = content
-
-            k -= 1
-
-        return solutions
-    elif len(syms) > matrix.rows:
-        # this system will have infinite number of solutions
-        # dependent on exactly len(syms) - i parameters
-        k, solutions = i - 1, {}
-
-        while k >= 0:
-            content = matrix[k, m]
-
-            # run back-substitution for variables
-            for j in range(k + 1, i):
-                content -= matrix[k, j]*solutions[syms[j]]
-
-            # run back-substitution for parameters
-            for j in range(i, m):
-                content -= matrix[k, j]*syms[j]
-
-            if do_simplify:
-                solutions[syms[k]] = simplify(content)
-            else:
-                solutions[syms[k]] = content
-
-            k -= 1
-
-        return solutions
+    if sol is S.EmptySet:
+        return None
+    elif isinstance(sol, FiniteSet):
+        assert len(sol) == 1
+        sol = sol.args[0]
+        return {sym:val for sym, val in zip(symbols, sol) if sym != val}
     else:
-        return []   # no solutions
+        raise RuntimeError("We should never get here!")
 
 
 def solve_undetermined_coeffs(equ, coeffs, sym, **flags):
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index 3f152dcf58..3b112d9748 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -2251,9 +2251,12 @@ def generate_derangements(perm):
     ========
     sympy.functions.combinatorial.factorials.subfactorial
     """
-    for p in multiset_permutations(perm):
-        if not any(i == j for i, j in zip(perm, p)):
-            yield p
+    p = multiset_permutations(perm)
+    indices = range(len(perm))
+    p0 = next(p)
+    for pi in p:
+        if all(pi[i] != p0[i] for i in indices):
+            yield pi
 
 
 def necklaces(n, k, free=False):


Test patch:
diff --git a/sympy/polys/tests/test_partfrac.py b/sympy/polys/tests/test_partfrac.py
index aa697ec448..0f7f033d99 100644
--- a/sympy/polys/tests/test_partfrac.py
+++ b/sympy/polys/tests/test_partfrac.py
@@ -85,6 +85,22 @@ def test_apart_extension():
 
     assert factor(together(apart(f)).expand()) == f
 
+    # https://github.com/sympy/sympy/issues/18531
+    from sympy.core import Mul
+    def mul2(expr):
+        # 2-arg mul hack...
+        return Mul(2, expr, evaluate=False)
+
+    f = ((x**2 + 1)**3/((x - 1)**2*(x + 1)**2*(-x**2 + 2*x + 1)*(x**2 + 2*x - 1)))
+    g = (1/mul2(x - sqrt(2) + 1)
+       - 1/mul2(x - sqrt(2) - 1)
+       + 1/mul2(x + 1 + sqrt(2))
+       - 1/mul2(x - 1 + sqrt(2))
+       + 1/mul2((x + 1)**2)
+       + 1/mul2((x - 1)**2))
+
+    assert apart(f, x, extension={sqrt(2)}) == g
+
 
 def test_apart_full():
     f = 1/(x**2 + 1)
diff --git a/sympy/simplify/tests/test_ratsimp.py b/sympy/simplify/tests/test_ratsimp.py
index 3dc6ee81b7..3f00f92313 100644
--- a/sympy/simplify/tests/test_ratsimp.py
+++ b/sympy/simplify/tests/test_ratsimp.py
@@ -44,7 +44,7 @@ def test_ratsimpmodprime():
     b = x - y
     F = [x*y**5 - x - y]
     assert ratsimpmodprime(a/b, F, x, y, order='lex') == \
-        (x**2 + x*y + x + y) / (x**2 - x*y)
+        (-x**2 - x*y - x - y) / (-x**2 + x*y)
 
     a = x + y**2 - 2
     b = x + y**2 - y - 1
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 37301d8109..a18637b851 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -326,6 +326,13 @@ def test_linear_system():
 
     assert solve([x + y + z + t, -z - t], x, y, z, t) == {x: -y, z: -t}
 
+    # https://github.com/sympy/sympy/issues/6420
+    M = Matrix([[0,    15.0, 10.0, 700.0],
+                [1,    1,    1,    100.0],
+                [0,    10.0, 5.0,  200.0],
+                [-5.0, 0,    0,    0    ]])
+
+    assert solve_linear_system(M, x, y, z) == {x: 0, y: -60.0, z: 160.0}
 
 def test_linear_system_function():
     a = Function('a')
@@ -1208,8 +1215,12 @@ def test_issue_5849():
         ans[0]) for ei in e] == [0, 0, I3 - I6, -I3 + I6, 0, 0, 0, 0, 0]
 
 
+# Should this work at all? Simpler examples fail e.g.:
+#    solve([x+y+z,x+y],[x,y])  ==  []
+# Here a solution only exists if I3 == I6 which is not generically true.
+@XFAIL
 def test_issue_5849_matrix():
-    '''Same as test_2750 but solved with the matrix solver.'''
+    '''Same as test_issue_5849 but solved with the matrix solver.'''
     I1, I2, I3, I4, I5, I6 = symbols('I1:7')
     dI1, dI4, dQ2, dQ4, Q2, Q4 = symbols('dI1,dI4,dQ2,dQ4,Q2,Q4')
 
@@ -1739,7 +1750,7 @@ def test_issue_5114_6611():
     ans = solve(list(eqs), list(v), simplify=False)
     # If time is taken to simplify then then 2617 below becomes
     # 1168 and the time is about 50 seconds instead of 2.
-    assert sum([s.count_ops() for s in ans.values()]) <= 2617
+    assert sum([s.count_ops() for s in ans.values()]) <= 3093
 
 
 def test_det_quick():
diff --git a/sympy/utilities/tests/test_iterables.py b/sympy/utilities/tests/test_iterables.py
index 729bcd8ccc..d7938ed032 100644
--- a/sympy/utilities/tests/test_iterables.py
+++ b/sympy/utilities/tests/test_iterables.py
@@ -543,7 +543,6 @@ def test_derangements():
         [2, 3, 0, 1], [2, 3, 1, 0], [3, 0, 1, 2], [3, 2, 0, 1], [3, 2, 1, 0]]
     assert list(generate_derangements([0, 1, 2, 2])) == [
         [2, 2, 0, 1], [2, 2, 1, 0]]
-    assert list(generate_derangements('ba')) == [list('ab')]
 
 
 def test_necklaces():