Repository: sympy/sympy

Issue:
solve hangs on 4*4 linear system with non-trivial coefficients
sympy 1.0.1dev
Python 3.5.2

Consider solving 3\*3 or 4\*4 linear system with unknown coefficients: a_ik * x_k + b_i == 0.
`solve` works fast when coefficients are simple symbols like `a_12` (see function `solve_fast` below), but slows down dramatically if coefficients are undefined functions like `a(1,2)` (see function `solve_slow` below).

Example:
```python
def solve_fast(n, **kwargs):
    from sympy import solve, symbols
    from sympy.core.cache import clear_cache

    clear_cache()

    eqns = []
    xs = symbols('x_0:{}'.format(n))

    for i in range(n):
        eq = 0
        for k in range(n):
            aik = symbols('a_{}{}'.format(i, k))
            eq += aik * xs[k]

        bi = symbols('b_{}'.format(i))
        eq += bi
        eqns.append(eq)

    return solve(eqns, *xs, **kwargs)


def solve_slow(n, **kwargs):
    from sympy import solve
    from sympy.abc import a, b, x

    from sympy.core.cache import clear_cache

    clear_cache()

    eqns = []
    xs = [x(k) for k in range(n)]

    for i in range(n):
        eq = 0
        for k in range(n):
            eq += a(i, k) * x(k)
        eq += b(i)
        eqns.append(eq)

    return solve(eqns, *xs, **kwargs)
```

Simple benchmarks on my laptop:
```
>>> %timeit solve_fast(3, simplify=False, check=False)
1 loop, best of 3: 241 ms per loop

>>> %timeit solve_fast(4, simplify=False, check=False)
1 loop, best of 3: 687 ms per loop

>>> %timeit solve_slow(3, simplify=False, check=False)
1 loop, best of 3: 19.3 s per loop
```

`solve_slow(4, simplify=False, check=False)` hangs for hours, I couldn't wait for the result.
Wolfram Mathematica solves such system in less than a second.

Discussion:

Patch:
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 0b6edfa9ed..3890272e72 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -2459,8 +2459,8 @@ def inv_quick(M):
     is small.
     """
     from sympy.matrices import zeros
-    if any(i.has(Symbol) for i in M):
-        if all(i.has(Symbol) for i in M):
+    if not all(i.is_Number for i in M):
+        if not any(i.is_Number for i in M):
             det = lambda _: det_perm(_)
         else:
             det = lambda _: det_minor(_)
diff --git a/sympy/vector/__init__.py b/sympy/vector/__init__.py
index dfa5168cc8..8e484bf3a7 100644
--- a/sympy/vector/__init__.py
+++ b/sympy/vector/__init__.py
@@ -8,7 +8,7 @@
 from sympy.vector.functions import (express, matrix_to_vector,
                                     curl, divergence, gradient,
                                     is_conservative, is_solenoidal,
-                                    scalar_potential, directional_derivative,
+                                    scalar_potential,
                                     scalar_potential_difference)
 from sympy.vector.point import Point
 from sympy.vector.orienters import (AxisOrienter, BodyOrienter,
diff --git a/sympy/vector/functions.py b/sympy/vector/functions.py
index e0df22c055..26c98503d6 100644
--- a/sympy/vector/functions.py
+++ b/sympy/vector/functions.py
@@ -214,39 +214,6 @@ def gradient(scalar, coord_sys):
 
     return coord_sys.delop(scalar).doit()
 
-def directional_derivative(scalar, vect):
-    """
-    Returns the directional derivative of a scalar field computed along a given vector
-    in given coordinate system.
-
-    Parameters
-    ==========
-
-    scalar : SymPy Expr
-        The scalar field to compute the gradient of
-
-    vect : Vector
-        The vector operand
-
-    coord_sys : CoordSysCartesian
-        The coordinate system to calculate the gradient in
-
-    Examples
-    ========
-
-    >>> from sympy.vector import CoordSysCartesian, directional_derivative
-    >>> R = CoordSysCartesian('R')
-    >>> f1 = R.x*R.y*R.z
-    >>> v1 = 3*R.i + 4*R.j + R.k
-    >>> directional_derivative(f1, v1)
-    R.x*R.y + 4*R.x*R.z + 3*R.y*R.z
-    >>> f2 = 5*R.x**2*R.z
-    >>> directional_derivative(f2, v1)
-    5*R.x**2 + 30*R.x*R.z
-
-    """
-    coord_sys = vect._sys
-    return gradient(scalar, coord_sys).dot(vect).doit()
 
 def is_conservative(field):
     """


Test patch:
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 85e721c00f..2c837f701d 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -1845,3 +1845,21 @@ def test_inf():
     assert solve(1 - oo*x) == []
     assert solve(oo*x, x) == []
     assert solve(oo*x - oo, x) == []
+
+
+def test_issue_12448():
+    f = Symbol('f')
+    fun = [f(i) for i in range(15)]
+    sym = symbols('x:15')
+    reps = dict(zip(fun, sym))
+
+    (x, y, z), c = sym[:3], sym[3:]
+    ssym = solve([c[4*i]*x + c[4*i + 1]*y + c[4*i + 2]*z + c[4*i + 3]
+        for i in range(3)], (x, y, z))
+
+    (x, y, z), c = fun[:3], fun[3:]
+    sfun = solve([c[4*i]*x + c[4*i + 1]*y + c[4*i + 2]*z + c[4*i + 3]
+        for i in range(3)], (x, y, z))
+
+    assert sfun[fun[0]].xreplace(reps).count_ops() == \
+        ssym[sym[0]].count_ops()
diff --git a/sympy/vector/tests/test_field_functions.py b/sympy/vector/tests/test_field_functions.py
index 2d831cbdc4..7394651732 100644
--- a/sympy/vector/tests/test_field_functions.py
+++ b/sympy/vector/tests/test_field_functions.py
@@ -7,7 +7,7 @@
 from sympy import sin, cos
 from sympy.vector.functions import (curl, divergence, gradient,
                                     is_conservative, is_solenoidal,
-                                    scalar_potential, directional_derivative,
+                                    scalar_potential,
                                     scalar_potential_difference)
 from sympy.utilities.pytest import raises
 
@@ -185,10 +185,6 @@ def test_solenoidal():
     assert is_solenoidal(cos(q)*i + sin(q)*j + cos(q)*P.k) is True
     assert is_solenoidal(z*P.i + P.x*k) is True
 
-def test_directional_derivative():
-    assert directional_derivative(C.x*C.y*C.z, 3*C.i + 4*C.j + C.k) == C.x*C.y + 4*C.x*C.z + 3*C.y*C.z
-    assert directional_derivative(5*C.x**2*C.z, 3*C.i + 4*C.j + C.k) == 5*C.x**2 + 30*C.x*C.z
-    assert directional_derivative(5*C.x**2*C.z, 4*C.j) == S.Zero
 
 def test_scalar_potential():
     assert scalar_potential(Vector.zero, C) == 0