Coverage for /usr/lib/python3/dist-packages/sympy/polys/polytools.py: 19%

2466 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1"""User-friendly public interface to polynomial functions. """ 

2 

3 

4from functools import wraps, reduce 

5from operator import mul 

6from typing import Optional 

7 

8from sympy.core import ( 

9 S, Expr, Add, Tuple 

10) 

11from sympy.core.basic import Basic 

12from sympy.core.decorators import _sympifyit 

13from sympy.core.exprtools import Factors, factor_nc, factor_terms 

14from sympy.core.evalf import ( 

15 pure_complex, evalf, fastlog, _evalf_with_bounded_error, quad_to_mpmath) 

16from sympy.core.function import Derivative 

17from sympy.core.mul import Mul, _keep_coeff 

18from sympy.core.numbers import ilcm, I, Integer, equal_valued 

19from sympy.core.relational import Relational, Equality 

20from sympy.core.sorting import ordered 

21from sympy.core.symbol import Dummy, Symbol 

22from sympy.core.sympify import sympify, _sympify 

23from sympy.core.traversal import preorder_traversal, bottom_up 

24from sympy.logic.boolalg import BooleanAtom 

25from sympy.polys import polyoptions as options 

26from sympy.polys.constructor import construct_domain 

27from sympy.polys.domains import FF, QQ, ZZ 

28from sympy.polys.domains.domainelement import DomainElement 

29from sympy.polys.fglmtools import matrix_fglm 

30from sympy.polys.groebnertools import groebner as _groebner 

31from sympy.polys.monomials import Monomial 

32from sympy.polys.orderings import monomial_key 

33from sympy.polys.polyclasses import DMP, DMF, ANP 

34from sympy.polys.polyerrors import ( 

35 OperationNotSupported, DomainError, 

36 CoercionFailed, UnificationFailed, 

37 GeneratorsNeeded, PolynomialError, 

38 MultivariatePolynomialError, 

39 ExactQuotientFailed, 

40 PolificationFailed, 

41 ComputationFailed, 

42 GeneratorsError, 

43) 

44from sympy.polys.polyutils import ( 

45 basic_from_dict, 

46 _sort_gens, 

47 _unify_gens, 

48 _dict_reorder, 

49 _dict_from_expr, 

50 _parallel_dict_from_expr, 

51) 

52from sympy.polys.rationaltools import together 

53from sympy.polys.rootisolation import dup_isolate_real_roots_list 

54from sympy.utilities import group, public, filldedent 

55from sympy.utilities.exceptions import sympy_deprecation_warning 

56from sympy.utilities.iterables import iterable, sift 

57 

58 

59# Required to avoid errors 

60import sympy.polys 

61 

62import mpmath 

63from mpmath.libmp.libhyper import NoConvergence 

64 

65 

66 

67def _polifyit(func): 

68 @wraps(func) 

69 def wrapper(f, g): 

70 g = _sympify(g) 

71 if isinstance(g, Poly): 

72 return func(f, g) 

73 elif isinstance(g, Expr): 

74 try: 

75 g = f.from_expr(g, *f.gens) 

76 except PolynomialError: 

77 if g.is_Matrix: 

78 return NotImplemented 

79 expr_method = getattr(f.as_expr(), func.__name__) 

80 result = expr_method(g) 

81 if result is not NotImplemented: 

82 sympy_deprecation_warning( 

83 """ 

84 Mixing Poly with non-polynomial expressions in binary 

85 operations is deprecated. Either explicitly convert 

86 the non-Poly operand to a Poly with as_poly() or 

87 convert the Poly to an Expr with as_expr(). 

88 """, 

89 deprecated_since_version="1.6", 

90 active_deprecations_target="deprecated-poly-nonpoly-binary-operations", 

91 ) 

92 return result 

93 else: 

94 return func(f, g) 

95 else: 

96 return NotImplemented 

97 return wrapper 

98 

99 

100 

101@public 

102class Poly(Basic): 

103 """ 

104 Generic class for representing and operating on polynomial expressions. 

105 

106 See :ref:`polys-docs` for general documentation. 

107 

108 Poly is a subclass of Basic rather than Expr but instances can be 

109 converted to Expr with the :py:meth:`~.Poly.as_expr` method. 

110 

111 .. deprecated:: 1.6 

112 

113 Combining Poly with non-Poly objects in binary operations is 

114 deprecated. Explicitly convert both objects to either Poly or Expr 

115 first. See :ref:`deprecated-poly-nonpoly-binary-operations`. 

116 

117 Examples 

118 ======== 

119 

120 >>> from sympy import Poly 

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

122 

123 Create a univariate polynomial: 

124 

125 >>> Poly(x*(x**2 + x - 1)**2) 

126 Poly(x**5 + 2*x**4 - x**3 - 2*x**2 + x, x, domain='ZZ') 

127 

128 Create a univariate polynomial with specific domain: 

129 

130 >>> from sympy import sqrt 

131 >>> Poly(x**2 + 2*x + sqrt(3), domain='R') 

132 Poly(1.0*x**2 + 2.0*x + 1.73205080756888, x, domain='RR') 

133 

134 Create a multivariate polynomial: 

135 

136 >>> Poly(y*x**2 + x*y + 1) 

137 Poly(x**2*y + x*y + 1, x, y, domain='ZZ') 

138 

139 Create a univariate polynomial, where y is a constant: 

140 

141 >>> Poly(y*x**2 + x*y + 1,x) 

142 Poly(y*x**2 + y*x + 1, x, domain='ZZ[y]') 

143 

144 You can evaluate the above polynomial as a function of y: 

145 

146 >>> Poly(y*x**2 + x*y + 1,x).eval(2) 

147 6*y + 1 

148 

149 See Also 

150 ======== 

151 

152 sympy.core.expr.Expr 

153 

154 """ 

155 

156 __slots__ = ('rep', 'gens') 

157 

158 is_commutative = True 

159 is_Poly = True 

160 _op_priority = 10.001 

161 

162 def __new__(cls, rep, *gens, **args): 

163 """Create a new polynomial instance out of something useful. """ 

164 opt = options.build_options(gens, args) 

165 

166 if 'order' in opt: 

167 raise NotImplementedError("'order' keyword is not implemented yet") 

168 

169 if isinstance(rep, (DMP, DMF, ANP, DomainElement)): 

170 return cls._from_domain_element(rep, opt) 

171 elif iterable(rep, exclude=str): 

172 if isinstance(rep, dict): 

173 return cls._from_dict(rep, opt) 

174 else: 

175 return cls._from_list(list(rep), opt) 

176 else: 

177 rep = sympify(rep) 

178 

179 if rep.is_Poly: 

180 return cls._from_poly(rep, opt) 

181 else: 

182 return cls._from_expr(rep, opt) 

183 

184 # Poly does not pass its args to Basic.__new__ to be stored in _args so we 

185 # have to emulate them here with an args property that derives from rep 

186 # and gens which are instance attributes. This also means we need to 

187 # define _hashable_content. The _hashable_content is rep and gens but args 

188 # uses expr instead of rep (expr is the Basic version of rep). Passing 

189 # expr in args means that Basic methods like subs should work. Using rep 

190 # otherwise means that Poly can remain more efficient than Basic by 

191 # avoiding creating a Basic instance just to be hashable. 

192 

193 @classmethod 

194 def new(cls, rep, *gens): 

195 """Construct :class:`Poly` instance from raw representation. """ 

196 if not isinstance(rep, DMP): 

197 raise PolynomialError( 

198 "invalid polynomial representation: %s" % rep) 

199 elif rep.lev != len(gens) - 1: 

200 raise PolynomialError("invalid arguments: %s, %s" % (rep, gens)) 

201 

202 obj = Basic.__new__(cls) 

203 obj.rep = rep 

204 obj.gens = gens 

205 

206 return obj 

207 

208 @property 

209 def expr(self): 

210 return basic_from_dict(self.rep.to_sympy_dict(), *self.gens) 

211 

212 @property 

213 def args(self): 

214 return (self.expr,) + self.gens 

215 

216 def _hashable_content(self): 

217 return (self.rep,) + self.gens 

218 

219 @classmethod 

220 def from_dict(cls, rep, *gens, **args): 

221 """Construct a polynomial from a ``dict``. """ 

222 opt = options.build_options(gens, args) 

223 return cls._from_dict(rep, opt) 

224 

225 @classmethod 

226 def from_list(cls, rep, *gens, **args): 

227 """Construct a polynomial from a ``list``. """ 

228 opt = options.build_options(gens, args) 

229 return cls._from_list(rep, opt) 

230 

231 @classmethod 

232 def from_poly(cls, rep, *gens, **args): 

233 """Construct a polynomial from a polynomial. """ 

234 opt = options.build_options(gens, args) 

235 return cls._from_poly(rep, opt) 

236 

237 @classmethod 

238 def from_expr(cls, rep, *gens, **args): 

239 """Construct a polynomial from an expression. """ 

240 opt = options.build_options(gens, args) 

241 return cls._from_expr(rep, opt) 

242 

243 @classmethod 

244 def _from_dict(cls, rep, opt): 

245 """Construct a polynomial from a ``dict``. """ 

246 gens = opt.gens 

247 

248 if not gens: 

249 raise GeneratorsNeeded( 

250 "Cannot initialize from 'dict' without generators") 

251 

252 level = len(gens) - 1 

253 domain = opt.domain 

254 

255 if domain is None: 

256 domain, rep = construct_domain(rep, opt=opt) 

257 else: 

258 for monom, coeff in rep.items(): 

259 rep[monom] = domain.convert(coeff) 

260 

261 return cls.new(DMP.from_dict(rep, level, domain), *gens) 

262 

263 @classmethod 

264 def _from_list(cls, rep, opt): 

265 """Construct a polynomial from a ``list``. """ 

266 gens = opt.gens 

267 

268 if not gens: 

269 raise GeneratorsNeeded( 

270 "Cannot initialize from 'list' without generators") 

271 elif len(gens) != 1: 

272 raise MultivariatePolynomialError( 

273 "'list' representation not supported") 

274 

275 level = len(gens) - 1 

276 domain = opt.domain 

277 

278 if domain is None: 

279 domain, rep = construct_domain(rep, opt=opt) 

280 else: 

281 rep = list(map(domain.convert, rep)) 

282 

283 return cls.new(DMP.from_list(rep, level, domain), *gens) 

284 

285 @classmethod 

286 def _from_poly(cls, rep, opt): 

287 """Construct a polynomial from a polynomial. """ 

288 if cls != rep.__class__: 

289 rep = cls.new(rep.rep, *rep.gens) 

290 

291 gens = opt.gens 

292 field = opt.field 

293 domain = opt.domain 

294 

295 if gens and rep.gens != gens: 

296 if set(rep.gens) != set(gens): 

297 return cls._from_expr(rep.as_expr(), opt) 

298 else: 

299 rep = rep.reorder(*gens) 

300 

301 if 'domain' in opt and domain: 

302 rep = rep.set_domain(domain) 

303 elif field is True: 

304 rep = rep.to_field() 

305 

306 return rep 

307 

308 @classmethod 

309 def _from_expr(cls, rep, opt): 

310 """Construct a polynomial from an expression. """ 

311 rep, opt = _dict_from_expr(rep, opt) 

312 return cls._from_dict(rep, opt) 

313 

314 @classmethod 

315 def _from_domain_element(cls, rep, opt): 

316 gens = opt.gens 

317 domain = opt.domain 

318 

319 level = len(gens) - 1 

320 rep = [domain.convert(rep)] 

321 

322 return cls.new(DMP.from_list(rep, level, domain), *gens) 

323 

324 def __hash__(self): 

325 return super().__hash__() 

326 

327 @property 

328 def free_symbols(self): 

329 """ 

330 Free symbols of a polynomial expression. 

331 

332 Examples 

333 ======== 

334 

335 >>> from sympy import Poly 

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

337 

338 >>> Poly(x**2 + 1).free_symbols 

339 {x} 

340 >>> Poly(x**2 + y).free_symbols 

341 {x, y} 

342 >>> Poly(x**2 + y, x).free_symbols 

343 {x, y} 

344 >>> Poly(x**2 + y, x, z).free_symbols 

345 {x, y} 

346 

347 """ 

348 symbols = set() 

349 gens = self.gens 

350 for i in range(len(gens)): 

351 for monom in self.monoms(): 

352 if monom[i]: 

353 symbols |= gens[i].free_symbols 

354 break 

355 

356 return symbols | self.free_symbols_in_domain 

357 

358 @property 

359 def free_symbols_in_domain(self): 

360 """ 

361 Free symbols of the domain of ``self``. 

362 

363 Examples 

364 ======== 

365 

366 >>> from sympy import Poly 

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

368 

369 >>> Poly(x**2 + 1).free_symbols_in_domain 

370 set() 

371 >>> Poly(x**2 + y).free_symbols_in_domain 

372 set() 

373 >>> Poly(x**2 + y, x).free_symbols_in_domain 

374 {y} 

375 

376 """ 

377 domain, symbols = self.rep.dom, set() 

378 

379 if domain.is_Composite: 

380 for gen in domain.symbols: 

381 symbols |= gen.free_symbols 

382 elif domain.is_EX: 

383 for coeff in self.coeffs(): 

384 symbols |= coeff.free_symbols 

385 

386 return symbols 

387 

388 @property 

389 def gen(self): 

390 """ 

391 Return the principal generator. 

392 

393 Examples 

394 ======== 

395 

396 >>> from sympy import Poly 

397 >>> from sympy.abc import x 

398 

399 >>> Poly(x**2 + 1, x).gen 

400 x 

401 

402 """ 

403 return self.gens[0] 

404 

405 @property 

406 def domain(self): 

407 """Get the ground domain of a :py:class:`~.Poly` 

408 

409 Returns 

410 ======= 

411 

412 :py:class:`~.Domain`: 

413 Ground domain of the :py:class:`~.Poly`. 

414 

415 Examples 

416 ======== 

417 

418 >>> from sympy import Poly, Symbol 

419 >>> x = Symbol('x') 

420 >>> p = Poly(x**2 + x) 

421 >>> p 

422 Poly(x**2 + x, x, domain='ZZ') 

423 >>> p.domain 

424 ZZ 

425 """ 

426 return self.get_domain() 

427 

428 @property 

429 def zero(self): 

430 """Return zero polynomial with ``self``'s properties. """ 

431 return self.new(self.rep.zero(self.rep.lev, self.rep.dom), *self.gens) 

432 

433 @property 

434 def one(self): 

435 """Return one polynomial with ``self``'s properties. """ 

436 return self.new(self.rep.one(self.rep.lev, self.rep.dom), *self.gens) 

437 

438 @property 

439 def unit(self): 

440 """Return unit polynomial with ``self``'s properties. """ 

441 return self.new(self.rep.unit(self.rep.lev, self.rep.dom), *self.gens) 

442 

443 def unify(f, g): 

444 """ 

445 Make ``f`` and ``g`` belong to the same domain. 

446 

447 Examples 

448 ======== 

449 

450 >>> from sympy import Poly 

451 >>> from sympy.abc import x 

452 

453 >>> f, g = Poly(x/2 + 1), Poly(2*x + 1) 

454 

455 >>> f 

456 Poly(1/2*x + 1, x, domain='QQ') 

457 >>> g 

458 Poly(2*x + 1, x, domain='ZZ') 

459 

460 >>> F, G = f.unify(g) 

461 

462 >>> F 

463 Poly(1/2*x + 1, x, domain='QQ') 

464 >>> G 

465 Poly(2*x + 1, x, domain='QQ') 

466 

467 """ 

468 _, per, F, G = f._unify(g) 

469 return per(F), per(G) 

470 

471 def _unify(f, g): 

472 g = sympify(g) 

473 

474 if not g.is_Poly: 

475 try: 

476 return f.rep.dom, f.per, f.rep, f.rep.per(f.rep.dom.from_sympy(g)) 

477 except CoercionFailed: 

478 raise UnificationFailed("Cannot unify %s with %s" % (f, g)) 

479 

480 if isinstance(f.rep, DMP) and isinstance(g.rep, DMP): 

481 gens = _unify_gens(f.gens, g.gens) 

482 

483 dom, lev = f.rep.dom.unify(g.rep.dom, gens), len(gens) - 1 

484 

485 if f.gens != gens: 

486 f_monoms, f_coeffs = _dict_reorder( 

487 f.rep.to_dict(), f.gens, gens) 

488 

489 if f.rep.dom != dom: 

490 f_coeffs = [dom.convert(c, f.rep.dom) for c in f_coeffs] 

491 

492 F = DMP(dict(list(zip(f_monoms, f_coeffs))), dom, lev) 

493 else: 

494 F = f.rep.convert(dom) 

495 

496 if g.gens != gens: 

497 g_monoms, g_coeffs = _dict_reorder( 

498 g.rep.to_dict(), g.gens, gens) 

499 

500 if g.rep.dom != dom: 

501 g_coeffs = [dom.convert(c, g.rep.dom) for c in g_coeffs] 

502 

503 G = DMP(dict(list(zip(g_monoms, g_coeffs))), dom, lev) 

504 else: 

505 G = g.rep.convert(dom) 

506 else: 

507 raise UnificationFailed("Cannot unify %s with %s" % (f, g)) 

508 

509 cls = f.__class__ 

510 

511 def per(rep, dom=dom, gens=gens, remove=None): 

512 if remove is not None: 

513 gens = gens[:remove] + gens[remove + 1:] 

514 

515 if not gens: 

516 return dom.to_sympy(rep) 

517 

518 return cls.new(rep, *gens) 

519 

520 return dom, per, F, G 

521 

522 def per(f, rep, gens=None, remove=None): 

523 """ 

524 Create a Poly out of the given representation. 

525 

526 Examples 

527 ======== 

528 

529 >>> from sympy import Poly, ZZ 

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

531 

532 >>> from sympy.polys.polyclasses import DMP 

533 

534 >>> a = Poly(x**2 + 1) 

535 

536 >>> a.per(DMP([ZZ(1), ZZ(1)], ZZ), gens=[y]) 

537 Poly(y + 1, y, domain='ZZ') 

538 

539 """ 

540 if gens is None: 

541 gens = f.gens 

542 

543 if remove is not None: 

544 gens = gens[:remove] + gens[remove + 1:] 

545 

546 if not gens: 

547 return f.rep.dom.to_sympy(rep) 

548 

549 return f.__class__.new(rep, *gens) 

550 

551 def set_domain(f, domain): 

552 """Set the ground domain of ``f``. """ 

553 opt = options.build_options(f.gens, {'domain': domain}) 

554 return f.per(f.rep.convert(opt.domain)) 

555 

556 def get_domain(f): 

557 """Get the ground domain of ``f``. """ 

558 return f.rep.dom 

559 

560 def set_modulus(f, modulus): 

561 """ 

562 Set the modulus of ``f``. 

563 

564 Examples 

565 ======== 

566 

567 >>> from sympy import Poly 

568 >>> from sympy.abc import x 

569 

570 >>> Poly(5*x**2 + 2*x - 1, x).set_modulus(2) 

571 Poly(x**2 + 1, x, modulus=2) 

572 

573 """ 

574 modulus = options.Modulus.preprocess(modulus) 

575 return f.set_domain(FF(modulus)) 

576 

577 def get_modulus(f): 

578 """ 

579 Get the modulus of ``f``. 

580 

581 Examples 

582 ======== 

583 

584 >>> from sympy import Poly 

585 >>> from sympy.abc import x 

586 

587 >>> Poly(x**2 + 1, modulus=2).get_modulus() 

588 2 

589 

590 """ 

591 domain = f.get_domain() 

592 

593 if domain.is_FiniteField: 

594 return Integer(domain.characteristic()) 

595 else: 

596 raise PolynomialError("not a polynomial over a Galois field") 

597 

598 def _eval_subs(f, old, new): 

599 """Internal implementation of :func:`subs`. """ 

600 if old in f.gens: 

601 if new.is_number: 

602 return f.eval(old, new) 

603 else: 

604 try: 

605 return f.replace(old, new) 

606 except PolynomialError: 

607 pass 

608 

609 return f.as_expr().subs(old, new) 

610 

611 def exclude(f): 

612 """ 

613 Remove unnecessary generators from ``f``. 

614 

615 Examples 

616 ======== 

617 

618 >>> from sympy import Poly 

619 >>> from sympy.abc import a, b, c, d, x 

620 

621 >>> Poly(a + x, a, b, c, d, x).exclude() 

622 Poly(a + x, a, x, domain='ZZ') 

623 

624 """ 

625 J, new = f.rep.exclude() 

626 gens = [gen for j, gen in enumerate(f.gens) if j not in J] 

627 

628 return f.per(new, gens=gens) 

629 

630 def replace(f, x, y=None, **_ignore): 

631 # XXX this does not match Basic's signature 

632 """ 

633 Replace ``x`` with ``y`` in generators list. 

634 

635 Examples 

636 ======== 

637 

638 >>> from sympy import Poly 

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

640 

641 >>> Poly(x**2 + 1, x).replace(x, y) 

642 Poly(y**2 + 1, y, domain='ZZ') 

643 

644 """ 

645 if y is None: 

646 if f.is_univariate: 

647 x, y = f.gen, x 

648 else: 

649 raise PolynomialError( 

650 "syntax supported only in univariate case") 

651 

652 if x == y or x not in f.gens: 

653 return f 

654 

655 if x in f.gens and y not in f.gens: 

656 dom = f.get_domain() 

657 

658 if not dom.is_Composite or y not in dom.symbols: 

659 gens = list(f.gens) 

660 gens[gens.index(x)] = y 

661 return f.per(f.rep, gens=gens) 

662 

663 raise PolynomialError("Cannot replace %s with %s in %s" % (x, y, f)) 

664 

665 def match(f, *args, **kwargs): 

666 """Match expression from Poly. See Basic.match()""" 

667 return f.as_expr().match(*args, **kwargs) 

668 

669 def reorder(f, *gens, **args): 

670 """ 

671 Efficiently apply new order of generators. 

672 

673 Examples 

674 ======== 

675 

676 >>> from sympy import Poly 

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

678 

679 >>> Poly(x**2 + x*y**2, x, y).reorder(y, x) 

680 Poly(y**2*x + x**2, y, x, domain='ZZ') 

681 

682 """ 

683 opt = options.Options((), args) 

684 

685 if not gens: 

686 gens = _sort_gens(f.gens, opt=opt) 

687 elif set(f.gens) != set(gens): 

688 raise PolynomialError( 

689 "generators list can differ only up to order of elements") 

690 

691 rep = dict(list(zip(*_dict_reorder(f.rep.to_dict(), f.gens, gens)))) 

692 

693 return f.per(DMP(rep, f.rep.dom, len(gens) - 1), gens=gens) 

694 

695 def ltrim(f, gen): 

696 """ 

697 Remove dummy generators from ``f`` that are to the left of 

698 specified ``gen`` in the generators as ordered. When ``gen`` 

699 is an integer, it refers to the generator located at that 

700 position within the tuple of generators of ``f``. 

701 

702 Examples 

703 ======== 

704 

705 >>> from sympy import Poly 

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

707 

708 >>> Poly(y**2 + y*z**2, x, y, z).ltrim(y) 

709 Poly(y**2 + y*z**2, y, z, domain='ZZ') 

710 >>> Poly(z, x, y, z).ltrim(-1) 

711 Poly(z, z, domain='ZZ') 

712 

713 """ 

714 rep = f.as_dict(native=True) 

715 j = f._gen_to_level(gen) 

716 

717 terms = {} 

718 

719 for monom, coeff in rep.items(): 

720 

721 if any(monom[:j]): 

722 # some generator is used in the portion to be trimmed 

723 raise PolynomialError("Cannot left trim %s" % f) 

724 

725 terms[monom[j:]] = coeff 

726 

727 gens = f.gens[j:] 

728 

729 return f.new(DMP.from_dict(terms, len(gens) - 1, f.rep.dom), *gens) 

730 

731 def has_only_gens(f, *gens): 

732 """ 

733 Return ``True`` if ``Poly(f, *gens)`` retains ground domain. 

734 

735 Examples 

736 ======== 

737 

738 >>> from sympy import Poly 

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

740 

741 >>> Poly(x*y + 1, x, y, z).has_only_gens(x, y) 

742 True 

743 >>> Poly(x*y + z, x, y, z).has_only_gens(x, y) 

744 False 

745 

746 """ 

747 indices = set() 

748 

749 for gen in gens: 

750 try: 

751 index = f.gens.index(gen) 

752 except ValueError: 

753 raise GeneratorsError( 

754 "%s doesn't have %s as generator" % (f, gen)) 

755 else: 

756 indices.add(index) 

757 

758 for monom in f.monoms(): 

759 for i, elt in enumerate(monom): 

760 if i not in indices and elt: 

761 return False 

762 

763 return True 

764 

765 def to_ring(f): 

766 """ 

767 Make the ground domain a ring. 

768 

769 Examples 

770 ======== 

771 

772 >>> from sympy import Poly, QQ 

773 >>> from sympy.abc import x 

774 

775 >>> Poly(x**2 + 1, domain=QQ).to_ring() 

776 Poly(x**2 + 1, x, domain='ZZ') 

777 

778 """ 

779 if hasattr(f.rep, 'to_ring'): 

780 result = f.rep.to_ring() 

781 else: # pragma: no cover 

782 raise OperationNotSupported(f, 'to_ring') 

783 

784 return f.per(result) 

785 

786 def to_field(f): 

787 """ 

788 Make the ground domain a field. 

789 

790 Examples 

791 ======== 

792 

793 >>> from sympy import Poly, ZZ 

794 >>> from sympy.abc import x 

795 

796 >>> Poly(x**2 + 1, x, domain=ZZ).to_field() 

797 Poly(x**2 + 1, x, domain='QQ') 

798 

799 """ 

800 if hasattr(f.rep, 'to_field'): 

801 result = f.rep.to_field() 

802 else: # pragma: no cover 

803 raise OperationNotSupported(f, 'to_field') 

804 

805 return f.per(result) 

806 

807 def to_exact(f): 

808 """ 

809 Make the ground domain exact. 

810 

811 Examples 

812 ======== 

813 

814 >>> from sympy import Poly, RR 

815 >>> from sympy.abc import x 

816 

817 >>> Poly(x**2 + 1.0, x, domain=RR).to_exact() 

818 Poly(x**2 + 1, x, domain='QQ') 

819 

820 """ 

821 if hasattr(f.rep, 'to_exact'): 

822 result = f.rep.to_exact() 

823 else: # pragma: no cover 

824 raise OperationNotSupported(f, 'to_exact') 

825 

826 return f.per(result) 

827 

828 def retract(f, field=None): 

829 """ 

830 Recalculate the ground domain of a polynomial. 

831 

832 Examples 

833 ======== 

834 

835 >>> from sympy import Poly 

836 >>> from sympy.abc import x 

837 

838 >>> f = Poly(x**2 + 1, x, domain='QQ[y]') 

839 >>> f 

840 Poly(x**2 + 1, x, domain='QQ[y]') 

841 

842 >>> f.retract() 

843 Poly(x**2 + 1, x, domain='ZZ') 

844 >>> f.retract(field=True) 

845 Poly(x**2 + 1, x, domain='QQ') 

846 

847 """ 

848 dom, rep = construct_domain(f.as_dict(zero=True), 

849 field=field, composite=f.domain.is_Composite or None) 

850 return f.from_dict(rep, f.gens, domain=dom) 

851 

852 def slice(f, x, m, n=None): 

853 """Take a continuous subsequence of terms of ``f``. """ 

854 if n is None: 

855 j, m, n = 0, x, m 

856 else: 

857 j = f._gen_to_level(x) 

858 

859 m, n = int(m), int(n) 

860 

861 if hasattr(f.rep, 'slice'): 

862 result = f.rep.slice(m, n, j) 

863 else: # pragma: no cover 

864 raise OperationNotSupported(f, 'slice') 

865 

866 return f.per(result) 

867 

868 def coeffs(f, order=None): 

869 """ 

870 Returns all non-zero coefficients from ``f`` in lex order. 

871 

872 Examples 

873 ======== 

874 

875 >>> from sympy import Poly 

876 >>> from sympy.abc import x 

877 

878 >>> Poly(x**3 + 2*x + 3, x).coeffs() 

879 [1, 2, 3] 

880 

881 See Also 

882 ======== 

883 all_coeffs 

884 coeff_monomial 

885 nth 

886 

887 """ 

888 return [f.rep.dom.to_sympy(c) for c in f.rep.coeffs(order=order)] 

889 

890 def monoms(f, order=None): 

891 """ 

892 Returns all non-zero monomials from ``f`` in lex order. 

893 

894 Examples 

895 ======== 

896 

897 >>> from sympy import Poly 

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

899 

900 >>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).monoms() 

901 [(2, 0), (1, 2), (1, 1), (0, 1)] 

902 

903 See Also 

904 ======== 

905 all_monoms 

906 

907 """ 

908 return f.rep.monoms(order=order) 

909 

910 def terms(f, order=None): 

911 """ 

912 Returns all non-zero terms from ``f`` in lex order. 

913 

914 Examples 

915 ======== 

916 

917 >>> from sympy import Poly 

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

919 

920 >>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).terms() 

921 [((2, 0), 1), ((1, 2), 2), ((1, 1), 1), ((0, 1), 3)] 

922 

923 See Also 

924 ======== 

925 all_terms 

926 

927 """ 

928 return [(m, f.rep.dom.to_sympy(c)) for m, c in f.rep.terms(order=order)] 

929 

930 def all_coeffs(f): 

931 """ 

932 Returns all coefficients from a univariate polynomial ``f``. 

933 

934 Examples 

935 ======== 

936 

937 >>> from sympy import Poly 

938 >>> from sympy.abc import x 

939 

940 >>> Poly(x**3 + 2*x - 1, x).all_coeffs() 

941 [1, 0, 2, -1] 

942 

943 """ 

944 return [f.rep.dom.to_sympy(c) for c in f.rep.all_coeffs()] 

945 

946 def all_monoms(f): 

947 """ 

948 Returns all monomials from a univariate polynomial ``f``. 

949 

950 Examples 

951 ======== 

952 

953 >>> from sympy import Poly 

954 >>> from sympy.abc import x 

955 

956 >>> Poly(x**3 + 2*x - 1, x).all_monoms() 

957 [(3,), (2,), (1,), (0,)] 

958 

959 See Also 

960 ======== 

961 all_terms 

962 

963 """ 

964 return f.rep.all_monoms() 

965 

966 def all_terms(f): 

967 """ 

968 Returns all terms from a univariate polynomial ``f``. 

969 

970 Examples 

971 ======== 

972 

973 >>> from sympy import Poly 

974 >>> from sympy.abc import x 

975 

976 >>> Poly(x**3 + 2*x - 1, x).all_terms() 

977 [((3,), 1), ((2,), 0), ((1,), 2), ((0,), -1)] 

978 

979 """ 

980 return [(m, f.rep.dom.to_sympy(c)) for m, c in f.rep.all_terms()] 

981 

982 def termwise(f, func, *gens, **args): 

983 """ 

984 Apply a function to all terms of ``f``. 

985 

986 Examples 

987 ======== 

988 

989 >>> from sympy import Poly 

990 >>> from sympy.abc import x 

991 

992 >>> def func(k, coeff): 

993 ... k = k[0] 

994 ... return coeff//10**(2-k) 

995 

996 >>> Poly(x**2 + 20*x + 400).termwise(func) 

997 Poly(x**2 + 2*x + 4, x, domain='ZZ') 

998 

999 """ 

1000 terms = {} 

1001 

1002 for monom, coeff in f.terms(): 

1003 result = func(monom, coeff) 

1004 

1005 if isinstance(result, tuple): 

1006 monom, coeff = result 

1007 else: 

1008 coeff = result 

1009 

1010 if coeff: 

1011 if monom not in terms: 

1012 terms[monom] = coeff 

1013 else: 

1014 raise PolynomialError( 

1015 "%s monomial was generated twice" % monom) 

1016 

1017 return f.from_dict(terms, *(gens or f.gens), **args) 

1018 

1019 def length(f): 

1020 """ 

1021 Returns the number of non-zero terms in ``f``. 

1022 

1023 Examples 

1024 ======== 

1025 

1026 >>> from sympy import Poly 

1027 >>> from sympy.abc import x 

1028 

1029 >>> Poly(x**2 + 2*x - 1).length() 

1030 3 

1031 

1032 """ 

1033 return len(f.as_dict()) 

1034 

1035 def as_dict(f, native=False, zero=False): 

1036 """ 

1037 Switch to a ``dict`` representation. 

1038 

1039 Examples 

1040 ======== 

1041 

1042 >>> from sympy import Poly 

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

1044 

1045 >>> Poly(x**2 + 2*x*y**2 - y, x, y).as_dict() 

1046 {(0, 1): -1, (1, 2): 2, (2, 0): 1} 

1047 

1048 """ 

1049 if native: 

1050 return f.rep.to_dict(zero=zero) 

1051 else: 

1052 return f.rep.to_sympy_dict(zero=zero) 

1053 

1054 def as_list(f, native=False): 

1055 """Switch to a ``list`` representation. """ 

1056 if native: 

1057 return f.rep.to_list() 

1058 else: 

1059 return f.rep.to_sympy_list() 

1060 

1061 def as_expr(f, *gens): 

1062 """ 

1063 Convert a Poly instance to an Expr instance. 

1064 

1065 Examples 

1066 ======== 

1067 

1068 >>> from sympy import Poly 

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

1070 

1071 >>> f = Poly(x**2 + 2*x*y**2 - y, x, y) 

1072 

1073 >>> f.as_expr() 

1074 x**2 + 2*x*y**2 - y 

1075 >>> f.as_expr({x: 5}) 

1076 10*y**2 - y + 25 

1077 >>> f.as_expr(5, 6) 

1078 379 

1079 

1080 """ 

1081 if not gens: 

1082 return f.expr 

1083 

1084 if len(gens) == 1 and isinstance(gens[0], dict): 

1085 mapping = gens[0] 

1086 gens = list(f.gens) 

1087 

1088 for gen, value in mapping.items(): 

1089 try: 

1090 index = gens.index(gen) 

1091 except ValueError: 

1092 raise GeneratorsError( 

1093 "%s doesn't have %s as generator" % (f, gen)) 

1094 else: 

1095 gens[index] = value 

1096 

1097 return basic_from_dict(f.rep.to_sympy_dict(), *gens) 

1098 

1099 def as_poly(self, *gens, **args): 

1100 """Converts ``self`` to a polynomial or returns ``None``. 

1101 

1102 >>> from sympy import sin 

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

1104 

1105 >>> print((x**2 + x*y).as_poly()) 

1106 Poly(x**2 + x*y, x, y, domain='ZZ') 

1107 

1108 >>> print((x**2 + x*y).as_poly(x, y)) 

1109 Poly(x**2 + x*y, x, y, domain='ZZ') 

1110 

1111 >>> print((x**2 + sin(y)).as_poly(x, y)) 

1112 None 

1113 

1114 """ 

1115 try: 

1116 poly = Poly(self, *gens, **args) 

1117 

1118 if not poly.is_Poly: 

1119 return None 

1120 else: 

1121 return poly 

1122 except PolynomialError: 

1123 return None 

1124 

1125 def lift(f): 

1126 """ 

1127 Convert algebraic coefficients to rationals. 

1128 

1129 Examples 

1130 ======== 

1131 

1132 >>> from sympy import Poly, I 

1133 >>> from sympy.abc import x 

1134 

1135 >>> Poly(x**2 + I*x + 1, x, extension=I).lift() 

1136 Poly(x**4 + 3*x**2 + 1, x, domain='QQ') 

1137 

1138 """ 

1139 if hasattr(f.rep, 'lift'): 

1140 result = f.rep.lift() 

1141 else: # pragma: no cover 

1142 raise OperationNotSupported(f, 'lift') 

1143 

1144 return f.per(result) 

1145 

1146 def deflate(f): 

1147 """ 

1148 Reduce degree of ``f`` by mapping ``x_i**m`` to ``y_i``. 

1149 

1150 Examples 

1151 ======== 

1152 

1153 >>> from sympy import Poly 

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

1155 

1156 >>> Poly(x**6*y**2 + x**3 + 1, x, y).deflate() 

1157 ((3, 2), Poly(x**2*y + x + 1, x, y, domain='ZZ')) 

1158 

1159 """ 

1160 if hasattr(f.rep, 'deflate'): 

1161 J, result = f.rep.deflate() 

1162 else: # pragma: no cover 

1163 raise OperationNotSupported(f, 'deflate') 

1164 

1165 return J, f.per(result) 

1166 

1167 def inject(f, front=False): 

1168 """ 

1169 Inject ground domain generators into ``f``. 

1170 

1171 Examples 

1172 ======== 

1173 

1174 >>> from sympy import Poly 

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

1176 

1177 >>> f = Poly(x**2*y + x*y**3 + x*y + 1, x) 

1178 

1179 >>> f.inject() 

1180 Poly(x**2*y + x*y**3 + x*y + 1, x, y, domain='ZZ') 

1181 >>> f.inject(front=True) 

1182 Poly(y**3*x + y*x**2 + y*x + 1, y, x, domain='ZZ') 

1183 

1184 """ 

1185 dom = f.rep.dom 

1186 

1187 if dom.is_Numerical: 

1188 return f 

1189 elif not dom.is_Poly: 

1190 raise DomainError("Cannot inject generators over %s" % dom) 

1191 

1192 if hasattr(f.rep, 'inject'): 

1193 result = f.rep.inject(front=front) 

1194 else: # pragma: no cover 

1195 raise OperationNotSupported(f, 'inject') 

1196 

1197 if front: 

1198 gens = dom.symbols + f.gens 

1199 else: 

1200 gens = f.gens + dom.symbols 

1201 

1202 return f.new(result, *gens) 

1203 

1204 def eject(f, *gens): 

1205 """ 

1206 Eject selected generators into the ground domain. 

1207 

1208 Examples 

1209 ======== 

1210 

1211 >>> from sympy import Poly 

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

1213 

1214 >>> f = Poly(x**2*y + x*y**3 + x*y + 1, x, y) 

1215 

1216 >>> f.eject(x) 

1217 Poly(x*y**3 + (x**2 + x)*y + 1, y, domain='ZZ[x]') 

1218 >>> f.eject(y) 

1219 Poly(y*x**2 + (y**3 + y)*x + 1, x, domain='ZZ[y]') 

1220 

1221 """ 

1222 dom = f.rep.dom 

1223 

1224 if not dom.is_Numerical: 

1225 raise DomainError("Cannot eject generators over %s" % dom) 

1226 

1227 k = len(gens) 

1228 

1229 if f.gens[:k] == gens: 

1230 _gens, front = f.gens[k:], True 

1231 elif f.gens[-k:] == gens: 

1232 _gens, front = f.gens[:-k], False 

1233 else: 

1234 raise NotImplementedError( 

1235 "can only eject front or back generators") 

1236 

1237 dom = dom.inject(*gens) 

1238 

1239 if hasattr(f.rep, 'eject'): 

1240 result = f.rep.eject(dom, front=front) 

1241 else: # pragma: no cover 

1242 raise OperationNotSupported(f, 'eject') 

1243 

1244 return f.new(result, *_gens) 

1245 

1246 def terms_gcd(f): 

1247 """ 

1248 Remove GCD of terms from the polynomial ``f``. 

1249 

1250 Examples 

1251 ======== 

1252 

1253 >>> from sympy import Poly 

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

1255 

1256 >>> Poly(x**6*y**2 + x**3*y, x, y).terms_gcd() 

1257 ((3, 1), Poly(x**3*y + 1, x, y, domain='ZZ')) 

1258 

1259 """ 

1260 if hasattr(f.rep, 'terms_gcd'): 

1261 J, result = f.rep.terms_gcd() 

1262 else: # pragma: no cover 

1263 raise OperationNotSupported(f, 'terms_gcd') 

1264 

1265 return J, f.per(result) 

1266 

1267 def add_ground(f, coeff): 

1268 """ 

1269 Add an element of the ground domain to ``f``. 

1270 

1271 Examples 

1272 ======== 

1273 

1274 >>> from sympy import Poly 

1275 >>> from sympy.abc import x 

1276 

1277 >>> Poly(x + 1).add_ground(2) 

1278 Poly(x + 3, x, domain='ZZ') 

1279 

1280 """ 

1281 if hasattr(f.rep, 'add_ground'): 

1282 result = f.rep.add_ground(coeff) 

1283 else: # pragma: no cover 

1284 raise OperationNotSupported(f, 'add_ground') 

1285 

1286 return f.per(result) 

1287 

1288 def sub_ground(f, coeff): 

1289 """ 

1290 Subtract an element of the ground domain from ``f``. 

1291 

1292 Examples 

1293 ======== 

1294 

1295 >>> from sympy import Poly 

1296 >>> from sympy.abc import x 

1297 

1298 >>> Poly(x + 1).sub_ground(2) 

1299 Poly(x - 1, x, domain='ZZ') 

1300 

1301 """ 

1302 if hasattr(f.rep, 'sub_ground'): 

1303 result = f.rep.sub_ground(coeff) 

1304 else: # pragma: no cover 

1305 raise OperationNotSupported(f, 'sub_ground') 

1306 

1307 return f.per(result) 

1308 

1309 def mul_ground(f, coeff): 

1310 """ 

1311 Multiply ``f`` by a an element of the ground domain. 

1312 

1313 Examples 

1314 ======== 

1315 

1316 >>> from sympy import Poly 

1317 >>> from sympy.abc import x 

1318 

1319 >>> Poly(x + 1).mul_ground(2) 

1320 Poly(2*x + 2, x, domain='ZZ') 

1321 

1322 """ 

1323 if hasattr(f.rep, 'mul_ground'): 

1324 result = f.rep.mul_ground(coeff) 

1325 else: # pragma: no cover 

1326 raise OperationNotSupported(f, 'mul_ground') 

1327 

1328 return f.per(result) 

1329 

1330 def quo_ground(f, coeff): 

1331 """ 

1332 Quotient of ``f`` by a an element of the ground domain. 

1333 

1334 Examples 

1335 ======== 

1336 

1337 >>> from sympy import Poly 

1338 >>> from sympy.abc import x 

1339 

1340 >>> Poly(2*x + 4).quo_ground(2) 

1341 Poly(x + 2, x, domain='ZZ') 

1342 

1343 >>> Poly(2*x + 3).quo_ground(2) 

1344 Poly(x + 1, x, domain='ZZ') 

1345 

1346 """ 

1347 if hasattr(f.rep, 'quo_ground'): 

1348 result = f.rep.quo_ground(coeff) 

1349 else: # pragma: no cover 

1350 raise OperationNotSupported(f, 'quo_ground') 

1351 

1352 return f.per(result) 

1353 

1354 def exquo_ground(f, coeff): 

1355 """ 

1356 Exact quotient of ``f`` by a an element of the ground domain. 

1357 

1358 Examples 

1359 ======== 

1360 

1361 >>> from sympy import Poly 

1362 >>> from sympy.abc import x 

1363 

1364 >>> Poly(2*x + 4).exquo_ground(2) 

1365 Poly(x + 2, x, domain='ZZ') 

1366 

1367 >>> Poly(2*x + 3).exquo_ground(2) 

1368 Traceback (most recent call last): 

1369 ... 

1370 ExactQuotientFailed: 2 does not divide 3 in ZZ 

1371 

1372 """ 

1373 if hasattr(f.rep, 'exquo_ground'): 

1374 result = f.rep.exquo_ground(coeff) 

1375 else: # pragma: no cover 

1376 raise OperationNotSupported(f, 'exquo_ground') 

1377 

1378 return f.per(result) 

1379 

1380 def abs(f): 

1381 """ 

1382 Make all coefficients in ``f`` positive. 

1383 

1384 Examples 

1385 ======== 

1386 

1387 >>> from sympy import Poly 

1388 >>> from sympy.abc import x 

1389 

1390 >>> Poly(x**2 - 1, x).abs() 

1391 Poly(x**2 + 1, x, domain='ZZ') 

1392 

1393 """ 

1394 if hasattr(f.rep, 'abs'): 

1395 result = f.rep.abs() 

1396 else: # pragma: no cover 

1397 raise OperationNotSupported(f, 'abs') 

1398 

1399 return f.per(result) 

1400 

1401 def neg(f): 

1402 """ 

1403 Negate all coefficients in ``f``. 

1404 

1405 Examples 

1406 ======== 

1407 

1408 >>> from sympy import Poly 

1409 >>> from sympy.abc import x 

1410 

1411 >>> Poly(x**2 - 1, x).neg() 

1412 Poly(-x**2 + 1, x, domain='ZZ') 

1413 

1414 >>> -Poly(x**2 - 1, x) 

1415 Poly(-x**2 + 1, x, domain='ZZ') 

1416 

1417 """ 

1418 if hasattr(f.rep, 'neg'): 

1419 result = f.rep.neg() 

1420 else: # pragma: no cover 

1421 raise OperationNotSupported(f, 'neg') 

1422 

1423 return f.per(result) 

1424 

1425 def add(f, g): 

1426 """ 

1427 Add two polynomials ``f`` and ``g``. 

1428 

1429 Examples 

1430 ======== 

1431 

1432 >>> from sympy import Poly 

1433 >>> from sympy.abc import x 

1434 

1435 >>> Poly(x**2 + 1, x).add(Poly(x - 2, x)) 

1436 Poly(x**2 + x - 1, x, domain='ZZ') 

1437 

1438 >>> Poly(x**2 + 1, x) + Poly(x - 2, x) 

1439 Poly(x**2 + x - 1, x, domain='ZZ') 

1440 

1441 """ 

1442 g = sympify(g) 

1443 

1444 if not g.is_Poly: 

1445 return f.add_ground(g) 

1446 

1447 _, per, F, G = f._unify(g) 

1448 

1449 if hasattr(f.rep, 'add'): 

1450 result = F.add(G) 

1451 else: # pragma: no cover 

1452 raise OperationNotSupported(f, 'add') 

1453 

1454 return per(result) 

1455 

1456 def sub(f, g): 

1457 """ 

1458 Subtract two polynomials ``f`` and ``g``. 

1459 

1460 Examples 

1461 ======== 

1462 

1463 >>> from sympy import Poly 

1464 >>> from sympy.abc import x 

1465 

1466 >>> Poly(x**2 + 1, x).sub(Poly(x - 2, x)) 

1467 Poly(x**2 - x + 3, x, domain='ZZ') 

1468 

1469 >>> Poly(x**2 + 1, x) - Poly(x - 2, x) 

1470 Poly(x**2 - x + 3, x, domain='ZZ') 

1471 

1472 """ 

1473 g = sympify(g) 

1474 

1475 if not g.is_Poly: 

1476 return f.sub_ground(g) 

1477 

1478 _, per, F, G = f._unify(g) 

1479 

1480 if hasattr(f.rep, 'sub'): 

1481 result = F.sub(G) 

1482 else: # pragma: no cover 

1483 raise OperationNotSupported(f, 'sub') 

1484 

1485 return per(result) 

1486 

1487 def mul(f, g): 

1488 """ 

1489 Multiply two polynomials ``f`` and ``g``. 

1490 

1491 Examples 

1492 ======== 

1493 

1494 >>> from sympy import Poly 

1495 >>> from sympy.abc import x 

1496 

1497 >>> Poly(x**2 + 1, x).mul(Poly(x - 2, x)) 

1498 Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ') 

1499 

1500 >>> Poly(x**2 + 1, x)*Poly(x - 2, x) 

1501 Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ') 

1502 

1503 """ 

1504 g = sympify(g) 

1505 

1506 if not g.is_Poly: 

1507 return f.mul_ground(g) 

1508 

1509 _, per, F, G = f._unify(g) 

1510 

1511 if hasattr(f.rep, 'mul'): 

1512 result = F.mul(G) 

1513 else: # pragma: no cover 

1514 raise OperationNotSupported(f, 'mul') 

1515 

1516 return per(result) 

1517 

1518 def sqr(f): 

1519 """ 

1520 Square a polynomial ``f``. 

1521 

1522 Examples 

1523 ======== 

1524 

1525 >>> from sympy import Poly 

1526 >>> from sympy.abc import x 

1527 

1528 >>> Poly(x - 2, x).sqr() 

1529 Poly(x**2 - 4*x + 4, x, domain='ZZ') 

1530 

1531 >>> Poly(x - 2, x)**2 

1532 Poly(x**2 - 4*x + 4, x, domain='ZZ') 

1533 

1534 """ 

1535 if hasattr(f.rep, 'sqr'): 

1536 result = f.rep.sqr() 

1537 else: # pragma: no cover 

1538 raise OperationNotSupported(f, 'sqr') 

1539 

1540 return f.per(result) 

1541 

1542 def pow(f, n): 

1543 """ 

1544 Raise ``f`` to a non-negative power ``n``. 

1545 

1546 Examples 

1547 ======== 

1548 

1549 >>> from sympy import Poly 

1550 >>> from sympy.abc import x 

1551 

1552 >>> Poly(x - 2, x).pow(3) 

1553 Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ') 

1554 

1555 >>> Poly(x - 2, x)**3 

1556 Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ') 

1557 

1558 """ 

1559 n = int(n) 

1560 

1561 if hasattr(f.rep, 'pow'): 

1562 result = f.rep.pow(n) 

1563 else: # pragma: no cover 

1564 raise OperationNotSupported(f, 'pow') 

1565 

1566 return f.per(result) 

1567 

1568 def pdiv(f, g): 

1569 """ 

1570 Polynomial pseudo-division of ``f`` by ``g``. 

1571 

1572 Examples 

1573 ======== 

1574 

1575 >>> from sympy import Poly 

1576 >>> from sympy.abc import x 

1577 

1578 >>> Poly(x**2 + 1, x).pdiv(Poly(2*x - 4, x)) 

1579 (Poly(2*x + 4, x, domain='ZZ'), Poly(20, x, domain='ZZ')) 

1580 

1581 """ 

1582 _, per, F, G = f._unify(g) 

1583 

1584 if hasattr(f.rep, 'pdiv'): 

1585 q, r = F.pdiv(G) 

1586 else: # pragma: no cover 

1587 raise OperationNotSupported(f, 'pdiv') 

1588 

1589 return per(q), per(r) 

1590 

1591 def prem(f, g): 

1592 """ 

1593 Polynomial pseudo-remainder of ``f`` by ``g``. 

1594 

1595 Caveat: The function prem(f, g, x) can be safely used to compute 

1596 in Z[x] _only_ subresultant polynomial remainder sequences (prs's). 

1597 

1598 To safely compute Euclidean and Sturmian prs's in Z[x] 

1599 employ anyone of the corresponding functions found in 

1600 the module sympy.polys.subresultants_qq_zz. The functions 

1601 in the module with suffix _pg compute prs's in Z[x] employing 

1602 rem(f, g, x), whereas the functions with suffix _amv 

1603 compute prs's in Z[x] employing rem_z(f, g, x). 

1604 

1605 The function rem_z(f, g, x) differs from prem(f, g, x) in that 

1606 to compute the remainder polynomials in Z[x] it premultiplies 

1607 the divident times the absolute value of the leading coefficient 

1608 of the divisor raised to the power degree(f, x) - degree(g, x) + 1. 

1609 

1610 

1611 Examples 

1612 ======== 

1613 

1614 >>> from sympy import Poly 

1615 >>> from sympy.abc import x 

1616 

1617 >>> Poly(x**2 + 1, x).prem(Poly(2*x - 4, x)) 

1618 Poly(20, x, domain='ZZ') 

1619 

1620 """ 

1621 _, per, F, G = f._unify(g) 

1622 

1623 if hasattr(f.rep, 'prem'): 

1624 result = F.prem(G) 

1625 else: # pragma: no cover 

1626 raise OperationNotSupported(f, 'prem') 

1627 

1628 return per(result) 

1629 

1630 def pquo(f, g): 

1631 """ 

1632 Polynomial pseudo-quotient of ``f`` by ``g``. 

1633 

1634 See the Caveat note in the function prem(f, g). 

1635 

1636 Examples 

1637 ======== 

1638 

1639 >>> from sympy import Poly 

1640 >>> from sympy.abc import x 

1641 

1642 >>> Poly(x**2 + 1, x).pquo(Poly(2*x - 4, x)) 

1643 Poly(2*x + 4, x, domain='ZZ') 

1644 

1645 >>> Poly(x**2 - 1, x).pquo(Poly(2*x - 2, x)) 

1646 Poly(2*x + 2, x, domain='ZZ') 

1647 

1648 """ 

1649 _, per, F, G = f._unify(g) 

1650 

1651 if hasattr(f.rep, 'pquo'): 

1652 result = F.pquo(G) 

1653 else: # pragma: no cover 

1654 raise OperationNotSupported(f, 'pquo') 

1655 

1656 return per(result) 

1657 

1658 def pexquo(f, g): 

1659 """ 

1660 Polynomial exact pseudo-quotient of ``f`` by ``g``. 

1661 

1662 Examples 

1663 ======== 

1664 

1665 >>> from sympy import Poly 

1666 >>> from sympy.abc import x 

1667 

1668 >>> Poly(x**2 - 1, x).pexquo(Poly(2*x - 2, x)) 

1669 Poly(2*x + 2, x, domain='ZZ') 

1670 

1671 >>> Poly(x**2 + 1, x).pexquo(Poly(2*x - 4, x)) 

1672 Traceback (most recent call last): 

1673 ... 

1674 ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1 

1675 

1676 """ 

1677 _, per, F, G = f._unify(g) 

1678 

1679 if hasattr(f.rep, 'pexquo'): 

1680 try: 

1681 result = F.pexquo(G) 

1682 except ExactQuotientFailed as exc: 

1683 raise exc.new(f.as_expr(), g.as_expr()) 

1684 else: # pragma: no cover 

1685 raise OperationNotSupported(f, 'pexquo') 

1686 

1687 return per(result) 

1688 

1689 def div(f, g, auto=True): 

1690 """ 

1691 Polynomial division with remainder of ``f`` by ``g``. 

1692 

1693 Examples 

1694 ======== 

1695 

1696 >>> from sympy import Poly 

1697 >>> from sympy.abc import x 

1698 

1699 >>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x)) 

1700 (Poly(1/2*x + 1, x, domain='QQ'), Poly(5, x, domain='QQ')) 

1701 

1702 >>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x), auto=False) 

1703 (Poly(0, x, domain='ZZ'), Poly(x**2 + 1, x, domain='ZZ')) 

1704 

1705 """ 

1706 dom, per, F, G = f._unify(g) 

1707 retract = False 

1708 

1709 if auto and dom.is_Ring and not dom.is_Field: 

1710 F, G = F.to_field(), G.to_field() 

1711 retract = True 

1712 

1713 if hasattr(f.rep, 'div'): 

1714 q, r = F.div(G) 

1715 else: # pragma: no cover 

1716 raise OperationNotSupported(f, 'div') 

1717 

1718 if retract: 

1719 try: 

1720 Q, R = q.to_ring(), r.to_ring() 

1721 except CoercionFailed: 

1722 pass 

1723 else: 

1724 q, r = Q, R 

1725 

1726 return per(q), per(r) 

1727 

1728 def rem(f, g, auto=True): 

1729 """ 

1730 Computes the polynomial remainder of ``f`` by ``g``. 

1731 

1732 Examples 

1733 ======== 

1734 

1735 >>> from sympy import Poly 

1736 >>> from sympy.abc import x 

1737 

1738 >>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x)) 

1739 Poly(5, x, domain='ZZ') 

1740 

1741 >>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x), auto=False) 

1742 Poly(x**2 + 1, x, domain='ZZ') 

1743 

1744 """ 

1745 dom, per, F, G = f._unify(g) 

1746 retract = False 

1747 

1748 if auto and dom.is_Ring and not dom.is_Field: 

1749 F, G = F.to_field(), G.to_field() 

1750 retract = True 

1751 

1752 if hasattr(f.rep, 'rem'): 

1753 r = F.rem(G) 

1754 else: # pragma: no cover 

1755 raise OperationNotSupported(f, 'rem') 

1756 

1757 if retract: 

1758 try: 

1759 r = r.to_ring() 

1760 except CoercionFailed: 

1761 pass 

1762 

1763 return per(r) 

1764 

1765 def quo(f, g, auto=True): 

1766 """ 

1767 Computes polynomial quotient of ``f`` by ``g``. 

1768 

1769 Examples 

1770 ======== 

1771 

1772 >>> from sympy import Poly 

1773 >>> from sympy.abc import x 

1774 

1775 >>> Poly(x**2 + 1, x).quo(Poly(2*x - 4, x)) 

1776 Poly(1/2*x + 1, x, domain='QQ') 

1777 

1778 >>> Poly(x**2 - 1, x).quo(Poly(x - 1, x)) 

1779 Poly(x + 1, x, domain='ZZ') 

1780 

1781 """ 

1782 dom, per, F, G = f._unify(g) 

1783 retract = False 

1784 

1785 if auto and dom.is_Ring and not dom.is_Field: 

1786 F, G = F.to_field(), G.to_field() 

1787 retract = True 

1788 

1789 if hasattr(f.rep, 'quo'): 

1790 q = F.quo(G) 

1791 else: # pragma: no cover 

1792 raise OperationNotSupported(f, 'quo') 

1793 

1794 if retract: 

1795 try: 

1796 q = q.to_ring() 

1797 except CoercionFailed: 

1798 pass 

1799 

1800 return per(q) 

1801 

1802 def exquo(f, g, auto=True): 

1803 """ 

1804 Computes polynomial exact quotient of ``f`` by ``g``. 

1805 

1806 Examples 

1807 ======== 

1808 

1809 >>> from sympy import Poly 

1810 >>> from sympy.abc import x 

1811 

1812 >>> Poly(x**2 - 1, x).exquo(Poly(x - 1, x)) 

1813 Poly(x + 1, x, domain='ZZ') 

1814 

1815 >>> Poly(x**2 + 1, x).exquo(Poly(2*x - 4, x)) 

1816 Traceback (most recent call last): 

1817 ... 

1818 ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1 

1819 

1820 """ 

1821 dom, per, F, G = f._unify(g) 

1822 retract = False 

1823 

1824 if auto and dom.is_Ring and not dom.is_Field: 

1825 F, G = F.to_field(), G.to_field() 

1826 retract = True 

1827 

1828 if hasattr(f.rep, 'exquo'): 

1829 try: 

1830 q = F.exquo(G) 

1831 except ExactQuotientFailed as exc: 

1832 raise exc.new(f.as_expr(), g.as_expr()) 

1833 else: # pragma: no cover 

1834 raise OperationNotSupported(f, 'exquo') 

1835 

1836 if retract: 

1837 try: 

1838 q = q.to_ring() 

1839 except CoercionFailed: 

1840 pass 

1841 

1842 return per(q) 

1843 

1844 def _gen_to_level(f, gen): 

1845 """Returns level associated with the given generator. """ 

1846 if isinstance(gen, int): 

1847 length = len(f.gens) 

1848 

1849 if -length <= gen < length: 

1850 if gen < 0: 

1851 return length + gen 

1852 else: 

1853 return gen 

1854 else: 

1855 raise PolynomialError("-%s <= gen < %s expected, got %s" % 

1856 (length, length, gen)) 

1857 else: 

1858 try: 

1859 return f.gens.index(sympify(gen)) 

1860 except ValueError: 

1861 raise PolynomialError( 

1862 "a valid generator expected, got %s" % gen) 

1863 

1864 def degree(f, gen=0): 

1865 """ 

1866 Returns degree of ``f`` in ``x_j``. 

1867 

1868 The degree of 0 is negative infinity. 

1869 

1870 Examples 

1871 ======== 

1872 

1873 >>> from sympy import Poly 

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

1875 

1876 >>> Poly(x**2 + y*x + 1, x, y).degree() 

1877 2 

1878 >>> Poly(x**2 + y*x + y, x, y).degree(y) 

1879 1 

1880 >>> Poly(0, x).degree() 

1881 -oo 

1882 

1883 """ 

1884 j = f._gen_to_level(gen) 

1885 

1886 if hasattr(f.rep, 'degree'): 

1887 return f.rep.degree(j) 

1888 else: # pragma: no cover 

1889 raise OperationNotSupported(f, 'degree') 

1890 

1891 def degree_list(f): 

1892 """ 

1893 Returns a list of degrees of ``f``. 

1894 

1895 Examples 

1896 ======== 

1897 

1898 >>> from sympy import Poly 

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

1900 

1901 >>> Poly(x**2 + y*x + 1, x, y).degree_list() 

1902 (2, 1) 

1903 

1904 """ 

1905 if hasattr(f.rep, 'degree_list'): 

1906 return f.rep.degree_list() 

1907 else: # pragma: no cover 

1908 raise OperationNotSupported(f, 'degree_list') 

1909 

1910 def total_degree(f): 

1911 """ 

1912 Returns the total degree of ``f``. 

1913 

1914 Examples 

1915 ======== 

1916 

1917 >>> from sympy import Poly 

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

1919 

1920 >>> Poly(x**2 + y*x + 1, x, y).total_degree() 

1921 2 

1922 >>> Poly(x + y**5, x, y).total_degree() 

1923 5 

1924 

1925 """ 

1926 if hasattr(f.rep, 'total_degree'): 

1927 return f.rep.total_degree() 

1928 else: # pragma: no cover 

1929 raise OperationNotSupported(f, 'total_degree') 

1930 

1931 def homogenize(f, s): 

1932 """ 

1933 Returns the homogeneous polynomial of ``f``. 

1934 

1935 A homogeneous polynomial is a polynomial whose all monomials with 

1936 non-zero coefficients have the same total degree. If you only 

1937 want to check if a polynomial is homogeneous, then use 

1938 :func:`Poly.is_homogeneous`. If you want not only to check if a 

1939 polynomial is homogeneous but also compute its homogeneous order, 

1940 then use :func:`Poly.homogeneous_order`. 

1941 

1942 Examples 

1943 ======== 

1944 

1945 >>> from sympy import Poly 

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

1947 

1948 >>> f = Poly(x**5 + 2*x**2*y**2 + 9*x*y**3) 

1949 >>> f.homogenize(z) 

1950 Poly(x**5 + 2*x**2*y**2*z + 9*x*y**3*z, x, y, z, domain='ZZ') 

1951 

1952 """ 

1953 if not isinstance(s, Symbol): 

1954 raise TypeError("``Symbol`` expected, got %s" % type(s)) 

1955 if s in f.gens: 

1956 i = f.gens.index(s) 

1957 gens = f.gens 

1958 else: 

1959 i = len(f.gens) 

1960 gens = f.gens + (s,) 

1961 if hasattr(f.rep, 'homogenize'): 

1962 return f.per(f.rep.homogenize(i), gens=gens) 

1963 raise OperationNotSupported(f, 'homogeneous_order') 

1964 

1965 def homogeneous_order(f): 

1966 """ 

1967 Returns the homogeneous order of ``f``. 

1968 

1969 A homogeneous polynomial is a polynomial whose all monomials with 

1970 non-zero coefficients have the same total degree. This degree is 

1971 the homogeneous order of ``f``. If you only want to check if a 

1972 polynomial is homogeneous, then use :func:`Poly.is_homogeneous`. 

1973 

1974 Examples 

1975 ======== 

1976 

1977 >>> from sympy import Poly 

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

1979 

1980 >>> f = Poly(x**5 + 2*x**3*y**2 + 9*x*y**4) 

1981 >>> f.homogeneous_order() 

1982 5 

1983 

1984 """ 

1985 if hasattr(f.rep, 'homogeneous_order'): 

1986 return f.rep.homogeneous_order() 

1987 else: # pragma: no cover 

1988 raise OperationNotSupported(f, 'homogeneous_order') 

1989 

1990 def LC(f, order=None): 

1991 """ 

1992 Returns the leading coefficient of ``f``. 

1993 

1994 Examples 

1995 ======== 

1996 

1997 >>> from sympy import Poly 

1998 >>> from sympy.abc import x 

1999 

2000 >>> Poly(4*x**3 + 2*x**2 + 3*x, x).LC() 

2001 4 

2002 

2003 """ 

2004 if order is not None: 

2005 return f.coeffs(order)[0] 

2006 

2007 if hasattr(f.rep, 'LC'): 

2008 result = f.rep.LC() 

2009 else: # pragma: no cover 

2010 raise OperationNotSupported(f, 'LC') 

2011 

2012 return f.rep.dom.to_sympy(result) 

2013 

2014 def TC(f): 

2015 """ 

2016 Returns the trailing coefficient of ``f``. 

2017 

2018 Examples 

2019 ======== 

2020 

2021 >>> from sympy import Poly 

2022 >>> from sympy.abc import x 

2023 

2024 >>> Poly(x**3 + 2*x**2 + 3*x, x).TC() 

2025 0 

2026 

2027 """ 

2028 if hasattr(f.rep, 'TC'): 

2029 result = f.rep.TC() 

2030 else: # pragma: no cover 

2031 raise OperationNotSupported(f, 'TC') 

2032 

2033 return f.rep.dom.to_sympy(result) 

2034 

2035 def EC(f, order=None): 

2036 """ 

2037 Returns the last non-zero coefficient of ``f``. 

2038 

2039 Examples 

2040 ======== 

2041 

2042 >>> from sympy import Poly 

2043 >>> from sympy.abc import x 

2044 

2045 >>> Poly(x**3 + 2*x**2 + 3*x, x).EC() 

2046 3 

2047 

2048 """ 

2049 if hasattr(f.rep, 'coeffs'): 

2050 return f.coeffs(order)[-1] 

2051 else: # pragma: no cover 

2052 raise OperationNotSupported(f, 'EC') 

2053 

2054 def coeff_monomial(f, monom): 

2055 """ 

2056 Returns the coefficient of ``monom`` in ``f`` if there, else None. 

2057 

2058 Examples 

2059 ======== 

2060 

2061 >>> from sympy import Poly, exp 

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

2063 

2064 >>> p = Poly(24*x*y*exp(8) + 23*x, x, y) 

2065 

2066 >>> p.coeff_monomial(x) 

2067 23 

2068 >>> p.coeff_monomial(y) 

2069 0 

2070 >>> p.coeff_monomial(x*y) 

2071 24*exp(8) 

2072 

2073 Note that ``Expr.coeff()`` behaves differently, collecting terms 

2074 if possible; the Poly must be converted to an Expr to use that 

2075 method, however: 

2076 

2077 >>> p.as_expr().coeff(x) 

2078 24*y*exp(8) + 23 

2079 >>> p.as_expr().coeff(y) 

2080 24*x*exp(8) 

2081 >>> p.as_expr().coeff(x*y) 

2082 24*exp(8) 

2083 

2084 See Also 

2085 ======== 

2086 nth: more efficient query using exponents of the monomial's generators 

2087 

2088 """ 

2089 return f.nth(*Monomial(monom, f.gens).exponents) 

2090 

2091 def nth(f, *N): 

2092 """ 

2093 Returns the ``n``-th coefficient of ``f`` where ``N`` are the 

2094 exponents of the generators in the term of interest. 

2095 

2096 Examples 

2097 ======== 

2098 

2099 >>> from sympy import Poly, sqrt 

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

2101 

2102 >>> Poly(x**3 + 2*x**2 + 3*x, x).nth(2) 

2103 2 

2104 >>> Poly(x**3 + 2*x*y**2 + y**2, x, y).nth(1, 2) 

2105 2 

2106 >>> Poly(4*sqrt(x)*y) 

2107 Poly(4*y*(sqrt(x)), y, sqrt(x), domain='ZZ') 

2108 >>> _.nth(1, 1) 

2109 4 

2110 

2111 See Also 

2112 ======== 

2113 coeff_monomial 

2114 

2115 """ 

2116 if hasattr(f.rep, 'nth'): 

2117 if len(N) != len(f.gens): 

2118 raise ValueError('exponent of each generator must be specified') 

2119 result = f.rep.nth(*list(map(int, N))) 

2120 else: # pragma: no cover 

2121 raise OperationNotSupported(f, 'nth') 

2122 

2123 return f.rep.dom.to_sympy(result) 

2124 

2125 def coeff(f, x, n=1, right=False): 

2126 # the semantics of coeff_monomial and Expr.coeff are different; 

2127 # if someone is working with a Poly, they should be aware of the 

2128 # differences and chose the method best suited for the query. 

2129 # Alternatively, a pure-polys method could be written here but 

2130 # at this time the ``right`` keyword would be ignored because Poly 

2131 # doesn't work with non-commutatives. 

2132 raise NotImplementedError( 

2133 'Either convert to Expr with `as_expr` method ' 

2134 'to use Expr\'s coeff method or else use the ' 

2135 '`coeff_monomial` method of Polys.') 

2136 

2137 def LM(f, order=None): 

2138 """ 

2139 Returns the leading monomial of ``f``. 

2140 

2141 The Leading monomial signifies the monomial having 

2142 the highest power of the principal generator in the 

2143 expression f. 

2144 

2145 Examples 

2146 ======== 

2147 

2148 >>> from sympy import Poly 

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

2150 

2151 >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LM() 

2152 x**2*y**0 

2153 

2154 """ 

2155 return Monomial(f.monoms(order)[0], f.gens) 

2156 

2157 def EM(f, order=None): 

2158 """ 

2159 Returns the last non-zero monomial of ``f``. 

2160 

2161 Examples 

2162 ======== 

2163 

2164 >>> from sympy import Poly 

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

2166 

2167 >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).EM() 

2168 x**0*y**1 

2169 

2170 """ 

2171 return Monomial(f.monoms(order)[-1], f.gens) 

2172 

2173 def LT(f, order=None): 

2174 """ 

2175 Returns the leading term of ``f``. 

2176 

2177 The Leading term signifies the term having 

2178 the highest power of the principal generator in the 

2179 expression f along with its coefficient. 

2180 

2181 Examples 

2182 ======== 

2183 

2184 >>> from sympy import Poly 

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

2186 

2187 >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LT() 

2188 (x**2*y**0, 4) 

2189 

2190 """ 

2191 monom, coeff = f.terms(order)[0] 

2192 return Monomial(monom, f.gens), coeff 

2193 

2194 def ET(f, order=None): 

2195 """ 

2196 Returns the last non-zero term of ``f``. 

2197 

2198 Examples 

2199 ======== 

2200 

2201 >>> from sympy import Poly 

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

2203 

2204 >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).ET() 

2205 (x**0*y**1, 3) 

2206 

2207 """ 

2208 monom, coeff = f.terms(order)[-1] 

2209 return Monomial(monom, f.gens), coeff 

2210 

2211 def max_norm(f): 

2212 """ 

2213 Returns maximum norm of ``f``. 

2214 

2215 Examples 

2216 ======== 

2217 

2218 >>> from sympy import Poly 

2219 >>> from sympy.abc import x 

2220 

2221 >>> Poly(-x**2 + 2*x - 3, x).max_norm() 

2222 3 

2223 

2224 """ 

2225 if hasattr(f.rep, 'max_norm'): 

2226 result = f.rep.max_norm() 

2227 else: # pragma: no cover 

2228 raise OperationNotSupported(f, 'max_norm') 

2229 

2230 return f.rep.dom.to_sympy(result) 

2231 

2232 def l1_norm(f): 

2233 """ 

2234 Returns l1 norm of ``f``. 

2235 

2236 Examples 

2237 ======== 

2238 

2239 >>> from sympy import Poly 

2240 >>> from sympy.abc import x 

2241 

2242 >>> Poly(-x**2 + 2*x - 3, x).l1_norm() 

2243 6 

2244 

2245 """ 

2246 if hasattr(f.rep, 'l1_norm'): 

2247 result = f.rep.l1_norm() 

2248 else: # pragma: no cover 

2249 raise OperationNotSupported(f, 'l1_norm') 

2250 

2251 return f.rep.dom.to_sympy(result) 

2252 

2253 def clear_denoms(self, convert=False): 

2254 """ 

2255 Clear denominators, but keep the ground domain. 

2256 

2257 Examples 

2258 ======== 

2259 

2260 >>> from sympy import Poly, S, QQ 

2261 >>> from sympy.abc import x 

2262 

2263 >>> f = Poly(x/2 + S(1)/3, x, domain=QQ) 

2264 

2265 >>> f.clear_denoms() 

2266 (6, Poly(3*x + 2, x, domain='QQ')) 

2267 >>> f.clear_denoms(convert=True) 

2268 (6, Poly(3*x + 2, x, domain='ZZ')) 

2269 

2270 """ 

2271 f = self 

2272 

2273 if not f.rep.dom.is_Field: 

2274 return S.One, f 

2275 

2276 dom = f.get_domain() 

2277 if dom.has_assoc_Ring: 

2278 dom = f.rep.dom.get_ring() 

2279 

2280 if hasattr(f.rep, 'clear_denoms'): 

2281 coeff, result = f.rep.clear_denoms() 

2282 else: # pragma: no cover 

2283 raise OperationNotSupported(f, 'clear_denoms') 

2284 

2285 coeff, f = dom.to_sympy(coeff), f.per(result) 

2286 

2287 if not convert or not dom.has_assoc_Ring: 

2288 return coeff, f 

2289 else: 

2290 return coeff, f.to_ring() 

2291 

2292 def rat_clear_denoms(self, g): 

2293 """ 

2294 Clear denominators in a rational function ``f/g``. 

2295 

2296 Examples 

2297 ======== 

2298 

2299 >>> from sympy import Poly 

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

2301 

2302 >>> f = Poly(x**2/y + 1, x) 

2303 >>> g = Poly(x**3 + y, x) 

2304 

2305 >>> p, q = f.rat_clear_denoms(g) 

2306 

2307 >>> p 

2308 Poly(x**2 + y, x, domain='ZZ[y]') 

2309 >>> q 

2310 Poly(y*x**3 + y**2, x, domain='ZZ[y]') 

2311 

2312 """ 

2313 f = self 

2314 

2315 dom, per, f, g = f._unify(g) 

2316 

2317 f = per(f) 

2318 g = per(g) 

2319 

2320 if not (dom.is_Field and dom.has_assoc_Ring): 

2321 return f, g 

2322 

2323 a, f = f.clear_denoms(convert=True) 

2324 b, g = g.clear_denoms(convert=True) 

2325 

2326 f = f.mul_ground(b) 

2327 g = g.mul_ground(a) 

2328 

2329 return f, g 

2330 

2331 def integrate(self, *specs, **args): 

2332 """ 

2333 Computes indefinite integral of ``f``. 

2334 

2335 Examples 

2336 ======== 

2337 

2338 >>> from sympy import Poly 

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

2340 

2341 >>> Poly(x**2 + 2*x + 1, x).integrate() 

2342 Poly(1/3*x**3 + x**2 + x, x, domain='QQ') 

2343 

2344 >>> Poly(x*y**2 + x, x, y).integrate((0, 1), (1, 0)) 

2345 Poly(1/2*x**2*y**2 + 1/2*x**2, x, y, domain='QQ') 

2346 

2347 """ 

2348 f = self 

2349 

2350 if args.get('auto', True) and f.rep.dom.is_Ring: 

2351 f = f.to_field() 

2352 

2353 if hasattr(f.rep, 'integrate'): 

2354 if not specs: 

2355 return f.per(f.rep.integrate(m=1)) 

2356 

2357 rep = f.rep 

2358 

2359 for spec in specs: 

2360 if isinstance(spec, tuple): 

2361 gen, m = spec 

2362 else: 

2363 gen, m = spec, 1 

2364 

2365 rep = rep.integrate(int(m), f._gen_to_level(gen)) 

2366 

2367 return f.per(rep) 

2368 else: # pragma: no cover 

2369 raise OperationNotSupported(f, 'integrate') 

2370 

2371 def diff(f, *specs, **kwargs): 

2372 """ 

2373 Computes partial derivative of ``f``. 

2374 

2375 Examples 

2376 ======== 

2377 

2378 >>> from sympy import Poly 

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

2380 

2381 >>> Poly(x**2 + 2*x + 1, x).diff() 

2382 Poly(2*x + 2, x, domain='ZZ') 

2383 

2384 >>> Poly(x*y**2 + x, x, y).diff((0, 0), (1, 1)) 

2385 Poly(2*x*y, x, y, domain='ZZ') 

2386 

2387 """ 

2388 if not kwargs.get('evaluate', True): 

2389 return Derivative(f, *specs, **kwargs) 

2390 

2391 if hasattr(f.rep, 'diff'): 

2392 if not specs: 

2393 return f.per(f.rep.diff(m=1)) 

2394 

2395 rep = f.rep 

2396 

2397 for spec in specs: 

2398 if isinstance(spec, tuple): 

2399 gen, m = spec 

2400 else: 

2401 gen, m = spec, 1 

2402 

2403 rep = rep.diff(int(m), f._gen_to_level(gen)) 

2404 

2405 return f.per(rep) 

2406 else: # pragma: no cover 

2407 raise OperationNotSupported(f, 'diff') 

2408 

2409 _eval_derivative = diff 

2410 

2411 def eval(self, x, a=None, auto=True): 

2412 """ 

2413 Evaluate ``f`` at ``a`` in the given variable. 

2414 

2415 Examples 

2416 ======== 

2417 

2418 >>> from sympy import Poly 

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

2420 

2421 >>> Poly(x**2 + 2*x + 3, x).eval(2) 

2422 11 

2423 

2424 >>> Poly(2*x*y + 3*x + y + 2, x, y).eval(x, 2) 

2425 Poly(5*y + 8, y, domain='ZZ') 

2426 

2427 >>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z) 

2428 

2429 >>> f.eval({x: 2}) 

2430 Poly(5*y + 2*z + 6, y, z, domain='ZZ') 

2431 >>> f.eval({x: 2, y: 5}) 

2432 Poly(2*z + 31, z, domain='ZZ') 

2433 >>> f.eval({x: 2, y: 5, z: 7}) 

2434 45 

2435 

2436 >>> f.eval((2, 5)) 

2437 Poly(2*z + 31, z, domain='ZZ') 

2438 >>> f(2, 5) 

2439 Poly(2*z + 31, z, domain='ZZ') 

2440 

2441 """ 

2442 f = self 

2443 

2444 if a is None: 

2445 if isinstance(x, dict): 

2446 mapping = x 

2447 

2448 for gen, value in mapping.items(): 

2449 f = f.eval(gen, value) 

2450 

2451 return f 

2452 elif isinstance(x, (tuple, list)): 

2453 values = x 

2454 

2455 if len(values) > len(f.gens): 

2456 raise ValueError("too many values provided") 

2457 

2458 for gen, value in zip(f.gens, values): 

2459 f = f.eval(gen, value) 

2460 

2461 return f 

2462 else: 

2463 j, a = 0, x 

2464 else: 

2465 j = f._gen_to_level(x) 

2466 

2467 if not hasattr(f.rep, 'eval'): # pragma: no cover 

2468 raise OperationNotSupported(f, 'eval') 

2469 

2470 try: 

2471 result = f.rep.eval(a, j) 

2472 except CoercionFailed: 

2473 if not auto: 

2474 raise DomainError("Cannot evaluate at %s in %s" % (a, f.rep.dom)) 

2475 else: 

2476 a_domain, [a] = construct_domain([a]) 

2477 new_domain = f.get_domain().unify_with_symbols(a_domain, f.gens) 

2478 

2479 f = f.set_domain(new_domain) 

2480 a = new_domain.convert(a, a_domain) 

2481 

2482 result = f.rep.eval(a, j) 

2483 

2484 return f.per(result, remove=j) 

2485 

2486 def __call__(f, *values): 

2487 """ 

2488 Evaluate ``f`` at the give values. 

2489 

2490 Examples 

2491 ======== 

2492 

2493 >>> from sympy import Poly 

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

2495 

2496 >>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z) 

2497 

2498 >>> f(2) 

2499 Poly(5*y + 2*z + 6, y, z, domain='ZZ') 

2500 >>> f(2, 5) 

2501 Poly(2*z + 31, z, domain='ZZ') 

2502 >>> f(2, 5, 7) 

2503 45 

2504 

2505 """ 

2506 return f.eval(values) 

2507 

2508 def half_gcdex(f, g, auto=True): 

2509 """ 

2510 Half extended Euclidean algorithm of ``f`` and ``g``. 

2511 

2512 Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``. 

2513 

2514 Examples 

2515 ======== 

2516 

2517 >>> from sympy import Poly 

2518 >>> from sympy.abc import x 

2519 

2520 >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15 

2521 >>> g = x**3 + x**2 - 4*x - 4 

2522 

2523 >>> Poly(f).half_gcdex(Poly(g)) 

2524 (Poly(-1/5*x + 3/5, x, domain='QQ'), Poly(x + 1, x, domain='QQ')) 

2525 

2526 """ 

2527 dom, per, F, G = f._unify(g) 

2528 

2529 if auto and dom.is_Ring: 

2530 F, G = F.to_field(), G.to_field() 

2531 

2532 if hasattr(f.rep, 'half_gcdex'): 

2533 s, h = F.half_gcdex(G) 

2534 else: # pragma: no cover 

2535 raise OperationNotSupported(f, 'half_gcdex') 

2536 

2537 return per(s), per(h) 

2538 

2539 def gcdex(f, g, auto=True): 

2540 """ 

2541 Extended Euclidean algorithm of ``f`` and ``g``. 

2542 

2543 Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``. 

2544 

2545 Examples 

2546 ======== 

2547 

2548 >>> from sympy import Poly 

2549 >>> from sympy.abc import x 

2550 

2551 >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15 

2552 >>> g = x**3 + x**2 - 4*x - 4 

2553 

2554 >>> Poly(f).gcdex(Poly(g)) 

2555 (Poly(-1/5*x + 3/5, x, domain='QQ'), 

2556 Poly(1/5*x**2 - 6/5*x + 2, x, domain='QQ'), 

2557 Poly(x + 1, x, domain='QQ')) 

2558 

2559 """ 

2560 dom, per, F, G = f._unify(g) 

2561 

2562 if auto and dom.is_Ring: 

2563 F, G = F.to_field(), G.to_field() 

2564 

2565 if hasattr(f.rep, 'gcdex'): 

2566 s, t, h = F.gcdex(G) 

2567 else: # pragma: no cover 

2568 raise OperationNotSupported(f, 'gcdex') 

2569 

2570 return per(s), per(t), per(h) 

2571 

2572 def invert(f, g, auto=True): 

2573 """ 

2574 Invert ``f`` modulo ``g`` when possible. 

2575 

2576 Examples 

2577 ======== 

2578 

2579 >>> from sympy import Poly 

2580 >>> from sympy.abc import x 

2581 

2582 >>> Poly(x**2 - 1, x).invert(Poly(2*x - 1, x)) 

2583 Poly(-4/3, x, domain='QQ') 

2584 

2585 >>> Poly(x**2 - 1, x).invert(Poly(x - 1, x)) 

2586 Traceback (most recent call last): 

2587 ... 

2588 NotInvertible: zero divisor 

2589 

2590 """ 

2591 dom, per, F, G = f._unify(g) 

2592 

2593 if auto and dom.is_Ring: 

2594 F, G = F.to_field(), G.to_field() 

2595 

2596 if hasattr(f.rep, 'invert'): 

2597 result = F.invert(G) 

2598 else: # pragma: no cover 

2599 raise OperationNotSupported(f, 'invert') 

2600 

2601 return per(result) 

2602 

2603 def revert(f, n): 

2604 """ 

2605 Compute ``f**(-1)`` mod ``x**n``. 

2606 

2607 Examples 

2608 ======== 

2609 

2610 >>> from sympy import Poly 

2611 >>> from sympy.abc import x 

2612 

2613 >>> Poly(1, x).revert(2) 

2614 Poly(1, x, domain='ZZ') 

2615 

2616 >>> Poly(1 + x, x).revert(1) 

2617 Poly(1, x, domain='ZZ') 

2618 

2619 >>> Poly(x**2 - 2, x).revert(2) 

2620 Traceback (most recent call last): 

2621 ... 

2622 NotReversible: only units are reversible in a ring 

2623 

2624 >>> Poly(1/x, x).revert(1) 

2625 Traceback (most recent call last): 

2626 ... 

2627 PolynomialError: 1/x contains an element of the generators set 

2628 

2629 """ 

2630 if hasattr(f.rep, 'revert'): 

2631 result = f.rep.revert(int(n)) 

2632 else: # pragma: no cover 

2633 raise OperationNotSupported(f, 'revert') 

2634 

2635 return f.per(result) 

2636 

2637 def subresultants(f, g): 

2638 """ 

2639 Computes the subresultant PRS of ``f`` and ``g``. 

2640 

2641 Examples 

2642 ======== 

2643 

2644 >>> from sympy import Poly 

2645 >>> from sympy.abc import x 

2646 

2647 >>> Poly(x**2 + 1, x).subresultants(Poly(x**2 - 1, x)) 

2648 [Poly(x**2 + 1, x, domain='ZZ'), 

2649 Poly(x**2 - 1, x, domain='ZZ'), 

2650 Poly(-2, x, domain='ZZ')] 

2651 

2652 """ 

2653 _, per, F, G = f._unify(g) 

2654 

2655 if hasattr(f.rep, 'subresultants'): 

2656 result = F.subresultants(G) 

2657 else: # pragma: no cover 

2658 raise OperationNotSupported(f, 'subresultants') 

2659 

2660 return list(map(per, result)) 

2661 

2662 def resultant(f, g, includePRS=False): 

2663 """ 

2664 Computes the resultant of ``f`` and ``g`` via PRS. 

2665 

2666 If includePRS=True, it includes the subresultant PRS in the result. 

2667 Because the PRS is used to calculate the resultant, this is more 

2668 efficient than calling :func:`subresultants` separately. 

2669 

2670 Examples 

2671 ======== 

2672 

2673 >>> from sympy import Poly 

2674 >>> from sympy.abc import x 

2675 

2676 >>> f = Poly(x**2 + 1, x) 

2677 

2678 >>> f.resultant(Poly(x**2 - 1, x)) 

2679 4 

2680 >>> f.resultant(Poly(x**2 - 1, x), includePRS=True) 

2681 (4, [Poly(x**2 + 1, x, domain='ZZ'), Poly(x**2 - 1, x, domain='ZZ'), 

2682 Poly(-2, x, domain='ZZ')]) 

2683 

2684 """ 

2685 _, per, F, G = f._unify(g) 

2686 

2687 if hasattr(f.rep, 'resultant'): 

2688 if includePRS: 

2689 result, R = F.resultant(G, includePRS=includePRS) 

2690 else: 

2691 result = F.resultant(G) 

2692 else: # pragma: no cover 

2693 raise OperationNotSupported(f, 'resultant') 

2694 

2695 if includePRS: 

2696 return (per(result, remove=0), list(map(per, R))) 

2697 return per(result, remove=0) 

2698 

2699 def discriminant(f): 

2700 """ 

2701 Computes the discriminant of ``f``. 

2702 

2703 Examples 

2704 ======== 

2705 

2706 >>> from sympy import Poly 

2707 >>> from sympy.abc import x 

2708 

2709 >>> Poly(x**2 + 2*x + 3, x).discriminant() 

2710 -8 

2711 

2712 """ 

2713 if hasattr(f.rep, 'discriminant'): 

2714 result = f.rep.discriminant() 

2715 else: # pragma: no cover 

2716 raise OperationNotSupported(f, 'discriminant') 

2717 

2718 return f.per(result, remove=0) 

2719 

2720 def dispersionset(f, g=None): 

2721 r"""Compute the *dispersion set* of two polynomials. 

2722 

2723 For two polynomials `f(x)` and `g(x)` with `\deg f > 0` 

2724 and `\deg g > 0` the dispersion set `\operatorname{J}(f, g)` is defined as: 

2725 

2726 .. math:: 

2727 \operatorname{J}(f, g) 

2728 & := \{a \in \mathbb{N}_0 | \gcd(f(x), g(x+a)) \neq 1\} \\ 

2729 & = \{a \in \mathbb{N}_0 | \deg \gcd(f(x), g(x+a)) \geq 1\} 

2730 

2731 For a single polynomial one defines `\operatorname{J}(f) := \operatorname{J}(f, f)`. 

2732 

2733 Examples 

2734 ======== 

2735 

2736 >>> from sympy import poly 

2737 >>> from sympy.polys.dispersion import dispersion, dispersionset 

2738 >>> from sympy.abc import x 

2739 

2740 Dispersion set and dispersion of a simple polynomial: 

2741 

2742 >>> fp = poly((x - 3)*(x + 3), x) 

2743 >>> sorted(dispersionset(fp)) 

2744 [0, 6] 

2745 >>> dispersion(fp) 

2746 6 

2747 

2748 Note that the definition of the dispersion is not symmetric: 

2749 

2750 >>> fp = poly(x**4 - 3*x**2 + 1, x) 

2751 >>> gp = fp.shift(-3) 

2752 >>> sorted(dispersionset(fp, gp)) 

2753 [2, 3, 4] 

2754 >>> dispersion(fp, gp) 

2755 4 

2756 >>> sorted(dispersionset(gp, fp)) 

2757 [] 

2758 >>> dispersion(gp, fp) 

2759 -oo 

2760 

2761 Computing the dispersion also works over field extensions: 

2762 

2763 >>> from sympy import sqrt 

2764 >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>') 

2765 >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>') 

2766 >>> sorted(dispersionset(fp, gp)) 

2767 [2] 

2768 >>> sorted(dispersionset(gp, fp)) 

2769 [1, 4] 

2770 

2771 We can even perform the computations for polynomials 

2772 having symbolic coefficients: 

2773 

2774 >>> from sympy.abc import a 

2775 >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x) 

2776 >>> sorted(dispersionset(fp)) 

2777 [0, 1] 

2778 

2779 See Also 

2780 ======== 

2781 

2782 dispersion 

2783 

2784 References 

2785 ========== 

2786 

2787 1. [ManWright94]_ 

2788 2. [Koepf98]_ 

2789 3. [Abramov71]_ 

2790 4. [Man93]_ 

2791 """ 

2792 from sympy.polys.dispersion import dispersionset 

2793 return dispersionset(f, g) 

2794 

2795 def dispersion(f, g=None): 

2796 r"""Compute the *dispersion* of polynomials. 

2797 

2798 For two polynomials `f(x)` and `g(x)` with `\deg f > 0` 

2799 and `\deg g > 0` the dispersion `\operatorname{dis}(f, g)` is defined as: 

2800 

2801 .. math:: 

2802 \operatorname{dis}(f, g) 

2803 & := \max\{ J(f,g) \cup \{0\} \} \\ 

2804 & = \max\{ \{a \in \mathbb{N} | \gcd(f(x), g(x+a)) \neq 1\} \cup \{0\} \} 

2805 

2806 and for a single polynomial `\operatorname{dis}(f) := \operatorname{dis}(f, f)`. 

2807 

2808 Examples 

2809 ======== 

2810 

2811 >>> from sympy import poly 

2812 >>> from sympy.polys.dispersion import dispersion, dispersionset 

2813 >>> from sympy.abc import x 

2814 

2815 Dispersion set and dispersion of a simple polynomial: 

2816 

2817 >>> fp = poly((x - 3)*(x + 3), x) 

2818 >>> sorted(dispersionset(fp)) 

2819 [0, 6] 

2820 >>> dispersion(fp) 

2821 6 

2822 

2823 Note that the definition of the dispersion is not symmetric: 

2824 

2825 >>> fp = poly(x**4 - 3*x**2 + 1, x) 

2826 >>> gp = fp.shift(-3) 

2827 >>> sorted(dispersionset(fp, gp)) 

2828 [2, 3, 4] 

2829 >>> dispersion(fp, gp) 

2830 4 

2831 >>> sorted(dispersionset(gp, fp)) 

2832 [] 

2833 >>> dispersion(gp, fp) 

2834 -oo 

2835 

2836 Computing the dispersion also works over field extensions: 

2837 

2838 >>> from sympy import sqrt 

2839 >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>') 

2840 >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>') 

2841 >>> sorted(dispersionset(fp, gp)) 

2842 [2] 

2843 >>> sorted(dispersionset(gp, fp)) 

2844 [1, 4] 

2845 

2846 We can even perform the computations for polynomials 

2847 having symbolic coefficients: 

2848 

2849 >>> from sympy.abc import a 

2850 >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x) 

2851 >>> sorted(dispersionset(fp)) 

2852 [0, 1] 

2853 

2854 See Also 

2855 ======== 

2856 

2857 dispersionset 

2858 

2859 References 

2860 ========== 

2861 

2862 1. [ManWright94]_ 

2863 2. [Koepf98]_ 

2864 3. [Abramov71]_ 

2865 4. [Man93]_ 

2866 """ 

2867 from sympy.polys.dispersion import dispersion 

2868 return dispersion(f, g) 

2869 

2870 def cofactors(f, g): 

2871 """ 

2872 Returns the GCD of ``f`` and ``g`` and their cofactors. 

2873 

2874 Returns polynomials ``(h, cff, cfg)`` such that ``h = gcd(f, g)``, and 

2875 ``cff = quo(f, h)`` and ``cfg = quo(g, h)`` are, so called, cofactors 

2876 of ``f`` and ``g``. 

2877 

2878 Examples 

2879 ======== 

2880 

2881 >>> from sympy import Poly 

2882 >>> from sympy.abc import x 

2883 

2884 >>> Poly(x**2 - 1, x).cofactors(Poly(x**2 - 3*x + 2, x)) 

2885 (Poly(x - 1, x, domain='ZZ'), 

2886 Poly(x + 1, x, domain='ZZ'), 

2887 Poly(x - 2, x, domain='ZZ')) 

2888 

2889 """ 

2890 _, per, F, G = f._unify(g) 

2891 

2892 if hasattr(f.rep, 'cofactors'): 

2893 h, cff, cfg = F.cofactors(G) 

2894 else: # pragma: no cover 

2895 raise OperationNotSupported(f, 'cofactors') 

2896 

2897 return per(h), per(cff), per(cfg) 

2898 

2899 def gcd(f, g): 

2900 """ 

2901 Returns the polynomial GCD of ``f`` and ``g``. 

2902 

2903 Examples 

2904 ======== 

2905 

2906 >>> from sympy import Poly 

2907 >>> from sympy.abc import x 

2908 

2909 >>> Poly(x**2 - 1, x).gcd(Poly(x**2 - 3*x + 2, x)) 

2910 Poly(x - 1, x, domain='ZZ') 

2911 

2912 """ 

2913 _, per, F, G = f._unify(g) 

2914 

2915 if hasattr(f.rep, 'gcd'): 

2916 result = F.gcd(G) 

2917 else: # pragma: no cover 

2918 raise OperationNotSupported(f, 'gcd') 

2919 

2920 return per(result) 

2921 

2922 def lcm(f, g): 

2923 """ 

2924 Returns polynomial LCM of ``f`` and ``g``. 

2925 

2926 Examples 

2927 ======== 

2928 

2929 >>> from sympy import Poly 

2930 >>> from sympy.abc import x 

2931 

2932 >>> Poly(x**2 - 1, x).lcm(Poly(x**2 - 3*x + 2, x)) 

2933 Poly(x**3 - 2*x**2 - x + 2, x, domain='ZZ') 

2934 

2935 """ 

2936 _, per, F, G = f._unify(g) 

2937 

2938 if hasattr(f.rep, 'lcm'): 

2939 result = F.lcm(G) 

2940 else: # pragma: no cover 

2941 raise OperationNotSupported(f, 'lcm') 

2942 

2943 return per(result) 

2944 

2945 def trunc(f, p): 

2946 """ 

2947 Reduce ``f`` modulo a constant ``p``. 

2948 

2949 Examples 

2950 ======== 

2951 

2952 >>> from sympy import Poly 

2953 >>> from sympy.abc import x 

2954 

2955 >>> Poly(2*x**3 + 3*x**2 + 5*x + 7, x).trunc(3) 

2956 Poly(-x**3 - x + 1, x, domain='ZZ') 

2957 

2958 """ 

2959 p = f.rep.dom.convert(p) 

2960 

2961 if hasattr(f.rep, 'trunc'): 

2962 result = f.rep.trunc(p) 

2963 else: # pragma: no cover 

2964 raise OperationNotSupported(f, 'trunc') 

2965 

2966 return f.per(result) 

2967 

2968 def monic(self, auto=True): 

2969 """ 

2970 Divides all coefficients by ``LC(f)``. 

2971 

2972 Examples 

2973 ======== 

2974 

2975 >>> from sympy import Poly, ZZ 

2976 >>> from sympy.abc import x 

2977 

2978 >>> Poly(3*x**2 + 6*x + 9, x, domain=ZZ).monic() 

2979 Poly(x**2 + 2*x + 3, x, domain='QQ') 

2980 

2981 >>> Poly(3*x**2 + 4*x + 2, x, domain=ZZ).monic() 

2982 Poly(x**2 + 4/3*x + 2/3, x, domain='QQ') 

2983 

2984 """ 

2985 f = self 

2986 

2987 if auto and f.rep.dom.is_Ring: 

2988 f = f.to_field() 

2989 

2990 if hasattr(f.rep, 'monic'): 

2991 result = f.rep.monic() 

2992 else: # pragma: no cover 

2993 raise OperationNotSupported(f, 'monic') 

2994 

2995 return f.per(result) 

2996 

2997 def content(f): 

2998 """ 

2999 Returns the GCD of polynomial coefficients. 

3000 

3001 Examples 

3002 ======== 

3003 

3004 >>> from sympy import Poly 

3005 >>> from sympy.abc import x 

3006 

3007 >>> Poly(6*x**2 + 8*x + 12, x).content() 

3008 2 

3009 

3010 """ 

3011 if hasattr(f.rep, 'content'): 

3012 result = f.rep.content() 

3013 else: # pragma: no cover 

3014 raise OperationNotSupported(f, 'content') 

3015 

3016 return f.rep.dom.to_sympy(result) 

3017 

3018 def primitive(f): 

3019 """ 

3020 Returns the content and a primitive form of ``f``. 

3021 

3022 Examples 

3023 ======== 

3024 

3025 >>> from sympy import Poly 

3026 >>> from sympy.abc import x 

3027 

3028 >>> Poly(2*x**2 + 8*x + 12, x).primitive() 

3029 (2, Poly(x**2 + 4*x + 6, x, domain='ZZ')) 

3030 

3031 """ 

3032 if hasattr(f.rep, 'primitive'): 

3033 cont, result = f.rep.primitive() 

3034 else: # pragma: no cover 

3035 raise OperationNotSupported(f, 'primitive') 

3036 

3037 return f.rep.dom.to_sympy(cont), f.per(result) 

3038 

3039 def compose(f, g): 

3040 """ 

3041 Computes the functional composition of ``f`` and ``g``. 

3042 

3043 Examples 

3044 ======== 

3045 

3046 >>> from sympy import Poly 

3047 >>> from sympy.abc import x 

3048 

3049 >>> Poly(x**2 + x, x).compose(Poly(x - 1, x)) 

3050 Poly(x**2 - x, x, domain='ZZ') 

3051 

3052 """ 

3053 _, per, F, G = f._unify(g) 

3054 

3055 if hasattr(f.rep, 'compose'): 

3056 result = F.compose(G) 

3057 else: # pragma: no cover 

3058 raise OperationNotSupported(f, 'compose') 

3059 

3060 return per(result) 

3061 

3062 def decompose(f): 

3063 """ 

3064 Computes a functional decomposition of ``f``. 

3065 

3066 Examples 

3067 ======== 

3068 

3069 >>> from sympy import Poly 

3070 >>> from sympy.abc import x 

3071 

3072 >>> Poly(x**4 + 2*x**3 - x - 1, x, domain='ZZ').decompose() 

3073 [Poly(x**2 - x - 1, x, domain='ZZ'), Poly(x**2 + x, x, domain='ZZ')] 

3074 

3075 """ 

3076 if hasattr(f.rep, 'decompose'): 

3077 result = f.rep.decompose() 

3078 else: # pragma: no cover 

3079 raise OperationNotSupported(f, 'decompose') 

3080 

3081 return list(map(f.per, result)) 

3082 

3083 def shift(f, a): 

3084 """ 

3085 Efficiently compute Taylor shift ``f(x + a)``. 

3086 

3087 Examples 

3088 ======== 

3089 

3090 >>> from sympy import Poly 

3091 >>> from sympy.abc import x 

3092 

3093 >>> Poly(x**2 - 2*x + 1, x).shift(2) 

3094 Poly(x**2 + 2*x + 1, x, domain='ZZ') 

3095 

3096 """ 

3097 if hasattr(f.rep, 'shift'): 

3098 result = f.rep.shift(a) 

3099 else: # pragma: no cover 

3100 raise OperationNotSupported(f, 'shift') 

3101 

3102 return f.per(result) 

3103 

3104 def transform(f, p, q): 

3105 """ 

3106 Efficiently evaluate the functional transformation ``q**n * f(p/q)``. 

3107 

3108 

3109 Examples 

3110 ======== 

3111 

3112 >>> from sympy import Poly 

3113 >>> from sympy.abc import x 

3114 

3115 >>> Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1, x), Poly(x - 1, x)) 

3116 Poly(4, x, domain='ZZ') 

3117 

3118 """ 

3119 P, Q = p.unify(q) 

3120 F, P = f.unify(P) 

3121 F, Q = F.unify(Q) 

3122 

3123 if hasattr(F.rep, 'transform'): 

3124 result = F.rep.transform(P.rep, Q.rep) 

3125 else: # pragma: no cover 

3126 raise OperationNotSupported(F, 'transform') 

3127 

3128 return F.per(result) 

3129 

3130 def sturm(self, auto=True): 

3131 """ 

3132 Computes the Sturm sequence of ``f``. 

3133 

3134 Examples 

3135 ======== 

3136 

3137 >>> from sympy import Poly 

3138 >>> from sympy.abc import x 

3139 

3140 >>> Poly(x**3 - 2*x**2 + x - 3, x).sturm() 

3141 [Poly(x**3 - 2*x**2 + x - 3, x, domain='QQ'), 

3142 Poly(3*x**2 - 4*x + 1, x, domain='QQ'), 

3143 Poly(2/9*x + 25/9, x, domain='QQ'), 

3144 Poly(-2079/4, x, domain='QQ')] 

3145 

3146 """ 

3147 f = self 

3148 

3149 if auto and f.rep.dom.is_Ring: 

3150 f = f.to_field() 

3151 

3152 if hasattr(f.rep, 'sturm'): 

3153 result = f.rep.sturm() 

3154 else: # pragma: no cover 

3155 raise OperationNotSupported(f, 'sturm') 

3156 

3157 return list(map(f.per, result)) 

3158 

3159 def gff_list(f): 

3160 """ 

3161 Computes greatest factorial factorization of ``f``. 

3162 

3163 Examples 

3164 ======== 

3165 

3166 >>> from sympy import Poly 

3167 >>> from sympy.abc import x 

3168 

3169 >>> f = x**5 + 2*x**4 - x**3 - 2*x**2 

3170 

3171 >>> Poly(f).gff_list() 

3172 [(Poly(x, x, domain='ZZ'), 1), (Poly(x + 2, x, domain='ZZ'), 4)] 

3173 

3174 """ 

3175 if hasattr(f.rep, 'gff_list'): 

3176 result = f.rep.gff_list() 

3177 else: # pragma: no cover 

3178 raise OperationNotSupported(f, 'gff_list') 

3179 

3180 return [(f.per(g), k) for g, k in result] 

3181 

3182 def norm(f): 

3183 """ 

3184 Computes the product, ``Norm(f)``, of the conjugates of 

3185 a polynomial ``f`` defined over a number field ``K``. 

3186 

3187 Examples 

3188 ======== 

3189 

3190 >>> from sympy import Poly, sqrt 

3191 >>> from sympy.abc import x 

3192 

3193 >>> a, b = sqrt(2), sqrt(3) 

3194 

3195 A polynomial over a quadratic extension. 

3196 Two conjugates x - a and x + a. 

3197 

3198 >>> f = Poly(x - a, x, extension=a) 

3199 >>> f.norm() 

3200 Poly(x**2 - 2, x, domain='QQ') 

3201 

3202 A polynomial over a quartic extension. 

3203 Four conjugates x - a, x - a, x + a and x + a. 

3204 

3205 >>> f = Poly(x - a, x, extension=(a, b)) 

3206 >>> f.norm() 

3207 Poly(x**4 - 4*x**2 + 4, x, domain='QQ') 

3208 

3209 """ 

3210 if hasattr(f.rep, 'norm'): 

3211 r = f.rep.norm() 

3212 else: # pragma: no cover 

3213 raise OperationNotSupported(f, 'norm') 

3214 

3215 return f.per(r) 

3216 

3217 def sqf_norm(f): 

3218 """ 

3219 Computes square-free norm of ``f``. 

3220 

3221 Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and 

3222 ``r(x) = Norm(g(x))`` is a square-free polynomial over ``K``, 

3223 where ``a`` is the algebraic extension of the ground domain. 

3224 

3225 Examples 

3226 ======== 

3227 

3228 >>> from sympy import Poly, sqrt 

3229 >>> from sympy.abc import x 

3230 

3231 >>> s, f, r = Poly(x**2 + 1, x, extension=[sqrt(3)]).sqf_norm() 

3232 

3233 >>> s 

3234 1 

3235 >>> f 

3236 Poly(x**2 - 2*sqrt(3)*x + 4, x, domain='QQ<sqrt(3)>') 

3237 >>> r 

3238 Poly(x**4 - 4*x**2 + 16, x, domain='QQ') 

3239 

3240 """ 

3241 if hasattr(f.rep, 'sqf_norm'): 

3242 s, g, r = f.rep.sqf_norm() 

3243 else: # pragma: no cover 

3244 raise OperationNotSupported(f, 'sqf_norm') 

3245 

3246 return s, f.per(g), f.per(r) 

3247 

3248 def sqf_part(f): 

3249 """ 

3250 Computes square-free part of ``f``. 

3251 

3252 Examples 

3253 ======== 

3254 

3255 >>> from sympy import Poly 

3256 >>> from sympy.abc import x 

3257 

3258 >>> Poly(x**3 - 3*x - 2, x).sqf_part() 

3259 Poly(x**2 - x - 2, x, domain='ZZ') 

3260 

3261 """ 

3262 if hasattr(f.rep, 'sqf_part'): 

3263 result = f.rep.sqf_part() 

3264 else: # pragma: no cover 

3265 raise OperationNotSupported(f, 'sqf_part') 

3266 

3267 return f.per(result) 

3268 

3269 def sqf_list(f, all=False): 

3270 """ 

3271 Returns a list of square-free factors of ``f``. 

3272 

3273 Examples 

3274 ======== 

3275 

3276 >>> from sympy import Poly 

3277 >>> from sympy.abc import x 

3278 

3279 >>> f = 2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16 

3280 

3281 >>> Poly(f).sqf_list() 

3282 (2, [(Poly(x + 1, x, domain='ZZ'), 2), 

3283 (Poly(x + 2, x, domain='ZZ'), 3)]) 

3284 

3285 >>> Poly(f).sqf_list(all=True) 

3286 (2, [(Poly(1, x, domain='ZZ'), 1), 

3287 (Poly(x + 1, x, domain='ZZ'), 2), 

3288 (Poly(x + 2, x, domain='ZZ'), 3)]) 

3289 

3290 """ 

3291 if hasattr(f.rep, 'sqf_list'): 

3292 coeff, factors = f.rep.sqf_list(all) 

3293 else: # pragma: no cover 

3294 raise OperationNotSupported(f, 'sqf_list') 

3295 

3296 return f.rep.dom.to_sympy(coeff), [(f.per(g), k) for g, k in factors] 

3297 

3298 def sqf_list_include(f, all=False): 

3299 """ 

3300 Returns a list of square-free factors of ``f``. 

3301 

3302 Examples 

3303 ======== 

3304 

3305 >>> from sympy import Poly, expand 

3306 >>> from sympy.abc import x 

3307 

3308 >>> f = expand(2*(x + 1)**3*x**4) 

3309 >>> f 

3310 2*x**7 + 6*x**6 + 6*x**5 + 2*x**4 

3311 

3312 >>> Poly(f).sqf_list_include() 

3313 [(Poly(2, x, domain='ZZ'), 1), 

3314 (Poly(x + 1, x, domain='ZZ'), 3), 

3315 (Poly(x, x, domain='ZZ'), 4)] 

3316 

3317 >>> Poly(f).sqf_list_include(all=True) 

3318 [(Poly(2, x, domain='ZZ'), 1), 

3319 (Poly(1, x, domain='ZZ'), 2), 

3320 (Poly(x + 1, x, domain='ZZ'), 3), 

3321 (Poly(x, x, domain='ZZ'), 4)] 

3322 

3323 """ 

3324 if hasattr(f.rep, 'sqf_list_include'): 

3325 factors = f.rep.sqf_list_include(all) 

3326 else: # pragma: no cover 

3327 raise OperationNotSupported(f, 'sqf_list_include') 

3328 

3329 return [(f.per(g), k) for g, k in factors] 

3330 

3331 def factor_list(f): 

3332 """ 

3333 Returns a list of irreducible factors of ``f``. 

3334 

3335 Examples 

3336 ======== 

3337 

3338 >>> from sympy import Poly 

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

3340 

3341 >>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y 

3342 

3343 >>> Poly(f).factor_list() 

3344 (2, [(Poly(x + y, x, y, domain='ZZ'), 1), 

3345 (Poly(x**2 + 1, x, y, domain='ZZ'), 2)]) 

3346 

3347 """ 

3348 if hasattr(f.rep, 'factor_list'): 

3349 try: 

3350 coeff, factors = f.rep.factor_list() 

3351 except DomainError: 

3352 return S.One, [(f, 1)] 

3353 else: # pragma: no cover 

3354 raise OperationNotSupported(f, 'factor_list') 

3355 

3356 return f.rep.dom.to_sympy(coeff), [(f.per(g), k) for g, k in factors] 

3357 

3358 def factor_list_include(f): 

3359 """ 

3360 Returns a list of irreducible factors of ``f``. 

3361 

3362 Examples 

3363 ======== 

3364 

3365 >>> from sympy import Poly 

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

3367 

3368 >>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y 

3369 

3370 >>> Poly(f).factor_list_include() 

3371 [(Poly(2*x + 2*y, x, y, domain='ZZ'), 1), 

3372 (Poly(x**2 + 1, x, y, domain='ZZ'), 2)] 

3373 

3374 """ 

3375 if hasattr(f.rep, 'factor_list_include'): 

3376 try: 

3377 factors = f.rep.factor_list_include() 

3378 except DomainError: 

3379 return [(f, 1)] 

3380 else: # pragma: no cover 

3381 raise OperationNotSupported(f, 'factor_list_include') 

3382 

3383 return [(f.per(g), k) for g, k in factors] 

3384 

3385 def intervals(f, all=False, eps=None, inf=None, sup=None, fast=False, sqf=False): 

3386 """ 

3387 Compute isolating intervals for roots of ``f``. 

3388 

3389 For real roots the Vincent-Akritas-Strzebonski (VAS) continued fractions method is used. 

3390 

3391 References 

3392 ========== 

3393 .. [#] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative Study of Two Real Root 

3394 Isolation Methods . Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005. 

3395 .. [#] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. Vigklas: Improving the 

3396 Performance of the Continued Fractions Method Using new Bounds of Positive Roots. Nonlinear 

3397 Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008. 

3398 

3399 Examples 

3400 ======== 

3401 

3402 >>> from sympy import Poly 

3403 >>> from sympy.abc import x 

3404 

3405 >>> Poly(x**2 - 3, x).intervals() 

3406 [((-2, -1), 1), ((1, 2), 1)] 

3407 >>> Poly(x**2 - 3, x).intervals(eps=1e-2) 

3408 [((-26/15, -19/11), 1), ((19/11, 26/15), 1)] 

3409 

3410 """ 

3411 if eps is not None: 

3412 eps = QQ.convert(eps) 

3413 

3414 if eps <= 0: 

3415 raise ValueError("'eps' must be a positive rational") 

3416 

3417 if inf is not None: 

3418 inf = QQ.convert(inf) 

3419 if sup is not None: 

3420 sup = QQ.convert(sup) 

3421 

3422 if hasattr(f.rep, 'intervals'): 

3423 result = f.rep.intervals( 

3424 all=all, eps=eps, inf=inf, sup=sup, fast=fast, sqf=sqf) 

3425 else: # pragma: no cover 

3426 raise OperationNotSupported(f, 'intervals') 

3427 

3428 if sqf: 

3429 def _real(interval): 

3430 s, t = interval 

3431 return (QQ.to_sympy(s), QQ.to_sympy(t)) 

3432 

3433 if not all: 

3434 return list(map(_real, result)) 

3435 

3436 def _complex(rectangle): 

3437 (u, v), (s, t) = rectangle 

3438 return (QQ.to_sympy(u) + I*QQ.to_sympy(v), 

3439 QQ.to_sympy(s) + I*QQ.to_sympy(t)) 

3440 

3441 real_part, complex_part = result 

3442 

3443 return list(map(_real, real_part)), list(map(_complex, complex_part)) 

3444 else: 

3445 def _real(interval): 

3446 (s, t), k = interval 

3447 return ((QQ.to_sympy(s), QQ.to_sympy(t)), k) 

3448 

3449 if not all: 

3450 return list(map(_real, result)) 

3451 

3452 def _complex(rectangle): 

3453 ((u, v), (s, t)), k = rectangle 

3454 return ((QQ.to_sympy(u) + I*QQ.to_sympy(v), 

3455 QQ.to_sympy(s) + I*QQ.to_sympy(t)), k) 

3456 

3457 real_part, complex_part = result 

3458 

3459 return list(map(_real, real_part)), list(map(_complex, complex_part)) 

3460 

3461 def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False): 

3462 """ 

3463 Refine an isolating interval of a root to the given precision. 

3464 

3465 Examples 

3466 ======== 

3467 

3468 >>> from sympy import Poly 

3469 >>> from sympy.abc import x 

3470 

3471 >>> Poly(x**2 - 3, x).refine_root(1, 2, eps=1e-2) 

3472 (19/11, 26/15) 

3473 

3474 """ 

3475 if check_sqf and not f.is_sqf: 

3476 raise PolynomialError("only square-free polynomials supported") 

3477 

3478 s, t = QQ.convert(s), QQ.convert(t) 

3479 

3480 if eps is not None: 

3481 eps = QQ.convert(eps) 

3482 

3483 if eps <= 0: 

3484 raise ValueError("'eps' must be a positive rational") 

3485 

3486 if steps is not None: 

3487 steps = int(steps) 

3488 elif eps is None: 

3489 steps = 1 

3490 

3491 if hasattr(f.rep, 'refine_root'): 

3492 S, T = f.rep.refine_root(s, t, eps=eps, steps=steps, fast=fast) 

3493 else: # pragma: no cover 

3494 raise OperationNotSupported(f, 'refine_root') 

3495 

3496 return QQ.to_sympy(S), QQ.to_sympy(T) 

3497 

3498 def count_roots(f, inf=None, sup=None): 

3499 """ 

3500 Return the number of roots of ``f`` in ``[inf, sup]`` interval. 

3501 

3502 Examples 

3503 ======== 

3504 

3505 >>> from sympy import Poly, I 

3506 >>> from sympy.abc import x 

3507 

3508 >>> Poly(x**4 - 4, x).count_roots(-3, 3) 

3509 2 

3510 >>> Poly(x**4 - 4, x).count_roots(0, 1 + 3*I) 

3511 1 

3512 

3513 """ 

3514 inf_real, sup_real = True, True 

3515 

3516 if inf is not None: 

3517 inf = sympify(inf) 

3518 

3519 if inf is S.NegativeInfinity: 

3520 inf = None 

3521 else: 

3522 re, im = inf.as_real_imag() 

3523 

3524 if not im: 

3525 inf = QQ.convert(inf) 

3526 else: 

3527 inf, inf_real = list(map(QQ.convert, (re, im))), False 

3528 

3529 if sup is not None: 

3530 sup = sympify(sup) 

3531 

3532 if sup is S.Infinity: 

3533 sup = None 

3534 else: 

3535 re, im = sup.as_real_imag() 

3536 

3537 if not im: 

3538 sup = QQ.convert(sup) 

3539 else: 

3540 sup, sup_real = list(map(QQ.convert, (re, im))), False 

3541 

3542 if inf_real and sup_real: 

3543 if hasattr(f.rep, 'count_real_roots'): 

3544 count = f.rep.count_real_roots(inf=inf, sup=sup) 

3545 else: # pragma: no cover 

3546 raise OperationNotSupported(f, 'count_real_roots') 

3547 else: 

3548 if inf_real and inf is not None: 

3549 inf = (inf, QQ.zero) 

3550 

3551 if sup_real and sup is not None: 

3552 sup = (sup, QQ.zero) 

3553 

3554 if hasattr(f.rep, 'count_complex_roots'): 

3555 count = f.rep.count_complex_roots(inf=inf, sup=sup) 

3556 else: # pragma: no cover 

3557 raise OperationNotSupported(f, 'count_complex_roots') 

3558 

3559 return Integer(count) 

3560 

3561 def root(f, index, radicals=True): 

3562 """ 

3563 Get an indexed root of a polynomial. 

3564 

3565 Examples 

3566 ======== 

3567 

3568 >>> from sympy import Poly 

3569 >>> from sympy.abc import x 

3570 

3571 >>> f = Poly(2*x**3 - 7*x**2 + 4*x + 4) 

3572 

3573 >>> f.root(0) 

3574 -1/2 

3575 >>> f.root(1) 

3576 2 

3577 >>> f.root(2) 

3578 2 

3579 >>> f.root(3) 

3580 Traceback (most recent call last): 

3581 ... 

3582 IndexError: root index out of [-3, 2] range, got 3 

3583 

3584 >>> Poly(x**5 + x + 1).root(0) 

3585 CRootOf(x**3 - x**2 + 1, 0) 

3586 

3587 """ 

3588 return sympy.polys.rootoftools.rootof(f, index, radicals=radicals) 

3589 

3590 def real_roots(f, multiple=True, radicals=True): 

3591 """ 

3592 Return a list of real roots with multiplicities. 

3593 

3594 Examples 

3595 ======== 

3596 

3597 >>> from sympy import Poly 

3598 >>> from sympy.abc import x 

3599 

3600 >>> Poly(2*x**3 - 7*x**2 + 4*x + 4).real_roots() 

3601 [-1/2, 2, 2] 

3602 >>> Poly(x**3 + x + 1).real_roots() 

3603 [CRootOf(x**3 + x + 1, 0)] 

3604 

3605 """ 

3606 reals = sympy.polys.rootoftools.CRootOf.real_roots(f, radicals=radicals) 

3607 

3608 if multiple: 

3609 return reals 

3610 else: 

3611 return group(reals, multiple=False) 

3612 

3613 def all_roots(f, multiple=True, radicals=True): 

3614 """ 

3615 Return a list of real and complex roots with multiplicities. 

3616 

3617 Examples 

3618 ======== 

3619 

3620 >>> from sympy import Poly 

3621 >>> from sympy.abc import x 

3622 

3623 >>> Poly(2*x**3 - 7*x**2 + 4*x + 4).all_roots() 

3624 [-1/2, 2, 2] 

3625 >>> Poly(x**3 + x + 1).all_roots() 

3626 [CRootOf(x**3 + x + 1, 0), 

3627 CRootOf(x**3 + x + 1, 1), 

3628 CRootOf(x**3 + x + 1, 2)] 

3629 

3630 """ 

3631 roots = sympy.polys.rootoftools.CRootOf.all_roots(f, radicals=radicals) 

3632 

3633 if multiple: 

3634 return roots 

3635 else: 

3636 return group(roots, multiple=False) 

3637 

3638 def nroots(f, n=15, maxsteps=50, cleanup=True): 

3639 """ 

3640 Compute numerical approximations of roots of ``f``. 

3641 

3642 Parameters 

3643 ========== 

3644 

3645 n ... the number of digits to calculate 

3646 maxsteps ... the maximum number of iterations to do 

3647 

3648 If the accuracy `n` cannot be reached in `maxsteps`, it will raise an 

3649 exception. You need to rerun with higher maxsteps. 

3650 

3651 Examples 

3652 ======== 

3653 

3654 >>> from sympy import Poly 

3655 >>> from sympy.abc import x 

3656 

3657 >>> Poly(x**2 - 3).nroots(n=15) 

3658 [-1.73205080756888, 1.73205080756888] 

3659 >>> Poly(x**2 - 3).nroots(n=30) 

3660 [-1.73205080756887729352744634151, 1.73205080756887729352744634151] 

3661 

3662 """ 

3663 if f.is_multivariate: 

3664 raise MultivariatePolynomialError( 

3665 "Cannot compute numerical roots of %s" % f) 

3666 

3667 if f.degree() <= 0: 

3668 return [] 

3669 

3670 # For integer and rational coefficients, convert them to integers only 

3671 # (for accuracy). Otherwise just try to convert the coefficients to 

3672 # mpmath.mpc and raise an exception if the conversion fails. 

3673 if f.rep.dom is ZZ: 

3674 coeffs = [int(coeff) for coeff in f.all_coeffs()] 

3675 elif f.rep.dom is QQ: 

3676 denoms = [coeff.q for coeff in f.all_coeffs()] 

3677 fac = ilcm(*denoms) 

3678 coeffs = [int(coeff*fac) for coeff in f.all_coeffs()] 

3679 else: 

3680 coeffs = [coeff.evalf(n=n).as_real_imag() 

3681 for coeff in f.all_coeffs()] 

3682 try: 

3683 coeffs = [mpmath.mpc(*coeff) for coeff in coeffs] 

3684 except TypeError: 

3685 raise DomainError("Numerical domain expected, got %s" % \ 

3686 f.rep.dom) 

3687 

3688 dps = mpmath.mp.dps 

3689 mpmath.mp.dps = n 

3690 

3691 from sympy.functions.elementary.complexes import sign 

3692 try: 

3693 # We need to add extra precision to guard against losing accuracy. 

3694 # 10 times the degree of the polynomial seems to work well. 

3695 roots = mpmath.polyroots(coeffs, maxsteps=maxsteps, 

3696 cleanup=cleanup, error=False, extraprec=f.degree()*10) 

3697 

3698 # Mpmath puts real roots first, then complex ones (as does all_roots) 

3699 # so we make sure this convention holds here, too. 

3700 roots = list(map(sympify, 

3701 sorted(roots, key=lambda r: (1 if r.imag else 0, r.real, abs(r.imag), sign(r.imag))))) 

3702 except NoConvergence: 

3703 try: 

3704 # If roots did not converge try again with more extra precision. 

3705 roots = mpmath.polyroots(coeffs, maxsteps=maxsteps, 

3706 cleanup=cleanup, error=False, extraprec=f.degree()*15) 

3707 roots = list(map(sympify, 

3708 sorted(roots, key=lambda r: (1 if r.imag else 0, r.real, abs(r.imag), sign(r.imag))))) 

3709 except NoConvergence: 

3710 raise NoConvergence( 

3711 'convergence to root failed; try n < %s or maxsteps > %s' % ( 

3712 n, maxsteps)) 

3713 finally: 

3714 mpmath.mp.dps = dps 

3715 

3716 return roots 

3717 

3718 def ground_roots(f): 

3719 """ 

3720 Compute roots of ``f`` by factorization in the ground domain. 

3721 

3722 Examples 

3723 ======== 

3724 

3725 >>> from sympy import Poly 

3726 >>> from sympy.abc import x 

3727 

3728 >>> Poly(x**6 - 4*x**4 + 4*x**3 - x**2).ground_roots() 

3729 {0: 2, 1: 2} 

3730 

3731 """ 

3732 if f.is_multivariate: 

3733 raise MultivariatePolynomialError( 

3734 "Cannot compute ground roots of %s" % f) 

3735 

3736 roots = {} 

3737 

3738 for factor, k in f.factor_list()[1]: 

3739 if factor.is_linear: 

3740 a, b = factor.all_coeffs() 

3741 roots[-b/a] = k 

3742 

3743 return roots 

3744 

3745 def nth_power_roots_poly(f, n): 

3746 """ 

3747 Construct a polynomial with n-th powers of roots of ``f``. 

3748 

3749 Examples 

3750 ======== 

3751 

3752 >>> from sympy import Poly 

3753 >>> from sympy.abc import x 

3754 

3755 >>> f = Poly(x**4 - x**2 + 1) 

3756 

3757 >>> f.nth_power_roots_poly(2) 

3758 Poly(x**4 - 2*x**3 + 3*x**2 - 2*x + 1, x, domain='ZZ') 

3759 >>> f.nth_power_roots_poly(3) 

3760 Poly(x**4 + 2*x**2 + 1, x, domain='ZZ') 

3761 >>> f.nth_power_roots_poly(4) 

3762 Poly(x**4 + 2*x**3 + 3*x**2 + 2*x + 1, x, domain='ZZ') 

3763 >>> f.nth_power_roots_poly(12) 

3764 Poly(x**4 - 4*x**3 + 6*x**2 - 4*x + 1, x, domain='ZZ') 

3765 

3766 """ 

3767 if f.is_multivariate: 

3768 raise MultivariatePolynomialError( 

3769 "must be a univariate polynomial") 

3770 

3771 N = sympify(n) 

3772 

3773 if N.is_Integer and N >= 1: 

3774 n = int(N) 

3775 else: 

3776 raise ValueError("'n' must an integer and n >= 1, got %s" % n) 

3777 

3778 x = f.gen 

3779 t = Dummy('t') 

3780 

3781 r = f.resultant(f.__class__.from_expr(x**n - t, x, t)) 

3782 

3783 return r.replace(t, x) 

3784 

3785 def same_root(f, a, b): 

3786 """ 

3787 Decide whether two roots of this polynomial are equal. 

3788 

3789 Examples 

3790 ======== 

3791 

3792 >>> from sympy import Poly, cyclotomic_poly, exp, I, pi 

3793 >>> f = Poly(cyclotomic_poly(5)) 

3794 >>> r0 = exp(2*I*pi/5) 

3795 >>> indices = [i for i, r in enumerate(f.all_roots()) if f.same_root(r, r0)] 

3796 >>> print(indices) 

3797 [3] 

3798 

3799 Raises 

3800 ====== 

3801 

3802 DomainError 

3803 If the domain of the polynomial is not :ref:`ZZ`, :ref:`QQ`, 

3804 :ref:`RR`, or :ref:`CC`. 

3805 MultivariatePolynomialError 

3806 If the polynomial is not univariate. 

3807 PolynomialError 

3808 If the polynomial is of degree < 2. 

3809 

3810 """ 

3811 if f.is_multivariate: 

3812 raise MultivariatePolynomialError( 

3813 "Must be a univariate polynomial") 

3814 

3815 dom_delta_sq = f.rep.mignotte_sep_bound_squared() 

3816 delta_sq = f.domain.get_field().to_sympy(dom_delta_sq) 

3817 # We have delta_sq = delta**2, where delta is a lower bound on the 

3818 # minimum separation between any two roots of this polynomial. 

3819 # Let eps = delta/3, and define eps_sq = eps**2 = delta**2/9. 

3820 eps_sq = delta_sq / 9 

3821 

3822 r, _, _, _ = evalf(1/eps_sq, 1, {}) 

3823 n = fastlog(r) 

3824 # Then 2^n > 1/eps**2. 

3825 m = (n // 2) + (n % 2) 

3826 # Then 2^(-m) < eps. 

3827 ev = lambda x: quad_to_mpmath(_evalf_with_bounded_error(x, m=m)) 

3828 

3829 # Then for any complex numbers a, b we will have 

3830 # |a - ev(a)| < eps and |b - ev(b)| < eps. 

3831 # So if |ev(a) - ev(b)|**2 < eps**2, then 

3832 # |ev(a) - ev(b)| < eps, hence |a - b| < 3*eps = delta. 

3833 A, B = ev(a), ev(b) 

3834 return (A.real - B.real)**2 + (A.imag - B.imag)**2 < eps_sq 

3835 

3836 def cancel(f, g, include=False): 

3837 """ 

3838 Cancel common factors in a rational function ``f/g``. 

3839 

3840 Examples 

3841 ======== 

3842 

3843 >>> from sympy import Poly 

3844 >>> from sympy.abc import x 

3845 

3846 >>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x)) 

3847 (1, Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ')) 

3848 

3849 >>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x), include=True) 

3850 (Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ')) 

3851 

3852 """ 

3853 dom, per, F, G = f._unify(g) 

3854 

3855 if hasattr(F, 'cancel'): 

3856 result = F.cancel(G, include=include) 

3857 else: # pragma: no cover 

3858 raise OperationNotSupported(f, 'cancel') 

3859 

3860 if not include: 

3861 if dom.has_assoc_Ring: 

3862 dom = dom.get_ring() 

3863 

3864 cp, cq, p, q = result 

3865 

3866 cp = dom.to_sympy(cp) 

3867 cq = dom.to_sympy(cq) 

3868 

3869 return cp/cq, per(p), per(q) 

3870 else: 

3871 return tuple(map(per, result)) 

3872 

3873 def make_monic_over_integers_by_scaling_roots(f): 

3874 """ 

3875 Turn any univariate polynomial over :ref:`QQ` or :ref:`ZZ` into a monic 

3876 polynomial over :ref:`ZZ`, by scaling the roots as necessary. 

3877 

3878 Explanation 

3879 =========== 

3880 

3881 This operation can be performed whether or not *f* is irreducible; when 

3882 it is, this can be understood as determining an algebraic integer 

3883 generating the same field as a root of *f*. 

3884 

3885 Examples 

3886 ======== 

3887 

3888 >>> from sympy import Poly, S 

3889 >>> from sympy.abc import x 

3890 >>> f = Poly(x**2/2 + S(1)/4 * x + S(1)/8, x, domain='QQ') 

3891 >>> f.make_monic_over_integers_by_scaling_roots() 

3892 (Poly(x**2 + 2*x + 4, x, domain='ZZ'), 4) 

3893 

3894 Returns 

3895 ======= 

3896 

3897 Pair ``(g, c)`` 

3898 g is the polynomial 

3899 

3900 c is the integer by which the roots had to be scaled 

3901 

3902 """ 

3903 if not f.is_univariate or f.domain not in [ZZ, QQ]: 

3904 raise ValueError('Polynomial must be univariate over ZZ or QQ.') 

3905 if f.is_monic and f.domain == ZZ: 

3906 return f, ZZ.one 

3907 else: 

3908 fm = f.monic() 

3909 c, _ = fm.clear_denoms() 

3910 return fm.transform(Poly(fm.gen), c).to_ring(), c 

3911 

3912 def galois_group(f, by_name=False, max_tries=30, randomize=False): 

3913 """ 

3914 Compute the Galois group of this polynomial. 

3915 

3916 Examples 

3917 ======== 

3918 

3919 >>> from sympy import Poly 

3920 >>> from sympy.abc import x 

3921 >>> f = Poly(x**4 - 2) 

3922 >>> G, _ = f.galois_group(by_name=True) 

3923 >>> print(G) 

3924 S4TransitiveSubgroups.D4 

3925 

3926 See Also 

3927 ======== 

3928 

3929 sympy.polys.numberfields.galoisgroups.galois_group 

3930 

3931 """ 

3932 from sympy.polys.numberfields.galoisgroups import ( 

3933 _galois_group_degree_3, _galois_group_degree_4_lookup, 

3934 _galois_group_degree_5_lookup_ext_factor, 

3935 _galois_group_degree_6_lookup, 

3936 ) 

3937 if (not f.is_univariate 

3938 or not f.is_irreducible 

3939 or f.domain not in [ZZ, QQ] 

3940 ): 

3941 raise ValueError('Polynomial must be irreducible and univariate over ZZ or QQ.') 

3942 gg = { 

3943 3: _galois_group_degree_3, 

3944 4: _galois_group_degree_4_lookup, 

3945 5: _galois_group_degree_5_lookup_ext_factor, 

3946 6: _galois_group_degree_6_lookup, 

3947 } 

3948 max_supported = max(gg.keys()) 

3949 n = f.degree() 

3950 if n > max_supported: 

3951 raise ValueError(f"Only polynomials up to degree {max_supported} are supported.") 

3952 elif n < 1: 

3953 raise ValueError("Constant polynomial has no Galois group.") 

3954 elif n == 1: 

3955 from sympy.combinatorics.galois import S1TransitiveSubgroups 

3956 name, alt = S1TransitiveSubgroups.S1, True 

3957 elif n == 2: 

3958 from sympy.combinatorics.galois import S2TransitiveSubgroups 

3959 name, alt = S2TransitiveSubgroups.S2, False 

3960 else: 

3961 g, _ = f.make_monic_over_integers_by_scaling_roots() 

3962 name, alt = gg[n](g, max_tries=max_tries, randomize=randomize) 

3963 G = name if by_name else name.get_perm_group() 

3964 return G, alt 

3965 

3966 @property 

3967 def is_zero(f): 

3968 """ 

3969 Returns ``True`` if ``f`` is a zero polynomial. 

3970 

3971 Examples 

3972 ======== 

3973 

3974 >>> from sympy import Poly 

3975 >>> from sympy.abc import x 

3976 

3977 >>> Poly(0, x).is_zero 

3978 True 

3979 >>> Poly(1, x).is_zero 

3980 False 

3981 

3982 """ 

3983 return f.rep.is_zero 

3984 

3985 @property 

3986 def is_one(f): 

3987 """ 

3988 Returns ``True`` if ``f`` is a unit polynomial. 

3989 

3990 Examples 

3991 ======== 

3992 

3993 >>> from sympy import Poly 

3994 >>> from sympy.abc import x 

3995 

3996 >>> Poly(0, x).is_one 

3997 False 

3998 >>> Poly(1, x).is_one 

3999 True 

4000 

4001 """ 

4002 return f.rep.is_one 

4003 

4004 @property 

4005 def is_sqf(f): 

4006 """ 

4007 Returns ``True`` if ``f`` is a square-free polynomial. 

4008 

4009 Examples 

4010 ======== 

4011 

4012 >>> from sympy import Poly 

4013 >>> from sympy.abc import x 

4014 

4015 >>> Poly(x**2 - 2*x + 1, x).is_sqf 

4016 False 

4017 >>> Poly(x**2 - 1, x).is_sqf 

4018 True 

4019 

4020 """ 

4021 return f.rep.is_sqf 

4022 

4023 @property 

4024 def is_monic(f): 

4025 """ 

4026 Returns ``True`` if the leading coefficient of ``f`` is one. 

4027 

4028 Examples 

4029 ======== 

4030 

4031 >>> from sympy import Poly 

4032 >>> from sympy.abc import x 

4033 

4034 >>> Poly(x + 2, x).is_monic 

4035 True 

4036 >>> Poly(2*x + 2, x).is_monic 

4037 False 

4038 

4039 """ 

4040 return f.rep.is_monic 

4041 

4042 @property 

4043 def is_primitive(f): 

4044 """ 

4045 Returns ``True`` if GCD of the coefficients of ``f`` is one. 

4046 

4047 Examples 

4048 ======== 

4049 

4050 >>> from sympy import Poly 

4051 >>> from sympy.abc import x 

4052 

4053 >>> Poly(2*x**2 + 6*x + 12, x).is_primitive 

4054 False 

4055 >>> Poly(x**2 + 3*x + 6, x).is_primitive 

4056 True 

4057 

4058 """ 

4059 return f.rep.is_primitive 

4060 

4061 @property 

4062 def is_ground(f): 

4063 """ 

4064 Returns ``True`` if ``f`` is an element of the ground domain. 

4065 

4066 Examples 

4067 ======== 

4068 

4069 >>> from sympy import Poly 

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

4071 

4072 >>> Poly(x, x).is_ground 

4073 False 

4074 >>> Poly(2, x).is_ground 

4075 True 

4076 >>> Poly(y, x).is_ground 

4077 True 

4078 

4079 """ 

4080 return f.rep.is_ground 

4081 

4082 @property 

4083 def is_linear(f): 

4084 """ 

4085 Returns ``True`` if ``f`` is linear in all its variables. 

4086 

4087 Examples 

4088 ======== 

4089 

4090 >>> from sympy import Poly 

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

4092 

4093 >>> Poly(x + y + 2, x, y).is_linear 

4094 True 

4095 >>> Poly(x*y + 2, x, y).is_linear 

4096 False 

4097 

4098 """ 

4099 return f.rep.is_linear 

4100 

4101 @property 

4102 def is_quadratic(f): 

4103 """ 

4104 Returns ``True`` if ``f`` is quadratic in all its variables. 

4105 

4106 Examples 

4107 ======== 

4108 

4109 >>> from sympy import Poly 

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

4111 

4112 >>> Poly(x*y + 2, x, y).is_quadratic 

4113 True 

4114 >>> Poly(x*y**2 + 2, x, y).is_quadratic 

4115 False 

4116 

4117 """ 

4118 return f.rep.is_quadratic 

4119 

4120 @property 

4121 def is_monomial(f): 

4122 """ 

4123 Returns ``True`` if ``f`` is zero or has only one term. 

4124 

4125 Examples 

4126 ======== 

4127 

4128 >>> from sympy import Poly 

4129 >>> from sympy.abc import x 

4130 

4131 >>> Poly(3*x**2, x).is_monomial 

4132 True 

4133 >>> Poly(3*x**2 + 1, x).is_monomial 

4134 False 

4135 

4136 """ 

4137 return f.rep.is_monomial 

4138 

4139 @property 

4140 def is_homogeneous(f): 

4141 """ 

4142 Returns ``True`` if ``f`` is a homogeneous polynomial. 

4143 

4144 A homogeneous polynomial is a polynomial whose all monomials with 

4145 non-zero coefficients have the same total degree. If you want not 

4146 only to check if a polynomial is homogeneous but also compute its 

4147 homogeneous order, then use :func:`Poly.homogeneous_order`. 

4148 

4149 Examples 

4150 ======== 

4151 

4152 >>> from sympy import Poly 

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

4154 

4155 >>> Poly(x**2 + x*y, x, y).is_homogeneous 

4156 True 

4157 >>> Poly(x**3 + x*y, x, y).is_homogeneous 

4158 False 

4159 

4160 """ 

4161 return f.rep.is_homogeneous 

4162 

4163 @property 

4164 def is_irreducible(f): 

4165 """ 

4166 Returns ``True`` if ``f`` has no factors over its domain. 

4167 

4168 Examples 

4169 ======== 

4170 

4171 >>> from sympy import Poly 

4172 >>> from sympy.abc import x 

4173 

4174 >>> Poly(x**2 + x + 1, x, modulus=2).is_irreducible 

4175 True 

4176 >>> Poly(x**2 + 1, x, modulus=2).is_irreducible 

4177 False 

4178 

4179 """ 

4180 return f.rep.is_irreducible 

4181 

4182 @property 

4183 def is_univariate(f): 

4184 """ 

4185 Returns ``True`` if ``f`` is a univariate polynomial. 

4186 

4187 Examples 

4188 ======== 

4189 

4190 >>> from sympy import Poly 

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

4192 

4193 >>> Poly(x**2 + x + 1, x).is_univariate 

4194 True 

4195 >>> Poly(x*y**2 + x*y + 1, x, y).is_univariate 

4196 False 

4197 >>> Poly(x*y**2 + x*y + 1, x).is_univariate 

4198 True 

4199 >>> Poly(x**2 + x + 1, x, y).is_univariate 

4200 False 

4201 

4202 """ 

4203 return len(f.gens) == 1 

4204 

4205 @property 

4206 def is_multivariate(f): 

4207 """ 

4208 Returns ``True`` if ``f`` is a multivariate polynomial. 

4209 

4210 Examples 

4211 ======== 

4212 

4213 >>> from sympy import Poly 

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

4215 

4216 >>> Poly(x**2 + x + 1, x).is_multivariate 

4217 False 

4218 >>> Poly(x*y**2 + x*y + 1, x, y).is_multivariate 

4219 True 

4220 >>> Poly(x*y**2 + x*y + 1, x).is_multivariate 

4221 False 

4222 >>> Poly(x**2 + x + 1, x, y).is_multivariate 

4223 True 

4224 

4225 """ 

4226 return len(f.gens) != 1 

4227 

4228 @property 

4229 def is_cyclotomic(f): 

4230 """ 

4231 Returns ``True`` if ``f`` is a cyclotomic polnomial. 

4232 

4233 Examples 

4234 ======== 

4235 

4236 >>> from sympy import Poly 

4237 >>> from sympy.abc import x 

4238 

4239 >>> f = x**16 + x**14 - x**10 + x**8 - x**6 + x**2 + 1 

4240 

4241 >>> Poly(f).is_cyclotomic 

4242 False 

4243 

4244 >>> g = x**16 + x**14 - x**10 - x**8 - x**6 + x**2 + 1 

4245 

4246 >>> Poly(g).is_cyclotomic 

4247 True 

4248 

4249 """ 

4250 return f.rep.is_cyclotomic 

4251 

4252 def __abs__(f): 

4253 return f.abs() 

4254 

4255 def __neg__(f): 

4256 return f.neg() 

4257 

4258 @_polifyit 

4259 def __add__(f, g): 

4260 return f.add(g) 

4261 

4262 @_polifyit 

4263 def __radd__(f, g): 

4264 return g.add(f) 

4265 

4266 @_polifyit 

4267 def __sub__(f, g): 

4268 return f.sub(g) 

4269 

4270 @_polifyit 

4271 def __rsub__(f, g): 

4272 return g.sub(f) 

4273 

4274 @_polifyit 

4275 def __mul__(f, g): 

4276 return f.mul(g) 

4277 

4278 @_polifyit 

4279 def __rmul__(f, g): 

4280 return g.mul(f) 

4281 

4282 @_sympifyit('n', NotImplemented) 

4283 def __pow__(f, n): 

4284 if n.is_Integer and n >= 0: 

4285 return f.pow(n) 

4286 else: 

4287 return NotImplemented 

4288 

4289 @_polifyit 

4290 def __divmod__(f, g): 

4291 return f.div(g) 

4292 

4293 @_polifyit 

4294 def __rdivmod__(f, g): 

4295 return g.div(f) 

4296 

4297 @_polifyit 

4298 def __mod__(f, g): 

4299 return f.rem(g) 

4300 

4301 @_polifyit 

4302 def __rmod__(f, g): 

4303 return g.rem(f) 

4304 

4305 @_polifyit 

4306 def __floordiv__(f, g): 

4307 return f.quo(g) 

4308 

4309 @_polifyit 

4310 def __rfloordiv__(f, g): 

4311 return g.quo(f) 

4312 

4313 @_sympifyit('g', NotImplemented) 

4314 def __truediv__(f, g): 

4315 return f.as_expr()/g.as_expr() 

4316 

4317 @_sympifyit('g', NotImplemented) 

4318 def __rtruediv__(f, g): 

4319 return g.as_expr()/f.as_expr() 

4320 

4321 @_sympifyit('other', NotImplemented) 

4322 def __eq__(self, other): 

4323 f, g = self, other 

4324 

4325 if not g.is_Poly: 

4326 try: 

4327 g = f.__class__(g, f.gens, domain=f.get_domain()) 

4328 except (PolynomialError, DomainError, CoercionFailed): 

4329 return False 

4330 

4331 if f.gens != g.gens: 

4332 return False 

4333 

4334 if f.rep.dom != g.rep.dom: 

4335 return False 

4336 

4337 return f.rep == g.rep 

4338 

4339 @_sympifyit('g', NotImplemented) 

4340 def __ne__(f, g): 

4341 return not f == g 

4342 

4343 def __bool__(f): 

4344 return not f.is_zero 

4345 

4346 def eq(f, g, strict=False): 

4347 if not strict: 

4348 return f == g 

4349 else: 

4350 return f._strict_eq(sympify(g)) 

4351 

4352 def ne(f, g, strict=False): 

4353 return not f.eq(g, strict=strict) 

4354 

4355 def _strict_eq(f, g): 

4356 return isinstance(g, f.__class__) and f.gens == g.gens and f.rep.eq(g.rep, strict=True) 

4357 

4358 

4359@public 

4360class PurePoly(Poly): 

4361 """Class for representing pure polynomials. """ 

4362 

4363 def _hashable_content(self): 

4364 """Allow SymPy to hash Poly instances. """ 

4365 return (self.rep,) 

4366 

4367 def __hash__(self): 

4368 return super().__hash__() 

4369 

4370 @property 

4371 def free_symbols(self): 

4372 """ 

4373 Free symbols of a polynomial. 

4374 

4375 Examples 

4376 ======== 

4377 

4378 >>> from sympy import PurePoly 

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

4380 

4381 >>> PurePoly(x**2 + 1).free_symbols 

4382 set() 

4383 >>> PurePoly(x**2 + y).free_symbols 

4384 set() 

4385 >>> PurePoly(x**2 + y, x).free_symbols 

4386 {y} 

4387 

4388 """ 

4389 return self.free_symbols_in_domain 

4390 

4391 @_sympifyit('other', NotImplemented) 

4392 def __eq__(self, other): 

4393 f, g = self, other 

4394 

4395 if not g.is_Poly: 

4396 try: 

4397 g = f.__class__(g, f.gens, domain=f.get_domain()) 

4398 except (PolynomialError, DomainError, CoercionFailed): 

4399 return False 

4400 

4401 if len(f.gens) != len(g.gens): 

4402 return False 

4403 

4404 if f.rep.dom != g.rep.dom: 

4405 try: 

4406 dom = f.rep.dom.unify(g.rep.dom, f.gens) 

4407 except UnificationFailed: 

4408 return False 

4409 

4410 f = f.set_domain(dom) 

4411 g = g.set_domain(dom) 

4412 

4413 return f.rep == g.rep 

4414 

4415 def _strict_eq(f, g): 

4416 return isinstance(g, f.__class__) and f.rep.eq(g.rep, strict=True) 

4417 

4418 def _unify(f, g): 

4419 g = sympify(g) 

4420 

4421 if not g.is_Poly: 

4422 try: 

4423 return f.rep.dom, f.per, f.rep, f.rep.per(f.rep.dom.from_sympy(g)) 

4424 except CoercionFailed: 

4425 raise UnificationFailed("Cannot unify %s with %s" % (f, g)) 

4426 

4427 if len(f.gens) != len(g.gens): 

4428 raise UnificationFailed("Cannot unify %s with %s" % (f, g)) 

4429 

4430 if not (isinstance(f.rep, DMP) and isinstance(g.rep, DMP)): 

4431 raise UnificationFailed("Cannot unify %s with %s" % (f, g)) 

4432 

4433 cls = f.__class__ 

4434 gens = f.gens 

4435 

4436 dom = f.rep.dom.unify(g.rep.dom, gens) 

4437 

4438 F = f.rep.convert(dom) 

4439 G = g.rep.convert(dom) 

4440 

4441 def per(rep, dom=dom, gens=gens, remove=None): 

4442 if remove is not None: 

4443 gens = gens[:remove] + gens[remove + 1:] 

4444 

4445 if not gens: 

4446 return dom.to_sympy(rep) 

4447 

4448 return cls.new(rep, *gens) 

4449 

4450 return dom, per, F, G 

4451 

4452 

4453@public 

4454def poly_from_expr(expr, *gens, **args): 

4455 """Construct a polynomial from an expression. """ 

4456 opt = options.build_options(gens, args) 

4457 return _poly_from_expr(expr, opt) 

4458 

4459 

4460def _poly_from_expr(expr, opt): 

4461 """Construct a polynomial from an expression. """ 

4462 orig, expr = expr, sympify(expr) 

4463 

4464 if not isinstance(expr, Basic): 

4465 raise PolificationFailed(opt, orig, expr) 

4466 elif expr.is_Poly: 

4467 poly = expr.__class__._from_poly(expr, opt) 

4468 

4469 opt.gens = poly.gens 

4470 opt.domain = poly.domain 

4471 

4472 if opt.polys is None: 

4473 opt.polys = True 

4474 

4475 return poly, opt 

4476 elif opt.expand: 

4477 expr = expr.expand() 

4478 

4479 rep, opt = _dict_from_expr(expr, opt) 

4480 if not opt.gens: 

4481 raise PolificationFailed(opt, orig, expr) 

4482 

4483 monoms, coeffs = list(zip(*list(rep.items()))) 

4484 domain = opt.domain 

4485 

4486 if domain is None: 

4487 opt.domain, coeffs = construct_domain(coeffs, opt=opt) 

4488 else: 

4489 coeffs = list(map(domain.from_sympy, coeffs)) 

4490 

4491 rep = dict(list(zip(monoms, coeffs))) 

4492 poly = Poly._from_dict(rep, opt) 

4493 

4494 if opt.polys is None: 

4495 opt.polys = False 

4496 

4497 return poly, opt 

4498 

4499 

4500@public 

4501def parallel_poly_from_expr(exprs, *gens, **args): 

4502 """Construct polynomials from expressions. """ 

4503 opt = options.build_options(gens, args) 

4504 return _parallel_poly_from_expr(exprs, opt) 

4505 

4506 

4507def _parallel_poly_from_expr(exprs, opt): 

4508 """Construct polynomials from expressions. """ 

4509 if len(exprs) == 2: 

4510 f, g = exprs 

4511 

4512 if isinstance(f, Poly) and isinstance(g, Poly): 

4513 f = f.__class__._from_poly(f, opt) 

4514 g = g.__class__._from_poly(g, opt) 

4515 

4516 f, g = f.unify(g) 

4517 

4518 opt.gens = f.gens 

4519 opt.domain = f.domain 

4520 

4521 if opt.polys is None: 

4522 opt.polys = True 

4523 

4524 return [f, g], opt 

4525 

4526 origs, exprs = list(exprs), [] 

4527 _exprs, _polys = [], [] 

4528 

4529 failed = False 

4530 

4531 for i, expr in enumerate(origs): 

4532 expr = sympify(expr) 

4533 

4534 if isinstance(expr, Basic): 

4535 if expr.is_Poly: 

4536 _polys.append(i) 

4537 else: 

4538 _exprs.append(i) 

4539 

4540 if opt.expand: 

4541 expr = expr.expand() 

4542 else: 

4543 failed = True 

4544 

4545 exprs.append(expr) 

4546 

4547 if failed: 

4548 raise PolificationFailed(opt, origs, exprs, True) 

4549 

4550 if _polys: 

4551 # XXX: this is a temporary solution 

4552 for i in _polys: 

4553 exprs[i] = exprs[i].as_expr() 

4554 

4555 reps, opt = _parallel_dict_from_expr(exprs, opt) 

4556 if not opt.gens: 

4557 raise PolificationFailed(opt, origs, exprs, True) 

4558 

4559 from sympy.functions.elementary.piecewise import Piecewise 

4560 for k in opt.gens: 

4561 if isinstance(k, Piecewise): 

4562 raise PolynomialError("Piecewise generators do not make sense") 

4563 

4564 coeffs_list, lengths = [], [] 

4565 

4566 all_monoms = [] 

4567 all_coeffs = [] 

4568 

4569 for rep in reps: 

4570 monoms, coeffs = list(zip(*list(rep.items()))) 

4571 

4572 coeffs_list.extend(coeffs) 

4573 all_monoms.append(monoms) 

4574 

4575 lengths.append(len(coeffs)) 

4576 

4577 domain = opt.domain 

4578 

4579 if domain is None: 

4580 opt.domain, coeffs_list = construct_domain(coeffs_list, opt=opt) 

4581 else: 

4582 coeffs_list = list(map(domain.from_sympy, coeffs_list)) 

4583 

4584 for k in lengths: 

4585 all_coeffs.append(coeffs_list[:k]) 

4586 coeffs_list = coeffs_list[k:] 

4587 

4588 polys = [] 

4589 

4590 for monoms, coeffs in zip(all_monoms, all_coeffs): 

4591 rep = dict(list(zip(monoms, coeffs))) 

4592 poly = Poly._from_dict(rep, opt) 

4593 polys.append(poly) 

4594 

4595 if opt.polys is None: 

4596 opt.polys = bool(_polys) 

4597 

4598 return polys, opt 

4599 

4600 

4601def _update_args(args, key, value): 

4602 """Add a new ``(key, value)`` pair to arguments ``dict``. """ 

4603 args = dict(args) 

4604 

4605 if key not in args: 

4606 args[key] = value 

4607 

4608 return args 

4609 

4610 

4611@public 

4612def degree(f, gen=0): 

4613 """ 

4614 Return the degree of ``f`` in the given variable. 

4615 

4616 The degree of 0 is negative infinity. 

4617 

4618 Examples 

4619 ======== 

4620 

4621 >>> from sympy import degree 

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

4623 

4624 >>> degree(x**2 + y*x + 1, gen=x) 

4625 2 

4626 >>> degree(x**2 + y*x + 1, gen=y) 

4627 1 

4628 >>> degree(0, x) 

4629 -oo 

4630 

4631 See also 

4632 ======== 

4633 

4634 sympy.polys.polytools.Poly.total_degree 

4635 degree_list 

4636 """ 

4637 

4638 f = sympify(f, strict=True) 

4639 gen_is_Num = sympify(gen, strict=True).is_Number 

4640 if f.is_Poly: 

4641 p = f 

4642 isNum = p.as_expr().is_Number 

4643 else: 

4644 isNum = f.is_Number 

4645 if not isNum: 

4646 if gen_is_Num: 

4647 p, _ = poly_from_expr(f) 

4648 else: 

4649 p, _ = poly_from_expr(f, gen) 

4650 

4651 if isNum: 

4652 return S.Zero if f else S.NegativeInfinity 

4653 

4654 if not gen_is_Num: 

4655 if f.is_Poly and gen not in p.gens: 

4656 # try recast without explicit gens 

4657 p, _ = poly_from_expr(f.as_expr()) 

4658 if gen not in p.gens: 

4659 return S.Zero 

4660 elif not f.is_Poly and len(f.free_symbols) > 1: 

4661 raise TypeError(filldedent(''' 

4662 A symbolic generator of interest is required for a multivariate 

4663 expression like func = %s, e.g. degree(func, gen = %s) instead of 

4664 degree(func, gen = %s). 

4665 ''' % (f, next(ordered(f.free_symbols)), gen))) 

4666 result = p.degree(gen) 

4667 return Integer(result) if isinstance(result, int) else S.NegativeInfinity 

4668 

4669 

4670@public 

4671def total_degree(f, *gens): 

4672 """ 

4673 Return the total_degree of ``f`` in the given variables. 

4674 

4675 Examples 

4676 ======== 

4677 >>> from sympy import total_degree, Poly 

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

4679 

4680 >>> total_degree(1) 

4681 0 

4682 >>> total_degree(x + x*y) 

4683 2 

4684 >>> total_degree(x + x*y, x) 

4685 1 

4686 

4687 If the expression is a Poly and no variables are given 

4688 then the generators of the Poly will be used: 

4689 

4690 >>> p = Poly(x + x*y, y) 

4691 >>> total_degree(p) 

4692 1 

4693 

4694 To deal with the underlying expression of the Poly, convert 

4695 it to an Expr: 

4696 

4697 >>> total_degree(p.as_expr()) 

4698 2 

4699 

4700 This is done automatically if any variables are given: 

4701 

4702 >>> total_degree(p, x) 

4703 1 

4704 

4705 See also 

4706 ======== 

4707 degree 

4708 """ 

4709 

4710 p = sympify(f) 

4711 if p.is_Poly: 

4712 p = p.as_expr() 

4713 if p.is_Number: 

4714 rv = 0 

4715 else: 

4716 if f.is_Poly: 

4717 gens = gens or f.gens 

4718 rv = Poly(p, gens).total_degree() 

4719 

4720 return Integer(rv) 

4721 

4722 

4723@public 

4724def degree_list(f, *gens, **args): 

4725 """ 

4726 Return a list of degrees of ``f`` in all variables. 

4727 

4728 Examples 

4729 ======== 

4730 

4731 >>> from sympy import degree_list 

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

4733 

4734 >>> degree_list(x**2 + y*x + 1) 

4735 (2, 1) 

4736 

4737 """ 

4738 options.allowed_flags(args, ['polys']) 

4739 

4740 try: 

4741 F, opt = poly_from_expr(f, *gens, **args) 

4742 except PolificationFailed as exc: 

4743 raise ComputationFailed('degree_list', 1, exc) 

4744 

4745 degrees = F.degree_list() 

4746 

4747 return tuple(map(Integer, degrees)) 

4748 

4749 

4750@public 

4751def LC(f, *gens, **args): 

4752 """ 

4753 Return the leading coefficient of ``f``. 

4754 

4755 Examples 

4756 ======== 

4757 

4758 >>> from sympy import LC 

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

4760 

4761 >>> LC(4*x**2 + 2*x*y**2 + x*y + 3*y) 

4762 4 

4763 

4764 """ 

4765 options.allowed_flags(args, ['polys']) 

4766 

4767 try: 

4768 F, opt = poly_from_expr(f, *gens, **args) 

4769 except PolificationFailed as exc: 

4770 raise ComputationFailed('LC', 1, exc) 

4771 

4772 return F.LC(order=opt.order) 

4773 

4774 

4775@public 

4776def LM(f, *gens, **args): 

4777 """ 

4778 Return the leading monomial of ``f``. 

4779 

4780 Examples 

4781 ======== 

4782 

4783 >>> from sympy import LM 

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

4785 

4786 >>> LM(4*x**2 + 2*x*y**2 + x*y + 3*y) 

4787 x**2 

4788 

4789 """ 

4790 options.allowed_flags(args, ['polys']) 

4791 

4792 try: 

4793 F, opt = poly_from_expr(f, *gens, **args) 

4794 except PolificationFailed as exc: 

4795 raise ComputationFailed('LM', 1, exc) 

4796 

4797 monom = F.LM(order=opt.order) 

4798 return monom.as_expr() 

4799 

4800 

4801@public 

4802def LT(f, *gens, **args): 

4803 """ 

4804 Return the leading term of ``f``. 

4805 

4806 Examples 

4807 ======== 

4808 

4809 >>> from sympy import LT 

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

4811 

4812 >>> LT(4*x**2 + 2*x*y**2 + x*y + 3*y) 

4813 4*x**2 

4814 

4815 """ 

4816 options.allowed_flags(args, ['polys']) 

4817 

4818 try: 

4819 F, opt = poly_from_expr(f, *gens, **args) 

4820 except PolificationFailed as exc: 

4821 raise ComputationFailed('LT', 1, exc) 

4822 

4823 monom, coeff = F.LT(order=opt.order) 

4824 return coeff*monom.as_expr() 

4825 

4826 

4827@public 

4828def pdiv(f, g, *gens, **args): 

4829 """ 

4830 Compute polynomial pseudo-division of ``f`` and ``g``. 

4831 

4832 Examples 

4833 ======== 

4834 

4835 >>> from sympy import pdiv 

4836 >>> from sympy.abc import x 

4837 

4838 >>> pdiv(x**2 + 1, 2*x - 4) 

4839 (2*x + 4, 20) 

4840 

4841 """ 

4842 options.allowed_flags(args, ['polys']) 

4843 

4844 try: 

4845 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

4846 except PolificationFailed as exc: 

4847 raise ComputationFailed('pdiv', 2, exc) 

4848 

4849 q, r = F.pdiv(G) 

4850 

4851 if not opt.polys: 

4852 return q.as_expr(), r.as_expr() 

4853 else: 

4854 return q, r 

4855 

4856 

4857@public 

4858def prem(f, g, *gens, **args): 

4859 """ 

4860 Compute polynomial pseudo-remainder of ``f`` and ``g``. 

4861 

4862 Examples 

4863 ======== 

4864 

4865 >>> from sympy import prem 

4866 >>> from sympy.abc import x 

4867 

4868 >>> prem(x**2 + 1, 2*x - 4) 

4869 20 

4870 

4871 """ 

4872 options.allowed_flags(args, ['polys']) 

4873 

4874 try: 

4875 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

4876 except PolificationFailed as exc: 

4877 raise ComputationFailed('prem', 2, exc) 

4878 

4879 r = F.prem(G) 

4880 

4881 if not opt.polys: 

4882 return r.as_expr() 

4883 else: 

4884 return r 

4885 

4886 

4887@public 

4888def pquo(f, g, *gens, **args): 

4889 """ 

4890 Compute polynomial pseudo-quotient of ``f`` and ``g``. 

4891 

4892 Examples 

4893 ======== 

4894 

4895 >>> from sympy import pquo 

4896 >>> from sympy.abc import x 

4897 

4898 >>> pquo(x**2 + 1, 2*x - 4) 

4899 2*x + 4 

4900 >>> pquo(x**2 - 1, 2*x - 1) 

4901 2*x + 1 

4902 

4903 """ 

4904 options.allowed_flags(args, ['polys']) 

4905 

4906 try: 

4907 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

4908 except PolificationFailed as exc: 

4909 raise ComputationFailed('pquo', 2, exc) 

4910 

4911 try: 

4912 q = F.pquo(G) 

4913 except ExactQuotientFailed: 

4914 raise ExactQuotientFailed(f, g) 

4915 

4916 if not opt.polys: 

4917 return q.as_expr() 

4918 else: 

4919 return q 

4920 

4921 

4922@public 

4923def pexquo(f, g, *gens, **args): 

4924 """ 

4925 Compute polynomial exact pseudo-quotient of ``f`` and ``g``. 

4926 

4927 Examples 

4928 ======== 

4929 

4930 >>> from sympy import pexquo 

4931 >>> from sympy.abc import x 

4932 

4933 >>> pexquo(x**2 - 1, 2*x - 2) 

4934 2*x + 2 

4935 

4936 >>> pexquo(x**2 + 1, 2*x - 4) 

4937 Traceback (most recent call last): 

4938 ... 

4939 ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1 

4940 

4941 """ 

4942 options.allowed_flags(args, ['polys']) 

4943 

4944 try: 

4945 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

4946 except PolificationFailed as exc: 

4947 raise ComputationFailed('pexquo', 2, exc) 

4948 

4949 q = F.pexquo(G) 

4950 

4951 if not opt.polys: 

4952 return q.as_expr() 

4953 else: 

4954 return q 

4955 

4956 

4957@public 

4958def div(f, g, *gens, **args): 

4959 """ 

4960 Compute polynomial division of ``f`` and ``g``. 

4961 

4962 Examples 

4963 ======== 

4964 

4965 >>> from sympy import div, ZZ, QQ 

4966 >>> from sympy.abc import x 

4967 

4968 >>> div(x**2 + 1, 2*x - 4, domain=ZZ) 

4969 (0, x**2 + 1) 

4970 >>> div(x**2 + 1, 2*x - 4, domain=QQ) 

4971 (x/2 + 1, 5) 

4972 

4973 """ 

4974 options.allowed_flags(args, ['auto', 'polys']) 

4975 

4976 try: 

4977 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

4978 except PolificationFailed as exc: 

4979 raise ComputationFailed('div', 2, exc) 

4980 

4981 q, r = F.div(G, auto=opt.auto) 

4982 

4983 if not opt.polys: 

4984 return q.as_expr(), r.as_expr() 

4985 else: 

4986 return q, r 

4987 

4988 

4989@public 

4990def rem(f, g, *gens, **args): 

4991 """ 

4992 Compute polynomial remainder of ``f`` and ``g``. 

4993 

4994 Examples 

4995 ======== 

4996 

4997 >>> from sympy import rem, ZZ, QQ 

4998 >>> from sympy.abc import x 

4999 

5000 >>> rem(x**2 + 1, 2*x - 4, domain=ZZ) 

5001 x**2 + 1 

5002 >>> rem(x**2 + 1, 2*x - 4, domain=QQ) 

5003 5 

5004 

5005 """ 

5006 options.allowed_flags(args, ['auto', 'polys']) 

5007 

5008 try: 

5009 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5010 except PolificationFailed as exc: 

5011 raise ComputationFailed('rem', 2, exc) 

5012 

5013 r = F.rem(G, auto=opt.auto) 

5014 

5015 if not opt.polys: 

5016 return r.as_expr() 

5017 else: 

5018 return r 

5019 

5020 

5021@public 

5022def quo(f, g, *gens, **args): 

5023 """ 

5024 Compute polynomial quotient of ``f`` and ``g``. 

5025 

5026 Examples 

5027 ======== 

5028 

5029 >>> from sympy import quo 

5030 >>> from sympy.abc import x 

5031 

5032 >>> quo(x**2 + 1, 2*x - 4) 

5033 x/2 + 1 

5034 >>> quo(x**2 - 1, x - 1) 

5035 x + 1 

5036 

5037 """ 

5038 options.allowed_flags(args, ['auto', 'polys']) 

5039 

5040 try: 

5041 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5042 except PolificationFailed as exc: 

5043 raise ComputationFailed('quo', 2, exc) 

5044 

5045 q = F.quo(G, auto=opt.auto) 

5046 

5047 if not opt.polys: 

5048 return q.as_expr() 

5049 else: 

5050 return q 

5051 

5052 

5053@public 

5054def exquo(f, g, *gens, **args): 

5055 """ 

5056 Compute polynomial exact quotient of ``f`` and ``g``. 

5057 

5058 Examples 

5059 ======== 

5060 

5061 >>> from sympy import exquo 

5062 >>> from sympy.abc import x 

5063 

5064 >>> exquo(x**2 - 1, x - 1) 

5065 x + 1 

5066 

5067 >>> exquo(x**2 + 1, 2*x - 4) 

5068 Traceback (most recent call last): 

5069 ... 

5070 ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1 

5071 

5072 """ 

5073 options.allowed_flags(args, ['auto', 'polys']) 

5074 

5075 try: 

5076 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5077 except PolificationFailed as exc: 

5078 raise ComputationFailed('exquo', 2, exc) 

5079 

5080 q = F.exquo(G, auto=opt.auto) 

5081 

5082 if not opt.polys: 

5083 return q.as_expr() 

5084 else: 

5085 return q 

5086 

5087 

5088@public 

5089def half_gcdex(f, g, *gens, **args): 

5090 """ 

5091 Half extended Euclidean algorithm of ``f`` and ``g``. 

5092 

5093 Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``. 

5094 

5095 Examples 

5096 ======== 

5097 

5098 >>> from sympy import half_gcdex 

5099 >>> from sympy.abc import x 

5100 

5101 >>> half_gcdex(x**4 - 2*x**3 - 6*x**2 + 12*x + 15, x**3 + x**2 - 4*x - 4) 

5102 (3/5 - x/5, x + 1) 

5103 

5104 """ 

5105 options.allowed_flags(args, ['auto', 'polys']) 

5106 

5107 try: 

5108 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5109 except PolificationFailed as exc: 

5110 domain, (a, b) = construct_domain(exc.exprs) 

5111 

5112 try: 

5113 s, h = domain.half_gcdex(a, b) 

5114 except NotImplementedError: 

5115 raise ComputationFailed('half_gcdex', 2, exc) 

5116 else: 

5117 return domain.to_sympy(s), domain.to_sympy(h) 

5118 

5119 s, h = F.half_gcdex(G, auto=opt.auto) 

5120 

5121 if not opt.polys: 

5122 return s.as_expr(), h.as_expr() 

5123 else: 

5124 return s, h 

5125 

5126 

5127@public 

5128def gcdex(f, g, *gens, **args): 

5129 """ 

5130 Extended Euclidean algorithm of ``f`` and ``g``. 

5131 

5132 Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``. 

5133 

5134 Examples 

5135 ======== 

5136 

5137 >>> from sympy import gcdex 

5138 >>> from sympy.abc import x 

5139 

5140 >>> gcdex(x**4 - 2*x**3 - 6*x**2 + 12*x + 15, x**3 + x**2 - 4*x - 4) 

5141 (3/5 - x/5, x**2/5 - 6*x/5 + 2, x + 1) 

5142 

5143 """ 

5144 options.allowed_flags(args, ['auto', 'polys']) 

5145 

5146 try: 

5147 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5148 except PolificationFailed as exc: 

5149 domain, (a, b) = construct_domain(exc.exprs) 

5150 

5151 try: 

5152 s, t, h = domain.gcdex(a, b) 

5153 except NotImplementedError: 

5154 raise ComputationFailed('gcdex', 2, exc) 

5155 else: 

5156 return domain.to_sympy(s), domain.to_sympy(t), domain.to_sympy(h) 

5157 

5158 s, t, h = F.gcdex(G, auto=opt.auto) 

5159 

5160 if not opt.polys: 

5161 return s.as_expr(), t.as_expr(), h.as_expr() 

5162 else: 

5163 return s, t, h 

5164 

5165 

5166@public 

5167def invert(f, g, *gens, **args): 

5168 """ 

5169 Invert ``f`` modulo ``g`` when possible. 

5170 

5171 Examples 

5172 ======== 

5173 

5174 >>> from sympy import invert, S, mod_inverse 

5175 >>> from sympy.abc import x 

5176 

5177 >>> invert(x**2 - 1, 2*x - 1) 

5178 -4/3 

5179 

5180 >>> invert(x**2 - 1, x - 1) 

5181 Traceback (most recent call last): 

5182 ... 

5183 NotInvertible: zero divisor 

5184 

5185 For more efficient inversion of Rationals, 

5186 use the :obj:`~.mod_inverse` function: 

5187 

5188 >>> mod_inverse(3, 5) 

5189 2 

5190 >>> (S(2)/5).invert(S(7)/3) 

5191 5/2 

5192 

5193 See Also 

5194 ======== 

5195 

5196 sympy.core.numbers.mod_inverse 

5197 

5198 """ 

5199 options.allowed_flags(args, ['auto', 'polys']) 

5200 

5201 try: 

5202 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5203 except PolificationFailed as exc: 

5204 domain, (a, b) = construct_domain(exc.exprs) 

5205 

5206 try: 

5207 return domain.to_sympy(domain.invert(a, b)) 

5208 except NotImplementedError: 

5209 raise ComputationFailed('invert', 2, exc) 

5210 

5211 h = F.invert(G, auto=opt.auto) 

5212 

5213 if not opt.polys: 

5214 return h.as_expr() 

5215 else: 

5216 return h 

5217 

5218 

5219@public 

5220def subresultants(f, g, *gens, **args): 

5221 """ 

5222 Compute subresultant PRS of ``f`` and ``g``. 

5223 

5224 Examples 

5225 ======== 

5226 

5227 >>> from sympy import subresultants 

5228 >>> from sympy.abc import x 

5229 

5230 >>> subresultants(x**2 + 1, x**2 - 1) 

5231 [x**2 + 1, x**2 - 1, -2] 

5232 

5233 """ 

5234 options.allowed_flags(args, ['polys']) 

5235 

5236 try: 

5237 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5238 except PolificationFailed as exc: 

5239 raise ComputationFailed('subresultants', 2, exc) 

5240 

5241 result = F.subresultants(G) 

5242 

5243 if not opt.polys: 

5244 return [r.as_expr() for r in result] 

5245 else: 

5246 return result 

5247 

5248 

5249@public 

5250def resultant(f, g, *gens, includePRS=False, **args): 

5251 """ 

5252 Compute resultant of ``f`` and ``g``. 

5253 

5254 Examples 

5255 ======== 

5256 

5257 >>> from sympy import resultant 

5258 >>> from sympy.abc import x 

5259 

5260 >>> resultant(x**2 + 1, x**2 - 1) 

5261 4 

5262 

5263 """ 

5264 options.allowed_flags(args, ['polys']) 

5265 

5266 try: 

5267 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5268 except PolificationFailed as exc: 

5269 raise ComputationFailed('resultant', 2, exc) 

5270 

5271 if includePRS: 

5272 result, R = F.resultant(G, includePRS=includePRS) 

5273 else: 

5274 result = F.resultant(G) 

5275 

5276 if not opt.polys: 

5277 if includePRS: 

5278 return result.as_expr(), [r.as_expr() for r in R] 

5279 return result.as_expr() 

5280 else: 

5281 if includePRS: 

5282 return result, R 

5283 return result 

5284 

5285 

5286@public 

5287def discriminant(f, *gens, **args): 

5288 """ 

5289 Compute discriminant of ``f``. 

5290 

5291 Examples 

5292 ======== 

5293 

5294 >>> from sympy import discriminant 

5295 >>> from sympy.abc import x 

5296 

5297 >>> discriminant(x**2 + 2*x + 3) 

5298 -8 

5299 

5300 """ 

5301 options.allowed_flags(args, ['polys']) 

5302 

5303 try: 

5304 F, opt = poly_from_expr(f, *gens, **args) 

5305 except PolificationFailed as exc: 

5306 raise ComputationFailed('discriminant', 1, exc) 

5307 

5308 result = F.discriminant() 

5309 

5310 if not opt.polys: 

5311 return result.as_expr() 

5312 else: 

5313 return result 

5314 

5315 

5316@public 

5317def cofactors(f, g, *gens, **args): 

5318 """ 

5319 Compute GCD and cofactors of ``f`` and ``g``. 

5320 

5321 Returns polynomials ``(h, cff, cfg)`` such that ``h = gcd(f, g)``, and 

5322 ``cff = quo(f, h)`` and ``cfg = quo(g, h)`` are, so called, cofactors 

5323 of ``f`` and ``g``. 

5324 

5325 Examples 

5326 ======== 

5327 

5328 >>> from sympy import cofactors 

5329 >>> from sympy.abc import x 

5330 

5331 >>> cofactors(x**2 - 1, x**2 - 3*x + 2) 

5332 (x - 1, x + 1, x - 2) 

5333 

5334 """ 

5335 options.allowed_flags(args, ['polys']) 

5336 

5337 try: 

5338 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5339 except PolificationFailed as exc: 

5340 domain, (a, b) = construct_domain(exc.exprs) 

5341 

5342 try: 

5343 h, cff, cfg = domain.cofactors(a, b) 

5344 except NotImplementedError: 

5345 raise ComputationFailed('cofactors', 2, exc) 

5346 else: 

5347 return domain.to_sympy(h), domain.to_sympy(cff), domain.to_sympy(cfg) 

5348 

5349 h, cff, cfg = F.cofactors(G) 

5350 

5351 if not opt.polys: 

5352 return h.as_expr(), cff.as_expr(), cfg.as_expr() 

5353 else: 

5354 return h, cff, cfg 

5355 

5356 

5357@public 

5358def gcd_list(seq, *gens, **args): 

5359 """ 

5360 Compute GCD of a list of polynomials. 

5361 

5362 Examples 

5363 ======== 

5364 

5365 >>> from sympy import gcd_list 

5366 >>> from sympy.abc import x 

5367 

5368 >>> gcd_list([x**3 - 1, x**2 - 1, x**2 - 3*x + 2]) 

5369 x - 1 

5370 

5371 """ 

5372 seq = sympify(seq) 

5373 

5374 def try_non_polynomial_gcd(seq): 

5375 if not gens and not args: 

5376 domain, numbers = construct_domain(seq) 

5377 

5378 if not numbers: 

5379 return domain.zero 

5380 elif domain.is_Numerical: 

5381 result, numbers = numbers[0], numbers[1:] 

5382 

5383 for number in numbers: 

5384 result = domain.gcd(result, number) 

5385 

5386 if domain.is_one(result): 

5387 break 

5388 

5389 return domain.to_sympy(result) 

5390 

5391 return None 

5392 

5393 result = try_non_polynomial_gcd(seq) 

5394 

5395 if result is not None: 

5396 return result 

5397 

5398 options.allowed_flags(args, ['polys']) 

5399 

5400 try: 

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

5402 

5403 # gcd for domain Q[irrational] (purely algebraic irrational) 

5404 if len(seq) > 1 and all(elt.is_algebraic and elt.is_irrational for elt in seq): 

5405 a = seq[-1] 

5406 lst = [ (a/elt).ratsimp() for elt in seq[:-1] ] 

5407 if all(frc.is_rational for frc in lst): 

5408 lc = 1 

5409 for frc in lst: 

5410 lc = lcm(lc, frc.as_numer_denom()[0]) 

5411 # abs ensures that the gcd is always non-negative 

5412 return abs(a/lc) 

5413 

5414 except PolificationFailed as exc: 

5415 result = try_non_polynomial_gcd(exc.exprs) 

5416 

5417 if result is not None: 

5418 return result 

5419 else: 

5420 raise ComputationFailed('gcd_list', len(seq), exc) 

5421 

5422 if not polys: 

5423 if not opt.polys: 

5424 return S.Zero 

5425 else: 

5426 return Poly(0, opt=opt) 

5427 

5428 result, polys = polys[0], polys[1:] 

5429 

5430 for poly in polys: 

5431 result = result.gcd(poly) 

5432 

5433 if result.is_one: 

5434 break 

5435 

5436 if not opt.polys: 

5437 return result.as_expr() 

5438 else: 

5439 return result 

5440 

5441 

5442@public 

5443def gcd(f, g=None, *gens, **args): 

5444 """ 

5445 Compute GCD of ``f`` and ``g``. 

5446 

5447 Examples 

5448 ======== 

5449 

5450 >>> from sympy import gcd 

5451 >>> from sympy.abc import x 

5452 

5453 >>> gcd(x**2 - 1, x**2 - 3*x + 2) 

5454 x - 1 

5455 

5456 """ 

5457 if hasattr(f, '__iter__'): 

5458 if g is not None: 

5459 gens = (g,) + gens 

5460 

5461 return gcd_list(f, *gens, **args) 

5462 elif g is None: 

5463 raise TypeError("gcd() takes 2 arguments or a sequence of arguments") 

5464 

5465 options.allowed_flags(args, ['polys']) 

5466 

5467 try: 

5468 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5469 

5470 # gcd for domain Q[irrational] (purely algebraic irrational) 

5471 a, b = map(sympify, (f, g)) 

5472 if a.is_algebraic and a.is_irrational and b.is_algebraic and b.is_irrational: 

5473 frc = (a/b).ratsimp() 

5474 if frc.is_rational: 

5475 # abs ensures that the returned gcd is always non-negative 

5476 return abs(a/frc.as_numer_denom()[0]) 

5477 

5478 except PolificationFailed as exc: 

5479 domain, (a, b) = construct_domain(exc.exprs) 

5480 

5481 try: 

5482 return domain.to_sympy(domain.gcd(a, b)) 

5483 except NotImplementedError: 

5484 raise ComputationFailed('gcd', 2, exc) 

5485 

5486 result = F.gcd(G) 

5487 

5488 if not opt.polys: 

5489 return result.as_expr() 

5490 else: 

5491 return result 

5492 

5493 

5494@public 

5495def lcm_list(seq, *gens, **args): 

5496 """ 

5497 Compute LCM of a list of polynomials. 

5498 

5499 Examples 

5500 ======== 

5501 

5502 >>> from sympy import lcm_list 

5503 >>> from sympy.abc import x 

5504 

5505 >>> lcm_list([x**3 - 1, x**2 - 1, x**2 - 3*x + 2]) 

5506 x**5 - x**4 - 2*x**3 - x**2 + x + 2 

5507 

5508 """ 

5509 seq = sympify(seq) 

5510 

5511 def try_non_polynomial_lcm(seq) -> Optional[Expr]: 

5512 if not gens and not args: 

5513 domain, numbers = construct_domain(seq) 

5514 

5515 if not numbers: 

5516 return domain.to_sympy(domain.one) 

5517 elif domain.is_Numerical: 

5518 result, numbers = numbers[0], numbers[1:] 

5519 

5520 for number in numbers: 

5521 result = domain.lcm(result, number) 

5522 

5523 return domain.to_sympy(result) 

5524 

5525 return None 

5526 

5527 result = try_non_polynomial_lcm(seq) 

5528 

5529 if result is not None: 

5530 return result 

5531 

5532 options.allowed_flags(args, ['polys']) 

5533 

5534 try: 

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

5536 

5537 # lcm for domain Q[irrational] (purely algebraic irrational) 

5538 if len(seq) > 1 and all(elt.is_algebraic and elt.is_irrational for elt in seq): 

5539 a = seq[-1] 

5540 lst = [ (a/elt).ratsimp() for elt in seq[:-1] ] 

5541 if all(frc.is_rational for frc in lst): 

5542 lc = 1 

5543 for frc in lst: 

5544 lc = lcm(lc, frc.as_numer_denom()[1]) 

5545 return a*lc 

5546 

5547 except PolificationFailed as exc: 

5548 result = try_non_polynomial_lcm(exc.exprs) 

5549 

5550 if result is not None: 

5551 return result 

5552 else: 

5553 raise ComputationFailed('lcm_list', len(seq), exc) 

5554 

5555 if not polys: 

5556 if not opt.polys: 

5557 return S.One 

5558 else: 

5559 return Poly(1, opt=opt) 

5560 

5561 result, polys = polys[0], polys[1:] 

5562 

5563 for poly in polys: 

5564 result = result.lcm(poly) 

5565 

5566 if not opt.polys: 

5567 return result.as_expr() 

5568 else: 

5569 return result 

5570 

5571 

5572@public 

5573def lcm(f, g=None, *gens, **args): 

5574 """ 

5575 Compute LCM of ``f`` and ``g``. 

5576 

5577 Examples 

5578 ======== 

5579 

5580 >>> from sympy import lcm 

5581 >>> from sympy.abc import x 

5582 

5583 >>> lcm(x**2 - 1, x**2 - 3*x + 2) 

5584 x**3 - 2*x**2 - x + 2 

5585 

5586 """ 

5587 if hasattr(f, '__iter__'): 

5588 if g is not None: 

5589 gens = (g,) + gens 

5590 

5591 return lcm_list(f, *gens, **args) 

5592 elif g is None: 

5593 raise TypeError("lcm() takes 2 arguments or a sequence of arguments") 

5594 

5595 options.allowed_flags(args, ['polys']) 

5596 

5597 try: 

5598 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5599 

5600 # lcm for domain Q[irrational] (purely algebraic irrational) 

5601 a, b = map(sympify, (f, g)) 

5602 if a.is_algebraic and a.is_irrational and b.is_algebraic and b.is_irrational: 

5603 frc = (a/b).ratsimp() 

5604 if frc.is_rational: 

5605 return a*frc.as_numer_denom()[1] 

5606 

5607 except PolificationFailed as exc: 

5608 domain, (a, b) = construct_domain(exc.exprs) 

5609 

5610 try: 

5611 return domain.to_sympy(domain.lcm(a, b)) 

5612 except NotImplementedError: 

5613 raise ComputationFailed('lcm', 2, exc) 

5614 

5615 result = F.lcm(G) 

5616 

5617 if not opt.polys: 

5618 return result.as_expr() 

5619 else: 

5620 return result 

5621 

5622 

5623@public 

5624def terms_gcd(f, *gens, **args): 

5625 """ 

5626 Remove GCD of terms from ``f``. 

5627 

5628 If the ``deep`` flag is True, then the arguments of ``f`` will have 

5629 terms_gcd applied to them. 

5630 

5631 If a fraction is factored out of ``f`` and ``f`` is an Add, then 

5632 an unevaluated Mul will be returned so that automatic simplification 

5633 does not redistribute it. The hint ``clear``, when set to False, can be 

5634 used to prevent such factoring when all coefficients are not fractions. 

5635 

5636 Examples 

5637 ======== 

5638 

5639 >>> from sympy import terms_gcd, cos 

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

5641 >>> terms_gcd(x**6*y**2 + x**3*y, x, y) 

5642 x**3*y*(x**3*y + 1) 

5643 

5644 The default action of polys routines is to expand the expression 

5645 given to them. terms_gcd follows this behavior: 

5646 

5647 >>> terms_gcd((3+3*x)*(x+x*y)) 

5648 3*x*(x*y + x + y + 1) 

5649 

5650 If this is not desired then the hint ``expand`` can be set to False. 

5651 In this case the expression will be treated as though it were comprised 

5652 of one or more terms: 

5653 

5654 >>> terms_gcd((3+3*x)*(x+x*y), expand=False) 

5655 (3*x + 3)*(x*y + x) 

5656 

5657 In order to traverse factors of a Mul or the arguments of other 

5658 functions, the ``deep`` hint can be used: 

5659 

5660 >>> terms_gcd((3 + 3*x)*(x + x*y), expand=False, deep=True) 

5661 3*x*(x + 1)*(y + 1) 

5662 >>> terms_gcd(cos(x + x*y), deep=True) 

5663 cos(x*(y + 1)) 

5664 

5665 Rationals are factored out by default: 

5666 

5667 >>> terms_gcd(x + y/2) 

5668 (2*x + y)/2 

5669 

5670 Only the y-term had a coefficient that was a fraction; if one 

5671 does not want to factor out the 1/2 in cases like this, the 

5672 flag ``clear`` can be set to False: 

5673 

5674 >>> terms_gcd(x + y/2, clear=False) 

5675 x + y/2 

5676 >>> terms_gcd(x*y/2 + y**2, clear=False) 

5677 y*(x/2 + y) 

5678 

5679 The ``clear`` flag is ignored if all coefficients are fractions: 

5680 

5681 >>> terms_gcd(x/3 + y/2, clear=False) 

5682 (2*x + 3*y)/6 

5683 

5684 See Also 

5685 ======== 

5686 sympy.core.exprtools.gcd_terms, sympy.core.exprtools.factor_terms 

5687 

5688 """ 

5689 

5690 orig = sympify(f) 

5691 

5692 if isinstance(f, Equality): 

5693 return Equality(*(terms_gcd(s, *gens, **args) for s in [f.lhs, f.rhs])) 

5694 elif isinstance(f, Relational): 

5695 raise TypeError("Inequalities cannot be used with terms_gcd. Found: %s" %(f,)) 

5696 

5697 if not isinstance(f, Expr) or f.is_Atom: 

5698 return orig 

5699 

5700 if args.get('deep', False): 

5701 new = f.func(*[terms_gcd(a, *gens, **args) for a in f.args]) 

5702 args.pop('deep') 

5703 args['expand'] = False 

5704 return terms_gcd(new, *gens, **args) 

5705 

5706 clear = args.pop('clear', True) 

5707 options.allowed_flags(args, ['polys']) 

5708 

5709 try: 

5710 F, opt = poly_from_expr(f, *gens, **args) 

5711 except PolificationFailed as exc: 

5712 return exc.expr 

5713 

5714 J, f = F.terms_gcd() 

5715 

5716 if opt.domain.is_Ring: 

5717 if opt.domain.is_Field: 

5718 denom, f = f.clear_denoms(convert=True) 

5719 

5720 coeff, f = f.primitive() 

5721 

5722 if opt.domain.is_Field: 

5723 coeff /= denom 

5724 else: 

5725 coeff = S.One 

5726 

5727 term = Mul(*[x**j for x, j in zip(f.gens, J)]) 

5728 if equal_valued(coeff, 1): 

5729 coeff = S.One 

5730 if term == 1: 

5731 return orig 

5732 

5733 if clear: 

5734 return _keep_coeff(coeff, term*f.as_expr()) 

5735 # base the clearing on the form of the original expression, not 

5736 # the (perhaps) Mul that we have now 

5737 coeff, f = _keep_coeff(coeff, f.as_expr(), clear=False).as_coeff_Mul() 

5738 return _keep_coeff(coeff, term*f, clear=False) 

5739 

5740 

5741@public 

5742def trunc(f, p, *gens, **args): 

5743 """ 

5744 Reduce ``f`` modulo a constant ``p``. 

5745 

5746 Examples 

5747 ======== 

5748 

5749 >>> from sympy import trunc 

5750 >>> from sympy.abc import x 

5751 

5752 >>> trunc(2*x**3 + 3*x**2 + 5*x + 7, 3) 

5753 -x**3 - x + 1 

5754 

5755 """ 

5756 options.allowed_flags(args, ['auto', 'polys']) 

5757 

5758 try: 

5759 F, opt = poly_from_expr(f, *gens, **args) 

5760 except PolificationFailed as exc: 

5761 raise ComputationFailed('trunc', 1, exc) 

5762 

5763 result = F.trunc(sympify(p)) 

5764 

5765 if not opt.polys: 

5766 return result.as_expr() 

5767 else: 

5768 return result 

5769 

5770 

5771@public 

5772def monic(f, *gens, **args): 

5773 """ 

5774 Divide all coefficients of ``f`` by ``LC(f)``. 

5775 

5776 Examples 

5777 ======== 

5778 

5779 >>> from sympy import monic 

5780 >>> from sympy.abc import x 

5781 

5782 >>> monic(3*x**2 + 4*x + 2) 

5783 x**2 + 4*x/3 + 2/3 

5784 

5785 """ 

5786 options.allowed_flags(args, ['auto', 'polys']) 

5787 

5788 try: 

5789 F, opt = poly_from_expr(f, *gens, **args) 

5790 except PolificationFailed as exc: 

5791 raise ComputationFailed('monic', 1, exc) 

5792 

5793 result = F.monic(auto=opt.auto) 

5794 

5795 if not opt.polys: 

5796 return result.as_expr() 

5797 else: 

5798 return result 

5799 

5800 

5801@public 

5802def content(f, *gens, **args): 

5803 """ 

5804 Compute GCD of coefficients of ``f``. 

5805 

5806 Examples 

5807 ======== 

5808 

5809 >>> from sympy import content 

5810 >>> from sympy.abc import x 

5811 

5812 >>> content(6*x**2 + 8*x + 12) 

5813 2 

5814 

5815 """ 

5816 options.allowed_flags(args, ['polys']) 

5817 

5818 try: 

5819 F, opt = poly_from_expr(f, *gens, **args) 

5820 except PolificationFailed as exc: 

5821 raise ComputationFailed('content', 1, exc) 

5822 

5823 return F.content() 

5824 

5825 

5826@public 

5827def primitive(f, *gens, **args): 

5828 """ 

5829 Compute content and the primitive form of ``f``. 

5830 

5831 Examples 

5832 ======== 

5833 

5834 >>> from sympy.polys.polytools import primitive 

5835 >>> from sympy.abc import x 

5836 

5837 >>> primitive(6*x**2 + 8*x + 12) 

5838 (2, 3*x**2 + 4*x + 6) 

5839 

5840 >>> eq = (2 + 2*x)*x + 2 

5841 

5842 Expansion is performed by default: 

5843 

5844 >>> primitive(eq) 

5845 (2, x**2 + x + 1) 

5846 

5847 Set ``expand`` to False to shut this off. Note that the 

5848 extraction will not be recursive; use the as_content_primitive method 

5849 for recursive, non-destructive Rational extraction. 

5850 

5851 >>> primitive(eq, expand=False) 

5852 (1, x*(2*x + 2) + 2) 

5853 

5854 >>> eq.as_content_primitive() 

5855 (2, x*(x + 1) + 1) 

5856 

5857 """ 

5858 options.allowed_flags(args, ['polys']) 

5859 

5860 try: 

5861 F, opt = poly_from_expr(f, *gens, **args) 

5862 except PolificationFailed as exc: 

5863 raise ComputationFailed('primitive', 1, exc) 

5864 

5865 cont, result = F.primitive() 

5866 if not opt.polys: 

5867 return cont, result.as_expr() 

5868 else: 

5869 return cont, result 

5870 

5871 

5872@public 

5873def compose(f, g, *gens, **args): 

5874 """ 

5875 Compute functional composition ``f(g)``. 

5876 

5877 Examples 

5878 ======== 

5879 

5880 >>> from sympy import compose 

5881 >>> from sympy.abc import x 

5882 

5883 >>> compose(x**2 + x, x - 1) 

5884 x**2 - x 

5885 

5886 """ 

5887 options.allowed_flags(args, ['polys']) 

5888 

5889 try: 

5890 (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) 

5891 except PolificationFailed as exc: 

5892 raise ComputationFailed('compose', 2, exc) 

5893 

5894 result = F.compose(G) 

5895 

5896 if not opt.polys: 

5897 return result.as_expr() 

5898 else: 

5899 return result 

5900 

5901 

5902@public 

5903def decompose(f, *gens, **args): 

5904 """ 

5905 Compute functional decomposition of ``f``. 

5906 

5907 Examples 

5908 ======== 

5909 

5910 >>> from sympy import decompose 

5911 >>> from sympy.abc import x 

5912 

5913 >>> decompose(x**4 + 2*x**3 - x - 1) 

5914 [x**2 - x - 1, x**2 + x] 

5915 

5916 """ 

5917 options.allowed_flags(args, ['polys']) 

5918 

5919 try: 

5920 F, opt = poly_from_expr(f, *gens, **args) 

5921 except PolificationFailed as exc: 

5922 raise ComputationFailed('decompose', 1, exc) 

5923 

5924 result = F.decompose() 

5925 

5926 if not opt.polys: 

5927 return [r.as_expr() for r in result] 

5928 else: 

5929 return result 

5930 

5931 

5932@public 

5933def sturm(f, *gens, **args): 

5934 """ 

5935 Compute Sturm sequence of ``f``. 

5936 

5937 Examples 

5938 ======== 

5939 

5940 >>> from sympy import sturm 

5941 >>> from sympy.abc import x 

5942 

5943 >>> sturm(x**3 - 2*x**2 + x - 3) 

5944 [x**3 - 2*x**2 + x - 3, 3*x**2 - 4*x + 1, 2*x/9 + 25/9, -2079/4] 

5945 

5946 """ 

5947 options.allowed_flags(args, ['auto', 'polys']) 

5948 

5949 try: 

5950 F, opt = poly_from_expr(f, *gens, **args) 

5951 except PolificationFailed as exc: 

5952 raise ComputationFailed('sturm', 1, exc) 

5953 

5954 result = F.sturm(auto=opt.auto) 

5955 

5956 if not opt.polys: 

5957 return [r.as_expr() for r in result] 

5958 else: 

5959 return result 

5960 

5961 

5962@public 

5963def gff_list(f, *gens, **args): 

5964 """ 

5965 Compute a list of greatest factorial factors of ``f``. 

5966 

5967 Note that the input to ff() and rf() should be Poly instances to use the 

5968 definitions here. 

5969 

5970 Examples 

5971 ======== 

5972 

5973 >>> from sympy import gff_list, ff, Poly 

5974 >>> from sympy.abc import x 

5975 

5976 >>> f = Poly(x**5 + 2*x**4 - x**3 - 2*x**2, x) 

5977 

5978 >>> gff_list(f) 

5979 [(Poly(x, x, domain='ZZ'), 1), (Poly(x + 2, x, domain='ZZ'), 4)] 

5980 

5981 >>> (ff(Poly(x), 1)*ff(Poly(x + 2), 4)) == f 

5982 True 

5983 

5984 >>> f = Poly(x**12 + 6*x**11 - 11*x**10 - 56*x**9 + 220*x**8 + 208*x**7 - \ 

5985 1401*x**6 + 1090*x**5 + 2715*x**4 - 6720*x**3 - 1092*x**2 + 5040*x, x) 

5986 

5987 >>> gff_list(f) 

5988 [(Poly(x**3 + 7, x, domain='ZZ'), 2), (Poly(x**2 + 5*x, x, domain='ZZ'), 3)] 

5989 

5990 >>> ff(Poly(x**3 + 7, x), 2)*ff(Poly(x**2 + 5*x, x), 3) == f 

5991 True 

5992 

5993 """ 

5994 options.allowed_flags(args, ['polys']) 

5995 

5996 try: 

5997 F, opt = poly_from_expr(f, *gens, **args) 

5998 except PolificationFailed as exc: 

5999 raise ComputationFailed('gff_list', 1, exc) 

6000 

6001 factors = F.gff_list() 

6002 

6003 if not opt.polys: 

6004 return [(g.as_expr(), k) for g, k in factors] 

6005 else: 

6006 return factors 

6007 

6008 

6009@public 

6010def gff(f, *gens, **args): 

6011 """Compute greatest factorial factorization of ``f``. """ 

6012 raise NotImplementedError('symbolic falling factorial') 

6013 

6014 

6015@public 

6016def sqf_norm(f, *gens, **args): 

6017 """ 

6018 Compute square-free norm of ``f``. 

6019 

6020 Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and 

6021 ``r(x) = Norm(g(x))`` is a square-free polynomial over ``K``, 

6022 where ``a`` is the algebraic extension of the ground domain. 

6023 

6024 Examples 

6025 ======== 

6026 

6027 >>> from sympy import sqf_norm, sqrt 

6028 >>> from sympy.abc import x 

6029 

6030 >>> sqf_norm(x**2 + 1, extension=[sqrt(3)]) 

6031 (1, x**2 - 2*sqrt(3)*x + 4, x**4 - 4*x**2 + 16) 

6032 

6033 """ 

6034 options.allowed_flags(args, ['polys']) 

6035 

6036 try: 

6037 F, opt = poly_from_expr(f, *gens, **args) 

6038 except PolificationFailed as exc: 

6039 raise ComputationFailed('sqf_norm', 1, exc) 

6040 

6041 s, g, r = F.sqf_norm() 

6042 

6043 if not opt.polys: 

6044 return Integer(s), g.as_expr(), r.as_expr() 

6045 else: 

6046 return Integer(s), g, r 

6047 

6048 

6049@public 

6050def sqf_part(f, *gens, **args): 

6051 """ 

6052 Compute square-free part of ``f``. 

6053 

6054 Examples 

6055 ======== 

6056 

6057 >>> from sympy import sqf_part 

6058 >>> from sympy.abc import x 

6059 

6060 >>> sqf_part(x**3 - 3*x - 2) 

6061 x**2 - x - 2 

6062 

6063 """ 

6064 options.allowed_flags(args, ['polys']) 

6065 

6066 try: 

6067 F, opt = poly_from_expr(f, *gens, **args) 

6068 except PolificationFailed as exc: 

6069 raise ComputationFailed('sqf_part', 1, exc) 

6070 

6071 result = F.sqf_part() 

6072 

6073 if not opt.polys: 

6074 return result.as_expr() 

6075 else: 

6076 return result 

6077 

6078 

6079def _sorted_factors(factors, method): 

6080 """Sort a list of ``(expr, exp)`` pairs. """ 

6081 if method == 'sqf': 

6082 def key(obj): 

6083 poly, exp = obj 

6084 rep = poly.rep.rep 

6085 return (exp, len(rep), len(poly.gens), str(poly.domain), rep) 

6086 else: 

6087 def key(obj): 

6088 poly, exp = obj 

6089 rep = poly.rep.rep 

6090 return (len(rep), len(poly.gens), exp, str(poly.domain), rep) 

6091 

6092 return sorted(factors, key=key) 

6093 

6094 

6095def _factors_product(factors): 

6096 """Multiply a list of ``(expr, exp)`` pairs. """ 

6097 return Mul(*[f.as_expr()**k for f, k in factors]) 

6098 

6099 

6100def _symbolic_factor_list(expr, opt, method): 

6101 """Helper function for :func:`_symbolic_factor`. """ 

6102 coeff, factors = S.One, [] 

6103 

6104 args = [i._eval_factor() if hasattr(i, '_eval_factor') else i 

6105 for i in Mul.make_args(expr)] 

6106 for arg in args: 

6107 if arg.is_Number or (isinstance(arg, Expr) and pure_complex(arg)): 

6108 coeff *= arg 

6109 continue 

6110 elif arg.is_Pow and arg.base != S.Exp1: 

6111 base, exp = arg.args 

6112 if base.is_Number and exp.is_Number: 

6113 coeff *= arg 

6114 continue 

6115 if base.is_Number: 

6116 factors.append((base, exp)) 

6117 continue 

6118 else: 

6119 base, exp = arg, S.One 

6120 

6121 try: 

6122 poly, _ = _poly_from_expr(base, opt) 

6123 except PolificationFailed as exc: 

6124 factors.append((exc.expr, exp)) 

6125 else: 

6126 func = getattr(poly, method + '_list') 

6127 

6128 _coeff, _factors = func() 

6129 if _coeff is not S.One: 

6130 if exp.is_Integer: 

6131 coeff *= _coeff**exp 

6132 elif _coeff.is_positive: 

6133 factors.append((_coeff, exp)) 

6134 else: 

6135 _factors.append((_coeff, S.One)) 

6136 

6137 if exp is S.One: 

6138 factors.extend(_factors) 

6139 elif exp.is_integer: 

6140 factors.extend([(f, k*exp) for f, k in _factors]) 

6141 else: 

6142 other = [] 

6143 

6144 for f, k in _factors: 

6145 if f.as_expr().is_positive: 

6146 factors.append((f, k*exp)) 

6147 else: 

6148 other.append((f, k)) 

6149 

6150 factors.append((_factors_product(other), exp)) 

6151 if method == 'sqf': 

6152 factors = [(reduce(mul, (f for f, _ in factors if _ == k)), k) 

6153 for k in {i for _, i in factors}] 

6154 

6155 return coeff, factors 

6156 

6157 

6158def _symbolic_factor(expr, opt, method): 

6159 """Helper function for :func:`_factor`. """ 

6160 if isinstance(expr, Expr): 

6161 if hasattr(expr,'_eval_factor'): 

6162 return expr._eval_factor() 

6163 coeff, factors = _symbolic_factor_list(together(expr, fraction=opt['fraction']), opt, method) 

6164 return _keep_coeff(coeff, _factors_product(factors)) 

6165 elif hasattr(expr, 'args'): 

6166 return expr.func(*[_symbolic_factor(arg, opt, method) for arg in expr.args]) 

6167 elif hasattr(expr, '__iter__'): 

6168 return expr.__class__([_symbolic_factor(arg, opt, method) for arg in expr]) 

6169 else: 

6170 return expr 

6171 

6172 

6173def _generic_factor_list(expr, gens, args, method): 

6174 """Helper function for :func:`sqf_list` and :func:`factor_list`. """ 

6175 options.allowed_flags(args, ['frac', 'polys']) 

6176 opt = options.build_options(gens, args) 

6177 

6178 expr = sympify(expr) 

6179 

6180 if isinstance(expr, (Expr, Poly)): 

6181 if isinstance(expr, Poly): 

6182 numer, denom = expr, 1 

6183 else: 

6184 numer, denom = together(expr).as_numer_denom() 

6185 

6186 cp, fp = _symbolic_factor_list(numer, opt, method) 

6187 cq, fq = _symbolic_factor_list(denom, opt, method) 

6188 

6189 if fq and not opt.frac: 

6190 raise PolynomialError("a polynomial expected, got %s" % expr) 

6191 

6192 _opt = opt.clone({"expand": True}) 

6193 

6194 for factors in (fp, fq): 

6195 for i, (f, k) in enumerate(factors): 

6196 if not f.is_Poly: 

6197 f, _ = _poly_from_expr(f, _opt) 

6198 factors[i] = (f, k) 

6199 

6200 fp = _sorted_factors(fp, method) 

6201 fq = _sorted_factors(fq, method) 

6202 

6203 if not opt.polys: 

6204 fp = [(f.as_expr(), k) for f, k in fp] 

6205 fq = [(f.as_expr(), k) for f, k in fq] 

6206 

6207 coeff = cp/cq 

6208 

6209 if not opt.frac: 

6210 return coeff, fp 

6211 else: 

6212 return coeff, fp, fq 

6213 else: 

6214 raise PolynomialError("a polynomial expected, got %s" % expr) 

6215 

6216 

6217def _generic_factor(expr, gens, args, method): 

6218 """Helper function for :func:`sqf` and :func:`factor`. """ 

6219 fraction = args.pop('fraction', True) 

6220 options.allowed_flags(args, []) 

6221 opt = options.build_options(gens, args) 

6222 opt['fraction'] = fraction 

6223 return _symbolic_factor(sympify(expr), opt, method) 

6224 

6225 

6226def to_rational_coeffs(f): 

6227 """ 

6228 try to transform a polynomial to have rational coefficients 

6229 

6230 try to find a transformation ``x = alpha*y`` 

6231 

6232 ``f(x) = lc*alpha**n * g(y)`` where ``g`` is a polynomial with 

6233 rational coefficients, ``lc`` the leading coefficient. 

6234 

6235 If this fails, try ``x = y + beta`` 

6236 ``f(x) = g(y)`` 

6237 

6238 Returns ``None`` if ``g`` not found; 

6239 ``(lc, alpha, None, g)`` in case of rescaling 

6240 ``(None, None, beta, g)`` in case of translation 

6241 

6242 Notes 

6243 ===== 

6244 

6245 Currently it transforms only polynomials without roots larger than 2. 

6246 

6247 Examples 

6248 ======== 

6249 

6250 >>> from sympy import sqrt, Poly, simplify 

6251 >>> from sympy.polys.polytools import to_rational_coeffs 

6252 >>> from sympy.abc import x 

6253 >>> p = Poly(((x**2-1)*(x-2)).subs({x:x*(1 + sqrt(2))}), x, domain='EX') 

6254 >>> lc, r, _, g = to_rational_coeffs(p) 

6255 >>> lc, r 

6256 (7 + 5*sqrt(2), 2 - 2*sqrt(2)) 

6257 >>> g 

6258 Poly(x**3 + x**2 - 1/4*x - 1/4, x, domain='QQ') 

6259 >>> r1 = simplify(1/r) 

6260 >>> Poly(lc*r**3*(g.as_expr()).subs({x:x*r1}), x, domain='EX') == p 

6261 True 

6262 

6263 """ 

6264 from sympy.simplify.simplify import simplify 

6265 

6266 def _try_rescale(f, f1=None): 

6267 """ 

6268 try rescaling ``x -> alpha*x`` to convert f to a polynomial 

6269 with rational coefficients. 

6270 Returns ``alpha, f``; if the rescaling is successful, 

6271 ``alpha`` is the rescaling factor, and ``f`` is the rescaled 

6272 polynomial; else ``alpha`` is ``None``. 

6273 """ 

6274 if not len(f.gens) == 1 or not (f.gens[0]).is_Atom: 

6275 return None, f 

6276 n = f.degree() 

6277 lc = f.LC() 

6278 f1 = f1 or f1.monic() 

6279 coeffs = f1.all_coeffs()[1:] 

6280 coeffs = [simplify(coeffx) for coeffx in coeffs] 

6281 if len(coeffs) > 1 and coeffs[-2]: 

6282 rescale1_x = simplify(coeffs[-2]/coeffs[-1]) 

6283 coeffs1 = [] 

6284 for i in range(len(coeffs)): 

6285 coeffx = simplify(coeffs[i]*rescale1_x**(i + 1)) 

6286 if not coeffx.is_rational: 

6287 break 

6288 coeffs1.append(coeffx) 

6289 else: 

6290 rescale_x = simplify(1/rescale1_x) 

6291 x = f.gens[0] 

6292 v = [x**n] 

6293 for i in range(1, n + 1): 

6294 v.append(coeffs1[i - 1]*x**(n - i)) 

6295 f = Add(*v) 

6296 f = Poly(f) 

6297 return lc, rescale_x, f 

6298 return None 

6299 

6300 def _try_translate(f, f1=None): 

6301 """ 

6302 try translating ``x -> x + alpha`` to convert f to a polynomial 

6303 with rational coefficients. 

6304 Returns ``alpha, f``; if the translating is successful, 

6305 ``alpha`` is the translating factor, and ``f`` is the shifted 

6306 polynomial; else ``alpha`` is ``None``. 

6307 """ 

6308 if not len(f.gens) == 1 or not (f.gens[0]).is_Atom: 

6309 return None, f 

6310 n = f.degree() 

6311 f1 = f1 or f1.monic() 

6312 coeffs = f1.all_coeffs()[1:] 

6313 c = simplify(coeffs[0]) 

6314 if c.is_Add and not c.is_rational: 

6315 rat, nonrat = sift(c.args, 

6316 lambda z: z.is_rational is True, binary=True) 

6317 alpha = -c.func(*nonrat)/n 

6318 f2 = f1.shift(alpha) 

6319 return alpha, f2 

6320 return None 

6321 

6322 def _has_square_roots(p): 

6323 """ 

6324 Return True if ``f`` is a sum with square roots but no other root 

6325 """ 

6326 coeffs = p.coeffs() 

6327 has_sq = False 

6328 for y in coeffs: 

6329 for x in Add.make_args(y): 

6330 f = Factors(x).factors 

6331 r = [wx.q for b, wx in f.items() if 

6332 b.is_number and wx.is_Rational and wx.q >= 2] 

6333 if not r: 

6334 continue 

6335 if min(r) == 2: 

6336 has_sq = True 

6337 if max(r) > 2: 

6338 return False 

6339 return has_sq 

6340 

6341 if f.get_domain().is_EX and _has_square_roots(f): 

6342 f1 = f.monic() 

6343 r = _try_rescale(f, f1) 

6344 if r: 

6345 return r[0], r[1], None, r[2] 

6346 else: 

6347 r = _try_translate(f, f1) 

6348 if r: 

6349 return None, None, r[0], r[1] 

6350 return None 

6351 

6352 

6353def _torational_factor_list(p, x): 

6354 """ 

6355 helper function to factor polynomial using to_rational_coeffs 

6356 

6357 Examples 

6358 ======== 

6359 

6360 >>> from sympy.polys.polytools import _torational_factor_list 

6361 >>> from sympy.abc import x 

6362 >>> from sympy import sqrt, expand, Mul 

6363 >>> p = expand(((x**2-1)*(x-2)).subs({x:x*(1 + sqrt(2))})) 

6364 >>> factors = _torational_factor_list(p, x); factors 

6365 (-2, [(-x*(1 + sqrt(2))/2 + 1, 1), (-x*(1 + sqrt(2)) - 1, 1), (-x*(1 + sqrt(2)) + 1, 1)]) 

6366 >>> expand(factors[0]*Mul(*[z[0] for z in factors[1]])) == p 

6367 True 

6368 >>> p = expand(((x**2-1)*(x-2)).subs({x:x + sqrt(2)})) 

6369 >>> factors = _torational_factor_list(p, x); factors 

6370 (1, [(x - 2 + sqrt(2), 1), (x - 1 + sqrt(2), 1), (x + 1 + sqrt(2), 1)]) 

6371 >>> expand(factors[0]*Mul(*[z[0] for z in factors[1]])) == p 

6372 True 

6373 

6374 """ 

6375 from sympy.simplify.simplify import simplify 

6376 p1 = Poly(p, x, domain='EX') 

6377 n = p1.degree() 

6378 res = to_rational_coeffs(p1) 

6379 if not res: 

6380 return None 

6381 lc, r, t, g = res 

6382 factors = factor_list(g.as_expr()) 

6383 if lc: 

6384 c = simplify(factors[0]*lc*r**n) 

6385 r1 = simplify(1/r) 

6386 a = [] 

6387 for z in factors[1:][0]: 

6388 a.append((simplify(z[0].subs({x: x*r1})), z[1])) 

6389 else: 

6390 c = factors[0] 

6391 a = [] 

6392 for z in factors[1:][0]: 

6393 a.append((z[0].subs({x: x - t}), z[1])) 

6394 return (c, a) 

6395 

6396 

6397@public 

6398def sqf_list(f, *gens, **args): 

6399 """ 

6400 Compute a list of square-free factors of ``f``. 

6401 

6402 Examples 

6403 ======== 

6404 

6405 >>> from sympy import sqf_list 

6406 >>> from sympy.abc import x 

6407 

6408 >>> sqf_list(2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16) 

6409 (2, [(x + 1, 2), (x + 2, 3)]) 

6410 

6411 """ 

6412 return _generic_factor_list(f, gens, args, method='sqf') 

6413 

6414 

6415@public 

6416def sqf(f, *gens, **args): 

6417 """ 

6418 Compute square-free factorization of ``f``. 

6419 

6420 Examples 

6421 ======== 

6422 

6423 >>> from sympy import sqf 

6424 >>> from sympy.abc import x 

6425 

6426 >>> sqf(2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16) 

6427 2*(x + 1)**2*(x + 2)**3 

6428 

6429 """ 

6430 return _generic_factor(f, gens, args, method='sqf') 

6431 

6432 

6433@public 

6434def factor_list(f, *gens, **args): 

6435 """ 

6436 Compute a list of irreducible factors of ``f``. 

6437 

6438 Examples 

6439 ======== 

6440 

6441 >>> from sympy import factor_list 

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

6443 

6444 >>> factor_list(2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y) 

6445 (2, [(x + y, 1), (x**2 + 1, 2)]) 

6446 

6447 """ 

6448 return _generic_factor_list(f, gens, args, method='factor') 

6449 

6450 

6451@public 

6452def factor(f, *gens, deep=False, **args): 

6453 """ 

6454 Compute the factorization of expression, ``f``, into irreducibles. (To 

6455 factor an integer into primes, use ``factorint``.) 

6456 

6457 There two modes implemented: symbolic and formal. If ``f`` is not an 

6458 instance of :class:`Poly` and generators are not specified, then the 

6459 former mode is used. Otherwise, the formal mode is used. 

6460 

6461 In symbolic mode, :func:`factor` will traverse the expression tree and 

6462 factor its components without any prior expansion, unless an instance 

6463 of :class:`~.Add` is encountered (in this case formal factorization is 

6464 used). This way :func:`factor` can handle large or symbolic exponents. 

6465 

6466 By default, the factorization is computed over the rationals. To factor 

6467 over other domain, e.g. an algebraic or finite field, use appropriate 

6468 options: ``extension``, ``modulus`` or ``domain``. 

6469 

6470 Examples 

6471 ======== 

6472 

6473 >>> from sympy import factor, sqrt, exp 

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

6475 

6476 >>> factor(2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y) 

6477 2*(x + y)*(x**2 + 1)**2 

6478 

6479 >>> factor(x**2 + 1) 

6480 x**2 + 1 

6481 >>> factor(x**2 + 1, modulus=2) 

6482 (x + 1)**2 

6483 >>> factor(x**2 + 1, gaussian=True) 

6484 (x - I)*(x + I) 

6485 

6486 >>> factor(x**2 - 2, extension=sqrt(2)) 

6487 (x - sqrt(2))*(x + sqrt(2)) 

6488 

6489 >>> factor((x**2 - 1)/(x**2 + 4*x + 4)) 

6490 (x - 1)*(x + 1)/(x + 2)**2 

6491 >>> factor((x**2 + 4*x + 4)**10000000*(x**2 + 1)) 

6492 (x + 2)**20000000*(x**2 + 1) 

6493 

6494 By default, factor deals with an expression as a whole: 

6495 

6496 >>> eq = 2**(x**2 + 2*x + 1) 

6497 >>> factor(eq) 

6498 2**(x**2 + 2*x + 1) 

6499 

6500 If the ``deep`` flag is True then subexpressions will 

6501 be factored: 

6502 

6503 >>> factor(eq, deep=True) 

6504 2**((x + 1)**2) 

6505 

6506 If the ``fraction`` flag is False then rational expressions 

6507 will not be combined. By default it is True. 

6508 

6509 >>> factor(5*x + 3*exp(2 - 7*x), deep=True) 

6510 (5*x*exp(7*x) + 3*exp(2))*exp(-7*x) 

6511 >>> factor(5*x + 3*exp(2 - 7*x), deep=True, fraction=False) 

6512 5*x + 3*exp(2)*exp(-7*x) 

6513 

6514 See Also 

6515 ======== 

6516 sympy.ntheory.factor_.factorint 

6517 

6518 """ 

6519 f = sympify(f) 

6520 if deep: 

6521 def _try_factor(expr): 

6522 """ 

6523 Factor, but avoid changing the expression when unable to. 

6524 """ 

6525 fac = factor(expr, *gens, **args) 

6526 if fac.is_Mul or fac.is_Pow: 

6527 return fac 

6528 return expr 

6529 

6530 f = bottom_up(f, _try_factor) 

6531 # clean up any subexpressions that may have been expanded 

6532 # while factoring out a larger expression 

6533 partials = {} 

6534 muladd = f.atoms(Mul, Add) 

6535 for p in muladd: 

6536 fac = factor(p, *gens, **args) 

6537 if (fac.is_Mul or fac.is_Pow) and fac != p: 

6538 partials[p] = fac 

6539 return f.xreplace(partials) 

6540 

6541 try: 

6542 return _generic_factor(f, gens, args, method='factor') 

6543 except PolynomialError as msg: 

6544 if not f.is_commutative: 

6545 return factor_nc(f) 

6546 else: 

6547 raise PolynomialError(msg) 

6548 

6549 

6550@public 

6551def intervals(F, all=False, eps=None, inf=None, sup=None, strict=False, fast=False, sqf=False): 

6552 """ 

6553 Compute isolating intervals for roots of ``f``. 

6554 

6555 Examples 

6556 ======== 

6557 

6558 >>> from sympy import intervals 

6559 >>> from sympy.abc import x 

6560 

6561 >>> intervals(x**2 - 3) 

6562 [((-2, -1), 1), ((1, 2), 1)] 

6563 >>> intervals(x**2 - 3, eps=1e-2) 

6564 [((-26/15, -19/11), 1), ((19/11, 26/15), 1)] 

6565 

6566 """ 

6567 if not hasattr(F, '__iter__'): 

6568 try: 

6569 F = Poly(F) 

6570 except GeneratorsNeeded: 

6571 return [] 

6572 

6573 return F.intervals(all=all, eps=eps, inf=inf, sup=sup, fast=fast, sqf=sqf) 

6574 else: 

6575 polys, opt = parallel_poly_from_expr(F, domain='QQ') 

6576 

6577 if len(opt.gens) > 1: 

6578 raise MultivariatePolynomialError 

6579 

6580 for i, poly in enumerate(polys): 

6581 polys[i] = poly.rep.rep 

6582 

6583 if eps is not None: 

6584 eps = opt.domain.convert(eps) 

6585 

6586 if eps <= 0: 

6587 raise ValueError("'eps' must be a positive rational") 

6588 

6589 if inf is not None: 

6590 inf = opt.domain.convert(inf) 

6591 if sup is not None: 

6592 sup = opt.domain.convert(sup) 

6593 

6594 intervals = dup_isolate_real_roots_list(polys, opt.domain, 

6595 eps=eps, inf=inf, sup=sup, strict=strict, fast=fast) 

6596 

6597 result = [] 

6598 

6599 for (s, t), indices in intervals: 

6600 s, t = opt.domain.to_sympy(s), opt.domain.to_sympy(t) 

6601 result.append(((s, t), indices)) 

6602 

6603 return result 

6604 

6605 

6606@public 

6607def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False): 

6608 """ 

6609 Refine an isolating interval of a root to the given precision. 

6610 

6611 Examples 

6612 ======== 

6613 

6614 >>> from sympy import refine_root 

6615 >>> from sympy.abc import x 

6616 

6617 >>> refine_root(x**2 - 3, 1, 2, eps=1e-2) 

6618 (19/11, 26/15) 

6619 

6620 """ 

6621 try: 

6622 F = Poly(f) 

6623 if not isinstance(f, Poly) and not F.gen.is_Symbol: 

6624 # root of sin(x) + 1 is -1 but when someone 

6625 # passes an Expr instead of Poly they may not expect 

6626 # that the generator will be sin(x), not x 

6627 raise PolynomialError("generator must be a Symbol") 

6628 except GeneratorsNeeded: 

6629 raise PolynomialError( 

6630 "Cannot refine a root of %s, not a polynomial" % f) 

6631 

6632 return F.refine_root(s, t, eps=eps, steps=steps, fast=fast, check_sqf=check_sqf) 

6633 

6634 

6635@public 

6636def count_roots(f, inf=None, sup=None): 

6637 """ 

6638 Return the number of roots of ``f`` in ``[inf, sup]`` interval. 

6639 

6640 If one of ``inf`` or ``sup`` is complex, it will return the number of roots 

6641 in the complex rectangle with corners at ``inf`` and ``sup``. 

6642 

6643 Examples 

6644 ======== 

6645 

6646 >>> from sympy import count_roots, I 

6647 >>> from sympy.abc import x 

6648 

6649 >>> count_roots(x**4 - 4, -3, 3) 

6650 2 

6651 >>> count_roots(x**4 - 4, 0, 1 + 3*I) 

6652 1 

6653 

6654 """ 

6655 try: 

6656 F = Poly(f, greedy=False) 

6657 if not isinstance(f, Poly) and not F.gen.is_Symbol: 

6658 # root of sin(x) + 1 is -1 but when someone 

6659 # passes an Expr instead of Poly they may not expect 

6660 # that the generator will be sin(x), not x 

6661 raise PolynomialError("generator must be a Symbol") 

6662 except GeneratorsNeeded: 

6663 raise PolynomialError("Cannot count roots of %s, not a polynomial" % f) 

6664 

6665 return F.count_roots(inf=inf, sup=sup) 

6666 

6667 

6668@public 

6669def real_roots(f, multiple=True): 

6670 """ 

6671 Return a list of real roots with multiplicities of ``f``. 

6672 

6673 Examples 

6674 ======== 

6675 

6676 >>> from sympy import real_roots 

6677 >>> from sympy.abc import x 

6678 

6679 >>> real_roots(2*x**3 - 7*x**2 + 4*x + 4) 

6680 [-1/2, 2, 2] 

6681 """ 

6682 try: 

6683 F = Poly(f, greedy=False) 

6684 if not isinstance(f, Poly) and not F.gen.is_Symbol: 

6685 # root of sin(x) + 1 is -1 but when someone 

6686 # passes an Expr instead of Poly they may not expect 

6687 # that the generator will be sin(x), not x 

6688 raise PolynomialError("generator must be a Symbol") 

6689 except GeneratorsNeeded: 

6690 raise PolynomialError( 

6691 "Cannot compute real roots of %s, not a polynomial" % f) 

6692 

6693 return F.real_roots(multiple=multiple) 

6694 

6695 

6696@public 

6697def nroots(f, n=15, maxsteps=50, cleanup=True): 

6698 """ 

6699 Compute numerical approximations of roots of ``f``. 

6700 

6701 Examples 

6702 ======== 

6703 

6704 >>> from sympy import nroots 

6705 >>> from sympy.abc import x 

6706 

6707 >>> nroots(x**2 - 3, n=15) 

6708 [-1.73205080756888, 1.73205080756888] 

6709 >>> nroots(x**2 - 3, n=30) 

6710 [-1.73205080756887729352744634151, 1.73205080756887729352744634151] 

6711 

6712 """ 

6713 try: 

6714 F = Poly(f, greedy=False) 

6715 if not isinstance(f, Poly) and not F.gen.is_Symbol: 

6716 # root of sin(x) + 1 is -1 but when someone 

6717 # passes an Expr instead of Poly they may not expect 

6718 # that the generator will be sin(x), not x 

6719 raise PolynomialError("generator must be a Symbol") 

6720 except GeneratorsNeeded: 

6721 raise PolynomialError( 

6722 "Cannot compute numerical roots of %s, not a polynomial" % f) 

6723 

6724 return F.nroots(n=n, maxsteps=maxsteps, cleanup=cleanup) 

6725 

6726 

6727@public 

6728def ground_roots(f, *gens, **args): 

6729 """ 

6730 Compute roots of ``f`` by factorization in the ground domain. 

6731 

6732 Examples 

6733 ======== 

6734 

6735 >>> from sympy import ground_roots 

6736 >>> from sympy.abc import x 

6737 

6738 >>> ground_roots(x**6 - 4*x**4 + 4*x**3 - x**2) 

6739 {0: 2, 1: 2} 

6740 

6741 """ 

6742 options.allowed_flags(args, []) 

6743 

6744 try: 

6745 F, opt = poly_from_expr(f, *gens, **args) 

6746 if not isinstance(f, Poly) and not F.gen.is_Symbol: 

6747 # root of sin(x) + 1 is -1 but when someone 

6748 # passes an Expr instead of Poly they may not expect 

6749 # that the generator will be sin(x), not x 

6750 raise PolynomialError("generator must be a Symbol") 

6751 except PolificationFailed as exc: 

6752 raise ComputationFailed('ground_roots', 1, exc) 

6753 

6754 return F.ground_roots() 

6755 

6756 

6757@public 

6758def nth_power_roots_poly(f, n, *gens, **args): 

6759 """ 

6760 Construct a polynomial with n-th powers of roots of ``f``. 

6761 

6762 Examples 

6763 ======== 

6764 

6765 >>> from sympy import nth_power_roots_poly, factor, roots 

6766 >>> from sympy.abc import x 

6767 

6768 >>> f = x**4 - x**2 + 1 

6769 >>> g = factor(nth_power_roots_poly(f, 2)) 

6770 

6771 >>> g 

6772 (x**2 - x + 1)**2 

6773 

6774 >>> R_f = [ (r**2).expand() for r in roots(f) ] 

6775 >>> R_g = roots(g).keys() 

6776 

6777 >>> set(R_f) == set(R_g) 

6778 True 

6779 

6780 """ 

6781 options.allowed_flags(args, []) 

6782 

6783 try: 

6784 F, opt = poly_from_expr(f, *gens, **args) 

6785 if not isinstance(f, Poly) and not F.gen.is_Symbol: 

6786 # root of sin(x) + 1 is -1 but when someone 

6787 # passes an Expr instead of Poly they may not expect 

6788 # that the generator will be sin(x), not x 

6789 raise PolynomialError("generator must be a Symbol") 

6790 except PolificationFailed as exc: 

6791 raise ComputationFailed('nth_power_roots_poly', 1, exc) 

6792 

6793 result = F.nth_power_roots_poly(n) 

6794 

6795 if not opt.polys: 

6796 return result.as_expr() 

6797 else: 

6798 return result 

6799 

6800 

6801@public 

6802def cancel(f, *gens, _signsimp=True, **args): 

6803 """ 

6804 Cancel common factors in a rational function ``f``. 

6805 

6806 Examples 

6807 ======== 

6808 

6809 >>> from sympy import cancel, sqrt, Symbol, together 

6810 >>> from sympy.abc import x 

6811 >>> A = Symbol('A', commutative=False) 

6812 

6813 >>> cancel((2*x**2 - 2)/(x**2 - 2*x + 1)) 

6814 (2*x + 2)/(x - 1) 

6815 >>> cancel((sqrt(3) + sqrt(15)*A)/(sqrt(2) + sqrt(10)*A)) 

6816 sqrt(6)/2 

6817 

6818 Note: due to automatic distribution of Rationals, a sum divided by an integer 

6819 will appear as a sum. To recover a rational form use `together` on the result: 

6820 

6821 >>> cancel(x/2 + 1) 

6822 x/2 + 1 

6823 >>> together(_) 

6824 (x + 2)/2 

6825 """ 

6826 from sympy.simplify.simplify import signsimp 

6827 from sympy.polys.rings import sring 

6828 options.allowed_flags(args, ['polys']) 

6829 

6830 f = sympify(f) 

6831 if _signsimp: 

6832 f = signsimp(f) 

6833 opt = {} 

6834 if 'polys' in args: 

6835 opt['polys'] = args['polys'] 

6836 

6837 if not isinstance(f, (tuple, Tuple)): 

6838 if f.is_Number or isinstance(f, Relational) or not isinstance(f, Expr): 

6839 return f 

6840 f = factor_terms(f, radical=True) 

6841 p, q = f.as_numer_denom() 

6842 

6843 elif len(f) == 2: 

6844 p, q = f 

6845 if isinstance(p, Poly) and isinstance(q, Poly): 

6846 opt['gens'] = p.gens 

6847 opt['domain'] = p.domain 

6848 opt['polys'] = opt.get('polys', True) 

6849 p, q = p.as_expr(), q.as_expr() 

6850 elif isinstance(f, Tuple): 

6851 return factor_terms(f) 

6852 else: 

6853 raise ValueError('unexpected argument: %s' % f) 

6854 

6855 from sympy.functions.elementary.piecewise import Piecewise 

6856 try: 

6857 if f.has(Piecewise): 

6858 raise PolynomialError() 

6859 R, (F, G) = sring((p, q), *gens, **args) 

6860 if not R.ngens: 

6861 if not isinstance(f, (tuple, Tuple)): 

6862 return f.expand() 

6863 else: 

6864 return S.One, p, q 

6865 except PolynomialError as msg: 

6866 if f.is_commutative and not f.has(Piecewise): 

6867 raise PolynomialError(msg) 

6868 # Handling of noncommutative and/or piecewise expressions 

6869 if f.is_Add or f.is_Mul: 

6870 c, nc = sift(f.args, lambda x: 

6871 x.is_commutative is True and not x.has(Piecewise), 

6872 binary=True) 

6873 nc = [cancel(i) for i in nc] 

6874 return f.func(cancel(f.func(*c)), *nc) 

6875 else: 

6876 reps = [] 

6877 pot = preorder_traversal(f) 

6878 next(pot) 

6879 for e in pot: 

6880 # XXX: This should really skip anything that's not Expr. 

6881 if isinstance(e, (tuple, Tuple, BooleanAtom)): 

6882 continue 

6883 try: 

6884 reps.append((e, cancel(e))) 

6885 pot.skip() # this was handled successfully 

6886 except NotImplementedError: 

6887 pass 

6888 return f.xreplace(dict(reps)) 

6889 

6890 c, (P, Q) = 1, F.cancel(G) 

6891 if opt.get('polys', False) and 'gens' not in opt: 

6892 opt['gens'] = R.symbols 

6893 

6894 if not isinstance(f, (tuple, Tuple)): 

6895 return c*(P.as_expr()/Q.as_expr()) 

6896 else: 

6897 P, Q = P.as_expr(), Q.as_expr() 

6898 if not opt.get('polys', False): 

6899 return c, P, Q 

6900 else: 

6901 return c, Poly(P, *gens, **opt), Poly(Q, *gens, **opt) 

6902 

6903 

6904@public 

6905def reduced(f, G, *gens, **args): 

6906 """ 

6907 Reduces a polynomial ``f`` modulo a set of polynomials ``G``. 

6908 

6909 Given a polynomial ``f`` and a set of polynomials ``G = (g_1, ..., g_n)``, 

6910 computes a set of quotients ``q = (q_1, ..., q_n)`` and the remainder ``r`` 

6911 such that ``f = q_1*g_1 + ... + q_n*g_n + r``, where ``r`` vanishes or ``r`` 

6912 is a completely reduced polynomial with respect to ``G``. 

6913 

6914 Examples 

6915 ======== 

6916 

6917 >>> from sympy import reduced 

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

6919 

6920 >>> reduced(2*x**4 + y**2 - x**2 + y**3, [x**3 - x, y**3 - y]) 

6921 ([2*x, 1], x**2 + y**2 + y) 

6922 

6923 """ 

6924 options.allowed_flags(args, ['polys', 'auto']) 

6925 

6926 try: 

6927 polys, opt = parallel_poly_from_expr([f] + list(G), *gens, **args) 

6928 except PolificationFailed as exc: 

6929 raise ComputationFailed('reduced', 0, exc) 

6930 

6931 domain = opt.domain 

6932 retract = False 

6933 

6934 if opt.auto and domain.is_Ring and not domain.is_Field: 

6935 opt = opt.clone({"domain": domain.get_field()}) 

6936 retract = True 

6937 

6938 from sympy.polys.rings import xring 

6939 _ring, _ = xring(opt.gens, opt.domain, opt.order) 

6940 

6941 for i, poly in enumerate(polys): 

6942 poly = poly.set_domain(opt.domain).rep.to_dict() 

6943 polys[i] = _ring.from_dict(poly) 

6944 

6945 Q, r = polys[0].div(polys[1:]) 

6946 

6947 Q = [Poly._from_dict(dict(q), opt) for q in Q] 

6948 r = Poly._from_dict(dict(r), opt) 

6949 

6950 if retract: 

6951 try: 

6952 _Q, _r = [q.to_ring() for q in Q], r.to_ring() 

6953 except CoercionFailed: 

6954 pass 

6955 else: 

6956 Q, r = _Q, _r 

6957 

6958 if not opt.polys: 

6959 return [q.as_expr() for q in Q], r.as_expr() 

6960 else: 

6961 return Q, r 

6962 

6963 

6964@public 

6965def groebner(F, *gens, **args): 

6966 """ 

6967 Computes the reduced Groebner basis for a set of polynomials. 

6968 

6969 Use the ``order`` argument to set the monomial ordering that will be 

6970 used to compute the basis. Allowed orders are ``lex``, ``grlex`` and 

6971 ``grevlex``. If no order is specified, it defaults to ``lex``. 

6972 

6973 For more information on Groebner bases, see the references and the docstring 

6974 of :func:`~.solve_poly_system`. 

6975 

6976 Examples 

6977 ======== 

6978 

6979 Example taken from [1]. 

6980 

6981 >>> from sympy import groebner 

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

6983 

6984 >>> F = [x*y - 2*y, 2*y**2 - x**2] 

6985 

6986 >>> groebner(F, x, y, order='lex') 

6987 GroebnerBasis([x**2 - 2*y**2, x*y - 2*y, y**3 - 2*y], x, y, 

6988 domain='ZZ', order='lex') 

6989 >>> groebner(F, x, y, order='grlex') 

6990 GroebnerBasis([y**3 - 2*y, x**2 - 2*y**2, x*y - 2*y], x, y, 

6991 domain='ZZ', order='grlex') 

6992 >>> groebner(F, x, y, order='grevlex') 

6993 GroebnerBasis([y**3 - 2*y, x**2 - 2*y**2, x*y - 2*y], x, y, 

6994 domain='ZZ', order='grevlex') 

6995 

6996 By default, an improved implementation of the Buchberger algorithm is 

6997 used. Optionally, an implementation of the F5B algorithm can be used. The 

6998 algorithm can be set using the ``method`` flag or with the 

6999 :func:`sympy.polys.polyconfig.setup` function. 

7000 

7001 >>> F = [x**2 - x - 1, (2*x - 1) * y - (x**10 - (1 - x)**10)] 

7002 

7003 >>> groebner(F, x, y, method='buchberger') 

7004 GroebnerBasis([x**2 - x - 1, y - 55], x, y, domain='ZZ', order='lex') 

7005 >>> groebner(F, x, y, method='f5b') 

7006 GroebnerBasis([x**2 - x - 1, y - 55], x, y, domain='ZZ', order='lex') 

7007 

7008 References 

7009 ========== 

7010 

7011 1. [Buchberger01]_ 

7012 2. [Cox97]_ 

7013 

7014 """ 

7015 return GroebnerBasis(F, *gens, **args) 

7016 

7017 

7018@public 

7019def is_zero_dimensional(F, *gens, **args): 

7020 """ 

7021 Checks if the ideal generated by a Groebner basis is zero-dimensional. 

7022 

7023 The algorithm checks if the set of monomials not divisible by the 

7024 leading monomial of any element of ``F`` is bounded. 

7025 

7026 References 

7027 ========== 

7028 

7029 David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties and 

7030 Algorithms, 3rd edition, p. 230 

7031 

7032 """ 

7033 return GroebnerBasis(F, *gens, **args).is_zero_dimensional 

7034 

7035 

7036@public 

7037class GroebnerBasis(Basic): 

7038 """Represents a reduced Groebner basis. """ 

7039 

7040 def __new__(cls, F, *gens, **args): 

7041 """Compute a reduced Groebner basis for a system of polynomials. """ 

7042 options.allowed_flags(args, ['polys', 'method']) 

7043 

7044 try: 

7045 polys, opt = parallel_poly_from_expr(F, *gens, **args) 

7046 except PolificationFailed as exc: 

7047 raise ComputationFailed('groebner', len(F), exc) 

7048 

7049 from sympy.polys.rings import PolyRing 

7050 ring = PolyRing(opt.gens, opt.domain, opt.order) 

7051 

7052 polys = [ring.from_dict(poly.rep.to_dict()) for poly in polys if poly] 

7053 

7054 G = _groebner(polys, ring, method=opt.method) 

7055 G = [Poly._from_dict(g, opt) for g in G] 

7056 

7057 return cls._new(G, opt) 

7058 

7059 @classmethod 

7060 def _new(cls, basis, options): 

7061 obj = Basic.__new__(cls) 

7062 

7063 obj._basis = tuple(basis) 

7064 obj._options = options 

7065 

7066 return obj 

7067 

7068 @property 

7069 def args(self): 

7070 basis = (p.as_expr() for p in self._basis) 

7071 return (Tuple(*basis), Tuple(*self._options.gens)) 

7072 

7073 @property 

7074 def exprs(self): 

7075 return [poly.as_expr() for poly in self._basis] 

7076 

7077 @property 

7078 def polys(self): 

7079 return list(self._basis) 

7080 

7081 @property 

7082 def gens(self): 

7083 return self._options.gens 

7084 

7085 @property 

7086 def domain(self): 

7087 return self._options.domain 

7088 

7089 @property 

7090 def order(self): 

7091 return self._options.order 

7092 

7093 def __len__(self): 

7094 return len(self._basis) 

7095 

7096 def __iter__(self): 

7097 if self._options.polys: 

7098 return iter(self.polys) 

7099 else: 

7100 return iter(self.exprs) 

7101 

7102 def __getitem__(self, item): 

7103 if self._options.polys: 

7104 basis = self.polys 

7105 else: 

7106 basis = self.exprs 

7107 

7108 return basis[item] 

7109 

7110 def __hash__(self): 

7111 return hash((self._basis, tuple(self._options.items()))) 

7112 

7113 def __eq__(self, other): 

7114 if isinstance(other, self.__class__): 

7115 return self._basis == other._basis and self._options == other._options 

7116 elif iterable(other): 

7117 return self.polys == list(other) or self.exprs == list(other) 

7118 else: 

7119 return False 

7120 

7121 def __ne__(self, other): 

7122 return not self == other 

7123 

7124 @property 

7125 def is_zero_dimensional(self): 

7126 """ 

7127 Checks if the ideal generated by a Groebner basis is zero-dimensional. 

7128 

7129 The algorithm checks if the set of monomials not divisible by the 

7130 leading monomial of any element of ``F`` is bounded. 

7131 

7132 References 

7133 ========== 

7134 

7135 David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties and 

7136 Algorithms, 3rd edition, p. 230 

7137 

7138 """ 

7139 def single_var(monomial): 

7140 return sum(map(bool, monomial)) == 1 

7141 

7142 exponents = Monomial([0]*len(self.gens)) 

7143 order = self._options.order 

7144 

7145 for poly in self.polys: 

7146 monomial = poly.LM(order=order) 

7147 

7148 if single_var(monomial): 

7149 exponents *= monomial 

7150 

7151 # If any element of the exponents vector is zero, then there's 

7152 # a variable for which there's no degree bound and the ideal 

7153 # generated by this Groebner basis isn't zero-dimensional. 

7154 return all(exponents) 

7155 

7156 def fglm(self, order): 

7157 """ 

7158 Convert a Groebner basis from one ordering to another. 

7159 

7160 The FGLM algorithm converts reduced Groebner bases of zero-dimensional 

7161 ideals from one ordering to another. This method is often used when it 

7162 is infeasible to compute a Groebner basis with respect to a particular 

7163 ordering directly. 

7164 

7165 Examples 

7166 ======== 

7167 

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

7169 >>> from sympy import groebner 

7170 

7171 >>> F = [x**2 - 3*y - x + 1, y**2 - 2*x + y - 1] 

7172 >>> G = groebner(F, x, y, order='grlex') 

7173 

7174 >>> list(G.fglm('lex')) 

7175 [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7] 

7176 >>> list(groebner(F, x, y, order='lex')) 

7177 [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7] 

7178 

7179 References 

7180 ========== 

7181 

7182 .. [1] J.C. Faugere, P. Gianni, D. Lazard, T. Mora (1994). Efficient 

7183 Computation of Zero-dimensional Groebner Bases by Change of 

7184 Ordering 

7185 

7186 """ 

7187 opt = self._options 

7188 

7189 src_order = opt.order 

7190 dst_order = monomial_key(order) 

7191 

7192 if src_order == dst_order: 

7193 return self 

7194 

7195 if not self.is_zero_dimensional: 

7196 raise NotImplementedError("Cannot convert Groebner bases of ideals with positive dimension") 

7197 

7198 polys = list(self._basis) 

7199 domain = opt.domain 

7200 

7201 opt = opt.clone({ 

7202 "domain": domain.get_field(), 

7203 "order": dst_order, 

7204 }) 

7205 

7206 from sympy.polys.rings import xring 

7207 _ring, _ = xring(opt.gens, opt.domain, src_order) 

7208 

7209 for i, poly in enumerate(polys): 

7210 poly = poly.set_domain(opt.domain).rep.to_dict() 

7211 polys[i] = _ring.from_dict(poly) 

7212 

7213 G = matrix_fglm(polys, _ring, dst_order) 

7214 G = [Poly._from_dict(dict(g), opt) for g in G] 

7215 

7216 if not domain.is_Field: 

7217 G = [g.clear_denoms(convert=True)[1] for g in G] 

7218 opt.domain = domain 

7219 

7220 return self._new(G, opt) 

7221 

7222 def reduce(self, expr, auto=True): 

7223 """ 

7224 Reduces a polynomial modulo a Groebner basis. 

7225 

7226 Given a polynomial ``f`` and a set of polynomials ``G = (g_1, ..., g_n)``, 

7227 computes a set of quotients ``q = (q_1, ..., q_n)`` and the remainder ``r`` 

7228 such that ``f = q_1*f_1 + ... + q_n*f_n + r``, where ``r`` vanishes or ``r`` 

7229 is a completely reduced polynomial with respect to ``G``. 

7230 

7231 Examples 

7232 ======== 

7233 

7234 >>> from sympy import groebner, expand 

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

7236 

7237 >>> f = 2*x**4 - x**2 + y**3 + y**2 

7238 >>> G = groebner([x**3 - x, y**3 - y]) 

7239 

7240 >>> G.reduce(f) 

7241 ([2*x, 1], x**2 + y**2 + y) 

7242 >>> Q, r = _ 

7243 

7244 >>> expand(sum(q*g for q, g in zip(Q, G)) + r) 

7245 2*x**4 - x**2 + y**3 + y**2 

7246 >>> _ == f 

7247 True 

7248 

7249 """ 

7250 poly = Poly._from_expr(expr, self._options) 

7251 polys = [poly] + list(self._basis) 

7252 

7253 opt = self._options 

7254 domain = opt.domain 

7255 

7256 retract = False 

7257 

7258 if auto and domain.is_Ring and not domain.is_Field: 

7259 opt = opt.clone({"domain": domain.get_field()}) 

7260 retract = True 

7261 

7262 from sympy.polys.rings import xring 

7263 _ring, _ = xring(opt.gens, opt.domain, opt.order) 

7264 

7265 for i, poly in enumerate(polys): 

7266 poly = poly.set_domain(opt.domain).rep.to_dict() 

7267 polys[i] = _ring.from_dict(poly) 

7268 

7269 Q, r = polys[0].div(polys[1:]) 

7270 

7271 Q = [Poly._from_dict(dict(q), opt) for q in Q] 

7272 r = Poly._from_dict(dict(r), opt) 

7273 

7274 if retract: 

7275 try: 

7276 _Q, _r = [q.to_ring() for q in Q], r.to_ring() 

7277 except CoercionFailed: 

7278 pass 

7279 else: 

7280 Q, r = _Q, _r 

7281 

7282 if not opt.polys: 

7283 return [q.as_expr() for q in Q], r.as_expr() 

7284 else: 

7285 return Q, r 

7286 

7287 def contains(self, poly): 

7288 """ 

7289 Check if ``poly`` belongs the ideal generated by ``self``. 

7290 

7291 Examples 

7292 ======== 

7293 

7294 >>> from sympy import groebner 

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

7296 

7297 >>> f = 2*x**3 + y**3 + 3*y 

7298 >>> G = groebner([x**2 + y**2 - 1, x*y - 2]) 

7299 

7300 >>> G.contains(f) 

7301 True 

7302 >>> G.contains(f + 1) 

7303 False 

7304 

7305 """ 

7306 return self.reduce(poly)[1] == 0 

7307 

7308 

7309@public 

7310def poly(expr, *gens, **args): 

7311 """ 

7312 Efficiently transform an expression into a polynomial. 

7313 

7314 Examples 

7315 ======== 

7316 

7317 >>> from sympy import poly 

7318 >>> from sympy.abc import x 

7319 

7320 >>> poly(x*(x**2 + x - 1)**2) 

7321 Poly(x**5 + 2*x**4 - x**3 - 2*x**2 + x, x, domain='ZZ') 

7322 

7323 """ 

7324 options.allowed_flags(args, []) 

7325 

7326 def _poly(expr, opt): 

7327 terms, poly_terms = [], [] 

7328 

7329 for term in Add.make_args(expr): 

7330 factors, poly_factors = [], [] 

7331 

7332 for factor in Mul.make_args(term): 

7333 if factor.is_Add: 

7334 poly_factors.append(_poly(factor, opt)) 

7335 elif factor.is_Pow and factor.base.is_Add and \ 

7336 factor.exp.is_Integer and factor.exp >= 0: 

7337 poly_factors.append( 

7338 _poly(factor.base, opt).pow(factor.exp)) 

7339 else: 

7340 factors.append(factor) 

7341 

7342 if not poly_factors: 

7343 terms.append(term) 

7344 else: 

7345 product = poly_factors[0] 

7346 

7347 for factor in poly_factors[1:]: 

7348 product = product.mul(factor) 

7349 

7350 if factors: 

7351 factor = Mul(*factors) 

7352 

7353 if factor.is_Number: 

7354 product = product.mul(factor) 

7355 else: 

7356 product = product.mul(Poly._from_expr(factor, opt)) 

7357 

7358 poly_terms.append(product) 

7359 

7360 if not poly_terms: 

7361 result = Poly._from_expr(expr, opt) 

7362 else: 

7363 result = poly_terms[0] 

7364 

7365 for term in poly_terms[1:]: 

7366 result = result.add(term) 

7367 

7368 if terms: 

7369 term = Add(*terms) 

7370 

7371 if term.is_Number: 

7372 result = result.add(term) 

7373 else: 

7374 result = result.add(Poly._from_expr(term, opt)) 

7375 

7376 return result.reorder(*opt.get('gens', ()), **args) 

7377 

7378 expr = sympify(expr) 

7379 

7380 if expr.is_Poly: 

7381 return Poly(expr, *gens, **args) 

7382 

7383 if 'expand' not in args: 

7384 args['expand'] = False 

7385 

7386 opt = options.build_options(gens, args) 

7387 

7388 return _poly(expr, opt) 

7389 

7390 

7391def named_poly(n, f, K, name, x, polys): 

7392 r"""Common interface to the low-level polynomial generating functions 

7393 in orthopolys and appellseqs. 

7394 

7395 Parameters 

7396 ========== 

7397 

7398 n : int 

7399 Index of the polynomial, which may or may not equal its degree. 

7400 f : callable 

7401 Low-level generating function to use. 

7402 K : Domain or None 

7403 Domain in which to perform the computations. If None, use the smallest 

7404 field containing the rationals and the extra parameters of x (see below). 

7405 name : str 

7406 Name of an arbitrary individual polynomial in the sequence generated 

7407 by f, only used in the error message for invalid n. 

7408 x : seq 

7409 The first element of this argument is the main variable of all 

7410 polynomials in this sequence. Any further elements are extra 

7411 parameters required by f. 

7412 polys : bool, optional 

7413 If True, return a Poly, otherwise (default) return an expression. 

7414 """ 

7415 if n < 0: 

7416 raise ValueError("Cannot generate %s of index %s" % (name, n)) 

7417 head, tail = x[0], x[1:] 

7418 if K is None: 

7419 K, tail = construct_domain(tail, field=True) 

7420 poly = DMP(f(int(n), *tail, K), K) 

7421 if head is None: 

7422 poly = PurePoly.new(poly, Dummy('x')) 

7423 else: 

7424 poly = Poly.new(poly, head) 

7425 return poly if polys else poly.as_expr()