Coverage for /usr/lib/python3/dist-packages/sympy/polys/rings.py: 25%

1"""Sparse polynomial rings. """

3from __future__ import annotations

4from typing import Any

6from operator import add, mul, lt, le, gt, ge

7from functools import reduce

8from types import GeneratorType

10from sympy.core.expr import Expr

11from sympy.core.numbers import igcd, oo

12from sympy.core.symbol import Symbol, symbols as _symbols

13from sympy.core.sympify import CantSympify, sympify

14from sympy.ntheory.multinomial import multinomial_coefficients

15from sympy.polys.compatibility import IPolys

16from sympy.polys.constructor import construct_domain

17from sympy.polys.densebasic import dmp_to_dict, dmp_from_dict

18from sympy.polys.domains.domainelement import DomainElement

19from sympy.polys.domains.polynomialring import PolynomialRing

20from sympy.polys.heuristicgcd import heugcd

21from sympy.polys.monomials import MonomialOps

22from sympy.polys.orderings import lex

23from sympy.polys.polyerrors import (

24 CoercionFailed, GeneratorsError,

25 ExactQuotientFailed, MultivariatePolynomialError)

26from sympy.polys.polyoptions import (Domain as DomainOpt,

27 Order as OrderOpt, build_options)

28from sympy.polys.polyutils import (expr_from_dict, _dict_reorder,

29 _parallel_dict_from_expr)

30from sympy.printing.defaults import DefaultPrinting

31from sympy.utilities import public, subsets

32from sympy.utilities.iterables import is_sequence

33from sympy.utilities.magic import pollute

35@public

36def ring(symbols, domain, order=lex):

37 """Construct a polynomial ring returning ``(ring, x_1, ..., x_n)``.

39 Parameters

40 ==========

42 symbols : str

43 Symbol/Expr or sequence of str, Symbol/Expr (non-empty)

44 domain : :class:`~.Domain` or coercible

45 order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``

47 Examples

48 ========

50 >>> from sympy.polys.rings import ring

51 >>> from sympy.polys.domains import ZZ

52 >>> from sympy.polys.orderings import lex

54 >>> R, x, y, z = ring("x,y,z", ZZ, lex)

55 >>> R

56 Polynomial ring in x, y, z over ZZ with lex order

57 >>> x + y + z

58 x + y + z

59 >>> type(_)

60 <class 'sympy.polys.rings.PolyElement'>

62 """

63 _ring = PolyRing(symbols, domain, order)

64 return (_ring,) + _ring.gens

66@public

67def xring(symbols, domain, order=lex):

68 """Construct a polynomial ring returning ``(ring, (x_1, ..., x_n))``.

70 Parameters

71 ==========

73 symbols : str

74 Symbol/Expr or sequence of str, Symbol/Expr (non-empty)

75 domain : :class:`~.Domain` or coercible

76 order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``

78 Examples

79 ========

81 >>> from sympy.polys.rings import xring

82 >>> from sympy.polys.domains import ZZ

83 >>> from sympy.polys.orderings import lex

85 >>> R, (x, y, z) = xring("x,y,z", ZZ, lex)

86 >>> R

87 Polynomial ring in x, y, z over ZZ with lex order

88 >>> x + y + z

89 x + y + z

90 >>> type(_)

91 <class 'sympy.polys.rings.PolyElement'>

93 """

94 _ring = PolyRing(symbols, domain, order)

95 return (_ring, _ring.gens)

97@public

98def vring(symbols, domain, order=lex):

99 """Construct a polynomial ring and inject ``x_1, ..., x_n`` into the global namespace.

100

101 Parameters

102 ==========

103

104 symbols : str

105 Symbol/Expr or sequence of str, Symbol/Expr (non-empty)

106 domain : :class:`~.Domain` or coercible

107 order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``

108

109 Examples

110 ========

111

112 >>> from sympy.polys.rings import vring

113 >>> from sympy.polys.domains import ZZ

114 >>> from sympy.polys.orderings import lex

115

116 >>> vring("x,y,z", ZZ, lex)

117 Polynomial ring in x, y, z over ZZ with lex order

118 >>> x + y + z # noqa:

119 x + y + z

120 >>> type(_)

121 <class 'sympy.polys.rings.PolyElement'>

122

123 """

124 _ring = PolyRing(symbols, domain, order)

125 pollute([ sym.name for sym in _ring.symbols ], _ring.gens)

126 return _ring

127

128@public

129def sring(exprs, *symbols, **options):

130 """Construct a ring deriving generators and domain from options and input expressions.

131

132 Parameters

133 ==========

134

135 exprs : :class:`~.Expr` or sequence of :class:`~.Expr` (sympifiable)

136 symbols : sequence of :class:`~.Symbol`/:class:`~.Expr`

137 options : keyword arguments understood by :class:`~.Options`

138

139 Examples

140 ========

141

142 >>> from sympy import sring, symbols

143

144 >>> x, y, z = symbols("x,y,z")

145 >>> R, f = sring(x + 2*y + 3*z)

146 >>> R

147 Polynomial ring in x, y, z over ZZ with lex order

148 >>> f

149 x + 2*y + 3*z

150 >>> type(_)

151 <class 'sympy.polys.rings.PolyElement'>

152

153 """

154 single = False

155

156 if not is_sequence(exprs):

157 exprs, single = [exprs], True

158

159 exprs = list(map(sympify, exprs))

160 opt = build_options(symbols, options)

161

162 # TODO: rewrite this so that it doesn't use expand() (see poly()).

163 reps, opt = _parallel_dict_from_expr(exprs, opt)

164

165 if opt.domain is None:

166 coeffs = sum([ list(rep.values()) for rep in reps ], [])

167

168 opt.domain, coeffs_dom = construct_domain(coeffs, opt=opt)

169

170 coeff_map = dict(zip(coeffs, coeffs_dom))

171 reps = [{m: coeff_map[c] for m, c in rep.items()} for rep in reps]

172

173 _ring = PolyRing(opt.gens, opt.domain, opt.order)

174 polys = list(map(_ring.from_dict, reps))

175

176 if single:

177 return (_ring, polys[0])

178 else:

179 return (_ring, polys)

180

181def _parse_symbols(symbols):

182 if isinstance(symbols, str):

183 return _symbols(symbols, seq=True) if symbols else ()

184 elif isinstance(symbols, Expr):

185 return (symbols,)

186 elif is_sequence(symbols):

187 if all(isinstance(s, str) for s in symbols):

188 return _symbols(symbols)

189 elif all(isinstance(s, Expr) for s in symbols):

190 return symbols

191

192 raise GeneratorsError("expected a string, Symbol or expression or a non-empty sequence of strings, Symbols or expressions")

193

194_ring_cache: dict[Any, Any] = {}

195

196class PolyRing(DefaultPrinting, IPolys):

197 """Multivariate distributed polynomial ring. """

198

199 def __new__(cls, symbols, domain, order=lex):

200 symbols = tuple(_parse_symbols(symbols))

201 ngens = len(symbols)

202 domain = DomainOpt.preprocess(domain)

203 order = OrderOpt.preprocess(order)

204

205 _hash_tuple = (cls.__name__, symbols, ngens, domain, order)

206 obj = _ring_cache.get(_hash_tuple)

207

208 if obj is None:

209 if domain.is_Composite and set(symbols) & set(domain.symbols):

210 raise GeneratorsError("polynomial ring and it's ground domain share generators")

211

212 obj = object.__new__(cls)

213 obj._hash_tuple = _hash_tuple

214 obj._hash = hash(_hash_tuple)

215 obj.dtype = type("PolyElement", (PolyElement,), {"ring": obj})

216 obj.symbols = symbols

217 obj.ngens = ngens

218 obj.domain = domain

219 obj.order = order

220

221 obj.zero_monom = (0,)*ngens

222 obj.gens = obj._gens()

223 obj._gens_set = set(obj.gens)

224

225 obj._one = [(obj.zero_monom, domain.one)]

226

227 if ngens:

228 # These expect monomials in at least one variable

229 codegen = MonomialOps(ngens)

230 obj.monomial_mul = codegen.mul()

231 obj.monomial_pow = codegen.pow()

232 obj.monomial_mulpow = codegen.mulpow()

233 obj.monomial_ldiv = codegen.ldiv()

234 obj.monomial_div = codegen.div()

235 obj.monomial_lcm = codegen.lcm()

236 obj.monomial_gcd = codegen.gcd()

237 else:

238 monunit = lambda a, b: ()

239 obj.monomial_mul = monunit

240 obj.monomial_pow = monunit

241 obj.monomial_mulpow = lambda a, b, c: ()

242 obj.monomial_ldiv = monunit

243 obj.monomial_div = monunit

244 obj.monomial_lcm = monunit

245 obj.monomial_gcd = monunit

246

247

248 if order is lex:

249 obj.leading_expv = max

250 else:

251 obj.leading_expv = lambda f: max(f, key=order)

252

253 for symbol, generator in zip(obj.symbols, obj.gens):

254 if isinstance(symbol, Symbol):

255 name = symbol.name

256

257 if not hasattr(obj, name):

258 setattr(obj, name, generator)

259

260 _ring_cache[_hash_tuple] = obj

261

262 return obj

263

264 def _gens(self):

265 """Return a list of polynomial generators. """

266 one = self.domain.one

267 _gens = []

268 for i in range(self.ngens):

269 expv = self.monomial_basis(i)

270 poly = self.zero

271 poly[expv] = one

272 _gens.append(poly)

273 return tuple(_gens)

274

275 def __getnewargs__(self):

276 return (self.symbols, self.domain, self.order)

277

278 def __getstate__(self):

279 state = self.__dict__.copy()

280 del state["leading_expv"]

281

282 for key, value in state.items():

283 if key.startswith("monomial_"):

284 del state[key]

285

286 return state

287

288 def __hash__(self):

289 return self._hash

290

291 def __eq__(self, other):

292 return isinstance(other, PolyRing) and \

293 (self.symbols, self.domain, self.ngens, self.order) == \

294 (other.symbols, other.domain, other.ngens, other.order)

295

296 def __ne__(self, other):

297 return not self == other

298

299 def clone(self, symbols=None, domain=None, order=None):

300 return self.__class__(symbols or self.symbols, domain or self.domain, order or self.order)

301

302 def monomial_basis(self, i):

303 """Return the ith-basis element. """

304 basis = [0]*self.ngens

305 basis[i] = 1

306 return tuple(basis)

307

308 @property

309 def zero(self):

310 return self.dtype()

311

312 @property

313 def one(self):

314 return self.dtype(self._one)

315

316 def domain_new(self, element, orig_domain=None):

317 return self.domain.convert(element, orig_domain)

318

319 def ground_new(self, coeff):

320 return self.term_new(self.zero_monom, coeff)

321

322 def term_new(self, monom, coeff):

323 coeff = self.domain_new(coeff)

324 poly = self.zero

325 if coeff:

326 poly[monom] = coeff

327 return poly

328

329 def ring_new(self, element):

330 if isinstance(element, PolyElement):

331 if self == element.ring:

332 return element

333 elif isinstance(self.domain, PolynomialRing) and self.domain.ring == element.ring:

334 return self.ground_new(element)

335 else:

336 raise NotImplementedError("conversion")

337 elif isinstance(element, str):

338 raise NotImplementedError("parsing")

339 elif isinstance(element, dict):

340 return self.from_dict(element)

341 elif isinstance(element, list):

342 try:

343 return self.from_terms(element)

344 except ValueError:

345 return self.from_list(element)

346 elif isinstance(element, Expr):

347 return self.from_expr(element)

348 else:

349 return self.ground_new(element)

350

351 __call__ = ring_new

352

353 def from_dict(self, element, orig_domain=None):

354 domain_new = self.domain_new

355 poly = self.zero

356

357 for monom, coeff in element.items():

358 coeff = domain_new(coeff, orig_domain)

359 if coeff:

360 poly[monom] = coeff

361

362 return poly

363

364 def from_terms(self, element, orig_domain=None):

365 return self.from_dict(dict(element), orig_domain)

366

367 def from_list(self, element):

368 return self.from_dict(dmp_to_dict(element, self.ngens-1, self.domain))

369

370 def _rebuild_expr(self, expr, mapping):

371 domain = self.domain

372

373 def _rebuild(expr):

374 generator = mapping.get(expr)

375

376 if generator is not None:

377 return generator

378 elif expr.is_Add:

379 return reduce(add, list(map(_rebuild, expr.args)))

380 elif expr.is_Mul:

381 return reduce(mul, list(map(_rebuild, expr.args)))

382 else:

383 # XXX: Use as_base_exp() to handle Pow(x, n) and also exp(n)

384 # XXX: E can be a generator e.g. sring([exp(2)]) -> ZZ[E]

385 base, exp = expr.as_base_exp()

386 if exp.is_Integer and exp > 1:

387 return _rebuild(base)**int(exp)

388 else:

389 return self.ground_new(domain.convert(expr))

390

391 return _rebuild(sympify(expr))

392

393 def from_expr(self, expr):

394 mapping = dict(list(zip(self.symbols, self.gens)))

395

396 try:

397 poly = self._rebuild_expr(expr, mapping)

398 except CoercionFailed:

399 raise ValueError("expected an expression convertible to a polynomial in %s, got %s" % (self, expr))

400 else:

401 return self.ring_new(poly)

402

403 def index(self, gen):

404 """Compute index of ``gen`` in ``self.gens``. """

405 if gen is None:

406 if self.ngens:

407 i = 0

408 else:

409 i = -1 # indicate impossible choice

410 elif isinstance(gen, int):

411 i = gen

412

413 if 0 <= i and i < self.ngens:

414 pass

415 elif -self.ngens <= i and i <= -1:

416 i = -i - 1

417 else:

418 raise ValueError("invalid generator index: %s" % gen)

419 elif isinstance(gen, self.dtype):

420 try:

421 i = self.gens.index(gen)

422 except ValueError:

423 raise ValueError("invalid generator: %s" % gen)

424 elif isinstance(gen, str):

425 try:

426 i = self.symbols.index(gen)

427 except ValueError:

428 raise ValueError("invalid generator: %s" % gen)

429 else:

430 raise ValueError("expected a polynomial generator, an integer, a string or None, got %s" % gen)

431

432 return i

433

434 def drop(self, *gens):

435 """Remove specified generators from this ring. """

436 indices = set(map(self.index, gens))

437 symbols = [ s for i, s in enumerate(self.symbols) if i not in indices ]

438

439 if not symbols:

440 return self.domain

441 else:

442 return self.clone(symbols=symbols)

443

444 def __getitem__(self, key):

445 symbols = self.symbols[key]

446

447 if not symbols:

448 return self.domain

449 else:

450 return self.clone(symbols=symbols)

451

452 def to_ground(self):

453 # TODO: should AlgebraicField be a Composite domain?

454 if self.domain.is_Composite or hasattr(self.domain, 'domain'):

455 return self.clone(domain=self.domain.domain)

456 else:

457 raise ValueError("%s is not a composite domain" % self.domain)

458

459 def to_domain(self):

460 return PolynomialRing(self)

461

462 def to_field(self):

463 from sympy.polys.fields import FracField

464 return FracField(self.symbols, self.domain, self.order)

465

466 @property

467 def is_univariate(self):

468 return len(self.gens) == 1

469

470 @property

471 def is_multivariate(self):

472 return len(self.gens) > 1

473

474 def add(self, *objs):

475 """

476 Add a sequence of polynomials or containers of polynomials.

477

478 Examples

479 ========

480

481 >>> from sympy.polys.rings import ring

482 >>> from sympy.polys.domains import ZZ

483

484 >>> R, x = ring("x", ZZ)

485 >>> R.add([ x**2 + 2*i + 3 for i in range(4) ])

486 4*x**2 + 24

487 >>> _.factor_list()

488 (4, [(x**2 + 6, 1)])

489

490 """

491 p = self.zero

492

493 for obj in objs:

494 if is_sequence(obj, include=GeneratorType):

495 p += self.add(*obj)

496 else:

497 p += obj

498

499 return p

500

501 def mul(self, *objs):

502 """

503 Multiply a sequence of polynomials or containers of polynomials.

504

505 Examples

506 ========

507

508 >>> from sympy.polys.rings import ring

509 >>> from sympy.polys.domains import ZZ

510

511 >>> R, x = ring("x", ZZ)

512 >>> R.mul([ x**2 + 2*i + 3 for i in range(4) ])

513 x**8 + 24*x**6 + 206*x**4 + 744*x**2 + 945

514 >>> _.factor_list()

515 (1, [(x**2 + 3, 1), (x**2 + 5, 1), (x**2 + 7, 1), (x**2 + 9, 1)])

516

517 """

518 p = self.one

519

520 for obj in objs:

521 if is_sequence(obj, include=GeneratorType):

522 p *= self.mul(*obj)

523 else:

524 p *= obj

525

526 return p

527

528 def drop_to_ground(self, *gens):

529 r"""

530 Remove specified generators from the ring and inject them into

531 its domain.

532 """

533 indices = set(map(self.index, gens))

534 symbols = [s for i, s in enumerate(self.symbols) if i not in indices]

535 gens = [gen for i, gen in enumerate(self.gens) if i not in indices]

536

537 if not symbols:

538 return self

539 else:

540 return self.clone(symbols=symbols, domain=self.drop(*gens))

541

542 def compose(self, other):

543 """Add the generators of ``other`` to ``self``"""

544 if self != other:

545 syms = set(self.symbols).union(set(other.symbols))

546 return self.clone(symbols=list(syms))

547 else:

548 return self

549

550 def add_gens(self, symbols):

551 """Add the elements of ``symbols`` as generators to ``self``"""

552 syms = set(self.symbols).union(set(symbols))

553 return self.clone(symbols=list(syms))

554

555 def symmetric_poly(self, n):

556 """

557 Return the elementary symmetric polynomial of degree *n* over

558 this ring's generators.

559 """

560 if n < 0 or n > self.ngens:

561 raise ValueError("Cannot generate symmetric polynomial of order %s for %s" % (n, self.gens))

562 elif not n:

563 return self.one

564 else:

565 poly = self.zero

566 for s in subsets(range(self.ngens), int(n)):

567 monom = tuple(int(i in s) for i in range(self.ngens))

568 poly += self.term_new(monom, self.domain.one)

569 return poly

570

571

572class PolyElement(DomainElement, DefaultPrinting, CantSympify, dict):

573 """Element of multivariate distributed polynomial ring. """

574

575 def new(self, init):

576 return self.__class__(init)

577

578 def parent(self):

579 return self.ring.to_domain()

580

581 def __getnewargs__(self):

582 return (self.ring, list(self.iterterms()))

583

584 _hash = None

585

586 def __hash__(self):

587 # XXX: This computes a hash of a dictionary, but currently we don't

588 # protect dictionary from being changed so any use site modifications

589 # will make hashing go wrong. Use this feature with caution until we

590 # figure out how to make a safe API without compromising speed of this

591 # low-level class.

592 _hash = self._hash

593 if _hash is None:

594 self._hash = _hash = hash((self.ring, frozenset(self.items())))

595 return _hash

596

597 def copy(self):

598 """Return a copy of polynomial self.

599

600 Polynomials are mutable; if one is interested in preserving

601 a polynomial, and one plans to use inplace operations, one

602 can copy the polynomial. This method makes a shallow copy.

603

604 Examples

605 ========

606

607 >>> from sympy.polys.domains import ZZ

608 >>> from sympy.polys.rings import ring

609

610 >>> R, x, y = ring('x, y', ZZ)

611 >>> p = (x + y)**2

612 >>> p1 = p.copy()

613 >>> p2 = p

614 >>> p[R.zero_monom] = 3

615 >>> p

616 x**2 + 2*x*y + y**2 + 3

617 >>> p1

618 x**2 + 2*x*y + y**2

619 >>> p2

620 x**2 + 2*x*y + y**2 + 3

621

622 """

623 return self.new(self)

624

625 def set_ring(self, new_ring):

626 if self.ring == new_ring:

627 return self

628 elif self.ring.symbols != new_ring.symbols:

629 terms = list(zip(*_dict_reorder(self, self.ring.symbols, new_ring.symbols)))

630 return new_ring.from_terms(terms, self.ring.domain)

631 else:

632 return new_ring.from_dict(self, self.ring.domain)

633

634 def as_expr(self, *symbols):

635 if not symbols:

636 symbols = self.ring.symbols

637 elif len(symbols) != self.ring.ngens:

638 raise ValueError(

639 "Wrong number of symbols, expected %s got %s" %

640 (self.ring.ngens, len(symbols))

641 )

642

643 return expr_from_dict(self.as_expr_dict(), *symbols)

644

645 def as_expr_dict(self):

646 to_sympy = self.ring.domain.to_sympy

647 return {monom: to_sympy(coeff) for monom, coeff in self.iterterms()}

648

649 def clear_denoms(self):

650 domain = self.ring.domain

651

652 if not domain.is_Field or not domain.has_assoc_Ring:

653 return domain.one, self

654

655 ground_ring = domain.get_ring()

656 common = ground_ring.one

657 lcm = ground_ring.lcm

658 denom = domain.denom

659

660 for coeff in self.values():

661 common = lcm(common, denom(coeff))

662

663 poly = self.new([ (k, v*common) for k, v in self.items() ])

664 return common, poly

665

666 def strip_zero(self):

667 """Eliminate monomials with zero coefficient. """

668 for k, v in list(self.items()):

669 if not v:

670 del self[k]

671

672 def __eq__(p1, p2):

673 """Equality test for polynomials.

674

675 Examples

676 ========

677

678 >>> from sympy.polys.domains import ZZ

679 >>> from sympy.polys.rings import ring

680

681 >>> _, x, y = ring('x, y', ZZ)

682 >>> p1 = (x + y)**2 + (x - y)**2

683 >>> p1 == 4*x*y

684 False

685 >>> p1 == 2*(x**2 + y**2)

686 True

687

688 """

689 if not p2:

690 return not p1

691 elif isinstance(p2, PolyElement) and p2.ring == p1.ring:

692 return dict.__eq__(p1, p2)

693 elif len(p1) > 1:

694 return False

695 else:

696 return p1.get(p1.ring.zero_monom) == p2

697

698 def __ne__(p1, p2):

699 return not p1 == p2

700

701 def almosteq(p1, p2, tolerance=None):

702 """Approximate equality test for polynomials. """

703 ring = p1.ring

704

705 if isinstance(p2, ring.dtype):

706 if set(p1.keys()) != set(p2.keys()):

707 return False

708

709 almosteq = ring.domain.almosteq

710

711 for k in p1.keys():

712 if not almosteq(p1[k], p2[k], tolerance):

713 return False

714 return True

715 elif len(p1) > 1:

716 return False

717 else:

718 try:

719 p2 = ring.domain.convert(p2)

720 except CoercionFailed:

721 return False

722 else:

723 return ring.domain.almosteq(p1.const(), p2, tolerance)

724

725 def sort_key(self):

726 return (len(self), self.terms())

727

728 def _cmp(p1, p2, op):

729 if isinstance(p2, p1.ring.dtype):

730 return op(p1.sort_key(), p2.sort_key())

731 else:

732 return NotImplemented

733

734 def __lt__(p1, p2):

735 return p1._cmp(p2, lt)

736 def __le__(p1, p2):

737 return p1._cmp(p2, le)

738 def __gt__(p1, p2):

739 return p1._cmp(p2, gt)

740 def __ge__(p1, p2):

741 return p1._cmp(p2, ge)

742

743 def _drop(self, gen):

744 ring = self.ring

745 i = ring.index(gen)

746

747 if ring.ngens == 1:

748 return i, ring.domain

749 else:

750 symbols = list(ring.symbols)

751 del symbols[i]

752 return i, ring.clone(symbols=symbols)

753

754 def drop(self, gen):

755 i, ring = self._drop(gen)

756

757 if self.ring.ngens == 1:

758 if self.is_ground:

759 return self.coeff(1)

760 else:

761 raise ValueError("Cannot drop %s" % gen)

762 else:

763 poly = ring.zero

764

765 for k, v in self.items():

766 if k[i] == 0:

767 K = list(k)

768 del K[i]

769 poly[tuple(K)] = v

770 else:

771 raise ValueError("Cannot drop %s" % gen)

772

773 return poly

774

775 def _drop_to_ground(self, gen):

776 ring = self.ring

777 i = ring.index(gen)

778

779 symbols = list(ring.symbols)

780 del symbols[i]

781 return i, ring.clone(symbols=symbols, domain=ring[i])

782

783 def drop_to_ground(self, gen):

784 if self.ring.ngens == 1:

785 raise ValueError("Cannot drop only generator to ground")

786

787 i, ring = self._drop_to_ground(gen)

788 poly = ring.zero

789 gen = ring.domain.gens[0]

790

791 for monom, coeff in self.iterterms():

792 mon = monom[:i] + monom[i+1:]

793 if mon not in poly:

794 poly[mon] = (gen**monom[i]).mul_ground(coeff)

795 else:

796 poly[mon] += (gen**monom[i]).mul_ground(coeff)

797

798 return poly

799

800 def to_dense(self):

801 return dmp_from_dict(self, self.ring.ngens-1, self.ring.domain)

802

803 def to_dict(self):

804 return dict(self)

805

806 def str(self, printer, precedence, exp_pattern, mul_symbol):

807 if not self:

808 return printer._print(self.ring.domain.zero)

809 prec_mul = precedence["Mul"]

810 prec_atom = precedence["Atom"]

811 ring = self.ring

812 symbols = ring.symbols

813 ngens = ring.ngens

814 zm = ring.zero_monom

815 sexpvs = []

816 for expv, coeff in self.terms():

817 negative = ring.domain.is_negative(coeff)

818 sign = " - " if negative else " + "

819 sexpvs.append(sign)

820 if expv == zm:

821 scoeff = printer._print(coeff)

822 if negative and scoeff.startswith("-"):

823 scoeff = scoeff[1:]

824 else:

825 if negative:

826 coeff = -coeff

827 if coeff != self.ring.domain.one:

828 scoeff = printer.parenthesize(coeff, prec_mul, strict=True)

829 else:

830 scoeff = ''

831 sexpv = []

832 for i in range(ngens):

833 exp = expv[i]

834 if not exp:

835 continue

836 symbol = printer.parenthesize(symbols[i], prec_atom, strict=True)

837 if exp != 1:

838 if exp != int(exp) or exp < 0:

839 sexp = printer.parenthesize(exp, prec_atom, strict=False)

840 else:

841 sexp = exp

842 sexpv.append(exp_pattern % (symbol, sexp))

843 else:

844 sexpv.append('%s' % symbol)

845 if scoeff:

846 sexpv = [scoeff] + sexpv

847 sexpvs.append(mul_symbol.join(sexpv))

848 if sexpvs[0] in [" + ", " - "]:

849 head = sexpvs.pop(0)

850 if head == " - ":

851 sexpvs.insert(0, "-")

852 return "".join(sexpvs)

853

854 @property

855 def is_generator(self):

856 return self in self.ring._gens_set

857

858 @property

859 def is_ground(self):

860 return not self or (len(self) == 1 and self.ring.zero_monom in self)

861

862 @property

863 def is_monomial(self):

864 return not self or (len(self) == 1 and self.LC == 1)

865

866 @property

867 def is_term(self):

868 return len(self) <= 1

869

870 @property

871 def is_negative(self):

872 return self.ring.domain.is_negative(self.LC)

873

874 @property

875 def is_positive(self):

876 return self.ring.domain.is_positive(self.LC)

877

878 @property

879 def is_nonnegative(self):

880 return self.ring.domain.is_nonnegative(self.LC)

881

882 @property

883 def is_nonpositive(self):

884 return self.ring.domain.is_nonpositive(self.LC)

885

886 @property

887 def is_zero(f):

888 return not f

889

890 @property

891 def is_one(f):

892 return f == f.ring.one

893

894 @property

895 def is_monic(f):

896 return f.ring.domain.is_one(f.LC)

897

898 @property

899 def is_primitive(f):

900 return f.ring.domain.is_one(f.content())

901

902 @property

903 def is_linear(f):

904 return all(sum(monom) <= 1 for monom in f.itermonoms())

905

906 @property

907 def is_quadratic(f):

908 return all(sum(monom) <= 2 for monom in f.itermonoms())

909

910 @property

911 def is_squarefree(f):

912 if not f.ring.ngens:

913 return True

914 return f.ring.dmp_sqf_p(f)

915

916 @property

917 def is_irreducible(f):

918 if not f.ring.ngens:

919 return True

920 return f.ring.dmp_irreducible_p(f)

921

922 @property

923 def is_cyclotomic(f):

924 if f.ring.is_univariate:

925 return f.ring.dup_cyclotomic_p(f)

926 else:

927 raise MultivariatePolynomialError("cyclotomic polynomial")

928

929 def __neg__(self):

930 return self.new([ (monom, -coeff) for monom, coeff in self.iterterms() ])

931

932 def __pos__(self):

933 return self

934

935 def __add__(p1, p2):

936 """Add two polynomials.

937

938 Examples

939 ========

940

941 >>> from sympy.polys.domains import ZZ

942 >>> from sympy.polys.rings import ring

943

944 >>> _, x, y = ring('x, y', ZZ)

945 >>> (x + y)**2 + (x - y)**2

946 2*x**2 + 2*y**2

947

948 """

949 if not p2:

950 return p1.copy()

951 ring = p1.ring

952 if isinstance(p2, ring.dtype):

953 p = p1.copy()

954 get = p.get

955 zero = ring.domain.zero

956 for k, v in p2.items():

957 v = get(k, zero) + v

958 if v:

959 p[k] = v

960 else:

961 del p[k]

962 return p

963 elif isinstance(p2, PolyElement):

964 if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:

965 pass

966 elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:

967 return p2.__radd__(p1)

968 else:

969 return NotImplemented

970

971 try:

972 cp2 = ring.domain_new(p2)

973 except CoercionFailed:

974 return NotImplemented

975 else:

976 p = p1.copy()

977 if not cp2:

978 return p

979 zm = ring.zero_monom

980 if zm not in p1.keys():

981 p[zm] = cp2

982 else:

983 if p2 == -p[zm]:

984 del p[zm]

985 else:

986 p[zm] += cp2

987 return p

988

989 def __radd__(p1, n):

990 p = p1.copy()

991 if not n:

992 return p

993 ring = p1.ring

994 try:

995 n = ring.domain_new(n)

996 except CoercionFailed:

997 return NotImplemented

998 else:

999 zm = ring.zero_monom

1000 if zm not in p1.keys():

1001 p[zm] = n

1002 else:

1003 if n == -p[zm]:

1004 del p[zm]

1005 else:

1006 p[zm] += n

1007 return p

1008

1009 def __sub__(p1, p2):

1010 """Subtract polynomial p2 from p1.

1011

1012 Examples

1013 ========

1014

1015 >>> from sympy.polys.domains import ZZ

1016 >>> from sympy.polys.rings import ring

1017

1018 >>> _, x, y = ring('x, y', ZZ)

1019 >>> p1 = x + y**2

1020 >>> p2 = x*y + y**2

1021 >>> p1 - p2

1022 -x*y + x

1023

1024 """

1025 if not p2:

1026 return p1.copy()

1027 ring = p1.ring

1028 if isinstance(p2, ring.dtype):

1029 p = p1.copy()

1030 get = p.get

1031 zero = ring.domain.zero

1032 for k, v in p2.items():

1033 v = get(k, zero) - v

1034 if v:

1035 p[k] = v

1036 else:

1037 del p[k]

1038 return p

1039 elif isinstance(p2, PolyElement):

1040 if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:

1041 pass

1042 elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:

1043 return p2.__rsub__(p1)

1044 else:

1045 return NotImplemented

1046

1047 try:

1048 p2 = ring.domain_new(p2)

1049 except CoercionFailed:

1050 return NotImplemented

1051 else:

1052 p = p1.copy()

1053 zm = ring.zero_monom

1054 if zm not in p1.keys():

1055 p[zm] = -p2

1056 else:

1057 if p2 == p[zm]:

1058 del p[zm]

1059 else:

1060 p[zm] -= p2

1061 return p

1062

1063 def __rsub__(p1, n):

1064 """n - p1 with n convertible to the coefficient domain.

1065

1066 Examples

1067 ========

1068

1069 >>> from sympy.polys.domains import ZZ

1070 >>> from sympy.polys.rings import ring

1071

1072 >>> _, x, y = ring('x, y', ZZ)

1073 >>> p = x + y

1074 >>> 4 - p

1075 -x - y + 4

1076

1077 """

1078 ring = p1.ring

1079 try:

1080 n = ring.domain_new(n)

1081 except CoercionFailed:

1082 return NotImplemented

1083 else:

1084 p = ring.zero

1085 for expv in p1:

1086 p[expv] = -p1[expv]

1087 p += n

1088 return p

1089

1090 def __mul__(p1, p2):

1091 """Multiply two polynomials.

1092

1093 Examples

1094 ========

1095

1096 >>> from sympy.polys.domains import QQ

1097 >>> from sympy.polys.rings import ring

1098

1099 >>> _, x, y = ring('x, y', QQ)

1100 >>> p1 = x + y

1101 >>> p2 = x - y

1102 >>> p1*p2

1103 x**2 - y**2

1104

1105 """

1106 ring = p1.ring

1107 p = ring.zero

1108 if not p1 or not p2:

1109 return p

1110 elif isinstance(p2, ring.dtype):

1111 get = p.get

1112 zero = ring.domain.zero

1113 monomial_mul = ring.monomial_mul

1114 p2it = list(p2.items())

1115 for exp1, v1 in p1.items():

1116 for exp2, v2 in p2it:

1117 exp = monomial_mul(exp1, exp2)

1118 p[exp] = get(exp, zero) + v1*v2

1119 p.strip_zero()

1120 return p

1121 elif isinstance(p2, PolyElement):

1122 if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:

1123 pass

1124 elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:

1125 return p2.__rmul__(p1)

1126 else:

1127 return NotImplemented

1128

1129 try:

1130 p2 = ring.domain_new(p2)

1131 except CoercionFailed:

1132 return NotImplemented

1133 else:

1134 for exp1, v1 in p1.items():

1135 v = v1*p2

1136 if v:

1137 p[exp1] = v

1138 return p

1139

1140 def __rmul__(p1, p2):

1141 """p2 * p1 with p2 in the coefficient domain of p1.

1142

1143 Examples

1144 ========

1145

1146 >>> from sympy.polys.domains import ZZ

1147 >>> from sympy.polys.rings import ring

1148

1149 >>> _, x, y = ring('x, y', ZZ)

1150 >>> p = x + y

1151 >>> 4 * p

1152 4*x + 4*y

1153

1154 """

1155 p = p1.ring.zero

1156 if not p2:

1157 return p

1158 try:

1159 p2 = p.ring.domain_new(p2)

1160 except CoercionFailed:

1161 return NotImplemented

1162 else:

1163 for exp1, v1 in p1.items():

1164 v = p2*v1

1165 if v:

1166 p[exp1] = v

1167 return p

1168

1169 def __pow__(self, n):

1170 """raise polynomial to power `n`

1171

1172 Examples

1173 ========

1174

1175 >>> from sympy.polys.domains import ZZ

1176 >>> from sympy.polys.rings import ring

1177

1178 >>> _, x, y = ring('x, y', ZZ)

1179 >>> p = x + y**2

1180 >>> p**3

1181 x**3 + 3*x**2*y**2 + 3*x*y**4 + y**6

1182

1183 """

1184 ring = self.ring

1185

1186 if not n:

1187 if self:

1188 return ring.one

1189 else:

1190 raise ValueError("0**0")

1191 elif len(self) == 1:

1192 monom, coeff = list(self.items())[0]

1193 p = ring.zero

1194 if coeff == ring.domain.one:

1195 p[ring.monomial_pow(monom, n)] = coeff

1196 else:

1197 p[ring.monomial_pow(monom, n)] = coeff**n

1198 return p

1199

1200 # For ring series, we need negative and rational exponent support only

1201 # with monomials.

1202 n = int(n)

1203 if n < 0:

1204 raise ValueError("Negative exponent")

1205

1206 elif n == 1:

1207 return self.copy()

1208 elif n == 2:

1209 return self.square()

1210 elif n == 3:

1211 return self*self.square()

1212 elif len(self) <= 5: # TODO: use an actual density measure

1213 return self._pow_multinomial(n)

1214 else:

1215 return self._pow_generic(n)

1216

1217 def _pow_generic(self, n):

1218 p = self.ring.one

1219 c = self

1220

1221 while True:

1222 if n & 1:

1223 p = p*c

1224 n -= 1

1225 if not n:

1226 break

1227

1228 c = c.square()

1229 n = n // 2

1230

1231 return p

1232

1233 def _pow_multinomial(self, n):

1234 multinomials = multinomial_coefficients(len(self), n).items()

1235 monomial_mulpow = self.ring.monomial_mulpow

1236 zero_monom = self.ring.zero_monom

1237 terms = self.items()

1238 zero = self.ring.domain.zero

1239 poly = self.ring.zero

1240

1241 for multinomial, multinomial_coeff in multinomials:

1242 product_monom = zero_monom

1243 product_coeff = multinomial_coeff

1244

1245 for exp, (monom, coeff) in zip(multinomial, terms):

1246 if exp:

1247 product_monom = monomial_mulpow(product_monom, monom, exp)

1248 product_coeff *= coeff**exp

1249

1250 monom = tuple(product_monom)

1251 coeff = product_coeff

1252

1253 coeff = poly.get(monom, zero) + coeff

1254

1255 if coeff:

1256 poly[monom] = coeff

1257 elif monom in poly:

1258 del poly[monom]

1259

1260 return poly

1261

1262 def square(self):

1263 """square of a polynomial

1264

1265 Examples

1266 ========

1267

1268 >>> from sympy.polys.rings import ring

1269 >>> from sympy.polys.domains import ZZ

1270

1271 >>> _, x, y = ring('x, y', ZZ)

1272 >>> p = x + y**2

1273 >>> p.square()

1274 x**2 + 2*x*y**2 + y**4

1275

1276 """

1277 ring = self.ring

1278 p = ring.zero

1279 get = p.get

1280 keys = list(self.keys())

1281 zero = ring.domain.zero

1282 monomial_mul = ring.monomial_mul

1283 for i in range(len(keys)):

1284 k1 = keys[i]

1285 pk = self[k1]

1286 for j in range(i):

1287 k2 = keys[j]

1288 exp = monomial_mul(k1, k2)

1289 p[exp] = get(exp, zero) + pk*self[k2]

1290 p = p.imul_num(2)

1291 get = p.get

1292 for k, v in self.items():

1293 k2 = monomial_mul(k, k)

1294 p[k2] = get(k2, zero) + v**2

1295 p.strip_zero()

1296 return p

1297

1298 def __divmod__(p1, p2):

1299 ring = p1.ring

1300

1301 if not p2:

1302 raise ZeroDivisionError("polynomial division")

1303 elif isinstance(p2, ring.dtype):

1304 return p1.div(p2)

1305 elif isinstance(p2, PolyElement):

1306 if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:

1307 pass

1308 elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:

1309 return p2.__rdivmod__(p1)

1310 else:

1311 return NotImplemented

1312

1313 try:

1314 p2 = ring.domain_new(p2)

1315 except CoercionFailed:

1316 return NotImplemented

1317 else:

1318 return (p1.quo_ground(p2), p1.rem_ground(p2))

1319

1320 def __rdivmod__(p1, p2):

1321 return NotImplemented

1322

1323 def __mod__(p1, p2):

1324 ring = p1.ring

1325

1326 if not p2:

1327 raise ZeroDivisionError("polynomial division")

1328 elif isinstance(p2, ring.dtype):

1329 return p1.rem(p2)

1330 elif isinstance(p2, PolyElement):

1331 if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:

1332 pass

1333 elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:

1334 return p2.__rmod__(p1)

1335 else:

1336 return NotImplemented

1337

1338 try:

1339 p2 = ring.domain_new(p2)

1340 except CoercionFailed:

1341 return NotImplemented

1342 else:

1343 return p1.rem_ground(p2)

1344

1345 def __rmod__(p1, p2):

1346 return NotImplemented

1347

1348 def __truediv__(p1, p2):

1349 ring = p1.ring

1350

1351 if not p2:

1352 raise ZeroDivisionError("polynomial division")

1353 elif isinstance(p2, ring.dtype):

1354 if p2.is_monomial:

1355 return p1*(p2**(-1))

1356 else:

1357 return p1.quo(p2)

1358 elif isinstance(p2, PolyElement):

1359 if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:

1360 pass

1361 elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:

1362 return p2.__rtruediv__(p1)

1363 else:

1364 return NotImplemented

1365

1366 try:

1367 p2 = ring.domain_new(p2)

1368 except CoercionFailed:

1369 return NotImplemented

1370 else:

1371 return p1.quo_ground(p2)

1372

1373 def __rtruediv__(p1, p2):

1374 return NotImplemented

1375

1376 __floordiv__ = __truediv__

1377 __rfloordiv__ = __rtruediv__

1378

1379 # TODO: use // (__floordiv__) for exquo()?

1380

1381 def _term_div(self):

1382 zm = self.ring.zero_monom

1383 domain = self.ring.domain

1384 domain_quo = domain.quo

1385 monomial_div = self.ring.monomial_div

1386

1387 if domain.is_Field:

1388 def term_div(a_lm_a_lc, b_lm_b_lc):

1389 a_lm, a_lc = a_lm_a_lc

1390 b_lm, b_lc = b_lm_b_lc

1391 if b_lm == zm: # apparently this is a very common case

1392 monom = a_lm

1393 else:

1394 monom = monomial_div(a_lm, b_lm)

1395 if monom is not None:

1396 return monom, domain_quo(a_lc, b_lc)

1397 else:

1398 return None

1399 else:

1400 def term_div(a_lm_a_lc, b_lm_b_lc):

1401 a_lm, a_lc = a_lm_a_lc

1402 b_lm, b_lc = b_lm_b_lc

1403 if b_lm == zm: # apparently this is a very common case

1404 monom = a_lm

1405 else:

1406 monom = monomial_div(a_lm, b_lm)

1407 if not (monom is None or a_lc % b_lc):

1408 return monom, domain_quo(a_lc, b_lc)

1409 else:

1410 return None

1411

1412 return term_div

1413

1414 def div(self, fv):

1415 """Division algorithm, see [CLO] p64.

1416

1417 fv array of polynomials

1418 return qv, r such that

1419 self = sum(fv[i]*qv[i]) + r

1420

1421 All polynomials are required not to be Laurent polynomials.

1422

1423 Examples

1424 ========

1425

1426 >>> from sympy.polys.rings import ring

1427 >>> from sympy.polys.domains import ZZ

1428

1429 >>> _, x, y = ring('x, y', ZZ)

1430 >>> f = x**3

1431 >>> f0 = x - y**2

1432 >>> f1 = x - y

1433 >>> qv, r = f.div((f0, f1))

1434 >>> qv[0]

1435 x**2 + x*y**2 + y**4

1436 >>> qv[1]

1437 0

1438 >>> r

1439 y**6

1440

1441 """

1442 ring = self.ring

1443 ret_single = False

1444 if isinstance(fv, PolyElement):

1445 ret_single = True

1446 fv = [fv]

1447 if not all(fv):

1448 raise ZeroDivisionError("polynomial division")

1449 if not self:

1450 if ret_single:

1451 return ring.zero, ring.zero

1452 else:

1453 return [], ring.zero

1454 for f in fv:

1455 if f.ring != ring:

1456 raise ValueError('self and f must have the same ring')

1457 s = len(fv)

1458 qv = [ring.zero for i in range(s)]

1459 p = self.copy()

1460 r = ring.zero

1461 term_div = self._term_div()

1462 expvs = [fx.leading_expv() for fx in fv]

1463 while p:

1464 i = 0

1465 divoccurred = 0

1466 while i < s and divoccurred == 0:

1467 expv = p.leading_expv()

1468 term = term_div((expv, p[expv]), (expvs[i], fv[i][expvs[i]]))

1469 if term is not None:

1470 expv1, c = term

1471 qv[i] = qv[i]._iadd_monom((expv1, c))

1472 p = p._iadd_poly_monom(fv[i], (expv1, -c))

1473 divoccurred = 1

1474 else:

1475 i += 1

1476 if not divoccurred:

1477 expv = p.leading_expv()

1478 r = r._iadd_monom((expv, p[expv]))

1479 del p[expv]

1480 if expv == ring.zero_monom:

1481 r += p

1482 if ret_single:

1483 if not qv:

1484 return ring.zero, r

1485 else:

1486 return qv[0], r

1487 else:

1488 return qv, r

1489

1490 def rem(self, G):

1491 f = self

1492 if isinstance(G, PolyElement):

1493 G = [G]

1494 if not all(G):

1495 raise ZeroDivisionError("polynomial division")

1496 ring = f.ring

1497 domain = ring.domain

1498 zero = domain.zero

1499 monomial_mul = ring.monomial_mul

1500 r = ring.zero

1501 term_div = f._term_div()

1502 ltf = f.LT

1503 f = f.copy()

1504 get = f.get

1505 while f:

1506 for g in G:

1507 tq = term_div(ltf, g.LT)

1508 if tq is not None:

1509 m, c = tq

1510 for mg, cg in g.iterterms():

1511 m1 = monomial_mul(mg, m)

1512 c1 = get(m1, zero) - c*cg

1513 if not c1:

1514 del f[m1]

1515 else:

1516 f[m1] = c1

1517 ltm = f.leading_expv()

1518 if ltm is not None:

1519 ltf = ltm, f[ltm]

1520

1521 break

1522 else:

1523 ltm, ltc = ltf

1524 if ltm in r:

1525 r[ltm] += ltc

1526 else:

1527 r[ltm] = ltc

1528 del f[ltm]

1529 ltm = f.leading_expv()

1530 if ltm is not None:

1531 ltf = ltm, f[ltm]

1532

1533 return r

1534

1535 def quo(f, G):

1536 return f.div(G)[0]

1537

1538 def exquo(f, G):

1539 q, r = f.div(G)

1540

1541 if not r:

1542 return q

1543 else:

1544 raise ExactQuotientFailed(f, G)

1545

1546 def _iadd_monom(self, mc):

1547 """add to self the monomial coeff*x0**i0*x1**i1*...

1548 unless self is a generator -- then just return the sum of the two.

1549

1550 mc is a tuple, (monom, coeff), where monomial is (i0, i1, ...)

1551

1552 Examples

1553 ========

1554

1555 >>> from sympy.polys.rings import ring

1556 >>> from sympy.polys.domains import ZZ

1557

1558 >>> _, x, y = ring('x, y', ZZ)

1559 >>> p = x**4 + 2*y

1560 >>> m = (1, 2)

1561 >>> p1 = p._iadd_monom((m, 5))

1562 >>> p1

1563 x**4 + 5*x*y**2 + 2*y

1564 >>> p1 is p

1565 True

1566 >>> p = x

1567 >>> p1 = p._iadd_monom((m, 5))

1568 >>> p1

1569 5*x*y**2 + x

1570 >>> p1 is p

1571 False

1572

1573 """

1574 if self in self.ring._gens_set:

1575 cpself = self.copy()

1576 else:

1577 cpself = self

1578 expv, coeff = mc

1579 c = cpself.get(expv)

1580 if c is None:

1581 cpself[expv] = coeff

1582 else:

1583 c += coeff

1584 if c:

1585 cpself[expv] = c

1586 else:

1587 del cpself[expv]

1588 return cpself

1589

1590 def _iadd_poly_monom(self, p2, mc):

1591 """add to self the product of (p)*(coeff*x0**i0*x1**i1*...)

1592 unless self is a generator -- then just return the sum of the two.

1593

1594 mc is a tuple, (monom, coeff), where monomial is (i0, i1, ...)

1595

1596 Examples

1597 ========

1598

1599 >>> from sympy.polys.rings import ring

1600 >>> from sympy.polys.domains import ZZ

1601

1602 >>> _, x, y, z = ring('x, y, z', ZZ)

1603 >>> p1 = x**4 + 2*y

1604 >>> p2 = y + z

1605 >>> m = (1, 2, 3)

1606 >>> p1 = p1._iadd_poly_monom(p2, (m, 3))

1607 >>> p1

1608 x**4 + 3*x*y**3*z**3 + 3*x*y**2*z**4 + 2*y

1609

1610 """

1611 p1 = self

1612 if p1 in p1.ring._gens_set:

1613 p1 = p1.copy()

1614 (m, c) = mc

1615 get = p1.get

1616 zero = p1.ring.domain.zero

1617 monomial_mul = p1.ring.monomial_mul

1618 for k, v in p2.items():

1619 ka = monomial_mul(k, m)

1620 coeff = get(ka, zero) + v*c

1621 if coeff:

1622 p1[ka] = coeff

1623 else:

1624 del p1[ka]

1625 return p1

1626

1627 def degree(f, x=None):

1628 """

1629 The leading degree in ``x`` or the main variable.

1630

1631 Note that the degree of 0 is negative infinity (the SymPy object -oo).

1632

1633 """

1634 i = f.ring.index(x)

1635

1636 if not f:

1637 return -oo

1638 elif i < 0:

1639 return 0

1640 else:

1641 return max([ monom[i] for monom in f.itermonoms() ])

1642

1643 def degrees(f):

1644 """

1645 A tuple containing leading degrees in all variables.

1646

1647 Note that the degree of 0 is negative infinity (the SymPy object -oo)

1648

1649 """

1650 if not f:

1651 return (-oo,)*f.ring.ngens

1652 else:

1653 return tuple(map(max, list(zip(*f.itermonoms()))))

1654

1655 def tail_degree(f, x=None):

1656 """

1657 The tail degree in ``x`` or the main variable.

1658

1659 Note that the degree of 0 is negative infinity (the SymPy object -oo)

1660

1661 """

1662 i = f.ring.index(x)

1663

1664 if not f:

1665 return -oo

1666 elif i < 0:

1667 return 0

1668 else:

1669 return min([ monom[i] for monom in f.itermonoms() ])

1670

1671 def tail_degrees(f):

1672 """

1673 A tuple containing tail degrees in all variables.

1674

1675 Note that the degree of 0 is negative infinity (the SymPy object -oo)

1676

1677 """

1678 if not f:

1679 return (-oo,)*f.ring.ngens

1680 else:

1681 return tuple(map(min, list(zip(*f.itermonoms()))))

1682

1683 def leading_expv(self):

1684 """Leading monomial tuple according to the monomial ordering.

1685

1686 Examples

1687 ========

1688

1689 >>> from sympy.polys.rings import ring

1690 >>> from sympy.polys.domains import ZZ

1691

1692 >>> _, x, y, z = ring('x, y, z', ZZ)

1693 >>> p = x**4 + x**3*y + x**2*z**2 + z**7

1694 >>> p.leading_expv()

1695 (4, 0, 0)

1696

1697 """

1698 if self:

1699 return self.ring.leading_expv(self)

1700 else:

1701 return None

1702

1703 def _get_coeff(self, expv):

1704 return self.get(expv, self.ring.domain.zero)

1705

1706 def coeff(self, element):

1707 """

1708 Returns the coefficient that stands next to the given monomial.

1709

1710 Parameters

1711 ==========

1712

1713 element : PolyElement (with ``is_monomial = True``) or 1

1714

1715 Examples

1716 ========

1717

1718 >>> from sympy.polys.rings import ring

1719 >>> from sympy.polys.domains import ZZ

1720

1721 >>> _, x, y, z = ring("x,y,z", ZZ)

1722 >>> f = 3*x**2*y - x*y*z + 7*z**3 + 23

1723

1724 >>> f.coeff(x**2*y)

1725 3

1726 >>> f.coeff(x*y)

1727 0

1728 >>> f.coeff(1)

1729 23

1730

1731 """

1732 if element == 1:

1733 return self._get_coeff(self.ring.zero_monom)

1734 elif isinstance(element, self.ring.dtype):

1735 terms = list(element.iterterms())

1736 if len(terms) == 1:

1737 monom, coeff = terms[0]

1738 if coeff == self.ring.domain.one:

1739 return self._get_coeff(monom)

1740

1741 raise ValueError("expected a monomial, got %s" % element)

1742

1743 def const(self):

1744 """Returns the constant coefficient. """

1745 return self._get_coeff(self.ring.zero_monom)

1746

1747 @property

1748 def LC(self):

1749 return self._get_coeff(self.leading_expv())

1750

1751 @property

1752 def LM(self):

1753 expv = self.leading_expv()

1754 if expv is None:

1755 return self.ring.zero_monom

1756 else:

1757 return expv

1758

1759 def leading_monom(self):

1760 """

1761 Leading monomial as a polynomial element.

1762

1763 Examples

1764 ========

1765

1766 >>> from sympy.polys.rings import ring

1767 >>> from sympy.polys.domains import ZZ

1768

1769 >>> _, x, y = ring('x, y', ZZ)

1770 >>> (3*x*y + y**2).leading_monom()

1771 x*y

1772

1773 """

1774 p = self.ring.zero

1775 expv = self.leading_expv()

1776 if expv:

1777 p[expv] = self.ring.domain.one

1778 return p

1779

1780 @property

1781 def LT(self):

1782 expv = self.leading_expv()

1783 if expv is None:

1784 return (self.ring.zero_monom, self.ring.domain.zero)

1785 else:

1786 return (expv, self._get_coeff(expv))

1787

1788 def leading_term(self):

1789 """Leading term as a polynomial element.

1790

1791 Examples

1792 ========

1793

1794 >>> from sympy.polys.rings import ring

1795 >>> from sympy.polys.domains import ZZ

1796

1797 >>> _, x, y = ring('x, y', ZZ)

1798 >>> (3*x*y + y**2).leading_term()

1799 3*x*y

1800

1801 """

1802 p = self.ring.zero

1803 expv = self.leading_expv()

1804 if expv is not None:

1805 p[expv] = self[expv]

1806 return p

1807

1808 def _sorted(self, seq, order):

1809 if order is None:

1810 order = self.ring.order

1811 else:

1812 order = OrderOpt.preprocess(order)

1813

1814 if order is lex:

1815 return sorted(seq, key=lambda monom: monom[0], reverse=True)

1816 else:

1817 return sorted(seq, key=lambda monom: order(monom[0]), reverse=True)

1818

1819 def coeffs(self, order=None):

1820 """Ordered list of polynomial coefficients.

1821

1822 Parameters

1823 ==========

1824

1825 order : :class:`~.MonomialOrder` or coercible, optional

1826

1827 Examples

1828 ========

1829

1830 >>> from sympy.polys.rings import ring

1831 >>> from sympy.polys.domains import ZZ

1832 >>> from sympy.polys.orderings import lex, grlex

1833

1834 >>> _, x, y = ring("x, y", ZZ, lex)

1835 >>> f = x*y**7 + 2*x**2*y**3

1836

1837 >>> f.coeffs()

1838 [2, 1]

1839 >>> f.coeffs(grlex)

1840 [1, 2]

1841

1842 """

1843 return [ coeff for _, coeff in self.terms(order) ]

1844

1845 def monoms(self, order=None):

1846 """Ordered list of polynomial monomials.

1847

1848 Parameters

1849 ==========

1850

1851 order : :class:`~.MonomialOrder` or coercible, optional

1852

1853 Examples

1854 ========

1855

1856 >>> from sympy.polys.rings import ring

1857 >>> from sympy.polys.domains import ZZ

1858 >>> from sympy.polys.orderings import lex, grlex

1859

1860 >>> _, x, y = ring("x, y", ZZ, lex)

1861 >>> f = x*y**7 + 2*x**2*y**3

1862

1863 >>> f.monoms()

1864 [(2, 3), (1, 7)]

1865 >>> f.monoms(grlex)

1866 [(1, 7), (2, 3)]

1867

1868 """

1869 return [ monom for monom, _ in self.terms(order) ]

1870

1871 def terms(self, order=None):

1872 """Ordered list of polynomial terms.

1873

1874 Parameters

1875 ==========

1876

1877 order : :class:`~.MonomialOrder` or coercible, optional

1878

1879 Examples

1880 ========

1881

1882 >>> from sympy.polys.rings import ring

1883 >>> from sympy.polys.domains import ZZ

1884 >>> from sympy.polys.orderings import lex, grlex

1885

1886 >>> _, x, y = ring("x, y", ZZ, lex)

1887 >>> f = x*y**7 + 2*x**2*y**3

1888

1889 >>> f.terms()

1890 [((2, 3), 2), ((1, 7), 1)]

1891 >>> f.terms(grlex)

1892 [((1, 7), 1), ((2, 3), 2)]

1893

1894 """

1895 return self._sorted(list(self.items()), order)

1896

1897 def itercoeffs(self):

1898 """Iterator over coefficients of a polynomial. """

1899 return iter(self.values())

1900

1901 def itermonoms(self):

1902 """Iterator over monomials of a polynomial. """

1903 return iter(self.keys())

1904

1905 def iterterms(self):

1906 """Iterator over terms of a polynomial. """

1907 return iter(self.items())

1908

1909 def listcoeffs(self):

1910 """Unordered list of polynomial coefficients. """

1911 return list(self.values())

1912

1913 def listmonoms(self):

1914 """Unordered list of polynomial monomials. """

1915 return list(self.keys())

1916

1917 def listterms(self):

1918 """Unordered list of polynomial terms. """

1919 return list(self.items())

1920

1921 def imul_num(p, c):

1922 """multiply inplace the polynomial p by an element in the

1923 coefficient ring, provided p is not one of the generators;

1924 else multiply not inplace

1925

1926 Examples

1927 ========

1928

1929 >>> from sympy.polys.rings import ring

1930 >>> from sympy.polys.domains import ZZ

1931

1932 >>> _, x, y = ring('x, y', ZZ)

1933 >>> p = x + y**2

1934 >>> p1 = p.imul_num(3)

1935 >>> p1

1936 3*x + 3*y**2

1937 >>> p1 is p

1938 True

1939 >>> p = x

1940 >>> p1 = p.imul_num(3)

1941 >>> p1

1942 3*x

1943 >>> p1 is p

1944 False

1945

1946 """

1947 if p in p.ring._gens_set:

1948 return p*c

1949 if not c:

1950 p.clear()

1951 return

1952 for exp in p:

1953 p[exp] *= c

1954 return p

1955

1956 def content(f):

1957 """Returns GCD of polynomial's coefficients. """

1958 domain = f.ring.domain

1959 cont = domain.zero

1960 gcd = domain.gcd

1961

1962 for coeff in f.itercoeffs():

1963 cont = gcd(cont, coeff)

1964

1965 return cont

1966

1967 def primitive(f):

1968 """Returns content and a primitive polynomial. """

1969 cont = f.content()

1970 return cont, f.quo_ground(cont)

1971

1972 def monic(f):

1973 """Divides all coefficients by the leading coefficient. """

1974 if not f:

1975 return f

1976 else:

1977 return f.quo_ground(f.LC)

1978

1979 def mul_ground(f, x):

1980 if not x:

1981 return f.ring.zero

1982

1983 terms = [ (monom, coeff*x) for monom, coeff in f.iterterms() ]

1984 return f.new(terms)

1985

1986 def mul_monom(f, monom):

1987 monomial_mul = f.ring.monomial_mul

1988 terms = [ (monomial_mul(f_monom, monom), f_coeff) for f_monom, f_coeff in f.items() ]

1989 return f.new(terms)

1990

1991 def mul_term(f, term):

1992 monom, coeff = term

1993

1994 if not f or not coeff:

1995 return f.ring.zero

1996 elif monom == f.ring.zero_monom:

1997 return f.mul_ground(coeff)

1998

1999 monomial_mul = f.ring.monomial_mul

2000 terms = [ (monomial_mul(f_monom, monom), f_coeff*coeff) for f_monom, f_coeff in f.items() ]

2001 return f.new(terms)

2002

2003 def quo_ground(f, x):

2004 domain = f.ring.domain

2005

2006 if not x:

2007 raise ZeroDivisionError('polynomial division')

2008 if not f or x == domain.one:

2009 return f

2010

2011 if domain.is_Field:

2012 quo = domain.quo

2013 terms = [ (monom, quo(coeff, x)) for monom, coeff in f.iterterms() ]

2014 else:

2015 terms = [ (monom, coeff // x) for monom, coeff in f.iterterms() if not (coeff % x) ]

2016

2017 return f.new(terms)

2018

2019 def quo_term(f, term):

2020 monom, coeff = term

2021

2022 if not coeff:

2023 raise ZeroDivisionError("polynomial division")

2024 elif not f:

2025 return f.ring.zero

2026 elif monom == f.ring.zero_monom:

2027 return f.quo_ground(coeff)

2028

2029 term_div = f._term_div()

2030

2031 terms = [ term_div(t, term) for t in f.iterterms() ]

2032 return f.new([ t for t in terms if t is not None ])

2033

2034 def trunc_ground(f, p):

2035 if f.ring.domain.is_ZZ:

2036 terms = []

2037

2038 for monom, coeff in f.iterterms():

2039 coeff = coeff % p

2040

2041 if coeff > p // 2:

2042 coeff = coeff - p

2043

2044 terms.append((monom, coeff))

2045 else:

2046 terms = [ (monom, coeff % p) for monom, coeff in f.iterterms() ]

2047

2048 poly = f.new(terms)

2049 poly.strip_zero()

2050 return poly

2051

2052 rem_ground = trunc_ground

2053

2054 def extract_ground(self, g):

2055 f = self

2056 fc = f.content()

2057 gc = g.content()

2058

2059 gcd = f.ring.domain.gcd(fc, gc)

2060

2061 f = f.quo_ground(gcd)

2062 g = g.quo_ground(gcd)

2063

2064 return gcd, f, g

2065

2066 def _norm(f, norm_func):

2067 if not f:

2068 return f.ring.domain.zero

2069 else:

2070 ground_abs = f.ring.domain.abs

2071 return norm_func([ ground_abs(coeff) for coeff in f.itercoeffs() ])

2072

2073 def max_norm(f):

2074 return f._norm(max)

2075

2076 def l1_norm(f):

2077 return f._norm(sum)

2078

2079 def deflate(f, *G):

2080 ring = f.ring

2081 polys = [f] + list(G)

2082

2083 J = [0]*ring.ngens

2084

2085 for p in polys:

2086 for monom in p.itermonoms():

2087 for i, m in enumerate(monom):

2088 J[i] = igcd(J[i], m)

2089

2090 for i, b in enumerate(J):

2091 if not b:

2092 J[i] = 1

2093

2094 J = tuple(J)

2095

2096 if all(b == 1 for b in J):

2097 return J, polys

2098

2099 H = []

2100

2101 for p in polys:

2102 h = ring.zero

2103

2104 for I, coeff in p.iterterms():

2105 N = [ i // j for i, j in zip(I, J) ]

2106 h[tuple(N)] = coeff

2107

2108 H.append(h)

2109

2110 return J, H

2111

2112 def inflate(f, J):

2113 poly = f.ring.zero

2114

2115 for I, coeff in f.iterterms():

2116 N = [ i*j for i, j in zip(I, J) ]

2117 poly[tuple(N)] = coeff

2118

2119 return poly

2120

2121 def lcm(self, g):

2122 f = self

2123 domain = f.ring.domain

2124

2125 if not domain.is_Field:

2126 fc, f = f.primitive()

2127 gc, g = g.primitive()

2128 c = domain.lcm(fc, gc)

2129

2130 h = (f*g).quo(f.gcd(g))

2131

2132 if not domain.is_Field:

2133 return h.mul_ground(c)

2134 else:

2135 return h.monic()

2136

2137 def gcd(f, g):

2138 return f.cofactors(g)[0]

2139

2140 def cofactors(f, g):

2141 if not f and not g:

2142 zero = f.ring.zero

2143 return zero, zero, zero

2144 elif not f:

2145 h, cff, cfg = f._gcd_zero(g)

2146 return h, cff, cfg

2147 elif not g:

2148 h, cfg, cff = g._gcd_zero(f)

2149 return h, cff, cfg

2150 elif len(f) == 1:

2151 h, cff, cfg = f._gcd_monom(g)

2152 return h, cff, cfg

2153 elif len(g) == 1:

2154 h, cfg, cff = g._gcd_monom(f)

2155 return h, cff, cfg

2156

2157 J, (f, g) = f.deflate(g)

2158 h, cff, cfg = f._gcd(g)

2159

2160 return (h.inflate(J), cff.inflate(J), cfg.inflate(J))

2161

2162 def _gcd_zero(f, g):

2163 one, zero = f.ring.one, f.ring.zero

2164 if g.is_nonnegative:

2165 return g, zero, one

2166 else:

2167 return -g, zero, -one

2168

2169 def _gcd_monom(f, g):

2170 ring = f.ring

2171 ground_gcd = ring.domain.gcd

2172 ground_quo = ring.domain.quo

2173 monomial_gcd = ring.monomial_gcd

2174 monomial_ldiv = ring.monomial_ldiv

2175 mf, cf = list(f.iterterms())[0]

2176 _mgcd, _cgcd = mf, cf

2177 for mg, cg in g.iterterms():

2178 _mgcd = monomial_gcd(_mgcd, mg)

2179 _cgcd = ground_gcd(_cgcd, cg)

2180 h = f.new([(_mgcd, _cgcd)])

2181 cff = f.new([(monomial_ldiv(mf, _mgcd), ground_quo(cf, _cgcd))])

2182 cfg = f.new([(monomial_ldiv(mg, _mgcd), ground_quo(cg, _cgcd)) for mg, cg in g.iterterms()])

2183 return h, cff, cfg

2184

2185 def _gcd(f, g):

2186 ring = f.ring

2187

2188 if ring.domain.is_QQ:

2189 return f._gcd_QQ(g)

2190 elif ring.domain.is_ZZ:

2191 return f._gcd_ZZ(g)

2192 else: # TODO: don't use dense representation (port PRS algorithms)

2193 return ring.dmp_inner_gcd(f, g)

2194

2195 def _gcd_ZZ(f, g):

2196 return heugcd(f, g)

2197

2198 def _gcd_QQ(self, g):

2199 f = self

2200 ring = f.ring

2201 new_ring = ring.clone(domain=ring.domain.get_ring())

2202

2203 cf, f = f.clear_denoms()

2204 cg, g = g.clear_denoms()

2205

2206 f = f.set_ring(new_ring)

2207 g = g.set_ring(new_ring)

2208

2209 h, cff, cfg = f._gcd_ZZ(g)

2210

2211 h = h.set_ring(ring)

2212 c, h = h.LC, h.monic()

2213

2214 cff = cff.set_ring(ring).mul_ground(ring.domain.quo(c, cf))

2215 cfg = cfg.set_ring(ring).mul_ground(ring.domain.quo(c, cg))

2216

2217 return h, cff, cfg

2218

2219 def cancel(self, g):

2220 """

2221 Cancel common factors in a rational function ``f/g``.

2222

2223 Examples

2224 ========

2225

2226 >>> from sympy.polys import ring, ZZ

2227 >>> R, x,y = ring("x,y", ZZ)

2228

2229 >>> (2*x**2 - 2).cancel(x**2 - 2*x + 1)

2230 (2*x + 2, x - 1)

2231

2232 """

2233 f = self

2234 ring = f.ring

2235

2236 if not f:

2237 return f, ring.one

2238

2239 domain = ring.domain

2240

2241 if not (domain.is_Field and domain.has_assoc_Ring):

2242 _, p, q = f.cofactors(g)

2243 else:

2244 new_ring = ring.clone(domain=domain.get_ring())

2245

2246 cq, f = f.clear_denoms()

2247 cp, g = g.clear_denoms()

2248

2249 f = f.set_ring(new_ring)

2250 g = g.set_ring(new_ring)

2251

2252 _, p, q = f.cofactors(g)

2253 _, cp, cq = new_ring.domain.cofactors(cp, cq)

2254

2255 p = p.set_ring(ring)

2256 q = q.set_ring(ring)

2257

2258 p = p.mul_ground(cp)

2259 q = q.mul_ground(cq)

2260

2261 # Make canonical with respect to sign or quadrant in the case of ZZ_I

2262 # or QQ_I. This ensures that the LC of the denominator is canonical by

2263 # multiplying top and bottom by a unit of the ring.

2264 u = q.canonical_unit()

2265 if u == domain.one:

2266 p, q = p, q

2267 elif u == -domain.one:

2268 p, q = -p, -q

2269 else:

2270 p = p.mul_ground(u)

2271 q = q.mul_ground(u)

2272

2273 return p, q

2274

2275 def canonical_unit(f):

2276 domain = f.ring.domain

2277 return domain.canonical_unit(f.LC)

2278

2279 def diff(f, x):

2280 """Computes partial derivative in ``x``.

2281

2282 Examples

2283 ========

2284

2285 >>> from sympy.polys.rings import ring

2286 >>> from sympy.polys.domains import ZZ

2287

2288 >>> _, x, y = ring("x,y", ZZ)

2289 >>> p = x + x**2*y**3

2290 >>> p.diff(x)

2291 2*x*y**3 + 1

2292

2293 """

2294 ring = f.ring

2295 i = ring.index(x)

2296 m = ring.monomial_basis(i)

2297 g = ring.zero

2298 for expv, coeff in f.iterterms():

2299 if expv[i]:

2300 e = ring.monomial_ldiv(expv, m)

2301 g[e] = ring.domain_new(coeff*expv[i])

2302 return g

2303

2304 def __call__(f, *values):

2305 if 0 < len(values) <= f.ring.ngens:

2306 return f.evaluate(list(zip(f.ring.gens, values)))

2307 else:

2308 raise ValueError("expected at least 1 and at most %s values, got %s" % (f.ring.ngens, len(values)))

2309

2310 def evaluate(self, x, a=None):

2311 f = self

2312

2313 if isinstance(x, list) and a is None:

2314 (X, a), x = x[0], x[1:]

2315 f = f.evaluate(X, a)

2316

2317 if not x:

2318 return f

2319 else:

2320 x = [ (Y.drop(X), a) for (Y, a) in x ]

2321 return f.evaluate(x)

2322

2323 ring = f.ring

2324 i = ring.index(x)

2325 a = ring.domain.convert(a)

2326

2327 if ring.ngens == 1:

2328 result = ring.domain.zero

2329

2330 for (n,), coeff in f.iterterms():

2331 result += coeff*a**n

2332

2333 return result

2334 else:

2335 poly = ring.drop(x).zero

2336

2337 for monom, coeff in f.iterterms():

2338 n, monom = monom[i], monom[:i] + monom[i+1:]

2339 coeff = coeff*a**n

2340

2341 if monom in poly:

2342 coeff = coeff + poly[monom]

2343

2344 if coeff:

2345 poly[monom] = coeff

2346 else:

2347 del poly[monom]

2348 else:

2349 if coeff:

2350 poly[monom] = coeff

2351

2352 return poly

2353

2354 def subs(self, x, a=None):

2355 f = self

2356

2357 if isinstance(x, list) and a is None:

2358 for X, a in x:

2359 f = f.subs(X, a)

2360 return f

2361

2362 ring = f.ring

2363 i = ring.index(x)

2364 a = ring.domain.convert(a)

2365

2366 if ring.ngens == 1:

2367 result = ring.domain.zero

2368

2369 for (n,), coeff in f.iterterms():

2370 result += coeff*a**n

2371

2372 return ring.ground_new(result)

2373 else:

2374 poly = ring.zero

2375

2376 for monom, coeff in f.iterterms():

2377 n, monom = monom[i], monom[:i] + (0,) + monom[i+1:]

2378 coeff = coeff*a**n

2379

2380 if monom in poly:

2381 coeff = coeff + poly[monom]

2382

2383 if coeff:

2384 poly[monom] = coeff

2385 else:

2386 del poly[monom]

2387 else:

2388 if coeff:

2389 poly[monom] = coeff

2390

2391 return poly

2392

2393 def symmetrize(self):

2394 r"""

2395 Rewrite *self* in terms of elementary symmetric polynomials.

2396

2397 Explanation

2398 ===========

2399

2400 If this :py:class:`~.PolyElement` belongs to a ring of $n$ variables,

2401 we can try to write it as a function of the elementary symmetric

2402 polynomials on $n$ variables. We compute a symmetric part, and a

2403 remainder for any part we were not able to symmetrize.

2404

2405 Examples

2406 ========

2407

2408 >>> from sympy.polys.rings import ring

2409 >>> from sympy.polys.domains import ZZ

2410 >>> R, x, y = ring("x,y", ZZ)

2411

2412 >>> f = x**2 + y**2

2413 >>> f.symmetrize()

2414 (x**2 - 2*y, 0, [(x, x + y), (y, x*y)])

2415

2416 >>> f = x**2 - y**2

2417 >>> f.symmetrize()

2418 (x**2 - 2*y, -2*y**2, [(x, x + y), (y, x*y)])

2419

2420 Returns

2421 =======

2422

2423 Triple ``(p, r, m)``

2424 ``p`` is a :py:class:`~.PolyElement` that represents our attempt

2425 to express *self* as a function of elementary symmetric

2426 polynomials. Each variable in ``p`` stands for one of the

2427 elementary symmetric polynomials. The correspondence is given

2428 by ``m``.

2429

2430 ``r`` is the remainder.

2431

2432 ``m`` is a list of pairs, giving the mapping from variables in

2433 ``p`` to elementary symmetric polynomials.

2434

2435 The triple satisfies the equation ``p.compose(m) + r == self``.

2436 If the remainder ``r`` is zero, *self* is symmetric. If it is

2437 nonzero, we were not able to represent *self* as symmetric.

2438

2439 See Also

2440 ========

2441

2442 sympy.polys.polyfuncs.symmetrize

2443

2444 References

2445 ==========

2446

2447 .. [1] Lauer, E. Algorithms for symmetrical polynomials, Proc. 1976

2448 ACM Symp. on Symbolic and Algebraic Computing, NY 242-247.

2449 https://dl.acm.org/doi/pdf/10.1145/800205.806342

2450

2451 """

2452 f = self.copy()

2453 ring = f.ring

2454 n = ring.ngens

2455

2456 if not n:

2457 return f, ring.zero, []

2458

2459 polys = [ring.symmetric_poly(i+1) for i in range(n)]

2460

2461 poly_powers = {}

2462 def get_poly_power(i, n):

2463 if (i, n) not in poly_powers:

2464 poly_powers[(i, n)] = polys[i]**n

2465 return poly_powers[(i, n)]

2466

2467 indices = list(range(n - 1))

2468 weights = list(range(n, 0, -1))

2469

2470 symmetric = ring.zero

2471

2472 while f:

2473 _height, _monom, _coeff = -1, None, None

2474

2475 for i, (monom, coeff) in enumerate(f.terms()):

2476 if all(monom[i] >= monom[i + 1] for i in indices):

2477 height = max([n*m for n, m in zip(weights, monom)])

2478

2479 if height > _height:

2480 _height, _monom, _coeff = height, monom, coeff

2481

2482 if _height != -1:

2483 monom, coeff = _monom, _coeff

2484 else:

2485 break

2486

2487 exponents = []

2488 for m1, m2 in zip(monom, monom[1:] + (0,)):

2489 exponents.append(m1 - m2)

2490

2491 symmetric += ring.term_new(tuple(exponents), coeff)

2492

2493 product = coeff

2494 for i, n in enumerate(exponents):

2495 product *= get_poly_power(i, n)

2496 f -= product

2497

2498 mapping = list(zip(ring.gens, polys))

2499

2500 return symmetric, f, mapping

2501

2502 def compose(f, x, a=None):

2503 ring = f.ring

2504 poly = ring.zero

2505 gens_map = dict(zip(ring.gens, range(ring.ngens)))

2506

2507 if a is not None:

2508 replacements = [(x, a)]

2509 else:

2510 if isinstance(x, list):

2511 replacements = list(x)

2512 elif isinstance(x, dict):

2513 replacements = sorted(x.items(), key=lambda k: gens_map[k[0]])

2514 else:

2515 raise ValueError("expected a generator, value pair a sequence of such pairs")

2516

2517 for k, (x, g) in enumerate(replacements):

2518 replacements[k] = (gens_map[x], ring.ring_new(g))

2519

2520 for monom, coeff in f.iterterms():

2521 monom = list(monom)

2522 subpoly = ring.one

2523

2524 for i, g in replacements:

2525 n, monom[i] = monom[i], 0

2526 if n:

2527 subpoly *= g**n

2528

2529 subpoly = subpoly.mul_term((tuple(monom), coeff))

2530 poly += subpoly

2531

2532 return poly

2533

2534 # TODO: following methods should point to polynomial

2535 # representation independent algorithm implementations.

2536

2537 def pdiv(f, g):

2538 return f.ring.dmp_pdiv(f, g)

2539

2540 def prem(f, g):

2541 return f.ring.dmp_prem(f, g)

2542

2543 def pquo(f, g):

2544 return f.ring.dmp_quo(f, g)

2545

2546 def pexquo(f, g):

2547 return f.ring.dmp_exquo(f, g)

2548

2549 def half_gcdex(f, g):

2550 return f.ring.dmp_half_gcdex(f, g)

2551

2552 def gcdex(f, g):

2553 return f.ring.dmp_gcdex(f, g)

2554

2555 def subresultants(f, g):

2556 return f.ring.dmp_subresultants(f, g)

2557

2558 def resultant(f, g):

2559 return f.ring.dmp_resultant(f, g)

2560

2561 def discriminant(f):

2562 return f.ring.dmp_discriminant(f)

2563

2564 def decompose(f):

2565 if f.ring.is_univariate:

2566 return f.ring.dup_decompose(f)

2567 else:

2568 raise MultivariatePolynomialError("polynomial decomposition")

2569

2570 def shift(f, a):

2571 if f.ring.is_univariate:

2572 return f.ring.dup_shift(f, a)

2573 else:

2574 raise MultivariatePolynomialError("polynomial shift")

2575

2576 def sturm(f):

2577 if f.ring.is_univariate:

2578 return f.ring.dup_sturm(f)

2579 else:

2580 raise MultivariatePolynomialError("sturm sequence")

2581

2582 def gff_list(f):

2583 return f.ring.dmp_gff_list(f)

2584

2585 def sqf_norm(f):

2586 return f.ring.dmp_sqf_norm(f)

2587

2588 def sqf_part(f):

2589 return f.ring.dmp_sqf_part(f)

2590

2591 def sqf_list(f, all=False):

2592 return f.ring.dmp_sqf_list(f, all=all)

2593

2594 def factor_list(f):

2595 return f.ring.dmp_factor_list(f)