forked from diofant/diofant
/
polysys.py
301 lines (217 loc) · 8.31 KB
/
polysys.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
"""Solvers of systems of polynomial equations. """
import collections
from diofant.core import S
from diofant.matrices import Matrix
from diofant.polys import Poly, groebner, sring
from diofant.polys.polytools import parallel_poly_from_expr
from diofant.polys.polyerrors import (ComputationFailed,
PolificationFailed, CoercionFailed)
from diofant.polys.solvers import solve_lin_sys
from diofant.simplify import rcollect, simplify
from diofant.utilities import default_sort_key, postfixes
__all__ = ('solve_linear_system', 'solve_poly_system')
class SolveFailed(Exception):
"""Raised when solver's conditions weren't met. """
def roots(p):
r = collections.defaultdict(int)
for v in p.all_roots():
r[v] += 1
return r
def solve_linear_system(system, *symbols, **flags):
r"""Solve system of linear equations.
Both under- and overdetermined systems are supported. The possible
number of solutions is zero, one or infinite.
Parameters
==========
system : Matrix
Nx(M+1) matrix, which means it has to be in augmented
form. This matrix will not be modified.
\*symbols : list
List of M Symbol's
Returns
=======
solution: dict or None
Respectively, this procedure will return None or
a dictionary with solutions. In the case of underdetermined
systems, all arbitrary parameters are skipped. This may
cause a situation in which an empty dictionary is returned.
In that case, all symbols can be assigned arbitrary values.
Examples
========
>>> from diofant.abc import x, y
Solve the following system::
x + 4 y == 2
-2 x + y == 14
>>> system = Matrix(( (1, 4, 2), (-2, 1, 14)))
>>> solve_linear_system(system, x, y)
{x: -6, y: 2}
A degenerate system returns an empty dictionary.
>>> system = Matrix(( (0,0,0), (0,0,0) ))
>>> solve_linear_system(system, x, y)
{}
See Also
========
diofant.matrices.matrices.MatrixBase.rref
"""
eqs = system*Matrix(symbols + (-1,))
domain, eqs = sring(eqs.transpose().tolist()[0], *symbols, field=True)
res = solve_lin_sys(eqs, domain)
if res is None:
return
for k in list(res.keys()):
s = domain.symbols[domain.index(k)]
res[s] = res[k].as_expr()
del res[k]
if flags.get('simplify', True):
res[s] = simplify(res[s])
return res
def solve_poly_system(seq, *gens, **args):
"""
Solve a system of polynomial equations.
Examples
========
>>> from diofant.abc import x, y
>>> solve_poly_system([x*y - 2*y, 2*y**2 - x**2], x, y)
[{x: 0, y: 0}, {x: 2, y: -sqrt(2)}, {x: 2, y: sqrt(2)}]
"""
try:
polys, opt = parallel_poly_from_expr(seq, *gens, **args)
except PolificationFailed as exc:
raise ComputationFailed('solve_poly_system', len(seq), exc)
if len(polys) == len(opt.gens) == 2:
f, g = polys
a, b = f.degree_list()
c, d = g.degree_list()
if a <= 2 and b <= 2 and c <= 2 and d <= 2:
try:
return solve_biquadratic(f, g, opt)
except SolveFailed:
pass
return solve_generic(polys, opt)
def solve_biquadratic(f, g, opt):
"""
Solve a system of two bivariate quadratic polynomial equations.
Examples
========
>>> from diofant.polys import Options, Poly
>>> from diofant.abc import x, y
>>> NewOption = Options((x, y), {'domain': 'ZZ'})
>>> a = Poly(y**2 - 4 + x, y, x, domain='ZZ')
>>> b = Poly(y*2 + 3*x - 7, y, x, domain='ZZ')
>>> solve_biquadratic(a, b, NewOption)
[{x: 1/3, y: 3}, {x: 41/27, y: 11/9}]
>>> a = Poly(y + x**2 - 3, y, x, domain='ZZ')
>>> b = Poly(-y + x - 4, y, x, domain='ZZ')
>>> solve_biquadratic(a, b, NewOption)
[{x: -sqrt(29)/2 + 7/2, y: -sqrt(29)/2 - 1/2},
{x: sqrt(29)/2 + 7/2, y: -1/2 + sqrt(29)/2}]
"""
G = groebner([f, g])
if len(G) == 1 and G[0].is_ground:
return
if len(G) != 2:
raise SolveFailed
p, q = G
x, y = opt.gens
p = Poly(p, x, expand=False)
q = q.ltrim(-1)
p_roots = [rcollect(expr, y) for expr in roots(p).keys()]
q_roots = list(roots(q).keys())
solutions = []
for q_root in q_roots:
for p_root in p_roots:
solution = {x: p_root.subs(y, q_root), y: q_root}
solutions.append(solution)
return sorted(solutions, key=default_sort_key)
def solve_generic(polys, opt):
"""
Solve a generic system of polynomial equations.
Returns all possible solutions over C[x_1, x_2, ..., x_m] of a
set F = { f_1, f_2, ..., f_n } of polynomial equations, using
Groebner basis approach. For now only zero-dimensional systems
are supported, which means F can have at most a finite number
of solutions.
The algorithm works by the fact that, supposing G is the basis
of F with respect to an elimination order (here lexicographic
order is used), G and F generate the same ideal, they have the
same set of solutions. By the elimination property, if G is a
reduced, zero-dimensional Groebner basis, then there exists an
univariate polynomial in G (in its last variable). This can be
solved by computing its roots. Substituting all computed roots
for the last (eliminated) variable in other elements of G, new
polynomial system is generated. Applying the above procedure
recursively, a finite number of solutions can be found.
References
==========
.. [Buchberger01] B. Buchberger, Groebner Bases: A Short
Introduction for Systems Theorists, In: R. Moreno-Diaz,
B. Buchberger, J.L. Freire, Proceedings of EUROCAST'01,
February, 2001
.. [Cox97] D. Cox, J. Little, D. O'Shea, Ideals, Varieties
and Algorithms, Springer, Second Edition, 1997, pp. 112
Examples
========
>>> from diofant.polys import Poly, Options
>>> from diofant.abc import x, y
>>> NewOption = Options((x, y), {'domain': 'ZZ'})
>>> a = Poly(x - y + 5, x, y, domain='ZZ')
>>> b = Poly(x + y - 3, x, y, domain='ZZ')
>>> solve_generic([a, b], NewOption)
[{x: -1, y: 4}]
>>> a = Poly(x - 2*y + 5, x, y, domain='ZZ')
>>> b = Poly(2*x - y - 3, x, y, domain='ZZ')
>>> solve_generic([a, b], NewOption)
[{x: 11/3, y: 13/3}]
>>> a = Poly(x**2 + y, x, y, domain='ZZ')
>>> b = Poly(x + y*4, x, y, domain='ZZ')
>>> solve_generic([a, b], NewOption)
[{x: 0, y: 0}, {x: 1/4, y: -1/16}]
"""
def _is_univariate(f):
"""Returns True if 'f' is univariate in its last variable. """
for monom in f.monoms():
if any(m > 0 for m in monom[:-1]):
return False
return True
def _subs_root(f, gen, zero):
"""Replace generator with a root so that the result is nice. """
p = f.as_expr({gen: zero})
if f.degree(gen) >= 2:
p = p.expand(deep=False)
return p
def _solve_reduced_system(system, gens):
"""Recursively solves reduced polynomial systems. """
basis = groebner(system, gens, polys=True)
if len(basis) == 1 and basis[0].is_ground:
return
univariate = list(filter(_is_univariate, basis))
if len(univariate) == 1:
f = univariate.pop()
else:
raise NotImplementedError("only zero-dimensional systems "
"supported (finite number of solutions)")
gens = f.gens
gen = gens[-1]
zeros = [k.doit() for k in roots(f.ltrim(gen)).keys()]
if len(basis) == 1:
return [{gen: zero} for zero in zeros]
solutions = []
for zero in zeros:
new_system = []
new_gens = gens[:-1]
for b in basis[:-1]:
eq = _subs_root(b, gen, zero)
if eq is not S.Zero:
new_system.append(eq)
for solution in _solve_reduced_system(new_system, new_gens):
solution[gen] = zero
solutions.append(solution)
return solutions
try:
result = _solve_reduced_system(polys, opt.gens)
except CoercionFailed: # pragma: no cover
raise NotImplementedError
if result is not None:
return sorted(result, key=default_sort_key)
else:
return