Repository: sympy/sympy
Pull request: add and use `has_xfree`
<!-- 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
fixes https://github.com/sympy/sympy/issues/23775 by removing auto-expansion by linsolve (while still handling Eq)
alternative to #23748

#### Brief description of what is fixed or changed

`has_free` is slower than needed for some cases since it does sub-expression matches:
```python
>>> eq = f(x+y+1)
>>> eq.has_free(f)
True
>>> eq.has_free(x + 1)
True
```

`has_xfree` is more strict and will only return True for full expression matches (similar to the way that `xreplace` makes replacements), so both of the above return False with `has_xfree`.

#### Other comments

`linsolve` and related functions that extract coefficients from linear systems fail when any term or factor depends on a provided symbol. Rather than providing a limited workaround for cases like treating `f(x)` as independent of `f(x).diff()`, those needing to solve such systems can replace all symbol-dependent objects with dummy symbols. The `recast_to_symbols` can do this:

```python
>>> from sympy.solvers.solveset import recast_to_symbols
>>> eqs = Tuple(f(x)+f(x).diff())
>>> e, k, d = recast_to_symbols(eqs, eqs.atoms(Derivative))
>>> linsolve(e, [f(x)]).xreplace(d)
{(-Derivative(f(x), x),)}
```
Alternatively (though hackish) one could temporarily define Derivative's `has_xfree` to return None
```python
>>> Derivative.has_xfree = lambda a, *b: None
>>> linsolve(eqs, [f(x)])
{(-Derivative(f(x), x),)}
>>> del Derivative.has_xfree
```
I tried to define `Derivative.bound_symbols` as `set(Derivative.variables)` so `iterfreeargs` would not return the function itself. This worked wrt `linsolve`, but many other tests failed. I am not sure of the implications, but it *seems* like the variables of differentation should be considered bound.
```diff
diff --git a/sympy/core/function.py b/sympy/core/function.py
index 551cf41fe8..0f0c6295f4 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -1480,6 +1480,10 @@ def __new__(cls, expr, *variables, **kwargs):
             expr = factor_terms(signsimp(expr))
         return expr

+    @property
+    def bound_symbols(cls):
+        return set(cls.variables)
+
     @property
     def canonical(cls):
         return cls.func(cls.expr,
```
then
```python
>>> from sympy.core.traversal import *
>>> list(iterfreeargs(f(x).diff(x,2)))
[Derivative(f(x), (x, 2)), 2]  --> in master: [Derivative(f(x), (x, 2)), f(x), (x, 2), x, x, 2]
>>> list(iterfreeargs(Integral(x,(x,1,2))))
[Integral(x, (x, 1, 2)), 1, 2]
>>> linsolve([f(x)+f(x).diff()], [f(x)])
{(-Derivative(f(x), x),)}
>>> integrate(Derivative(f(x), x), x)
x*Derivative(f(x), x)  --> in master: f(x)
>>> integrate(x**3*Derivative(f(x), (x, 4)))
x**4*Derivative(f(x), (x, 4))/4 --> in master: x**3*Derivative(f(x), (x, 3)) - 3*x**2*Derivative(f(x), (x, 2)) + 6*x*Derivative(f(x), x) - 6*f(x)
```

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

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 -->
* solvers
  * `linsolve` no longer fully expands expressions so any nonlinear terms that would cancel after expansion now raise
  * `linear_coeffs` is more efficient and can return results as a dictionary (which is more economical for sparse equations)
  * `linear_eq_to_matrix` can now handle system expressed in terms of functions
* core
  * `Basic.has_xfree` has been added. Like `xreplace`, it targets only full expression matching
<!-- END RELEASE NOTES -->


Problem statement:
linsolve can return NaN
The routine to identify coefficients of a linear system is fast and can quickly distinguish linear from non-linear system, but the obtained coefficients are not expanded and should probably be so before using them to compute a solution for the system: the coefficient of `x` in the following is 0
```python
>>> linsolve([x*(y+1)*(y-1)-x*y**2+x-2],[x])
{(NaN,)}
>>> linsolve([0],[x])
{(x,)}
```
This expansion/check should not be done in the coefficient extraction routine, however.

Hints from issue discussion:
The coefficient should be expanded automatically when constructing the domain. It fails because the wrong domain is used because of a bug in `construct_domain`:
```python
In [13]: construct_domain([-y**2 + (y-1)*(y+1) + y])
Out[13]: (ZZ[y], [y - 1])

In [14]: construct_domain([-y**2 + (y-1)*(y+1) + 1])
Out[14]: (EX, [EX(-y**2 + (y - 1)*(y + 1) + 1)])
```
This seems to fix it:
```diff
diff --git a/sympy/polys/constructor.py b/sympy/polys/constructor.py
index dad69c5..f6b0990 100644
--- a/sympy/polys/constructor.py
+++ b/sympy/polys/constructor.py
@@ -139,8 +139,6 @@ def _construct_composite(coeffs, opt):
         denoms.append(denom)
 
     polys, gens = parallel_dict_from_basic(numers + denoms)  # XXX: sorting
-    if not gens:
-        return None
 
     if opt.composite is None:
         if any(gen.is_number and gen.is_algebraic for gen in gens):
diff --git a/sympy/polys/matrices/linsolve.py b/sympy/polys/matrices/linsolve.py
index 75ae26d..3cbe078 100644
--- a/sympy/polys/matrices/linsolve.py
+++ b/sympy/polys/matrices/linsolve.py
@@ -124,9 +124,15 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     sym2index = dict(zip(syms, range(nsyms)))
     eqsdict = []
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
-        eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
+        eqdict = {}
+        for s, c in eq.items():
+            ec = elem_map[c]
+            if ec:
+                eqdict[sym2index[s]] = ec
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            erhs = elem_map[rhs]
+            if erhs:
+                eqdict[nsyms] = - erhs
         if eqdict:
             eqsdict.append(eqdict)
     sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)

```

Commit message: don't integrate 0 wrt t

It doesn't hurt, but it changes the form of the solution when a Sparse or Dense matrix is used because the Sparse automatically doesn't provide zeros (so they aren't integrated) while Dense does.

Patch to review:
diff --git a/sympy/core/basic.py b/sympy/core/basic.py
index 14640c32d9d3..1fcf7fa9092f 100644
--- a/sympy/core/basic.py
+++ b/sympy/core/basic.py
@@ -1258,6 +1258,33 @@ def has(self, *patterns):
         """
         return self._has(iterargs, *patterns)
 
+    def has_xfree(self, s):
+        """return True if self has any of the patterns in s as a
+        free argument, else False. This is like `Basic.has_free`
+        but this will only report exact argument matches.
+
+        Examples
+        ========
+
+        >>> from sympy import Function
+        >>> from sympy.abc import x, y
+        >>> f = Function('f')
+        >>> f(x).has_xfree({f})
+        False
+        >>> f(x).has_xfree({f(x)})
+        True
+        >>> f(x + 1).has_xfree({x})
+        True
+        >>> f(x + 1).has_xfree({x + 1})
+        True
+        >>> f(x + y + 1).has_xfree({x + 1})
+        False
+        """
+        # protect O(1) containment check by requiring:
+        if not type(s) in (dict, set):
+            raise ValueError('expecting set or dict argument')
+        return any(a in s for a in iterfreeargs(self))
+
     @cacheit
     def has_free(self, *patterns):
         """return True if self has object(s) ``x`` as a free expression
@@ -1284,8 +1311,26 @@ def has_free(self, *patterns):
         True
         >>> (x + y + 1).has_free(y + 1)
         True
-
         """
+        if not patterns:
+            return False
+        p0 = patterns[0]
+        if len(patterns) == 1 and iterable(p0) and not isinstance(p0, Basic):
+            # Basic can contain iterables (though not non-Basic, ideally)
+            # but don't encourage mixed passing patterns
+            raise ValueError(filldedent('''
+                Expecting 1 or more Basic args, not a single
+                non-Basic iterable. Don't forget to unpack
+                iterables: `eq.has_free(*patterns)`'''))
+        # try quick test first
+        try:
+            s = set(patterns)
+            rv = self.has_xfree(s)
+            if rv:
+                return rv
+        except TypeError:
+            pass
+        # now try matching through slower _has
         return self._has(iterfreeargs, *patterns)
 
     def _has(self, iterargs, *patterns):
diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py
index d94201d65dc2..51376c91c7a3 100644
--- a/sympy/core/tests/test_expr.py
+++ b/sympy/core/tests/test_expr.py
@@ -31,6 +31,7 @@
 from sympy.polys.polytools import factor, cancel, Poly
 from sympy.polys.rationaltools import together
 from sympy.series.order import O
+from sympy.sets.sets import FiniteSet
 from sympy.simplify.combsimp import combsimp
 from sympy.simplify.gammasimp import gammasimp
 from sympy.simplify.powsimp import powsimp
@@ -1624,12 +1625,31 @@ def test_has_free():
     assert Integral(f(x), (f(x), 1, y)).has_free(y)
     assert not Integral(f(x), (f(x), 1, y)).has_free(x)
     assert not Integral(f(x), (f(x), 1, y)).has_free(f(x))
+    # simple extraction
+    assert (x + 1 + y).has_free(x + 1)
+    assert not (x + 2 + y).has_free(x + 1)
+    assert (2 + 3*x*y).has_free(3*x)
+    raises(ValueError, lambda: x.has_free({x, y}))
+    s = FiniteSet(1, 2)
+    assert Piecewise((s, x > 3), (4, True)).has_free(s)
+    assert not Piecewise((1, x > 3), (4, True)).has_free(s)
+    # can't make set of these, but fallback will handle
+    assert not x.has_free(y, [])
+
+
+def test_has_xfree():
+    assert (x + 1).has_xfree({x})
+    assert ((x + 1)**2).has_xfree({x + 1})
+    assert not (x + y + 1).has_xfree({x + 1})
+    raises(ValueError, lambda: x.has_xfree(x))
+    raises(ValueError, lambda: x.has_xfree([x]))
 
 
 def test_issue_5300():
     x = Symbol('x', commutative=False)
     assert x*sqrt(2)/sqrt(6) == x*sqrt(3)/3
 
+
 def test_floordiv():
     from sympy.functions.elementary.integers import floor
     assert x // y == floor(x / y)
diff --git a/sympy/core/traversal.py b/sympy/core/traversal.py
index 980ff4a03381..0615c221a538 100644
--- a/sympy/core/traversal.py
+++ b/sympy/core/traversal.py
@@ -41,6 +41,7 @@ def iterfreeargs(expr, _first=True):
 
     Examples
     ========
+
     >>> from sympy import Integral, Function
     >>> from sympy.abc import x
     >>> f = Function('f')
@@ -66,8 +67,6 @@ def iterfreeargs(expr, _first=True):
             pass  # for cases like f being an arg
 
 
-
-
 class preorder_traversal:
     """
     Do a pre-order traversal of a tree.
diff --git a/sympy/polys/matrices/linsolve.py b/sympy/polys/matrices/linsolve.py
index 75ae26d669b2..ebf37e461e78 100644
--- a/sympy/polys/matrices/linsolve.py
+++ b/sympy/polys/matrices/linsolve.py
@@ -42,8 +42,11 @@
     sdm_nullspace_from_rref
 )
 
+from sympy.utilities.misc import filldedent
+
 
 def _linsolve(eqs, syms):
+
     """Solve a linear system of equations.
 
     Examples
@@ -69,8 +72,8 @@ def _linsolve(eqs, syms):
     nsyms = len(syms)
 
     # Convert to sparse augmented matrix (len(eqs) x (nsyms+1))
-    eqsdict, rhs = _linear_eq_to_dict(eqs, syms)
-    Aaug = sympy_dict_to_dm(eqsdict, rhs, syms)
+    eqsdict, const = _linear_eq_to_dict(eqs, syms)
+    Aaug = sympy_dict_to_dm(eqsdict, const, syms)
     K = Aaug.domain
 
     # sdm_irref has issues with float matrices. This uses the ddm_rref()
@@ -126,53 +129,55 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
         eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            eqdict[nsyms] = -elem_map[rhs]
         if eqdict:
             eqsdict.append(eqdict)
-    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)
+    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K)
     return sdm_aug
 
 
-def _expand_eqs_deprecated(eqs):
-    """Use expand to cancel nonlinear terms.
+def _linear_eq_to_dict(eqs, syms):
+    """Convert a system Expr/Eq equations into dict form, returning
+    the coefficient dictionaries and a list of syms-independent terms
+    from each expression in ``eqs```.
 
-    This approach matches previous behaviour of linsolve but should be
-    deprecated.
+    Examples
+    ========
+
+    >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict
+    >>> from sympy.abc import x
+    >>> _linear_eq_to_dict([2*x + 3], {x})
+    ([{x: 2}], [3])
     """
-    def expand_eq(eq):
-        if eq.is_Equality:
-            eq = eq.lhs - eq.rhs
-        return eq.expand()
+    coeffs = []
+    ind = []
+    symset = set(syms)
+    for i, e in enumerate(eqs):
+        c, d = _lin_eq2dict(e, symset)
+        coeffs.append(d)
+        ind.append(c)
+    return coeffs, ind
 
-    return [expand_eq(eq) for eq in eqs]
 
+def _lin_eq2dict(a, symset):
+    """return (c, d) where c is the sym-independent part of ``a`` and
+    ``d`` is an efficiently calculated dictionary mapping symbols to
+    their coefficients. A PolyNonlinearError is raised if non-linearity
+    is detected.
 
-def _linear_eq_to_dict(eqs, syms):
-    """Convert a system Expr/Eq equations into dict form"""
-    try:
-        return _linear_eq_to_dict_inner(eqs, syms)
-    except PolyNonlinearError:
-        # XXX: This should be deprecated:
-        eqs = _expand_eqs_deprecated(eqs)
-        return _linear_eq_to_dict_inner(eqs, syms)
-
-
-def _linear_eq_to_dict_inner(eqs, syms):
-    """Convert a system Expr/Eq equations into dict form"""
-    syms = set(syms)
-    eqsdict, eqs_rhs = [], []
-    for eq in eqs:
-        rhs, eqdict = _lin_eq2dict(eq, syms)
-        eqsdict.append(eqdict)
-        eqs_rhs.append(rhs)
-    return eqsdict, eqs_rhs
+    The values in the dictionary will be non-zero.
 
+    Examples
+    ========
 
-def _lin_eq2dict(a, symset):
-    """Efficiently convert a linear equation to a dict of coefficients"""
+    >>> from sympy.polys.matrices.linsolve import _lin_eq2dict
+    >>> from sympy.abc import x, y
+    >>> _lin_eq2dict(x + 2*y + 3, {x, y})
+    (3, {x: 1, y: 2})
+    """
     if a in symset:
         return S.Zero, {a: S.One}
-    elif a.is_Add:
+    if a.is_Add:
         terms_list = defaultdict(list)
         coeff_list = []
         for ai in a.args:
@@ -183,7 +188,7 @@ def _lin_eq2dict(a, symset):
         coeff = Add(*coeff_list)
         terms = {sym: Add(*coeffs) for sym, coeffs in terms_list.items()}
         return coeff, terms
-    elif a.is_Mul:
+    if a.is_Mul:
         terms = terms_coeff = None
         coeff_list = []
         for ai in a.args:
@@ -194,16 +199,30 @@ def _lin_eq2dict(a, symset):
                 terms = ti
                 terms_coeff = ci
             else:
-                raise PolyNonlinearError
-        coeff = Mul(*coeff_list)
+                # since ti is not null and we already have
+                # a term, this is a cross term
+                raise PolyNonlinearError(filldedent('''
+                    nonlinear cross-term: %s''' % a))
+        coeff = Mul._from_args(coeff_list)
         if terms is None:
             return coeff, {}
         else:
             terms = {sym: coeff * c for sym, c in terms.items()}
             return  coeff * terms_coeff, terms
-    elif a.is_Equality:
-        return _lin_eq2dict(a.lhs - a.rhs, symset)
-    elif not a.has_free(*symset):
+    if a.is_Equality:
+        (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+            for ai in a.args]
+        # there were no nonlinear errors so now
+        # cancellation is allowed
+        coeff -= cR
+        for k, v in tR.items():
+            if k in terms:
+                terms[k] -= v
+            else:
+                terms[k] = v
+        # don't store coefficients of 0, however
+        terms = {k: v for k, v in terms.items() if v}
+        return coeff, terms
+    if not a.has_xfree(symset):
         return a, {}
-    else:
-        raise PolyNonlinearError
+    raise PolyNonlinearError('nonlinear term: %s' % a)
diff --git a/sympy/polys/matrices/tests/test_linsolve.py b/sympy/polys/matrices/tests/test_linsolve.py
index e82f3f6af664..9d8cd7eb9feb 100644
--- a/sympy/polys/matrices/tests/test_linsolve.py
+++ b/sympy/polys/matrices/tests/test_linsolve.py
@@ -103,6 +103,9 @@ def all_close(sol1, sol2, eps=1e-15):
 
 
 def test__linsolve_deprecated():
-    assert _linsolve([Eq(x**2, x**2+y)], [x, y]) == {x:x, y:S.Zero}
-    assert _linsolve([(x+y)**2-x**2], [x]) == {x:-y/2}
-    assert _linsolve([Eq((x+y)**2, x**2)], [x]) == {x:-y/2}
+    raises(PolyNonlinearError, lambda:
+        _linsolve([Eq(x**2, x**2 + y)], [x, y]))
+    raises(PolyNonlinearError, lambda:
+        _linsolve([(x + y)**2 - x**2], [x]))
+    raises(PolyNonlinearError, lambda:
+        _linsolve([Eq((x + y)**2, x**2)], [x]))
diff --git a/sympy/solvers/ode/systems.py b/sympy/solvers/ode/systems.py
index 82bdea435e8d..98263983e6f7 100644
--- a/sympy/solvers/ode/systems.py
+++ b/sympy/solvers/ode/systems.py
@@ -197,7 +197,6 @@ def simpcoeff(coeff, wrt2):
         rep = {}
 
     sol = [Eq(s.lhs, simprhs(s.rhs, rep, wrt1, wrt2)) for s in sol]
-
     return sol
 
 
@@ -468,13 +467,7 @@ def linear_ode_to_matrix(eqs, funcs, t, order):
 
     for o in range(order, -1, -1):
         # Work from the highest derivative down
-        funcs_deriv = [func.diff(t, o) for func in funcs]
-
-        # linear_eq_to_matrix expects a proper symbol so substitute e.g.
-        # Derivative(x(t), t) for a Dummy.
-        rep = {func_deriv: Dummy() for func_deriv in funcs_deriv}
-        eqs = [eq.subs(rep) for eq in eqs]
-        syms = [rep[func_deriv] for func_deriv in funcs_deriv]
+        syms = [func.diff(t, o) for func in funcs]
 
         # Ai is the matrix for X(t).diff(t, o)
         # eqs is minus the remainder of the equations.
@@ -947,7 +940,7 @@ def linodesolve(A, t, b=None, B=None, type="auto", doit=False,
     # constants = numbered_symbols(prefix='C', cls=Dummy, start=const_idx+1)
     Cvect = Matrix(list(Dummy() for _ in range(n)))
 
-    if any(type == typ for typ in ["type2", "type4", "type6"]) and b is None:
+    if b is None and any(type == typ for typ in ["type2", "type4", "type6"]):
         b = zeros(n, 1)
 
     is_transformed = tau is not None
@@ -973,6 +966,7 @@ def linodesolve(A, t, b=None, B=None, type="auto", doit=False,
                 A = system_info['A']
                 b = system_info['b']
 
+    intx_wrtt = lambda x: Integral(x, t) if x else 0
     if type in ("type1", "type2", "type5", "type6"):
         P, J = matrix_exp_jordan_form(A, t)
         P = simplify(P)
@@ -981,8 +975,7 @@ def linodesolve(A, t, b=None, B=None, type="auto", doit=False,
             sol_vector = P * (J * Cvect)
         else:
             Jinv = J.subs(t, -t)
-            sol_vector = P * J * ((Jinv * P.inv() * b).applyfunc(lambda x: Integral(x, t)) + Cvect)
-
+            sol_vector = P * J * ((Jinv * P.inv() * b).applyfunc(intx_wrtt) + Cvect)
     else:
         if B is None:
             B, _ = _is_commutative_anti_derivative(A, t)
@@ -990,7 +983,7 @@ def linodesolve(A, t, b=None, B=None, type="auto", doit=False,
         if type == "type3":
             sol_vector = B.exp() * Cvect
         else:
-            sol_vector = B.exp() * (((-B).exp() * b).applyfunc(lambda x: Integral(x, t)) + Cvect)
+            sol_vector = B.exp() * (((-B).exp() * b).applyfunc(intx_wrtt) + Cvect)
 
     if is_transformed:
         sol_vector = sol_vector.subs(t, tau)
diff --git a/sympy/solvers/ode/tests/test_systems.py b/sympy/solvers/ode/tests/test_systems.py
index dd08d644f7a2..2a2b1e155c6d 100644
--- a/sympy/solvers/ode/tests/test_systems.py
+++ b/sympy/solvers/ode/tests/test_systems.py
@@ -1071,8 +1071,8 @@ def test_sysode_linear_neq_order1_type2():
     eqs6 = [Eq(Derivative(f(x), x), -9*f(x) - 4*g(x)),
             Eq(Derivative(g(x), x), -4*g(x)),
             Eq(Derivative(h(x), x), h(x) + exp(x))]
-    sol6 = [Eq(f(x), C1*exp(-4*x)*Rational(-4, 5) + C2*exp(-9*x)),
-            Eq(g(x), C1*exp(-4*x)),
+    sol6 = [Eq(f(x), C2*exp(-4*x)*Rational(-4, 5) + C1*exp(-9*x)),
+            Eq(g(x), C2*exp(-4*x)),
             Eq(h(x), C3*exp(x) + x*exp(x))]
     assert dsolve(eqs6) == sol6
     assert checksysodesol(eqs6, sol6) == (True, [0, 0, 0])
@@ -1647,11 +1647,11 @@ def test_higher_order_to_first_order_9():
 
     eqs9 = [f(x) + g(x) - 2*exp(I*x) + 2*Derivative(f(x), x) + Derivative(f(x), (x, 2)),
             f(x) + g(x) - 2*exp(I*x) + 2*Derivative(g(x), x) + Derivative(g(x), (x, 2))]
-    sol9 =  [Eq(f(x), -C1 + C2*exp(-2*x)/2 - (C3/2 - C4/2)*exp(-x)*cos(x)
-                    + (C3/2 + C4/2)*exp(-x)*sin(x) + 2*((1 - 2*I)*exp(I*x)*sin(x)**2/5)
+    sol9 =  [Eq(f(x), -C1 + C4*exp(-2*x)/2 - (C2/2 - C3/2)*exp(-x)*cos(x)
+                    + (C2/2 + C3/2)*exp(-x)*sin(x) + 2*((1 - 2*I)*exp(I*x)*sin(x)**2/5)
                     + 2*((1 - 2*I)*exp(I*x)*cos(x)**2/5)),
-            Eq(g(x), C1 - C2*exp(-2*x)/2 - (C3/2 - C4/2)*exp(-x)*cos(x)
-                    + (C3/2 + C4/2)*exp(-x)*sin(x) + 2*((1 - 2*I)*exp(I*x)*sin(x)**2/5)
+            Eq(g(x), C1 - C4*exp(-2*x)/2 - (C2/2 - C3/2)*exp(-x)*cos(x)
+                    + (C2/2 + C3/2)*exp(-x)*sin(x) + 2*((1 - 2*I)*exp(I*x)*sin(x)**2/5)
                     + 2*((1 - 2*I)*exp(I*x)*cos(x)**2/5))]
     assert dsolve(eqs9) == sol9
     assert checksysodesol(eqs9, sol9) == (True, [0, 0])
@@ -1977,11 +1977,8 @@ def test_linodesolve():
 
     # non-homogeneous term assumed to be 0
     sol1 = [-C1*exp(-t/2 + sqrt(5)*t/2)/2 + sqrt(5)*C1*exp(-t/2 + sqrt(5)*t/2)/2 - sqrt(5)*C2*exp(-sqrt(5)*t/2
-                - t/2)/2 - C2*exp(-sqrt(5)*t/2 - t/2)/2 - exp(-t/2 + sqrt(5)*t/2)*Integral(0, t)/2 +
-                sqrt(5)*exp(-t/2 + sqrt(5)*t/2)*Integral(0, t)/2 - sqrt(5)*exp(-sqrt(5)*t/2 - t/2)*Integral(0, t)/2
-                - exp(-sqrt(5)*t/2 - t/2)*Integral(0, t)/2,
-            C1*exp(-t/2 + sqrt(5)*t/2) + C2*exp(-sqrt(5)*t/2 - t/2)
-                + exp(-t/2 + sqrt(5)*t/2)*Integral(0, t) + exp(-sqrt(5)*t/2 - t/2)*Integral(0, t)]
+                - t/2)/2 - C2*exp(-sqrt(5)*t/2 - t/2)/2,
+            C1*exp(-t/2 + sqrt(5)*t/2) + C2*exp(-sqrt(5)*t/2 - t/2)]
     assert constant_renumber(linodesolve(A, t, type="type2"), variables=[t]) == sol1
 
     # Testing the Errors
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index b0c49825a988..c19e88868803 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -12,7 +12,7 @@
 """
 from sympy.core.sympify import sympify
 from sympy.core import (S, Pow, Dummy, pi, Expr, Wild, Mul, Equality,
-                        Add)
+                        Add, Basic)
 from sympy.core.containers import Tuple
 from sympy.core.function import (Lambda, expand_complex, AppliedUndef,
                                 expand_log, _mexpand, expand_trig, nfloat)
@@ -23,7 +23,7 @@
 from sympy.core.sorting import default_sort_key, ordered
 from sympy.core.symbol import Symbol, _uniquely_named_symbol
 from sympy.core.sympify import _sympify
-from sympy.core.traversal import iterfreeargs
+from sympy.polys.matrices.linsolve import _linear_eq_to_dict, _lin_eq2dict
 from sympy.polys.polyroots import UnsolvableFactorError
 from sympy.simplify.simplify import simplify, fraction, trigsimp, nsimplify
 from sympy.simplify import powdenest, logcombine
@@ -38,7 +38,7 @@
 from sympy.sets import (FiniteSet, imageset, Interval, Intersection,
                         Union, ConditionSet, ImageSet, Complement, Contains)
 from sympy.sets.sets import Set, ProductSet
-from sympy.matrices import Matrix, MatrixBase
+from sympy.matrices import SparseMatrix, MatrixBase
 from sympy.ntheory import totient
 from sympy.ntheory.factor_ import divisors
 from sympy.ntheory.residue_ntheory import discrete_log, nthroot_mod
@@ -54,11 +54,10 @@
 from sympy.solvers.polysys import solve_poly_system
 from sympy.utilities import filldedent
 from sympy.utilities.iterables import (numbered_symbols, has_dups,
-                                       is_sequence)
+                                       is_sequence, iterable)
 from sympy.calculus.util import periodicity, continuous_domain, function_range
 
 from types import GeneratorType
-from collections import defaultdict
 
 
 class NonlinearError(ValueError):
@@ -2411,7 +2410,7 @@ def solvify(f, symbol, domain):
 ###############################################################################
 
 
-def linear_coeffs(eq, *syms, **_kw):
+def linear_coeffs(eq, *syms, dict=False):
     """Return a list whose elements are the coefficients of the
     corresponding symbols in the sum of terms in  ``eq``.
     The additive constant is returned as the last element of the
@@ -2422,69 +2421,82 @@ def linear_coeffs(eq, *syms, **_kw):
 
     NonlinearError
         The equation contains a nonlinear term
+    ValueError
+        duplicate or unordered symbols are passed
+
+    Parameters
+    ==========
+
+    dict - (default False) when True, return coefficients as a
+        dictionary with coefficients keyed to syms that were present;
+        key 1 gives the constant term
 
     Examples
     ========
 
     >>> from sympy.solvers.solveset import linear_coeffs
     >>> from sympy.abc import x, y, z
-
     >>> linear_coeffs(3*x + 2*y - 1, x, y)
     [3, 2, -1]
 
     It is not necessary to expand the expression:
 
-    >>> linear_coeffs(x + y*(z*(x*3 + 2) + 3), x)
-    [3*y*z + 1, y*(2*z + 3)]
+        >>> linear_coeffs(x + y*(z*(x*3 + 2) + 3), x)
+        [3*y*z + 1, y*(2*z + 3)]
 
-    But if there are nonlinear or cross terms -- even if they would
-    cancel after simplification -- an error is raised so the situation
-    does not pass silently past the caller's attention:
+    When nonlinear is detected, an error will be raised:
 
-    >>> eq = 1/x*(x - 1) + 1/x
-    >>> linear_coeffs(eq.expand(), x)
-    [0, 1]
-    >>> linear_coeffs(eq, x)
-    Traceback (most recent call last):
-    ...
-    NonlinearError: nonlinear term encountered: 1/x
+        * even if they would cancel after expansion (so the
+        situation does not pass silently past the caller's
+        attention)
 
-    >>> linear_coeffs(x*(y + 1) - x*y, x, y)
-    Traceback (most recent call last):
-    ...
-    NonlinearError: nonlinear term encountered: x*(y + 1)
+        >>> eq = 1/x*(x - 1) + 1/x
+        >>> linear_coeffs(eq.expand(), x)
+        [0, 1]
+        >>> linear_coeffs(eq, x)
+        Traceback (most recent call last):
+        ...
+        NonlinearError:
+        nonlinear in given generators
+
+        * when there are cross terms
+
+        >>> linear_coeffs(x*(y + 1), x, y)
+        Traceback (most recent call last):
+        ...
+        NonlinearError:
+        symbol-dependent cross-terms encountered
+
+        * when there are terms that contain an expression
+        dependent on the symbols that is not linear
+
+        >>> linear_coeffs(x**2, x)
+        Traceback (most recent call last):
+        ...
+        NonlinearError:
+        nonlinear in given generators
     """
-    d = defaultdict(list)
     eq = _sympify(eq)
+    if len(syms) == 1 and iterable(syms[0]) and not isinstance(syms[0], Basic):
+        raise ValueError('expecting unpacked symbols, *syms')
     symset = set(syms)
     if len(symset) != len(syms):
         raise ValueError('duplicate symbols given')
-    has = set(iterfreeargs(eq)) & symset
-    if not has:
-        return [S.Zero]*len(syms) + [eq]
-    c, terms = eq.as_coeff_add(*has)
-    d[0].extend(Add.make_args(c))
-    for t in terms:
-        m, f = t.as_coeff_mul(*has)
-        if len(f) != 1:
-            break
-        f = f[0]
-        if f in symset:
-            d[f].append(m)
-        elif f.is_Add:
-            d1 = linear_coeffs(f, *has, **{'dict': True})
-            d[0].append(m*d1.pop(0))
-            for xf, vf in d1.items():
-                d[xf].append(m*vf)
-        else:
-            break
-    else:
-        for k, v in d.items():
-            d[k] = Add(*v)
-        if not _kw:
-            return [d.get(s, S.Zero) for s in syms]+ [d[0]]
-        return d  # default is still list but this won't matter
-    raise NonlinearError('nonlinear term encountered: %s' % t)
+    try:
+        c, d = _lin_eq2dict(eq, symset)
+    except PolyNonlinearError as err:
+        raise NonlinearError(str(err))
+    if dict:
+        if c:
+            d[S.One] = c
+        return d
+    rv = [S.Zero]*(len(syms) + 1)
+    rv[-1] = c
+    for i, k in enumerate(syms):
+        if k not in d:
+            continue
+        rv[i] = d[k]
+    return rv
 
 
 def linear_eq_to_matrix(equations, *symbols):
@@ -2531,39 +2543,39 @@ def linear_eq_to_matrix(equations, *symbols):
     The coefficients (numerical or symbolic) of the symbols will
     be returned as matrices:
 
-    >>> eqns = [c*x + z - 1 - c, y + z, x - y]
-    >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
-    >>> A
-    Matrix([
-    [c,  0, 1],
-    [0,  1, 1],
-    [1, -1, 0]])
-    >>> b
-    Matrix([
-    [c + 1],
-    [    0],
-    [    0]])
+        >>> eqns = [c*x + z - 1 - c, y + z, x - y]
+        >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
+        >>> A
+        Matrix([
+        [c,  0, 1],
+        [0,  1, 1],
+        [1, -1, 0]])
+        >>> b
+        Matrix([
+        [c + 1],
+        [    0],
+        [    0]])
 
     This routine does not simplify expressions and will raise an error
     if nonlinearity is encountered:
 
-    >>> eqns = [
-    ...     (x**2 - 3*x)/(x - 3) - 3,
-    ...     y**2 - 3*y - y*(y - 4) + x - 4]
-    >>> linear_eq_to_matrix(eqns, [x, y])
-    Traceback (most recent call last):
-    ...
-    NonlinearError:
-    The term (x**2 - 3*x)/(x - 3) is nonlinear in {x, y}
+            >>> eqns = [
+            ...     (x**2 - 3*x)/(x - 3) - 3,
+            ...     y**2 - 3*y - y*(y - 4) + x - 4]
+            >>> linear_eq_to_matrix(eqns, [x, y])
+            Traceback (most recent call last):
+            ...
+            NonlinearError:
+            symbol-dependent term can be ignored using `strict=False`
 
-    Simplifying these equations will discard the removable singularity
-    in the first, reveal the linear structure of the second:
+        Simplifying these equations will discard the removable singularity
+        in the first and reveal the linear structure of the second:
 
-    >>> [e.simplify() for e in eqns]
-    [x - 3, x + y - 4]
+            >>> [e.simplify() for e in eqns]
+            [x - 3, x + y - 4]
 
-    Any such simplification needed to eliminate nonlinear terms must
-    be done before calling this routine.
+        Any such simplification needed to eliminate nonlinear terms must
+        be done *before* calling this routine.
     """
     if not symbols:
         raise ValueError(filldedent('''
@@ -2574,12 +2586,6 @@ def linear_eq_to_matrix(equations, *symbols):
     if hasattr(symbols[0], '__iter__'):
         symbols = symbols[0]
 
-    for i in symbols:
-        if not isinstance(i, Symbol):
-            raise ValueError(filldedent('''
-            Expecting a Symbol but got %s
-            ''' % i))
-
     if has_dups(symbols):
         raise ValueError('Symbols must be unique')
 
@@ -2594,14 +2600,19 @@ def linear_eq_to_matrix(equations, *symbols):
             Eq or Matrix.
             '''))
 
-    A, b = [], []
-    for i, f in enumerate(equations):
-        if isinstance(f, Equality):
-            f = f.rewrite(Add, evaluate=False)
-        coeff_list = linear_coeffs(f, *symbols)
-        b.append(-coeff_list.pop())
-        A.append(coeff_list)
-    A, b = map(Matrix, (A, b))
+    # construct the dictionaries
+    try:
+        eq, c = _linear_eq_to_dict(equations, symbols)
+    except PolyNonlinearError as err:
+        raise NonlinearError(str(err))
+    # prepare output matrices
+    n, m = shape = len(eq), len(symbols)
+    ix = dict(zip(symbols, range(m)))
+    dat = {(row, ix[k]): d[k] for row, d in enumerate(eq) for k in d}
+    rhs = [-i for i in c]
+    del c
+    A = SparseMatrix(*shape, dat)
+    b = SparseMatrix(n, 1, rhs)
     return A, b
 
 
@@ -2760,16 +2771,23 @@ def linsolve(system, *symbols):
     >>> linsolve([], x)
     EmptySet
 
-    * An error is raised if, after expansion, any nonlinearity
-      is detected:
+    * An error is raised if any nonlinearity is detected, even
+    if it could be removed with expansion
+
+    >>> linsolve([x*(1/x - 1)], x)
+    Traceback (most recent call last):
+    ...
+    NonlinearError: nonlinear term: 1/x
+
+    >>> linsolve([x*(y + 1)], x, y)
+    Traceback (most recent call last):
+    ...
+    NonlinearError: nonlinear cross-term: x*(y + 1)
 
-    >>> linsolve([x*(1/x - 1), (y - 1)**2 - y**2 + 1], x, y)
-    {(1, 1)}
     >>> linsolve([x**2 - 1], x)
     Traceback (most recent call last):
     ...
-    NonlinearError:
-    nonlinear term encountered: x**2
+    NonlinearError: nonlinear term: x**2
     """
     if not system:
         return S.EmptySet
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 87def61c7632..64676fcc3e71 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1362,6 +1362,10 @@ def test_abs_invert_solvify():
 
 
 def test_linear_eq_to_matrix():
+    assert linear_eq_to_matrix(0, x) == (Matrix([[0]]), Matrix([[0]]))
+    assert linear_eq_to_matrix(1, x) == (Matrix([[0]]), Matrix([[-1]]))
+
+    # integer coefficients
     eqns1 = [2*x + y - 2*z - 3, x - y - z, x + y + 3*z - 12]
     eqns2 = [Eq(3*x + 2*y - z, 1), Eq(2*x - 2*y + 4*z, -2), -2*x + y - 2*z]
 
@@ -1379,17 +1383,24 @@ def test_linear_eq_to_matrix():
     assert A == Matrix([[a*b, b, c], [d + e, f, g], [i, j, k]])
     assert B == Matrix([[d], [h], [l]])
 
-    # raise ValueError if
+    # raise Errors if
     # 1) no symbols are given
     raises(ValueError, lambda: linear_eq_to_matrix(eqns3))
     # 2) there are duplicates
     raises(ValueError, lambda: linear_eq_to_matrix(eqns3, [x, x, y]))
-    # 3) there are non-symbols
-    raises(ValueError, lambda: linear_eq_to_matrix(eqns3, [x, 1/a, y]))
-    # 4) a nonlinear term is detected in the original expression
+    # 3) a nonlinear term is detected in the original expression
     raises(NonlinearError, lambda: linear_eq_to_matrix(Eq(1/x + x, 1/x), [x]))
+    raises(NonlinearError, lambda: linear_eq_to_matrix([x**2], [x]))
+    raises(NonlinearError, lambda: linear_eq_to_matrix([x*y], [x, y]))
+    # 4) Eq being used to represent equations autoevaluates
+    # (use unevaluated Eq instead)
+    raises(ValueError, lambda: linear_eq_to_matrix(Eq(x, x), x))
+    raises(ValueError, lambda: linear_eq_to_matrix(Eq(x, x + 1), x))
+
+
+    # if non-symbols are passed, the user is responsible for interpreting
+    assert linear_eq_to_matrix([x], [1/x]) == (Matrix([[0]]), Matrix([[-x]]))
 
-    assert linear_eq_to_matrix(1, x) == (Matrix([[0]]), Matrix([[-1]]))
     # issue 15195
     assert linear_eq_to_matrix(x + y*(z*(3*x + 2) + 3), x) == (
         Matrix([[3*y*z + 1]]), Matrix([[-y*(2*z + 3)]]))
@@ -1502,12 +1513,11 @@ def test_linsolve():
     assert linsolve(Eqns, x, y) == {
             (kilo*newton*Rational(-28, 3), kN*Rational(4, 3))}
 
-    # linsolve fully expands expressions, so removable singularities
-    # and other nonlinearity does not raise an error
+    # linsolve does not allow expansion (real or implemented)
+    # to remove singularities, but it will cancel linear terms
     assert linsolve([Eq(x, x + y)], [x, y]) == {(x, 0)}
-    assert linsolve([Eq(1/x, 1/x + y)], [x, y]) == {(x, 0)}
-    assert linsolve([Eq(y/x, y/x + y)], [x, y]) == {(x, 0)}
-    assert linsolve([Eq(x*(x + 1), x**2 + y)], [x, y]) == {(y, y)}
+    raises(NonlinearError, lambda:
+        linsolve([Eq(x**2, x**2 + y)], [x, y]))
 
     # corner cases
     #
@@ -2785,8 +2795,11 @@ def test_linear_coeffs():
         linear_coeffs(x, x, x))
     assert linear_coeffs(a*(x + y), x, y) == [a, a, 0]
     assert linear_coeffs(1.0, x, y) == [0, 0, 1.0]
+    # don't include coefficients of 0
+    assert linear_coeffs(Eq(x, x + y), x, y, dict=True) == {y: 1}
+    assert linear_coeffs(0, x, y, dict=True) == {}
+
 
-# modular tests
 def test_is_modular():
     assert _is_modular(y, x) is False
     assert _is_modular(Mod(x, 3) - 1, x) is True


Code review comments:
diff_hunk: @@ -2594,14 +2600,19 @@ def linear_eq_to_matrix(equations, *symbols):
             Eq or Matrix.
             '''))
 
-    A, b = [], []
-    for i, f in enumerate(equations):
-        if isinstance(f, Equality):
-            f = f.rewrite(Add, evaluate=False)
-        coeff_list = linear_coeffs(f, *symbols)
-        b.append(-coeff_list.pop())
-        A.append(coeff_list)
-    A, b = map(Matrix, (A, b))
+    # construct the dictionaries
+    try:
+        eq, c = _linear_eq_to_dict(equations, symbols)
+    except PolyNonlinearError as err:
+        raise NonlinearError(str(err))
+    # prepare output matrices
+    n, m = shape = len(eq), len(symbols)
+    ix = dict(zip(symbols, range(m)))
+    dat = {(row, ix[k]): d[k] for row, d in enumerate(eq) for k in d}
+    rhs = [-i for i in c]
+    del c
+    A = SparseMatrix(*shape, dat)
+    b = SparseMatrix(n, 1, rhs)
path: sympy/solvers/solveset.py
text: @user1:
I pretty much want to abolish SparseMatrix so it would be better not to add new uses of it.

diff_hunk: @@ -1258,6 +1258,33 @@ def has(self, *patterns):
         """
         return self._has(iterargs, *patterns)
 
+    def has_xfree(self, s):
path: sympy/core/basic.py
text: @user1:
Maybe this should be type-hinted as `s: set[Basic]`.

@author:
I don't follow you. `s` should be a set, not Basic.

@user1:
It should be a `set` of `Basic` which is what `set[Basic]` means in typing syntax.

diff_hunk: @@ -194,16 +199,30 @@ def _lin_eq2dict(a, symset):
                 terms = ti
                 terms_coeff = ci
             else:
-                raise PolyNonlinearError
-        coeff = Mul(*coeff_list)
+                # since ti is not null and we already have
+                # a term, this is a cross term
+                raise PolyNonlinearError(filldedent('''
+                    nonlinear cross-term: %s''' % a))
+        coeff = Mul._from_args(coeff_list)
         if terms is None:
             return coeff, {}
         else:
             terms = {sym: coeff * c for sym, c in terms.items()}
             return  coeff * terms_coeff, terms
-    elif a.is_Equality:
-        return _lin_eq2dict(a.lhs - a.rhs, symset)
-    elif not a.has_free(*symset):
+    if a.is_Equality:
+        (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+            for ai in a.args]
path: sympy/polys/matrices/linsolve.py
text: @user1:
Maybe equalities should be handled at a higher level before calling this function.

@author:
I had already considered that. It is here because otherwise anyone using any of the linear coefficient extracting routines would have to duplicate this code to carefully make sure there is no linearity and then allow cancellation of linear terms.

@user1:
It could be handled by `_linear_eq_to_dict` for example.

The checking for zeros probably needs to be handled more broadly than just for `Equality`.

diff_hunk: @@ -126,53 +129,70 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
         eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            eqdict[nsyms] = -elem_map[rhs]
         if eqdict:
             eqsdict.append(eqdict)
-    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)
+    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K)
     return sdm_aug
 
 
-def _expand_eqs_deprecated(eqs):
-    """Use expand to cancel nonlinear terms.
+def _linear_eq_to_dict(eqs, syms):
+    """Convert a system Expr/Eq equations into dict form, returning
+    the coefficient dictionaries and a list of syms-independent terms
+    from each expression in ``eqs```.
 
-    This approach matches previous behaviour of linsolve but should be
-    deprecated.
+    Examples
+    ========
+
+    >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict
+    >>> from sympy.abc import x
+    >>> _linear_eq_to_dict([2*x + 3], {x})
+    ([{x: 2}], [3])
     """
-    def expand_eq(eq):
-        if eq.is_Equality:
-            eq = eq.lhs - eq.rhs
-        return eq.expand()
+    coeffs = []
+    ind = []
+    symset = set(syms)
+    for i, e in enumerate(eqs):
+        if e.is_Equality:
+            (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+                for ai in e.args]
path: sympy/polys/matrices/linsolve.py
text: @user1:
A loop to avoid repeating a line just twice is not really a good application of DRY. It's a lot easier to understand the straight-forward repetitive code:
```python
coeff, terms = _lin_eq2dict(e.lhs, symset)
cR, tR = _lin_eq2dict(e.rhs, symset)
```

diff_hunk: @@ -126,53 +129,70 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
         eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            eqdict[nsyms] = -elem_map[rhs]
         if eqdict:
             eqsdict.append(eqdict)
-    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)
+    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K)
     return sdm_aug
 
 
-def _expand_eqs_deprecated(eqs):
-    """Use expand to cancel nonlinear terms.
+def _linear_eq_to_dict(eqs, syms):
+    """Convert a system Expr/Eq equations into dict form, returning
+    the coefficient dictionaries and a list of syms-independent terms
+    from each expression in ``eqs```.
 
-    This approach matches previous behaviour of linsolve but should be
-    deprecated.
+    Examples
+    ========
+
+    >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict
+    >>> from sympy.abc import x
+    >>> _linear_eq_to_dict([2*x + 3], {x})
+    ([{x: 2}], [3])
     """
-    def expand_eq(eq):
-        if eq.is_Equality:
-            eq = eq.lhs - eq.rhs
-        return eq.expand()
+    coeffs = []
+    ind = []
+    symset = set(syms)
+    for i, e in enumerate(eqs):
+        if e.is_Equality:
+            (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+                for ai in e.args]
+            # there were no nonlinear errors so now
+            # cancellation is allowed
+            coeff -= cR
+            for k, v in tR.items():
+                if k in terms:
+                    terms[k] -= v
+                else:
+                    terms[k] = v
path: sympy/polys/matrices/linsolve.py
text: @user1:
Should this not be `-v`?

@author:
Hmm.. that should be tested.

diff_hunk: @@ -126,53 +129,70 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
         eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            eqdict[nsyms] = -elem_map[rhs]
         if eqdict:
             eqsdict.append(eqdict)
-    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)
+    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K)
     return sdm_aug
 
 
-def _expand_eqs_deprecated(eqs):
-    """Use expand to cancel nonlinear terms.
+def _linear_eq_to_dict(eqs, syms):
+    """Convert a system Expr/Eq equations into dict form, returning
+    the coefficient dictionaries and a list of syms-independent terms
+    from each expression in ``eqs```.
 
-    This approach matches previous behaviour of linsolve but should be
-    deprecated.
+    Examples
+    ========
+
+    >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict
+    >>> from sympy.abc import x
+    >>> _linear_eq_to_dict([2*x + 3], {x})
+    ([{x: 2}], [3])
     """
-    def expand_eq(eq):
-        if eq.is_Equality:
-            eq = eq.lhs - eq.rhs
-        return eq.expand()
+    coeffs = []
+    ind = []
+    symset = set(syms)
+    for i, e in enumerate(eqs):
+        if e.is_Equality:
+            (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+                for ai in e.args]
+            # there were no nonlinear errors so now
+            # cancellation is allowed
+            coeff -= cR
+            for k, v in tR.items():
+                if k in terms:
+                    terms[k] -= v
+                else:
+                    terms[k] = v
+            # don't store coefficients of 0, however
+            terms = {k: v for k, v in terms.items() if v}
+            c, d = coeff, terms
+        else:
+            c, d = _lin_eq2dict(e, symset)
+        coeffs.append(d)
+        ind.append(c)
+    return coeffs, ind
 
-    return [expand_eq(eq) for eq in eqs]
 
+def _lin_eq2dict(a, symset):
+    """return (c, d) where c is the sym-independent part of ``a`` and
+    ``d`` is an efficiently calculated dictionary mapping symbols to
+    their coefficients. A PolyNonlinearError is raised if non-linearity
+    is detected.
 
-def _linear_eq_to_dict(eqs, syms):
-    """Convert a system Expr/Eq equations into dict form"""
-    try:
-        return _linear_eq_to_dict_inner(eqs, syms)
-    except PolyNonlinearError:
-        # XXX: This should be deprecated:
-        eqs = _expand_eqs_deprecated(eqs)
-        return _linear_eq_to_dict_inner(eqs, syms)
-
-
-def _linear_eq_to_dict_inner(eqs, syms):
-    """Convert a system Expr/Eq equations into dict form"""
-    syms = set(syms)
-    eqsdict, eqs_rhs = [], []
-    for eq in eqs:
-        rhs, eqdict = _lin_eq2dict(eq, syms)
-        eqsdict.append(eqdict)
-        eqs_rhs.append(rhs)
-    return eqsdict, eqs_rhs
+    The values in the dictionary will be non-zero.
 
+    Examples
+    ========
 
-def _lin_eq2dict(a, symset):
-    """Efficiently convert a linear equation to a dict of coefficients"""
+    >>> from sympy.polys.matrices.linsolve import _lin_eq2dict
+    >>> from sympy.abc import x, y
+    >>> _lin_eq2dict(x + 2*y + 3, {x, y})
+    (3, {x: 1, y: 2})
+    """
     if a in symset:
         return S.Zero, {a: S.One}
-    elif a.is_Add:
+    if a.is_Add:
path: sympy/polys/matrices/linsolve.py
text: @user1:
Why have you changed these `elif` to `if`? It's much clearer to use `elif` to show that all cases are mutually exclusive without needing to scan through to see that every control path either returns or raises.

@author:
Preference. If every block clearly ends with a return (or starts with an elif) coverage is complete.

diff_hunk: @@ -126,53 +129,70 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
         eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            eqdict[nsyms] = -elem_map[rhs]
         if eqdict:
             eqsdict.append(eqdict)
-    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)
+    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K)
     return sdm_aug
 
 
-def _expand_eqs_deprecated(eqs):
-    """Use expand to cancel nonlinear terms.
+def _linear_eq_to_dict(eqs, syms):
+    """Convert a system Expr/Eq equations into dict form, returning
+    the coefficient dictionaries and a list of syms-independent terms
+    from each expression in ``eqs```.
 
-    This approach matches previous behaviour of linsolve but should be
-    deprecated.
+    Examples
+    ========
+
+    >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict
+    >>> from sympy.abc import x
+    >>> _linear_eq_to_dict([2*x + 3], {x})
+    ([{x: 2}], [3])
     """
-    def expand_eq(eq):
-        if eq.is_Equality:
-            eq = eq.lhs - eq.rhs
-        return eq.expand()
+    coeffs = []
+    ind = []
+    symset = set(syms)
+    for i, e in enumerate(eqs):
+        if e.is_Equality:
+            (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+                for ai in e.args]
path: sympy/polys/matrices/linsolve.py
text: @author:
```suggestion
            coeff, terms = _lin_eq2dict(e.lhs, symset)
            cR, tR = _lin_eq2dict(e.rhs, symset)
```

diff_hunk: @@ -126,53 +129,70 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
         eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            eqdict[nsyms] = -elem_map[rhs]
         if eqdict:
             eqsdict.append(eqdict)
-    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)
+    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K)
     return sdm_aug
 
 
-def _expand_eqs_deprecated(eqs):
-    """Use expand to cancel nonlinear terms.
+def _linear_eq_to_dict(eqs, syms):
+    """Convert a system Expr/Eq equations into dict form, returning
+    the coefficient dictionaries and a list of syms-independent terms
+    from each expression in ``eqs```.
 
-    This approach matches previous behaviour of linsolve but should be
-    deprecated.
+    Examples
+    ========
+
+    >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict
+    >>> from sympy.abc import x
+    >>> _linear_eq_to_dict([2*x + 3], {x})
+    ([{x: 2}], [3])
     """
-    def expand_eq(eq):
-        if eq.is_Equality:
-            eq = eq.lhs - eq.rhs
-        return eq.expand()
+    coeffs = []
+    ind = []
+    symset = set(syms)
+    for i, e in enumerate(eqs):
+        if e.is_Equality:
+            (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+                for ai in e.args]
+            # there were no nonlinear errors so now
+            # cancellation is allowed
+            coeff -= cR
+            for k, v in tR.items():
+                if k in terms:
+                    terms[k] -= v
+                else:
+                    terms[k] = v
path: sympy/polys/matrices/linsolve.py
text: @author:
```suggestion
                    terms[k] = -v
```

diff_hunk: @@ -126,53 +129,70 @@ def sympy_dict_to_dm(eqs_coeffs, eqs_rhs, syms):
     for eq, rhs in zip(eqs_coeffs, eqs_rhs):
         eqdict = {sym2index[s]: elem_map[c] for s, c in eq.items()}
         if rhs:
-            eqdict[nsyms] = - elem_map[rhs]
+            eqdict[nsyms] = -elem_map[rhs]
         if eqdict:
             eqsdict.append(eqdict)
-    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms+1), K)
+    sdm_aug = SDM(enumerate(eqsdict), (neqs, nsyms + 1), K)
     return sdm_aug
 
 
-def _expand_eqs_deprecated(eqs):
-    """Use expand to cancel nonlinear terms.
+def _linear_eq_to_dict(eqs, syms):
+    """Convert a system Expr/Eq equations into dict form, returning
+    the coefficient dictionaries and a list of syms-independent terms
+    from each expression in ``eqs```.
 
-    This approach matches previous behaviour of linsolve but should be
-    deprecated.
+    Examples
+    ========
+
+    >>> from sympy.polys.matrices.linsolve import _linear_eq_to_dict
+    >>> from sympy.abc import x
+    >>> _linear_eq_to_dict([2*x + 3], {x})
+    ([{x: 2}], [3])
     """
-    def expand_eq(eq):
-        if eq.is_Equality:
-            eq = eq.lhs - eq.rhs
-        return eq.expand()
+    coeffs = []
+    ind = []
+    symset = set(syms)
+    for i, e in enumerate(eqs):
+        if e.is_Equality:
+            (coeff, terms), (cR, tR) = [_lin_eq2dict(ai, symset)
+                for ai in e.args]
+            # there were no nonlinear errors so now
+            # cancellation is allowed
+            coeff -= cR
+            for k, v in tR.items():
+                if k in terms:
+                    terms[k] -= v
+                else:
+                    terms[k] = v
+            # don't store coefficients of 0, however
+            terms = {k: v for k, v in terms.items() if v}
+            c, d = coeff, terms
+        else:
+            c, d = _lin_eq2dict(e, symset)
+        coeffs.append(d)
+        ind.append(c)
+    return coeffs, ind
 
-    return [expand_eq(eq) for eq in eqs]
 
+def _lin_eq2dict(a, symset):
+    """return (c, d) where c is the sym-independent part of ``a`` and
+    ``d`` is an efficiently calculated dictionary mapping symbols to
+    their coefficients. A PolyNonlinearError is raised if non-linearity
+    is detected.
 
-def _linear_eq_to_dict(eqs, syms):
-    """Convert a system Expr/Eq equations into dict form"""
-    try:
-        return _linear_eq_to_dict_inner(eqs, syms)
-    except PolyNonlinearError:
-        # XXX: This should be deprecated:
-        eqs = _expand_eqs_deprecated(eqs)
-        return _linear_eq_to_dict_inner(eqs, syms)
-
-
-def _linear_eq_to_dict_inner(eqs, syms):
-    """Convert a system Expr/Eq equations into dict form"""
-    syms = set(syms)
-    eqsdict, eqs_rhs = [], []
-    for eq in eqs:
-        rhs, eqdict = _lin_eq2dict(eq, syms)
-        eqsdict.append(eqdict)
-        eqs_rhs.append(rhs)
-    return eqsdict, eqs_rhs
+    The values in the dictionary will be non-zero.
 
+    Examples
+    ========
 
-def _lin_eq2dict(a, symset):
-    """Efficiently convert a linear equation to a dict of coefficients"""
+    >>> from sympy.polys.matrices.linsolve import _lin_eq2dict
+    >>> from sympy.abc import x, y
+    >>> _lin_eq2dict(x + 2*y + 3, {x, y})
+    (3, {x: 1, y: 2})
+    """
     if a in symset:
         return S.Zero, {a: S.One}
-    elif a.is_Add:
+    if a.is_Add:
path: sympy/polys/matrices/linsolve.py
text: @author:
```suggestion
    elif a.is_Add:
```

diff_hunk: @@ -183,7 +203,7 @@ def _lin_eq2dict(a, symset):
         coeff = Add(*coeff_list)
         terms = {sym: Add(*coeffs) for sym, coeffs in terms_list.items()}
         return coeff, terms
-    elif a.is_Mul:
+    if a.is_Mul:
path: sympy/polys/matrices/linsolve.py
text: @author:
```suggestion
    elif a.is_Mul:
```

diff_hunk: @@ -194,16 +214,16 @@ def _lin_eq2dict(a, symset):
                 terms = ti
                 terms_coeff = ci
             else:
-                raise PolyNonlinearError
-        coeff = Mul(*coeff_list)
+                # since ti is not null and we already have
+                # a term, this is a cross term
+                raise PolyNonlinearError(filldedent('''
+                    nonlinear cross-term: %s''' % a))
+        coeff = Mul._from_args(coeff_list)
         if terms is None:
             return coeff, {}
         else:
             terms = {sym: coeff * c for sym, c in terms.items()}
             return  coeff * terms_coeff, terms
-    elif a.is_Equality:
-        return _lin_eq2dict(a.lhs - a.rhs, symset)
-    elif not a.has_free(*symset):
+    if not a.has_xfree(symset):
path: sympy/polys/matrices/linsolve.py
text: @author:
```suggestion
    elif not a.has_xfree(symset):
```

diff_hunk: @@ -194,16 +214,16 @@ def _lin_eq2dict(a, symset):
                 terms = ti
                 terms_coeff = ci
             else:
-                raise PolyNonlinearError
-        coeff = Mul(*coeff_list)
+                # since ti is not null and we already have
+                # a term, this is a cross term
+                raise PolyNonlinearError(filldedent('''
+                    nonlinear cross-term: %s''' % a))
+        coeff = Mul._from_args(coeff_list)
         if terms is None:
             return coeff, {}
         else:
             terms = {sym: coeff * c for sym, c in terms.items()}
             return  coeff * terms_coeff, terms
-    elif a.is_Equality:
-        return _lin_eq2dict(a.lhs - a.rhs, symset)
-    elif not a.has_free(*symset):
+    if not a.has_xfree(symset):
         return a, {}
-    else:
-        raise PolyNonlinearError
+    raise PolyNonlinearError('nonlinear term: %s' % a)
path: sympy/polys/matrices/linsolve.py
text: @author:
```suggestion
    else:
        raise PolyNonlinearError('nonlinear term: %s' % a)
```

diff_hunk: @@ -1624,12 +1625,31 @@ def test_has_free():
     assert Integral(f(x), (f(x), 1, y)).has_free(y)
     assert not Integral(f(x), (f(x), 1, y)).has_free(x)
     assert not Integral(f(x), (f(x), 1, y)).has_free(f(x))
+    # simple extraction
+    assert (x + 1 + y).has_free(x + 1)
+    assert not (x + 2 + y).has_free(x + 1)
+    assert (2 + 3*x*y).has_free(3*x)
+    raises(ValueError, lambda: x.has_free({x, y}))
+    s = FiniteSet(1, 2)
+    assert Piecewise((s, x > 3), (4, True)).has_free(s)
+    assert not Piecewise((1, x > 3), (4, True)).has_free(s)
+    # can't make set of these, but fallback will handle
+    assert not x.has_free(y, [])
path: sympy/core/tests/test_expr.py
text: @user1:
Shouldn't this just give a TypeError?

@author:
I'm not sure that it is yet a valid assumption that all Basic args are hashable. If it is, then an error can be raised.

@user1:
Basic args should be Basic and therefore should be hashable. If they aren't then that's a bug and it should be fixed. We shouldn't try to work around bugs like that because it just makes them less visible. There shouldn't be any code that attempts to allow for Basic to have non-Basic args.

diff_hunk: @@ -1258,6 +1258,33 @@ def has(self, *patterns):
         """
         return self._has(iterargs, *patterns)
 
+    def has_xfree(self, s: set[Basic]):
+        """return True if self has any of the patterns in s as a
+        free argument, else False. This is like `Basic.has_free`
+        but this will only report exact argument matches.
+
+        Examples
+        ========
+
+        >>> from sympy import Function
+        >>> from sympy.abc import x, y
+        >>> f = Function('f')
+        >>> f(x).has_xfree({f})
+        False
+        >>> f(x).has_xfree({f(x)})
+        True
+        >>> f(x + 1).has_xfree({x})
+        True
+        >>> f(x + 1).has_xfree({x + 1})
+        True
+        >>> f(x + y + 1).has_xfree({x + 1})
+        False
+        """
+        # protect O(1) containment check by requiring:
+        if not type(s) in (dict, set):
+            raise ValueError('expecting set or dict argument')
path: sympy/core/basic.py
text: @user1:
I guess TypeError is more appropriate. The type hint in the signature suggests that only `set` is allowed. I suggest only allowing `set` here. I don't see a strong benefit in anything else.

@author:
I put the dict there because you indicated that dict or set would have O(1) lookup. And there is no reason to force the user to convert dict to set if it can be used directly. I can make up a case where someone would want to pass a dict, but I'm not sure *I* would do so in practice. But I have appreciated the flexibility of Python to just "do the right thing" if no rules are broken...and to enforce as few rules as necessary.

diff_hunk: @@ -1258,6 +1258,33 @@ def has(self, *patterns):
         """
         return self._has(iterargs, *patterns)
 
+    def has_xfree(self, s: set[Basic]):
path: sympy/core/basic.py
text: @user1:
If we're type-hinting then we can also add the return type `-> bool`.

diff_hunk: @@ -1284,8 +1311,27 @@ def has_free(self, *patterns):
         True
         >>> (x + y + 1).has_free(y + 1)
         True
-
         """
+        if not patterns:
+            return False
+        p0 = patterns[0]
+        if len(patterns) == 1 and iterable(p0) and not isinstance(p0, Basic):
+            # Basic can contain iterables (though not non-Basic, ideally)
+            # but don't encourage mixed passing patterns
+            raise ValueError(filldedent('''
+                Expecting 1 or more Basic args, not a single
+                non-Basic iterable. Don't forget to unpack
+                iterables: `eq.has_free(*patterns)`'''))
+        # try quick test first
+        try:
+            s = set(patterns)
+        except TypeError:
+            pass  # patterns had a non-hashable element
path: sympy/core/basic.py
text: @user1:
As I've said before I would prefer not to catch `TypeError`. There are too many places that do this. A `TypeError` should almost always mean a bug and so the error should not generally be caught.

@author:
Please make a suggestion since there is no other error to catch and the scope has been limited to a single obvious case of why it might be raised. 

The reason this is here is because if we can make a set of the items and this test passes then it is a quick exit. And if not (either unhashable items or pattern not found) then we go to the next step where unhashable things are thrown out. In case of unhashable I don't want to let this raise for the user because the next step will ignore them (or could its own TypeError).

@user1:
Why we can't we just let this TypeError bubble up?

@author:
I guess we can see what fails...will make the change

diff_hunk: @@ -1284,8 +1311,27 @@ def has_free(self, *patterns):
         True
         >>> (x + y + 1).has_free(y + 1)
         True
-
         """
+        if not patterns:
+            return False
+        p0 = patterns[0]
+        if len(patterns) == 1 and iterable(p0) and not isinstance(p0, Basic):
+            # Basic can contain iterables (though not non-Basic, ideally)
+            # but don't encourage mixed passing patterns
+            raise ValueError(filldedent('''
path: sympy/core/basic.py
text: @user1:
`TypeError` would make more sense here.

diff_hunk: @@ -1258,6 +1258,33 @@ def has(self, *patterns):
         """
         return self._has(iterargs, *patterns)
 
+    def has_xfree(self, s: set[Basic]):
path: sympy/core/basic.py
text: @author:
```suggestion
    def has_xfree(self, s: set[Basic]) -> bool:
```

diff_hunk: @@ -1284,8 +1311,26 @@ def has_free(self, *patterns):
         True
         >>> (x + y + 1).has_free(y + 1)
         True
-
         """
+        if not patterns:
+            return False
+        p0 = patterns[0]
+        if len(patterns) == 1 and iterable(p0) and not isinstance(p0, Basic):
+            # Basic can contain iterables (though not non-Basic, ideally)
+            # but don't encourage mixed passing patterns
+            raise ValueError(filldedent('''
+                Expecting 1 or more Basic args, not a single
+                non-Basic iterable. Don't forget to unpack
+                iterables: `eq.has_free(*patterns)`'''))
+        # try quick test first
+        try:
+            s = set(patterns)
+            rv = self.has_xfree(s)
+            if rv:
+                return rv
+        except TypeError:
path: sympy/core/basic.py
text: @user1:
It's usually a bad idea to catch `TypeError`. Why is it being caught here?

@author:
In case of an un unhashable pattern which is not handled until we get to `_has`.

@user1:
What is an example of an unhashable pattern? Would they not be hashable if sympified?

@user1:
I would still prefer not to catch `TypeError`. It is not a good class of exceptions to catch. At least if we are catching it then the `try` block needs to be restricted in scope to a single operation that is the one where `TypeError` is expected and there should be comments explaining precisely where the `TypeError` can be expected to come from.

diff_hunk: @@ -1284,8 +1311,27 @@ def has_free(self, *patterns):
         True
         >>> (x + y + 1).has_free(y + 1)
         True
-
         """
+        if not patterns:
+            return False
+        p0 = patterns[0]
+        if len(patterns) == 1 and iterable(p0) and not isinstance(p0, Basic):
+            # Basic can contain iterables (though not non-Basic, ideally)
+            # but don't encourage mixed passing patterns
+            raise ValueError(filldedent('''
path: sympy/core/basic.py
text: @author:
```suggestion
            raise TypeError(filldedent('''
```

diff_hunk: @@ -1624,12 +1625,31 @@ def test_has_free():
     assert Integral(f(x), (f(x), 1, y)).has_free(y)
     assert not Integral(f(x), (f(x), 1, y)).has_free(x)
     assert not Integral(f(x), (f(x), 1, y)).has_free(f(x))
+    # simple extraction
+    assert (x + 1 + y).has_free(x + 1)
+    assert not (x + 2 + y).has_free(x + 1)
+    assert (2 + 3*x*y).has_free(3*x)
+    raises(ValueError, lambda: x.has_free({x, y}))
path: sympy/core/tests/test_expr.py
text: @author:
```suggestion
    raises(TypeError, lambda: x.has_free({x, y}))
```

diff_hunk: @@ -1624,12 +1625,31 @@ def test_has_free():
     assert Integral(f(x), (f(x), 1, y)).has_free(y)
     assert not Integral(f(x), (f(x), 1, y)).has_free(x)
     assert not Integral(f(x), (f(x), 1, y)).has_free(f(x))
+    # simple extraction
+    assert (x + 1 + y).has_free(x + 1)
+    assert not (x + 2 + y).has_free(x + 1)
+    assert (2 + 3*x*y).has_free(3*x)
+    raises(TypeError, lambda: x.has_free({x, y}))
+    s = FiniteSet(1, 2)
+    assert Piecewise((s, x > 3), (4, True)).has_free(s)
+    assert not Piecewise((1, x > 3), (4, True)).has_free(s)
+    # can't make set of these, but fallback will handle
+    assert not x.has_free(y, [])
path: sympy/core/tests/test_expr.py
text: @author:
```suggestion
    raises(TypeError, lambda: x.has_free(y, []))
```

diff_hunk: @@ -1284,8 +1311,27 @@ def has_free(self, *patterns):
         True
         >>> (x + y + 1).has_free(y + 1)
         True
-
         """
+        if not patterns:
+            return False
+        p0 = patterns[0]
+        if len(patterns) == 1 and iterable(p0) and not isinstance(p0, Basic):
+            # Basic can contain iterables (though not non-Basic, ideally)
+            # but don't encourage mixed passing patterns
+            raise TypeError(filldedent('''
+                Expecting 1 or more Basic args, not a single
+                non-Basic iterable. Don't forget to unpack
+                iterables: `eq.has_free(*patterns)`'''))
+        # try quick test first
+        try:
+            s = set(patterns)
+        except TypeError:
+            pass  # patterns had a non-hashable element
path: sympy/core/basic.py
text: @author:
```suggestion
        s = set(patterns)
```

diff_hunk: @@ -1258,6 +1258,33 @@ def has(self, *patterns):
         """
         return self._has(iterargs, *patterns)
 
+    def has_xfree(self, s: set[Basic]):
+        """return True if self has any of the patterns in s as a
+        free argument, else False. This is like `Basic.has_free`
+        but this will only report exact argument matches.
+
+        Examples
+        ========
+
+        >>> from sympy import Function
+        >>> from sympy.abc import x, y
+        >>> f = Function('f')
+        >>> f(x).has_xfree({f})
+        False
+        >>> f(x).has_xfree({f(x)})
+        True
+        >>> f(x + 1).has_xfree({x})
+        True
+        >>> f(x + 1).has_xfree({x + 1})
+        True
+        >>> f(x + y + 1).has_xfree({x + 1})
+        False
+        """
+        # protect O(1) containment check by requiring:
+        if not type(s) in (dict, set):
+            raise ValueError('expecting set or dict argument')
path: sympy/core/basic.py
text: @author:
```suggestion
        if type(s) is not set:
            raise TypeError('expecting set argument')
```

diff_hunk: @@ -2785,8 +2797,11 @@ def test_linear_coeffs():
         linear_coeffs(x, x, x))
     assert linear_coeffs(a*(x + y), x, y) == [a, a, 0]
     assert linear_coeffs(1.0, x, y) == [0, 0, 1.0]
+    # don't include coefficients of 0
+    assert linear_coeffs(Eq(x, x + y), x, y, dict=True) == {y: 1}
path: sympy/solvers/tests/test_solveset.py
text: @author:
```suggestion
    assert linear_coeffs(Eq(x, x + y), x, y, dict=True) == {y: -1}
```

diff_hunk: @@ -1624,12 +1625,31 @@ def test_has_free():
     assert Integral(f(x), (f(x), 1, y)).has_free(y)
     assert not Integral(f(x), (f(x), 1, y)).has_free(x)
     assert not Integral(f(x), (f(x), 1, y)).has_free(f(x))
+    # simple extraction
+    assert (x + 1 + y).has_free(x + 1)
+    assert not (x + 2 + y).has_free(x + 1)
+    assert (2 + 3*x*y).has_free(3*x)
+    raises(TypeError, lambda: x.has_free({x, y}))
+    s = FiniteSet(1, 2)
+    assert Piecewise((s, x > 3), (4, True)).has_free(s)
+    assert not Piecewise((1, x > 3), (4, True)).has_free(s)
+    # can't make set of these, but fallback will handle
+    raises(TypeError, lambda: x.has_free(y, []))
+
+
+def test_has_xfree():
+    assert (x + 1).has_xfree({x})
+    assert ((x + 1)**2).has_xfree({x + 1})
+    assert not (x + y + 1).has_xfree({x + 1})
+    raises(ValueError, lambda: x.has_xfree(x))
+    raises(ValueError, lambda: x.has_xfree([x]))
path: sympy/core/tests/test_expr.py
text: @author:
```suggestion
    raises(TypeError, lambda: x.has_xfree(x))
    raises(TypeError, lambda: x.has_xfree([x]))
```

diff_hunk: @@ -2760,16 +2771,23 @@ def linsolve(system, *symbols):
     >>> linsolve([], x)
     EmptySet
 
-    * An error is raised if, after expansion, any nonlinearity
-      is detected:
+    * An error is raised if any nonlinearity is detected, even
+    if it could be removed with expansion
path: sympy/solvers/solveset.py
text: @author:
```suggestion
    * An error is raised if any nonlinearity is detected, even
      if it could be removed with expansion
```

diff_hunk: @@ -1284,8 +1311,24 @@ def has_free(self, *patterns):
         True
         >>> (x + y + 1).has_free(y + 1)
         True
-
         """
+        if not patterns:
+            return False
+        p0 = patterns[0]
+        if len(patterns) == 1 and iterable(p0) and not isinstance(p0, Basic):
+            # Basic can contain iterables (though not non-Basic, ideally)
+            # but don't encourage mixed passing patterns
+            raise TypeError(filldedent('''
+                Expecting 1 or more Basic args, not a single
+                non-Basic iterable. Don't forget to unpack
+                iterables: `eq.has_free(*patterns)`'''))
+        # try quick test first
+        s = set(patterns)
+        else:
+            rv = self.has_xfree(s)
+            if rv:
+                return rv
path: sympy/core/basic.py
text: @author:
```suggestion
        rv = self.has_xfree(s)
        if rv:
            return rv
```

diff_hunk: @@ -2760,16 +2771,23 @@ def linsolve(system, *symbols):
     >>> linsolve([], x)
     EmptySet
 
-    * An error is raised if, after expansion, any nonlinearity
-      is detected:
+    * An error is raised if any nonlinearity is detected, even
+    if it could be removed with expansion
path: sympy/solvers/solveset.py
text: @user1:
I'm guessing that the indentation here is the sphinx problem:
```
/home/runner/work/sympy/sympy/sympy/solvers/solveset.py:docstring of sympy.solvers.solveset.linsolve:169: WARNING: Bullet list ends without a blank line; unexpected unindent.
```

@author:
I used the editor to search for leading space followed by `*`, looking for this case. Not sure how I missed it! Thanks.