Coverage for /usr/lib/python3/dist-packages/sympy/polys/solvers.py: 14%
145 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Low-level linear systems solver. """
4from sympy.utilities.exceptions import sympy_deprecation_warning
5from sympy.utilities.iterables import connected_components
7from sympy.core.sympify import sympify
8from sympy.core.numbers import Integer, Rational
9from sympy.matrices.dense import MutableDenseMatrix
10from sympy.polys.domains import ZZ, QQ
12from sympy.polys.domains import EX
13from sympy.polys.rings import sring
14from sympy.polys.polyerrors import NotInvertible
15from sympy.polys.domainmatrix import DomainMatrix
18class PolyNonlinearError(Exception):
19 """Raised by solve_lin_sys for nonlinear equations"""
20 pass
23class RawMatrix(MutableDenseMatrix):
24 """
25 .. deprecated:: 1.9
27 This class fundamentally is broken by design. Use ``DomainMatrix`` if
28 you want a matrix over the polys domains or ``Matrix`` for a matrix
29 with ``Expr`` elements. The ``RawMatrix`` class will be removed/broken
30 in future in order to reestablish the invariant that the elements of a
31 Matrix should be of type ``Expr``.
33 """
34 _sympify = staticmethod(lambda x: x) # type: ignore
36 def __init__(self, *args, **kwargs):
37 sympy_deprecation_warning(
38 """
39 The RawMatrix class is deprecated. Use either DomainMatrix or
40 Matrix instead.
41 """,
42 deprecated_since_version="1.9",
43 active_deprecations_target="deprecated-rawmatrix",
44 )
46 domain = ZZ
47 for i in range(self.rows):
48 for j in range(self.cols):
49 val = self[i,j]
50 if getattr(val, 'is_Poly', False):
51 K = val.domain[val.gens]
52 val_sympy = val.as_expr()
53 elif hasattr(val, 'parent'):
54 K = val.parent()
55 val_sympy = K.to_sympy(val)
56 elif isinstance(val, (int, Integer)):
57 K = ZZ
58 val_sympy = sympify(val)
59 elif isinstance(val, Rational):
60 K = QQ
61 val_sympy = val
62 else:
63 for K in ZZ, QQ:
64 if K.of_type(val):
65 val_sympy = K.to_sympy(val)
66 break
67 else:
68 raise TypeError
69 domain = domain.unify(K)
70 self[i,j] = val_sympy
71 self.ring = domain
74def eqs_to_matrix(eqs_coeffs, eqs_rhs, gens, domain):
75 """Get matrix from linear equations in dict format.
77 Explanation
78 ===========
80 Get the matrix representation of a system of linear equations represented
81 as dicts with low-level DomainElement coefficients. This is an
82 *internal* function that is used by solve_lin_sys.
84 Parameters
85 ==========
87 eqs_coeffs: list[dict[Symbol, DomainElement]]
88 The left hand sides of the equations as dicts mapping from symbols to
89 coefficients where the coefficients are instances of
90 DomainElement.
91 eqs_rhs: list[DomainElements]
92 The right hand sides of the equations as instances of
93 DomainElement.
94 gens: list[Symbol]
95 The unknowns in the system of equations.
96 domain: Domain
97 The domain for coefficients of both lhs and rhs.
99 Returns
100 =======
102 The augmented matrix representation of the system as a DomainMatrix.
104 Examples
105 ========
107 >>> from sympy import symbols, ZZ
108 >>> from sympy.polys.solvers import eqs_to_matrix
109 >>> x, y = symbols('x, y')
110 >>> eqs_coeff = [{x:ZZ(1), y:ZZ(1)}, {x:ZZ(1), y:ZZ(-1)}]
111 >>> eqs_rhs = [ZZ(0), ZZ(-1)]
112 >>> eqs_to_matrix(eqs_coeff, eqs_rhs, [x, y], ZZ)
113 DomainMatrix([[1, 1, 0], [1, -1, 1]], (2, 3), ZZ)
115 See also
116 ========
118 solve_lin_sys: Uses :func:`~eqs_to_matrix` internally
119 """
120 sym2index = {x: n for n, x in enumerate(gens)}
121 nrows = len(eqs_coeffs)
122 ncols = len(gens) + 1
123 rows = [[domain.zero] * ncols for _ in range(nrows)]
124 for row, eq_coeff, eq_rhs in zip(rows, eqs_coeffs, eqs_rhs):
125 for sym, coeff in eq_coeff.items():
126 row[sym2index[sym]] = domain.convert(coeff)
127 row[-1] = -domain.convert(eq_rhs)
129 return DomainMatrix(rows, (nrows, ncols), domain)
132def sympy_eqs_to_ring(eqs, symbols):
133 """Convert a system of equations from Expr to a PolyRing
135 Explanation
136 ===========
138 High-level functions like ``solve`` expect Expr as inputs but can use
139 ``solve_lin_sys`` internally. This function converts equations from
140 ``Expr`` to the low-level poly types used by the ``solve_lin_sys``
141 function.
143 Parameters
144 ==========
146 eqs: List of Expr
147 A list of equations as Expr instances
148 symbols: List of Symbol
149 A list of the symbols that are the unknowns in the system of
150 equations.
152 Returns
153 =======
155 Tuple[List[PolyElement], Ring]: The equations as PolyElement instances
156 and the ring of polynomials within which each equation is represented.
158 Examples
159 ========
161 >>> from sympy import symbols
162 >>> from sympy.polys.solvers import sympy_eqs_to_ring
163 >>> a, x, y = symbols('a, x, y')
164 >>> eqs = [x-y, x+a*y]
165 >>> eqs_ring, ring = sympy_eqs_to_ring(eqs, [x, y])
166 >>> eqs_ring
167 [x - y, x + a*y]
168 >>> type(eqs_ring[0])
169 <class 'sympy.polys.rings.PolyElement'>
170 >>> ring
171 ZZ(a)[x,y]
173 With the equations in this form they can be passed to ``solve_lin_sys``:
175 >>> from sympy.polys.solvers import solve_lin_sys
176 >>> solve_lin_sys(eqs_ring, ring)
177 {y: 0, x: 0}
178 """
179 try:
180 K, eqs_K = sring(eqs, symbols, field=True, extension=True)
181 except NotInvertible:
182 # https://github.com/sympy/sympy/issues/18874
183 K, eqs_K = sring(eqs, symbols, domain=EX)
184 return eqs_K, K.to_domain()
187def solve_lin_sys(eqs, ring, _raw=True):
188 """Solve a system of linear equations from a PolynomialRing
190 Explanation
191 ===========
193 Solves a system of linear equations given as PolyElement instances of a
194 PolynomialRing. The basic arithmetic is carried out using instance of
195 DomainElement which is more efficient than :class:`~sympy.core.expr.Expr`
196 for the most common inputs.
198 While this is a public function it is intended primarily for internal use
199 so its interface is not necessarily convenient. Users are suggested to use
200 the :func:`sympy.solvers.solveset.linsolve` function (which uses this
201 function internally) instead.
203 Parameters
204 ==========
206 eqs: list[PolyElement]
207 The linear equations to be solved as elements of a
208 PolynomialRing (assumed equal to zero).
209 ring: PolynomialRing
210 The polynomial ring from which eqs are drawn. The generators of this
211 ring are the unknowns to be solved for and the domain of the ring is
212 the domain of the coefficients of the system of equations.
213 _raw: bool
214 If *_raw* is False, the keys and values in the returned dictionary
215 will be of type Expr (and the unit of the field will be removed from
216 the keys) otherwise the low-level polys types will be returned, e.g.
217 PolyElement: PythonRational.
219 Returns
220 =======
222 ``None`` if the system has no solution.
224 dict[Symbol, Expr] if _raw=False
226 dict[Symbol, DomainElement] if _raw=True.
228 Examples
229 ========
231 >>> from sympy import symbols
232 >>> from sympy.polys.solvers import solve_lin_sys, sympy_eqs_to_ring
233 >>> x, y = symbols('x, y')
234 >>> eqs = [x - y, x + y - 2]
235 >>> eqs_ring, ring = sympy_eqs_to_ring(eqs, [x, y])
236 >>> solve_lin_sys(eqs_ring, ring)
237 {y: 1, x: 1}
239 Passing ``_raw=False`` returns the same result except that the keys are
240 ``Expr`` rather than low-level poly types.
242 >>> solve_lin_sys(eqs_ring, ring, _raw=False)
243 {x: 1, y: 1}
245 See also
246 ========
248 sympy_eqs_to_ring: prepares the inputs to ``solve_lin_sys``.
249 linsolve: ``linsolve`` uses ``solve_lin_sys`` internally.
250 sympy.solvers.solvers.solve: ``solve`` uses ``solve_lin_sys`` internally.
251 """
252 as_expr = not _raw
254 assert ring.domain.is_Field
256 eqs_dict = [dict(eq) for eq in eqs]
258 one_monom = ring.one.monoms()[0]
259 zero = ring.domain.zero
261 eqs_rhs = []
262 eqs_coeffs = []
263 for eq_dict in eqs_dict:
264 eq_rhs = eq_dict.pop(one_monom, zero)
265 eq_coeffs = {}
266 for monom, coeff in eq_dict.items():
267 if sum(monom) != 1:
268 msg = "Nonlinear term encountered in solve_lin_sys"
269 raise PolyNonlinearError(msg)
270 eq_coeffs[ring.gens[monom.index(1)]] = coeff
271 if not eq_coeffs:
272 if not eq_rhs:
273 continue
274 else:
275 return None
276 eqs_rhs.append(eq_rhs)
277 eqs_coeffs.append(eq_coeffs)
279 result = _solve_lin_sys(eqs_coeffs, eqs_rhs, ring)
281 if result is not None and as_expr:
283 def to_sympy(x):
284 as_expr = getattr(x, 'as_expr', None)
285 if as_expr:
286 return as_expr()
287 else:
288 return ring.domain.to_sympy(x)
290 tresult = {to_sympy(sym): to_sympy(val) for sym, val in result.items()}
292 # Remove 1.0x
293 result = {}
294 for k, v in tresult.items():
295 if k.is_Mul:
296 c, s = k.as_coeff_Mul()
297 result[s] = v/c
298 else:
299 result[k] = v
301 return result
304def _solve_lin_sys(eqs_coeffs, eqs_rhs, ring):
305 """Solve a linear system from dict of PolynomialRing coefficients
307 Explanation
308 ===========
310 This is an **internal** function used by :func:`solve_lin_sys` after the
311 equations have been preprocessed. The role of this function is to split
312 the system into connected components and pass those to
313 :func:`_solve_lin_sys_component`.
315 Examples
316 ========
318 Setup a system for $x-y=0$ and $x+y=2$ and solve:
320 >>> from sympy import symbols, sring
321 >>> from sympy.polys.solvers import _solve_lin_sys
322 >>> x, y = symbols('x, y')
323 >>> R, (xr, yr) = sring([x, y], [x, y])
324 >>> eqs = [{xr:R.one, yr:-R.one}, {xr:R.one, yr:R.one}]
325 >>> eqs_rhs = [R.zero, -2*R.one]
326 >>> _solve_lin_sys(eqs, eqs_rhs, R)
327 {y: 1, x: 1}
329 See also
330 ========
332 solve_lin_sys: This function is used internally by :func:`solve_lin_sys`.
333 """
334 V = ring.gens
335 E = []
336 for eq_coeffs in eqs_coeffs:
337 syms = list(eq_coeffs)
338 E.extend(zip(syms[:-1], syms[1:]))
339 G = V, E
341 components = connected_components(G)
343 sym2comp = {}
344 for n, component in enumerate(components):
345 for sym in component:
346 sym2comp[sym] = n
348 subsystems = [([], []) for _ in range(len(components))]
349 for eq_coeff, eq_rhs in zip(eqs_coeffs, eqs_rhs):
350 sym = next(iter(eq_coeff), None)
351 sub_coeff, sub_rhs = subsystems[sym2comp[sym]]
352 sub_coeff.append(eq_coeff)
353 sub_rhs.append(eq_rhs)
355 sol = {}
356 for subsystem in subsystems:
357 subsol = _solve_lin_sys_component(subsystem[0], subsystem[1], ring)
358 if subsol is None:
359 return None
360 sol.update(subsol)
362 return sol
365def _solve_lin_sys_component(eqs_coeffs, eqs_rhs, ring):
366 """Solve a linear system from dict of PolynomialRing coefficients
368 Explanation
369 ===========
371 This is an **internal** function used by :func:`solve_lin_sys` after the
372 equations have been preprocessed. After :func:`_solve_lin_sys` splits the
373 system into connected components this function is called for each
374 component. The system of equations is solved using Gauss-Jordan
375 elimination with division followed by back-substitution.
377 Examples
378 ========
380 Setup a system for $x-y=0$ and $x+y=2$ and solve:
382 >>> from sympy import symbols, sring
383 >>> from sympy.polys.solvers import _solve_lin_sys_component
384 >>> x, y = symbols('x, y')
385 >>> R, (xr, yr) = sring([x, y], [x, y])
386 >>> eqs = [{xr:R.one, yr:-R.one}, {xr:R.one, yr:R.one}]
387 >>> eqs_rhs = [R.zero, -2*R.one]
388 >>> _solve_lin_sys_component(eqs, eqs_rhs, R)
389 {y: 1, x: 1}
391 See also
392 ========
394 solve_lin_sys: This function is used internally by :func:`solve_lin_sys`.
395 """
397 # transform from equations to matrix form
398 matrix = eqs_to_matrix(eqs_coeffs, eqs_rhs, ring.gens, ring.domain)
400 # convert to a field for rref
401 if not matrix.domain.is_Field:
402 matrix = matrix.to_field()
404 # solve by row-reduction
405 echelon, pivots = matrix.rref()
407 # construct the returnable form of the solutions
408 keys = ring.gens
410 if pivots and pivots[-1] == len(keys):
411 return None
413 if len(pivots) == len(keys):
414 sol = []
415 for s in [row[-1] for row in echelon.rep.to_ddm()]:
416 a = s
417 sol.append(a)
418 sols = dict(zip(keys, sol))
419 else:
420 sols = {}
421 g = ring.gens
422 # Extract ground domain coefficients and convert to the ring:
423 if hasattr(ring, 'ring'):
424 convert = ring.ring.ground_new
425 else:
426 convert = ring.ground_new
427 echelon = echelon.rep.to_ddm()
428 vals_set = {v for row in echelon for v in row}
429 vals_map = {v: convert(v) for v in vals_set}
430 echelon = [[vals_map[eij] for eij in ei] for ei in echelon]
431 for i, p in enumerate(pivots):
432 v = echelon[i][-1] - sum(echelon[i][j]*g[j] for j in range(p+1, len(g)) if echelon[i][j])
433 sols[keys[p]] = v
435 return sols