Coverage for /usr/lib/python3/dist-packages/sympy/solvers/polysys.py: 11%

1"""Solvers of systems of polynomial equations. """ 

2import itertools 

3 

4from sympy.core import S 

5from sympy.core.sorting import default_sort_key 

6from sympy.polys import Poly, groebner, roots 

7from sympy.polys.polytools import parallel_poly_from_expr 

8from sympy.polys.polyerrors import (ComputationFailed, 

9 PolificationFailed, CoercionFailed) 

10from sympy.simplify import rcollect 

11from sympy.utilities import postfixes 

12from sympy.utilities.misc import filldedent 

13 

14 

15class SolveFailed(Exception): 

16 """Raised when solver's conditions were not met. """ 

17 

18 

19def solve_poly_system(seq, *gens, strict=False, **args): 

20 """ 

21 Return a list of solutions for the system of polynomial equations 

22 or else None. 

23 

24 Parameters 

25 ========== 

26 

27 seq: a list/tuple/set 

28 Listing all the equations that are needed to be solved 

29 gens: generators 

30 generators of the equations in seq for which we want the 

31 solutions 

32 strict: a boolean (default is False) 

33 if strict is True, NotImplementedError will be raised if 

34 the solution is known to be incomplete (which can occur if 

35 not all solutions are expressible in radicals) 

36 args: Keyword arguments 

37 Special options for solving the equations. 

38 

39 

40 Returns 

41 ======= 

42 

43 List[Tuple] 

44 a list of tuples with elements being solutions for the 

45 symbols in the order they were passed as gens 

46 None 

47 None is returned when the computed basis contains only the ground. 

48 

49 Examples 

50 ======== 

51 

52 >>> from sympy import solve_poly_system 

53 >>> from sympy.abc import x, y 

54 

55 >>> solve_poly_system([x*y - 2*y, 2*y**2 - x**2], x, y) 

56 [(0, 0), (2, -sqrt(2)), (2, sqrt(2))] 

57 

58 >>> solve_poly_system([x**5 - x + y**3, y**2 - 1], x, y, strict=True) 

59 Traceback (most recent call last): 

60 ... 

61 UnsolvableFactorError 

62 

63 """ 

64 try: 

65 polys, opt = parallel_poly_from_expr(seq, *gens, **args) 

66 except PolificationFailed as exc: 

67 raise ComputationFailed('solve_poly_system', len(seq), exc) 

68 

69 if len(polys) == len(opt.gens) == 2: 

70 f, g = polys 

71 

72 if all(i <= 2 for i in f.degree_list() + g.degree_list()): 

73 try: 

74 return solve_biquadratic(f, g, opt) 

75 except SolveFailed: 

76 pass 

77 

78 return solve_generic(polys, opt, strict=strict) 

79 

80 

81def solve_biquadratic(f, g, opt): 

82 """Solve a system of two bivariate quadratic polynomial equations. 

83 

84 Parameters 

85 ========== 

86 

87 f: a single Expr or Poly 

88 First equation 

89 g: a single Expr or Poly 

90 Second Equation 

91 opt: an Options object 

92 For specifying keyword arguments and generators 

93 

94 Returns 

95 ======= 

96 

97 List[Tuple] 

98 a list of tuples with elements being solutions for the 

99 symbols in the order they were passed as gens 

100 None 

101 None is returned when the computed basis contains only the ground. 

102 

103 Examples 

104 ======== 

105 

106 >>> from sympy import Options, Poly 

107 >>> from sympy.abc import x, y 

108 >>> from sympy.solvers.polysys import solve_biquadratic 

109 >>> NewOption = Options((x, y), {'domain': 'ZZ'}) 

110 

111 >>> a = Poly(y**2 - 4 + x, y, x, domain='ZZ') 

112 >>> b = Poly(y*2 + 3*x - 7, y, x, domain='ZZ') 

113 >>> solve_biquadratic(a, b, NewOption) 

114 [(1/3, 3), (41/27, 11/9)] 

115 

116 >>> a = Poly(y + x**2 - 3, y, x, domain='ZZ') 

117 >>> b = Poly(-y + x - 4, y, x, domain='ZZ') 

118 >>> solve_biquadratic(a, b, NewOption) 

119 [(7/2 - sqrt(29)/2, -sqrt(29)/2 - 1/2), (sqrt(29)/2 + 7/2, -1/2 + \ 

120 sqrt(29)/2)] 

121 """ 

122 G = groebner([f, g]) 

123 

124 if len(G) == 1 and G[0].is_ground: 

125 return None 

126 

127 if len(G) != 2: 

128 raise SolveFailed 

129 

130 x, y = opt.gens 

131 p, q = G 

132 if not p.gcd(q).is_ground: 

133 # not 0-dimensional 

134 raise SolveFailed 

135 

136 p = Poly(p, x, expand=False) 

137 p_roots = [rcollect(expr, y) for expr in roots(p).keys()] 

138 

139 q = q.ltrim(-1) 

140 q_roots = list(roots(q).keys()) 

141 

142 solutions = [(p_root.subs(y, q_root), q_root) for q_root, p_root in 

143 itertools.product(q_roots, p_roots)] 

144 

145 return sorted(solutions, key=default_sort_key) 

146 

147 

148def solve_generic(polys, opt, strict=False): 

149 """ 

150 Solve a generic system of polynomial equations. 

151 

152 Returns all possible solutions over C[x_1, x_2, ..., x_m] of a 

153 set F = { f_1, f_2, ..., f_n } of polynomial equations, using 

154 Groebner basis approach. For now only zero-dimensional systems 

155 are supported, which means F can have at most a finite number 

156 of solutions. If the basis contains only the ground, None is 

157 returned. 

158 

159 The algorithm works by the fact that, supposing G is the basis 

160 of F with respect to an elimination order (here lexicographic 

161 order is used), G and F generate the same ideal, they have the 

162 same set of solutions. By the elimination property, if G is a 

163 reduced, zero-dimensional Groebner basis, then there exists an 

164 univariate polynomial in G (in its last variable). This can be 

165 solved by computing its roots. Substituting all computed roots 

166 for the last (eliminated) variable in other elements of G, new 

167 polynomial system is generated. Applying the above procedure 

168 recursively, a finite number of solutions can be found. 

169 

170 The ability of finding all solutions by this procedure depends 

171 on the root finding algorithms. If no solutions were found, it 

172 means only that roots() failed, but the system is solvable. To 

173 overcome this difficulty use numerical algorithms instead. 

174 

175 Parameters 

176 ========== 

177 

178 polys: a list/tuple/set 

179 Listing all the polynomial equations that are needed to be solved 

180 opt: an Options object 

181 For specifying keyword arguments and generators 

182 strict: a boolean 

183 If strict is True, NotImplementedError will be raised if the solution 

184 is known to be incomplete 

185 

186 Returns 

187 ======= 

188 

189 List[Tuple] 

190 a list of tuples with elements being solutions for the 

191 symbols in the order they were passed as gens 

192 None 

193 None is returned when the computed basis contains only the ground. 

194 

195 References 

196 ========== 

197 

198 .. [Buchberger01] B. Buchberger, Groebner Bases: A Short 

199 Introduction for Systems Theorists, In: R. Moreno-Diaz, 

200 B. Buchberger, J.L. Freire, Proceedings of EUROCAST'01, 

201 February, 2001 

202 

203 .. [Cox97] D. Cox, J. Little, D. O'Shea, Ideals, Varieties 

204 and Algorithms, Springer, Second Edition, 1997, pp. 112 

205 

206 Raises 

207 ======== 

208 

209 NotImplementedError 

210 If the system is not zero-dimensional (does not have a finite 

211 number of solutions) 

212 

213 UnsolvableFactorError 

214 If ``strict`` is True and not all solution components are 

215 expressible in radicals 

216 

217 Examples 

218 ======== 

219 

220 >>> from sympy import Poly, Options 

221 >>> from sympy.solvers.polysys import solve_generic 

222 >>> from sympy.abc import x, y 

223 >>> NewOption = Options((x, y), {'domain': 'ZZ'}) 

224 

225 >>> a = Poly(x - y + 5, x, y, domain='ZZ') 

226 >>> b = Poly(x + y - 3, x, y, domain='ZZ') 

227 >>> solve_generic([a, b], NewOption) 

228 [(-1, 4)] 

229 

230 >>> a = Poly(x - 2*y + 5, x, y, domain='ZZ') 

231 >>> b = Poly(2*x - y - 3, x, y, domain='ZZ') 

232 >>> solve_generic([a, b], NewOption) 

233 [(11/3, 13/3)] 

234 

235 >>> a = Poly(x**2 + y, x, y, domain='ZZ') 

236 >>> b = Poly(x + y*4, x, y, domain='ZZ') 

237 >>> solve_generic([a, b], NewOption) 

238 [(0, 0), (1/4, -1/16)] 

239 

240 >>> a = Poly(x**5 - x + y**3, x, y, domain='ZZ') 

241 >>> b = Poly(y**2 - 1, x, y, domain='ZZ') 

242 >>> solve_generic([a, b], NewOption, strict=True) 

243 Traceback (most recent call last): 

244 ... 

245 UnsolvableFactorError 

246 

247 """ 

248 def _is_univariate(f): 

249 """Returns True if 'f' is univariate in its last variable. """ 

250 for monom in f.monoms(): 

251 if any(monom[:-1]): 

252 return False 

253 

254 return True 

255 

256 def _subs_root(f, gen, zero): 

257 """Replace generator with a root so that the result is nice. """ 

258 p = f.as_expr({gen: zero}) 

259 

260 if f.degree(gen) >= 2: 

261 p = p.expand(deep=False) 

262 

263 return p 

264 

265 def _solve_reduced_system(system, gens, entry=False): 

266 """Recursively solves reduced polynomial systems. """ 

267 if len(system) == len(gens) == 1: 

268 # the below line will produce UnsolvableFactorError if 

269 # strict=True and the solution from `roots` is incomplete 

270 zeros = list(roots(system[0], gens[-1], strict=strict).keys()) 

271 return [(zero,) for zero in zeros] 

272 

273 basis = groebner(system, gens, polys=True) 

274 

275 if len(basis) == 1 and basis[0].is_ground: 

276 if not entry: 

277 return [] 

278 else: 

279 return None 

280 

281 univariate = list(filter(_is_univariate, basis)) 

282 

283 if len(basis) < len(gens): 

284 raise NotImplementedError(filldedent(''' 

285 only zero-dimensional systems supported 

286 (finite number of solutions) 

287 ''')) 

288 

289 if len(univariate) == 1: 

290 f = univariate.pop() 

291 else: 

292 raise NotImplementedError(filldedent(''' 

293 only zero-dimensional systems supported 

294 (finite number of solutions) 

295 ''')) 

296 

297 gens = f.gens 

298 gen = gens[-1] 

299 

300 # the below line will produce UnsolvableFactorError if 

301 # strict=True and the solution from `roots` is incomplete 

302 zeros = list(roots(f.ltrim(gen), strict=strict).keys()) 

303 

304 if not zeros: 

305 return [] 

306 

307 if len(basis) == 1: 

308 return [(zero,) for zero in zeros] 

309 

310 solutions = [] 

311 

312 for zero in zeros: 

313 new_system = [] 

314 new_gens = gens[:-1] 

315 

316 for b in basis[:-1]: 

317 eq = _subs_root(b, gen, zero) 

318 

319 if eq is not S.Zero: 

320 new_system.append(eq) 

321 

322 for solution in _solve_reduced_system(new_system, new_gens): 

323 solutions.append(solution + (zero,)) 

324 

325 if solutions and len(solutions[0]) != len(gens): 

326 raise NotImplementedError(filldedent(''' 

327 only zero-dimensional systems supported 

328 (finite number of solutions) 

329 ''')) 

330 return solutions 

331 

332 try: 

333 result = _solve_reduced_system(polys, opt.gens, entry=True) 

334 except CoercionFailed: 

335 raise NotImplementedError 

336 

337 if result is not None: 

338 return sorted(result, key=default_sort_key) 

339 

340 

341def solve_triangulated(polys, *gens, **args): 

342 """ 

343 Solve a polynomial system using Gianni-Kalkbrenner algorithm. 

344 

345 The algorithm proceeds by computing one Groebner basis in the ground 

346 domain and then by iteratively computing polynomial factorizations in 

347 appropriately constructed algebraic extensions of the ground domain. 

348 

349 Parameters 

350 ========== 

351 

352 polys: a list/tuple/set 

353 Listing all the equations that are needed to be solved 

354 gens: generators 

355 generators of the equations in polys for which we want the 

356 solutions 

357 args: Keyword arguments 

358 Special options for solving the equations 

359 

360 Returns 

361 ======= 

362 

363 List[Tuple] 

364 A List of tuples. Solutions for symbols that satisfy the 

365 equations listed in polys 

366 

367 Examples 

368 ======== 

369 

370 >>> from sympy import solve_triangulated 

371 >>> from sympy.abc import x, y, z 

372 

373 >>> F = [x**2 + y + z - 1, x + y**2 + z - 1, x + y + z**2 - 1] 

374 

375 >>> solve_triangulated(F, x, y, z) 

376 [(0, 0, 1), (0, 1, 0), (1, 0, 0)] 

377 

378 References 

379 ========== 

380 

381 1. Patrizia Gianni, Teo Mora, Algebraic Solution of System of 

382 Polynomial Equations using Groebner Bases, AAECC-5 on Applied Algebra, 

383 Algebraic Algorithms and Error-Correcting Codes, LNCS 356 247--257, 1989 

384 

385 """ 

386 G = groebner(polys, gens, polys=True) 

387 G = list(reversed(G)) 

388 

389 domain = args.get('domain') 

390 

391 if domain is not None: 

392 for i, g in enumerate(G): 

393 G[i] = g.set_domain(domain) 

394 

395 f, G = G[0].ltrim(-1), G[1:] 

396 dom = f.get_domain() 

397 

398 zeros = f.ground_roots() 

399 solutions = set() 

400 

401 for zero in zeros: 

402 solutions.add(((zero,), dom)) 

403 

404 var_seq = reversed(gens[:-1]) 

405 vars_seq = postfixes(gens[1:]) 

406 

407 for var, vars in zip(var_seq, vars_seq): 

408 _solutions = set() 

409 

410 for values, dom in solutions: 

411 H, mapping = [], list(zip(vars, values)) 

412 

413 for g in G: 

414 _vars = (var,) + vars 

415 

416 if g.has_only_gens(*_vars) and g.degree(var) != 0: 

417 h = g.ltrim(var).eval(dict(mapping)) 

418 

419 if g.degree(var) == h.degree(): 

420 H.append(h) 

421 

422 p = min(H, key=lambda h: h.degree()) 

423 zeros = p.ground_roots() 

424 

425 for zero in zeros: 

426 if not zero.is_Rational: 

427 dom_zero = dom.algebraic_field(zero) 

428 else: 

429 dom_zero = dom 

430 

431 _solutions.add(((zero,) + values, dom_zero)) 

432 

433 solutions = _solutions 

434 

435 solutions = list(solutions) 

436 

437 for i, (solution, _) in enumerate(solutions): 

438 solutions[i] = solution 

439 

440 return sorted(solutions, key=default_sort_key)