1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467 |
- """Sparse polynomial rings. """
- from typing import Any, Dict as tDict
- from operator import add, mul, lt, le, gt, ge
- from functools import reduce
- from types import GeneratorType
- from sympy.core.expr import Expr
- from sympy.core.numbers import igcd, oo
- from sympy.core.symbol import Symbol, symbols as _symbols
- from sympy.core.sympify import CantSympify, sympify
- from sympy.ntheory.multinomial import multinomial_coefficients
- from sympy.polys.compatibility import IPolys
- from sympy.polys.constructor import construct_domain
- from sympy.polys.densebasic import dmp_to_dict, dmp_from_dict
- from sympy.polys.domains.domainelement import DomainElement
- from sympy.polys.domains.polynomialring import PolynomialRing
- from sympy.polys.heuristicgcd import heugcd
- from sympy.polys.monomials import MonomialOps
- from sympy.polys.orderings import lex
- from sympy.polys.polyerrors import (
- CoercionFailed, GeneratorsError,
- ExactQuotientFailed, MultivariatePolynomialError)
- from sympy.polys.polyoptions import (Domain as DomainOpt,
- Order as OrderOpt, build_options)
- from sympy.polys.polyutils import (expr_from_dict, _dict_reorder,
- _parallel_dict_from_expr)
- from sympy.printing.defaults import DefaultPrinting
- from sympy.utilities import public
- from sympy.utilities.iterables import is_sequence
- from sympy.utilities.magic import pollute
- @public
- def ring(symbols, domain, order=lex):
- """Construct a polynomial ring returning ``(ring, x_1, ..., x_n)``.
- Parameters
- ==========
- symbols : str
- Symbol/Expr or sequence of str, Symbol/Expr (non-empty)
- domain : :class:`~.Domain` or coercible
- order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.orderings import lex
- >>> R, x, y, z = ring("x,y,z", ZZ, lex)
- >>> R
- Polynomial ring in x, y, z over ZZ with lex order
- >>> x + y + z
- x + y + z
- >>> type(_)
- <class 'sympy.polys.rings.PolyElement'>
- """
- _ring = PolyRing(symbols, domain, order)
- return (_ring,) + _ring.gens
- @public
- def xring(symbols, domain, order=lex):
- """Construct a polynomial ring returning ``(ring, (x_1, ..., x_n))``.
- Parameters
- ==========
- symbols : str
- Symbol/Expr or sequence of str, Symbol/Expr (non-empty)
- domain : :class:`~.Domain` or coercible
- order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``
- Examples
- ========
- >>> from sympy.polys.rings import xring
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.orderings import lex
- >>> R, (x, y, z) = xring("x,y,z", ZZ, lex)
- >>> R
- Polynomial ring in x, y, z over ZZ with lex order
- >>> x + y + z
- x + y + z
- >>> type(_)
- <class 'sympy.polys.rings.PolyElement'>
- """
- _ring = PolyRing(symbols, domain, order)
- return (_ring, _ring.gens)
- @public
- def vring(symbols, domain, order=lex):
- """Construct a polynomial ring and inject ``x_1, ..., x_n`` into the global namespace.
- Parameters
- ==========
- symbols : str
- Symbol/Expr or sequence of str, Symbol/Expr (non-empty)
- domain : :class:`~.Domain` or coercible
- order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``
- Examples
- ========
- >>> from sympy.polys.rings import vring
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.orderings import lex
- >>> vring("x,y,z", ZZ, lex)
- Polynomial ring in x, y, z over ZZ with lex order
- >>> x + y + z # noqa:
- x + y + z
- >>> type(_)
- <class 'sympy.polys.rings.PolyElement'>
- """
- _ring = PolyRing(symbols, domain, order)
- pollute([ sym.name for sym in _ring.symbols ], _ring.gens)
- return _ring
- @public
- def sring(exprs, *symbols, **options):
- """Construct a ring deriving generators and domain from options and input expressions.
- Parameters
- ==========
- exprs : :class:`~.Expr` or sequence of :class:`~.Expr` (sympifiable)
- symbols : sequence of :class:`~.Symbol`/:class:`~.Expr`
- options : keyword arguments understood by :class:`~.Options`
- Examples
- ========
- >>> from sympy import sring, symbols
- >>> x, y, z = symbols("x,y,z")
- >>> R, f = sring(x + 2*y + 3*z)
- >>> R
- Polynomial ring in x, y, z over ZZ with lex order
- >>> f
- x + 2*y + 3*z
- >>> type(_)
- <class 'sympy.polys.rings.PolyElement'>
- """
- single = False
- if not is_sequence(exprs):
- exprs, single = [exprs], True
- exprs = list(map(sympify, exprs))
- opt = build_options(symbols, options)
- # TODO: rewrite this so that it doesn't use expand() (see poly()).
- reps, opt = _parallel_dict_from_expr(exprs, opt)
- if opt.domain is None:
- coeffs = sum([ list(rep.values()) for rep in reps ], [])
- opt.domain, coeffs_dom = construct_domain(coeffs, opt=opt)
- coeff_map = dict(zip(coeffs, coeffs_dom))
- reps = [{m: coeff_map[c] for m, c in rep.items()} for rep in reps]
- _ring = PolyRing(opt.gens, opt.domain, opt.order)
- polys = list(map(_ring.from_dict, reps))
- if single:
- return (_ring, polys[0])
- else:
- return (_ring, polys)
- def _parse_symbols(symbols):
- if isinstance(symbols, str):
- return _symbols(symbols, seq=True) if symbols else ()
- elif isinstance(symbols, Expr):
- return (symbols,)
- elif is_sequence(symbols):
- if all(isinstance(s, str) for s in symbols):
- return _symbols(symbols)
- elif all(isinstance(s, Expr) for s in symbols):
- return symbols
- raise GeneratorsError("expected a string, Symbol or expression or a non-empty sequence of strings, Symbols or expressions")
- _ring_cache = {} # type: tDict[Any, Any]
- class PolyRing(DefaultPrinting, IPolys):
- """Multivariate distributed polynomial ring. """
- def __new__(cls, symbols, domain, order=lex):
- symbols = tuple(_parse_symbols(symbols))
- ngens = len(symbols)
- domain = DomainOpt.preprocess(domain)
- order = OrderOpt.preprocess(order)
- _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
- obj = _ring_cache.get(_hash_tuple)
- if obj is None:
- if domain.is_Composite and set(symbols) & set(domain.symbols):
- raise GeneratorsError("polynomial ring and it's ground domain share generators")
- obj = object.__new__(cls)
- obj._hash_tuple = _hash_tuple
- obj._hash = hash(_hash_tuple)
- obj.dtype = type("PolyElement", (PolyElement,), {"ring": obj})
- obj.symbols = symbols
- obj.ngens = ngens
- obj.domain = domain
- obj.order = order
- obj.zero_monom = (0,)*ngens
- obj.gens = obj._gens()
- obj._gens_set = set(obj.gens)
- obj._one = [(obj.zero_monom, domain.one)]
- if ngens:
- # These expect monomials in at least one variable
- codegen = MonomialOps(ngens)
- obj.monomial_mul = codegen.mul()
- obj.monomial_pow = codegen.pow()
- obj.monomial_mulpow = codegen.mulpow()
- obj.monomial_ldiv = codegen.ldiv()
- obj.monomial_div = codegen.div()
- obj.monomial_lcm = codegen.lcm()
- obj.monomial_gcd = codegen.gcd()
- else:
- monunit = lambda a, b: ()
- obj.monomial_mul = monunit
- obj.monomial_pow = monunit
- obj.monomial_mulpow = lambda a, b, c: ()
- obj.monomial_ldiv = monunit
- obj.monomial_div = monunit
- obj.monomial_lcm = monunit
- obj.monomial_gcd = monunit
- if order is lex:
- obj.leading_expv = max
- else:
- obj.leading_expv = lambda f: max(f, key=order)
- for symbol, generator in zip(obj.symbols, obj.gens):
- if isinstance(symbol, Symbol):
- name = symbol.name
- if not hasattr(obj, name):
- setattr(obj, name, generator)
- _ring_cache[_hash_tuple] = obj
- return obj
- def _gens(self):
- """Return a list of polynomial generators. """
- one = self.domain.one
- _gens = []
- for i in range(self.ngens):
- expv = self.monomial_basis(i)
- poly = self.zero
- poly[expv] = one
- _gens.append(poly)
- return tuple(_gens)
- def __getnewargs__(self):
- return (self.symbols, self.domain, self.order)
- def __getstate__(self):
- state = self.__dict__.copy()
- del state["leading_expv"]
- for key, value in state.items():
- if key.startswith("monomial_"):
- del state[key]
- return state
- def __hash__(self):
- return self._hash
- def __eq__(self, other):
- return isinstance(other, PolyRing) and \
- (self.symbols, self.domain, self.ngens, self.order) == \
- (other.symbols, other.domain, other.ngens, other.order)
- def __ne__(self, other):
- return not self == other
- def clone(self, symbols=None, domain=None, order=None):
- return self.__class__(symbols or self.symbols, domain or self.domain, order or self.order)
- def monomial_basis(self, i):
- """Return the ith-basis element. """
- basis = [0]*self.ngens
- basis[i] = 1
- return tuple(basis)
- @property
- def zero(self):
- return self.dtype()
- @property
- def one(self):
- return self.dtype(self._one)
- def domain_new(self, element, orig_domain=None):
- return self.domain.convert(element, orig_domain)
- def ground_new(self, coeff):
- return self.term_new(self.zero_monom, coeff)
- def term_new(self, monom, coeff):
- coeff = self.domain_new(coeff)
- poly = self.zero
- if coeff:
- poly[monom] = coeff
- return poly
- def ring_new(self, element):
- if isinstance(element, PolyElement):
- if self == element.ring:
- return element
- elif isinstance(self.domain, PolynomialRing) and self.domain.ring == element.ring:
- return self.ground_new(element)
- else:
- raise NotImplementedError("conversion")
- elif isinstance(element, str):
- raise NotImplementedError("parsing")
- elif isinstance(element, dict):
- return self.from_dict(element)
- elif isinstance(element, list):
- try:
- return self.from_terms(element)
- except ValueError:
- return self.from_list(element)
- elif isinstance(element, Expr):
- return self.from_expr(element)
- else:
- return self.ground_new(element)
- __call__ = ring_new
- def from_dict(self, element, orig_domain=None):
- domain_new = self.domain_new
- poly = self.zero
- for monom, coeff in element.items():
- coeff = domain_new(coeff, orig_domain)
- if coeff:
- poly[monom] = coeff
- return poly
- def from_terms(self, element, orig_domain=None):
- return self.from_dict(dict(element), orig_domain)
- def from_list(self, element):
- return self.from_dict(dmp_to_dict(element, self.ngens-1, self.domain))
- def _rebuild_expr(self, expr, mapping):
- domain = self.domain
- def _rebuild(expr):
- generator = mapping.get(expr)
- if generator is not None:
- return generator
- elif expr.is_Add:
- return reduce(add, list(map(_rebuild, expr.args)))
- elif expr.is_Mul:
- return reduce(mul, list(map(_rebuild, expr.args)))
- else:
- # XXX: Use as_base_exp() to handle Pow(x, n) and also exp(n)
- # XXX: E can be a generator e.g. sring([exp(2)]) -> ZZ[E]
- base, exp = expr.as_base_exp()
- if exp.is_Integer and exp > 1:
- return _rebuild(base)**int(exp)
- else:
- return self.ground_new(domain.convert(expr))
- return _rebuild(sympify(expr))
- def from_expr(self, expr):
- mapping = dict(list(zip(self.symbols, self.gens)))
- try:
- poly = self._rebuild_expr(expr, mapping)
- except CoercionFailed:
- raise ValueError("expected an expression convertible to a polynomial in %s, got %s" % (self, expr))
- else:
- return self.ring_new(poly)
- def index(self, gen):
- """Compute index of ``gen`` in ``self.gens``. """
- if gen is None:
- if self.ngens:
- i = 0
- else:
- i = -1 # indicate impossible choice
- elif isinstance(gen, int):
- i = gen
- if 0 <= i and i < self.ngens:
- pass
- elif -self.ngens <= i and i <= -1:
- i = -i - 1
- else:
- raise ValueError("invalid generator index: %s" % gen)
- elif isinstance(gen, self.dtype):
- try:
- i = self.gens.index(gen)
- except ValueError:
- raise ValueError("invalid generator: %s" % gen)
- elif isinstance(gen, str):
- try:
- i = self.symbols.index(gen)
- except ValueError:
- raise ValueError("invalid generator: %s" % gen)
- else:
- raise ValueError("expected a polynomial generator, an integer, a string or None, got %s" % gen)
- return i
- def drop(self, *gens):
- """Remove specified generators from this ring. """
- indices = set(map(self.index, gens))
- symbols = [ s for i, s in enumerate(self.symbols) if i not in indices ]
- if not symbols:
- return self.domain
- else:
- return self.clone(symbols=symbols)
- def __getitem__(self, key):
- symbols = self.symbols[key]
- if not symbols:
- return self.domain
- else:
- return self.clone(symbols=symbols)
- def to_ground(self):
- # TODO: should AlgebraicField be a Composite domain?
- if self.domain.is_Composite or hasattr(self.domain, 'domain'):
- return self.clone(domain=self.domain.domain)
- else:
- raise ValueError("%s is not a composite domain" % self.domain)
- def to_domain(self):
- return PolynomialRing(self)
- def to_field(self):
- from sympy.polys.fields import FracField
- return FracField(self.symbols, self.domain, self.order)
- @property
- def is_univariate(self):
- return len(self.gens) == 1
- @property
- def is_multivariate(self):
- return len(self.gens) > 1
- def add(self, *objs):
- """
- Add a sequence of polynomials or containers of polynomials.
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> R, x = ring("x", ZZ)
- >>> R.add([ x**2 + 2*i + 3 for i in range(4) ])
- 4*x**2 + 24
- >>> _.factor_list()
- (4, [(x**2 + 6, 1)])
- """
- p = self.zero
- for obj in objs:
- if is_sequence(obj, include=GeneratorType):
- p += self.add(*obj)
- else:
- p += obj
- return p
- def mul(self, *objs):
- """
- Multiply a sequence of polynomials or containers of polynomials.
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> R, x = ring("x", ZZ)
- >>> R.mul([ x**2 + 2*i + 3 for i in range(4) ])
- x**8 + 24*x**6 + 206*x**4 + 744*x**2 + 945
- >>> _.factor_list()
- (1, [(x**2 + 3, 1), (x**2 + 5, 1), (x**2 + 7, 1), (x**2 + 9, 1)])
- """
- p = self.one
- for obj in objs:
- if is_sequence(obj, include=GeneratorType):
- p *= self.mul(*obj)
- else:
- p *= obj
- return p
- def drop_to_ground(self, *gens):
- r"""
- Remove specified generators from the ring and inject them into
- its domain.
- """
- indices = set(map(self.index, gens))
- symbols = [s for i, s in enumerate(self.symbols) if i not in indices]
- gens = [gen for i, gen in enumerate(self.gens) if i not in indices]
- if not symbols:
- return self
- else:
- return self.clone(symbols=symbols, domain=self.drop(*gens))
- def compose(self, other):
- """Add the generators of ``other`` to ``self``"""
- if self != other:
- syms = set(self.symbols).union(set(other.symbols))
- return self.clone(symbols=list(syms))
- else:
- return self
- def add_gens(self, symbols):
- """Add the elements of ``symbols`` as generators to ``self``"""
- syms = set(self.symbols).union(set(symbols))
- return self.clone(symbols=list(syms))
- class PolyElement(DomainElement, DefaultPrinting, CantSympify, dict):
- """Element of multivariate distributed polynomial ring. """
- def new(self, init):
- return self.__class__(init)
- def parent(self):
- return self.ring.to_domain()
- def __getnewargs__(self):
- return (self.ring, list(self.iterterms()))
- _hash = None
- def __hash__(self):
- # XXX: This computes a hash of a dictionary, but currently we don't
- # protect dictionary from being changed so any use site modifications
- # will make hashing go wrong. Use this feature with caution until we
- # figure out how to make a safe API without compromising speed of this
- # low-level class.
- _hash = self._hash
- if _hash is None:
- self._hash = _hash = hash((self.ring, frozenset(self.items())))
- return _hash
- def copy(self):
- """Return a copy of polynomial self.
- Polynomials are mutable; if one is interested in preserving
- a polynomial, and one plans to use inplace operations, one
- can copy the polynomial. This method makes a shallow copy.
- Examples
- ========
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.rings import ring
- >>> R, x, y = ring('x, y', ZZ)
- >>> p = (x + y)**2
- >>> p1 = p.copy()
- >>> p2 = p
- >>> p[R.zero_monom] = 3
- >>> p
- x**2 + 2*x*y + y**2 + 3
- >>> p1
- x**2 + 2*x*y + y**2
- >>> p2
- x**2 + 2*x*y + y**2 + 3
- """
- return self.new(self)
- def set_ring(self, new_ring):
- if self.ring == new_ring:
- return self
- elif self.ring.symbols != new_ring.symbols:
- terms = list(zip(*_dict_reorder(self, self.ring.symbols, new_ring.symbols)))
- return new_ring.from_terms(terms, self.ring.domain)
- else:
- return new_ring.from_dict(self, self.ring.domain)
- def as_expr(self, *symbols):
- if symbols and len(symbols) != self.ring.ngens:
- raise ValueError("not enough symbols, expected %s got %s" % (self.ring.ngens, len(symbols)))
- else:
- symbols = self.ring.symbols
- return expr_from_dict(self.as_expr_dict(), *symbols)
- def as_expr_dict(self):
- to_sympy = self.ring.domain.to_sympy
- return {monom: to_sympy(coeff) for monom, coeff in self.iterterms()}
- def clear_denoms(self):
- domain = self.ring.domain
- if not domain.is_Field or not domain.has_assoc_Ring:
- return domain.one, self
- ground_ring = domain.get_ring()
- common = ground_ring.one
- lcm = ground_ring.lcm
- denom = domain.denom
- for coeff in self.values():
- common = lcm(common, denom(coeff))
- poly = self.new([ (k, v*common) for k, v in self.items() ])
- return common, poly
- def strip_zero(self):
- """Eliminate monomials with zero coefficient. """
- for k, v in list(self.items()):
- if not v:
- del self[k]
- def __eq__(p1, p2):
- """Equality test for polynomials.
- Examples
- ========
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.rings import ring
- >>> _, x, y = ring('x, y', ZZ)
- >>> p1 = (x + y)**2 + (x - y)**2
- >>> p1 == 4*x*y
- False
- >>> p1 == 2*(x**2 + y**2)
- True
- """
- if not p2:
- return not p1
- elif isinstance(p2, PolyElement) and p2.ring == p1.ring:
- return dict.__eq__(p1, p2)
- elif len(p1) > 1:
- return False
- else:
- return p1.get(p1.ring.zero_monom) == p2
- def __ne__(p1, p2):
- return not p1 == p2
- def almosteq(p1, p2, tolerance=None):
- """Approximate equality test for polynomials. """
- ring = p1.ring
- if isinstance(p2, ring.dtype):
- if set(p1.keys()) != set(p2.keys()):
- return False
- almosteq = ring.domain.almosteq
- for k in p1.keys():
- if not almosteq(p1[k], p2[k], tolerance):
- return False
- return True
- elif len(p1) > 1:
- return False
- else:
- try:
- p2 = ring.domain.convert(p2)
- except CoercionFailed:
- return False
- else:
- return ring.domain.almosteq(p1.const(), p2, tolerance)
- def sort_key(self):
- return (len(self), self.terms())
- def _cmp(p1, p2, op):
- if isinstance(p2, p1.ring.dtype):
- return op(p1.sort_key(), p2.sort_key())
- else:
- return NotImplemented
- def __lt__(p1, p2):
- return p1._cmp(p2, lt)
- def __le__(p1, p2):
- return p1._cmp(p2, le)
- def __gt__(p1, p2):
- return p1._cmp(p2, gt)
- def __ge__(p1, p2):
- return p1._cmp(p2, ge)
- def _drop(self, gen):
- ring = self.ring
- i = ring.index(gen)
- if ring.ngens == 1:
- return i, ring.domain
- else:
- symbols = list(ring.symbols)
- del symbols[i]
- return i, ring.clone(symbols=symbols)
- def drop(self, gen):
- i, ring = self._drop(gen)
- if self.ring.ngens == 1:
- if self.is_ground:
- return self.coeff(1)
- else:
- raise ValueError("Cannot drop %s" % gen)
- else:
- poly = ring.zero
- for k, v in self.items():
- if k[i] == 0:
- K = list(k)
- del K[i]
- poly[tuple(K)] = v
- else:
- raise ValueError("Cannot drop %s" % gen)
- return poly
- def _drop_to_ground(self, gen):
- ring = self.ring
- i = ring.index(gen)
- symbols = list(ring.symbols)
- del symbols[i]
- return i, ring.clone(symbols=symbols, domain=ring[i])
- def drop_to_ground(self, gen):
- if self.ring.ngens == 1:
- raise ValueError("Cannot drop only generator to ground")
- i, ring = self._drop_to_ground(gen)
- poly = ring.zero
- gen = ring.domain.gens[0]
- for monom, coeff in self.iterterms():
- mon = monom[:i] + monom[i+1:]
- if mon not in poly:
- poly[mon] = (gen**monom[i]).mul_ground(coeff)
- else:
- poly[mon] += (gen**monom[i]).mul_ground(coeff)
- return poly
- def to_dense(self):
- return dmp_from_dict(self, self.ring.ngens-1, self.ring.domain)
- def to_dict(self):
- return dict(self)
- def str(self, printer, precedence, exp_pattern, mul_symbol):
- if not self:
- return printer._print(self.ring.domain.zero)
- prec_mul = precedence["Mul"]
- prec_atom = precedence["Atom"]
- ring = self.ring
- symbols = ring.symbols
- ngens = ring.ngens
- zm = ring.zero_monom
- sexpvs = []
- for expv, coeff in self.terms():
- negative = ring.domain.is_negative(coeff)
- sign = " - " if negative else " + "
- sexpvs.append(sign)
- if expv == zm:
- scoeff = printer._print(coeff)
- if negative and scoeff.startswith("-"):
- scoeff = scoeff[1:]
- else:
- if negative:
- coeff = -coeff
- if coeff != self.ring.one:
- scoeff = printer.parenthesize(coeff, prec_mul, strict=True)
- else:
- scoeff = ''
- sexpv = []
- for i in range(ngens):
- exp = expv[i]
- if not exp:
- continue
- symbol = printer.parenthesize(symbols[i], prec_atom, strict=True)
- if exp != 1:
- if exp != int(exp) or exp < 0:
- sexp = printer.parenthesize(exp, prec_atom, strict=False)
- else:
- sexp = exp
- sexpv.append(exp_pattern % (symbol, sexp))
- else:
- sexpv.append('%s' % symbol)
- if scoeff:
- sexpv = [scoeff] + sexpv
- sexpvs.append(mul_symbol.join(sexpv))
- if sexpvs[0] in [" + ", " - "]:
- head = sexpvs.pop(0)
- if head == " - ":
- sexpvs.insert(0, "-")
- return "".join(sexpvs)
- @property
- def is_generator(self):
- return self in self.ring._gens_set
- @property
- def is_ground(self):
- return not self or (len(self) == 1 and self.ring.zero_monom in self)
- @property
- def is_monomial(self):
- return not self or (len(self) == 1 and self.LC == 1)
- @property
- def is_term(self):
- return len(self) <= 1
- @property
- def is_negative(self):
- return self.ring.domain.is_negative(self.LC)
- @property
- def is_positive(self):
- return self.ring.domain.is_positive(self.LC)
- @property
- def is_nonnegative(self):
- return self.ring.domain.is_nonnegative(self.LC)
- @property
- def is_nonpositive(self):
- return self.ring.domain.is_nonpositive(self.LC)
- @property
- def is_zero(f):
- return not f
- @property
- def is_one(f):
- return f == f.ring.one
- @property
- def is_monic(f):
- return f.ring.domain.is_one(f.LC)
- @property
- def is_primitive(f):
- return f.ring.domain.is_one(f.content())
- @property
- def is_linear(f):
- return all(sum(monom) <= 1 for monom in f.itermonoms())
- @property
- def is_quadratic(f):
- return all(sum(monom) <= 2 for monom in f.itermonoms())
- @property
- def is_squarefree(f):
- if not f.ring.ngens:
- return True
- return f.ring.dmp_sqf_p(f)
- @property
- def is_irreducible(f):
- if not f.ring.ngens:
- return True
- return f.ring.dmp_irreducible_p(f)
- @property
- def is_cyclotomic(f):
- if f.ring.is_univariate:
- return f.ring.dup_cyclotomic_p(f)
- else:
- raise MultivariatePolynomialError("cyclotomic polynomial")
- def __neg__(self):
- return self.new([ (monom, -coeff) for monom, coeff in self.iterterms() ])
- def __pos__(self):
- return self
- def __add__(p1, p2):
- """Add two polynomials.
- Examples
- ========
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.rings import ring
- >>> _, x, y = ring('x, y', ZZ)
- >>> (x + y)**2 + (x - y)**2
- 2*x**2 + 2*y**2
- """
- if not p2:
- return p1.copy()
- ring = p1.ring
- if isinstance(p2, ring.dtype):
- p = p1.copy()
- get = p.get
- zero = ring.domain.zero
- for k, v in p2.items():
- v = get(k, zero) + v
- if v:
- p[k] = v
- else:
- del p[k]
- return p
- elif isinstance(p2, PolyElement):
- if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
- pass
- elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
- return p2.__radd__(p1)
- else:
- return NotImplemented
- try:
- cp2 = ring.domain_new(p2)
- except CoercionFailed:
- return NotImplemented
- else:
- p = p1.copy()
- if not cp2:
- return p
- zm = ring.zero_monom
- if zm not in p1.keys():
- p[zm] = cp2
- else:
- if p2 == -p[zm]:
- del p[zm]
- else:
- p[zm] += cp2
- return p
- def __radd__(p1, n):
- p = p1.copy()
- if not n:
- return p
- ring = p1.ring
- try:
- n = ring.domain_new(n)
- except CoercionFailed:
- return NotImplemented
- else:
- zm = ring.zero_monom
- if zm not in p1.keys():
- p[zm] = n
- else:
- if n == -p[zm]:
- del p[zm]
- else:
- p[zm] += n
- return p
- def __sub__(p1, p2):
- """Subtract polynomial p2 from p1.
- Examples
- ========
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.rings import ring
- >>> _, x, y = ring('x, y', ZZ)
- >>> p1 = x + y**2
- >>> p2 = x*y + y**2
- >>> p1 - p2
- -x*y + x
- """
- if not p2:
- return p1.copy()
- ring = p1.ring
- if isinstance(p2, ring.dtype):
- p = p1.copy()
- get = p.get
- zero = ring.domain.zero
- for k, v in p2.items():
- v = get(k, zero) - v
- if v:
- p[k] = v
- else:
- del p[k]
- return p
- elif isinstance(p2, PolyElement):
- if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
- pass
- elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
- return p2.__rsub__(p1)
- else:
- return NotImplemented
- try:
- p2 = ring.domain_new(p2)
- except CoercionFailed:
- return NotImplemented
- else:
- p = p1.copy()
- zm = ring.zero_monom
- if zm not in p1.keys():
- p[zm] = -p2
- else:
- if p2 == p[zm]:
- del p[zm]
- else:
- p[zm] -= p2
- return p
- def __rsub__(p1, n):
- """n - p1 with n convertible to the coefficient domain.
- Examples
- ========
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.rings import ring
- >>> _, x, y = ring('x, y', ZZ)
- >>> p = x + y
- >>> 4 - p
- -x - y + 4
- """
- ring = p1.ring
- try:
- n = ring.domain_new(n)
- except CoercionFailed:
- return NotImplemented
- else:
- p = ring.zero
- for expv in p1:
- p[expv] = -p1[expv]
- p += n
- return p
- def __mul__(p1, p2):
- """Multiply two polynomials.
- Examples
- ========
- >>> from sympy.polys.domains import QQ
- >>> from sympy.polys.rings import ring
- >>> _, x, y = ring('x, y', QQ)
- >>> p1 = x + y
- >>> p2 = x - y
- >>> p1*p2
- x**2 - y**2
- """
- ring = p1.ring
- p = ring.zero
- if not p1 or not p2:
- return p
- elif isinstance(p2, ring.dtype):
- get = p.get
- zero = ring.domain.zero
- monomial_mul = ring.monomial_mul
- p2it = list(p2.items())
- for exp1, v1 in p1.items():
- for exp2, v2 in p2it:
- exp = monomial_mul(exp1, exp2)
- p[exp] = get(exp, zero) + v1*v2
- p.strip_zero()
- return p
- elif isinstance(p2, PolyElement):
- if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
- pass
- elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
- return p2.__rmul__(p1)
- else:
- return NotImplemented
- try:
- p2 = ring.domain_new(p2)
- except CoercionFailed:
- return NotImplemented
- else:
- for exp1, v1 in p1.items():
- v = v1*p2
- if v:
- p[exp1] = v
- return p
- def __rmul__(p1, p2):
- """p2 * p1 with p2 in the coefficient domain of p1.
- Examples
- ========
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.rings import ring
- >>> _, x, y = ring('x, y', ZZ)
- >>> p = x + y
- >>> 4 * p
- 4*x + 4*y
- """
- p = p1.ring.zero
- if not p2:
- return p
- try:
- p2 = p.ring.domain_new(p2)
- except CoercionFailed:
- return NotImplemented
- else:
- for exp1, v1 in p1.items():
- v = p2*v1
- if v:
- p[exp1] = v
- return p
- def __pow__(self, n):
- """raise polynomial to power `n`
- Examples
- ========
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.rings import ring
- >>> _, x, y = ring('x, y', ZZ)
- >>> p = x + y**2
- >>> p**3
- x**3 + 3*x**2*y**2 + 3*x*y**4 + y**6
- """
- ring = self.ring
- if not n:
- if self:
- return ring.one
- else:
- raise ValueError("0**0")
- elif len(self) == 1:
- monom, coeff = list(self.items())[0]
- p = ring.zero
- if coeff == 1:
- p[ring.monomial_pow(monom, n)] = coeff
- else:
- p[ring.monomial_pow(monom, n)] = coeff**n
- return p
- # For ring series, we need negative and rational exponent support only
- # with monomials.
- n = int(n)
- if n < 0:
- raise ValueError("Negative exponent")
- elif n == 1:
- return self.copy()
- elif n == 2:
- return self.square()
- elif n == 3:
- return self*self.square()
- elif len(self) <= 5: # TODO: use an actual density measure
- return self._pow_multinomial(n)
- else:
- return self._pow_generic(n)
- def _pow_generic(self, n):
- p = self.ring.one
- c = self
- while True:
- if n & 1:
- p = p*c
- n -= 1
- if not n:
- break
- c = c.square()
- n = n // 2
- return p
- def _pow_multinomial(self, n):
- multinomials = multinomial_coefficients(len(self), n).items()
- monomial_mulpow = self.ring.monomial_mulpow
- zero_monom = self.ring.zero_monom
- terms = self.items()
- zero = self.ring.domain.zero
- poly = self.ring.zero
- for multinomial, multinomial_coeff in multinomials:
- product_monom = zero_monom
- product_coeff = multinomial_coeff
- for exp, (monom, coeff) in zip(multinomial, terms):
- if exp:
- product_monom = monomial_mulpow(product_monom, monom, exp)
- product_coeff *= coeff**exp
- monom = tuple(product_monom)
- coeff = product_coeff
- coeff = poly.get(monom, zero) + coeff
- if coeff:
- poly[monom] = coeff
- elif monom in poly:
- del poly[monom]
- return poly
- def square(self):
- """square of a polynomial
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y = ring('x, y', ZZ)
- >>> p = x + y**2
- >>> p.square()
- x**2 + 2*x*y**2 + y**4
- """
- ring = self.ring
- p = ring.zero
- get = p.get
- keys = list(self.keys())
- zero = ring.domain.zero
- monomial_mul = ring.monomial_mul
- for i in range(len(keys)):
- k1 = keys[i]
- pk = self[k1]
- for j in range(i):
- k2 = keys[j]
- exp = monomial_mul(k1, k2)
- p[exp] = get(exp, zero) + pk*self[k2]
- p = p.imul_num(2)
- get = p.get
- for k, v in self.items():
- k2 = monomial_mul(k, k)
- p[k2] = get(k2, zero) + v**2
- p.strip_zero()
- return p
- def __divmod__(p1, p2):
- ring = p1.ring
- if not p2:
- raise ZeroDivisionError("polynomial division")
- elif isinstance(p2, ring.dtype):
- return p1.div(p2)
- elif isinstance(p2, PolyElement):
- if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
- pass
- elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
- return p2.__rdivmod__(p1)
- else:
- return NotImplemented
- try:
- p2 = ring.domain_new(p2)
- except CoercionFailed:
- return NotImplemented
- else:
- return (p1.quo_ground(p2), p1.rem_ground(p2))
- def __rdivmod__(p1, p2):
- return NotImplemented
- def __mod__(p1, p2):
- ring = p1.ring
- if not p2:
- raise ZeroDivisionError("polynomial division")
- elif isinstance(p2, ring.dtype):
- return p1.rem(p2)
- elif isinstance(p2, PolyElement):
- if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
- pass
- elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
- return p2.__rmod__(p1)
- else:
- return NotImplemented
- try:
- p2 = ring.domain_new(p2)
- except CoercionFailed:
- return NotImplemented
- else:
- return p1.rem_ground(p2)
- def __rmod__(p1, p2):
- return NotImplemented
- def __truediv__(p1, p2):
- ring = p1.ring
- if not p2:
- raise ZeroDivisionError("polynomial division")
- elif isinstance(p2, ring.dtype):
- if p2.is_monomial:
- return p1*(p2**(-1))
- else:
- return p1.quo(p2)
- elif isinstance(p2, PolyElement):
- if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
- pass
- elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
- return p2.__rtruediv__(p1)
- else:
- return NotImplemented
- try:
- p2 = ring.domain_new(p2)
- except CoercionFailed:
- return NotImplemented
- else:
- return p1.quo_ground(p2)
- def __rtruediv__(p1, p2):
- return NotImplemented
- __floordiv__ = __truediv__
- __rfloordiv__ = __rtruediv__
- # TODO: use // (__floordiv__) for exquo()?
- def _term_div(self):
- zm = self.ring.zero_monom
- domain = self.ring.domain
- domain_quo = domain.quo
- monomial_div = self.ring.monomial_div
- if domain.is_Field:
- def term_div(a_lm_a_lc, b_lm_b_lc):
- a_lm, a_lc = a_lm_a_lc
- b_lm, b_lc = b_lm_b_lc
- if b_lm == zm: # apparently this is a very common case
- monom = a_lm
- else:
- monom = monomial_div(a_lm, b_lm)
- if monom is not None:
- return monom, domain_quo(a_lc, b_lc)
- else:
- return None
- else:
- def term_div(a_lm_a_lc, b_lm_b_lc):
- a_lm, a_lc = a_lm_a_lc
- b_lm, b_lc = b_lm_b_lc
- if b_lm == zm: # apparently this is a very common case
- monom = a_lm
- else:
- monom = monomial_div(a_lm, b_lm)
- if not (monom is None or a_lc % b_lc):
- return monom, domain_quo(a_lc, b_lc)
- else:
- return None
- return term_div
- def div(self, fv):
- """Division algorithm, see [CLO] p64.
- fv array of polynomials
- return qv, r such that
- self = sum(fv[i]*qv[i]) + r
- All polynomials are required not to be Laurent polynomials.
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y = ring('x, y', ZZ)
- >>> f = x**3
- >>> f0 = x - y**2
- >>> f1 = x - y
- >>> qv, r = f.div((f0, f1))
- >>> qv[0]
- x**2 + x*y**2 + y**4
- >>> qv[1]
- 0
- >>> r
- y**6
- """
- ring = self.ring
- ret_single = False
- if isinstance(fv, PolyElement):
- ret_single = True
- fv = [fv]
- if not all(fv):
- raise ZeroDivisionError("polynomial division")
- if not self:
- if ret_single:
- return ring.zero, ring.zero
- else:
- return [], ring.zero
- for f in fv:
- if f.ring != ring:
- raise ValueError('self and f must have the same ring')
- s = len(fv)
- qv = [ring.zero for i in range(s)]
- p = self.copy()
- r = ring.zero
- term_div = self._term_div()
- expvs = [fx.leading_expv() for fx in fv]
- while p:
- i = 0
- divoccurred = 0
- while i < s and divoccurred == 0:
- expv = p.leading_expv()
- term = term_div((expv, p[expv]), (expvs[i], fv[i][expvs[i]]))
- if term is not None:
- expv1, c = term
- qv[i] = qv[i]._iadd_monom((expv1, c))
- p = p._iadd_poly_monom(fv[i], (expv1, -c))
- divoccurred = 1
- else:
- i += 1
- if not divoccurred:
- expv = p.leading_expv()
- r = r._iadd_monom((expv, p[expv]))
- del p[expv]
- if expv == ring.zero_monom:
- r += p
- if ret_single:
- if not qv:
- return ring.zero, r
- else:
- return qv[0], r
- else:
- return qv, r
- def rem(self, G):
- f = self
- if isinstance(G, PolyElement):
- G = [G]
- if not all(G):
- raise ZeroDivisionError("polynomial division")
- ring = f.ring
- domain = ring.domain
- zero = domain.zero
- monomial_mul = ring.monomial_mul
- r = ring.zero
- term_div = f._term_div()
- ltf = f.LT
- f = f.copy()
- get = f.get
- while f:
- for g in G:
- tq = term_div(ltf, g.LT)
- if tq is not None:
- m, c = tq
- for mg, cg in g.iterterms():
- m1 = monomial_mul(mg, m)
- c1 = get(m1, zero) - c*cg
- if not c1:
- del f[m1]
- else:
- f[m1] = c1
- ltm = f.leading_expv()
- if ltm is not None:
- ltf = ltm, f[ltm]
- break
- else:
- ltm, ltc = ltf
- if ltm in r:
- r[ltm] += ltc
- else:
- r[ltm] = ltc
- del f[ltm]
- ltm = f.leading_expv()
- if ltm is not None:
- ltf = ltm, f[ltm]
- return r
- def quo(f, G):
- return f.div(G)[0]
- def exquo(f, G):
- q, r = f.div(G)
- if not r:
- return q
- else:
- raise ExactQuotientFailed(f, G)
- def _iadd_monom(self, mc):
- """add to self the monomial coeff*x0**i0*x1**i1*...
- unless self is a generator -- then just return the sum of the two.
- mc is a tuple, (monom, coeff), where monomial is (i0, i1, ...)
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y = ring('x, y', ZZ)
- >>> p = x**4 + 2*y
- >>> m = (1, 2)
- >>> p1 = p._iadd_monom((m, 5))
- >>> p1
- x**4 + 5*x*y**2 + 2*y
- >>> p1 is p
- True
- >>> p = x
- >>> p1 = p._iadd_monom((m, 5))
- >>> p1
- 5*x*y**2 + x
- >>> p1 is p
- False
- """
- if self in self.ring._gens_set:
- cpself = self.copy()
- else:
- cpself = self
- expv, coeff = mc
- c = cpself.get(expv)
- if c is None:
- cpself[expv] = coeff
- else:
- c += coeff
- if c:
- cpself[expv] = c
- else:
- del cpself[expv]
- return cpself
- def _iadd_poly_monom(self, p2, mc):
- """add to self the product of (p)*(coeff*x0**i0*x1**i1*...)
- unless self is a generator -- then just return the sum of the two.
- mc is a tuple, (monom, coeff), where monomial is (i0, i1, ...)
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y, z = ring('x, y, z', ZZ)
- >>> p1 = x**4 + 2*y
- >>> p2 = y + z
- >>> m = (1, 2, 3)
- >>> p1 = p1._iadd_poly_monom(p2, (m, 3))
- >>> p1
- x**4 + 3*x*y**3*z**3 + 3*x*y**2*z**4 + 2*y
- """
- p1 = self
- if p1 in p1.ring._gens_set:
- p1 = p1.copy()
- (m, c) = mc
- get = p1.get
- zero = p1.ring.domain.zero
- monomial_mul = p1.ring.monomial_mul
- for k, v in p2.items():
- ka = monomial_mul(k, m)
- coeff = get(ka, zero) + v*c
- if coeff:
- p1[ka] = coeff
- else:
- del p1[ka]
- return p1
- def degree(f, x=None):
- """
- The leading degree in ``x`` or the main variable.
- Note that the degree of 0 is negative infinity (the SymPy object -oo).
- """
- i = f.ring.index(x)
- if not f:
- return -oo
- elif i < 0:
- return 0
- else:
- return max([ monom[i] for monom in f.itermonoms() ])
- def degrees(f):
- """
- A tuple containing leading degrees in all variables.
- Note that the degree of 0 is negative infinity (the SymPy object -oo)
- """
- if not f:
- return (-oo,)*f.ring.ngens
- else:
- return tuple(map(max, list(zip(*f.itermonoms()))))
- def tail_degree(f, x=None):
- """
- The tail degree in ``x`` or the main variable.
- Note that the degree of 0 is negative infinity (the SymPy object -oo)
- """
- i = f.ring.index(x)
- if not f:
- return -oo
- elif i < 0:
- return 0
- else:
- return min([ monom[i] for monom in f.itermonoms() ])
- def tail_degrees(f):
- """
- A tuple containing tail degrees in all variables.
- Note that the degree of 0 is negative infinity (the SymPy object -oo)
- """
- if not f:
- return (-oo,)*f.ring.ngens
- else:
- return tuple(map(min, list(zip(*f.itermonoms()))))
- def leading_expv(self):
- """Leading monomial tuple according to the monomial ordering.
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y, z = ring('x, y, z', ZZ)
- >>> p = x**4 + x**3*y + x**2*z**2 + z**7
- >>> p.leading_expv()
- (4, 0, 0)
- """
- if self:
- return self.ring.leading_expv(self)
- else:
- return None
- def _get_coeff(self, expv):
- return self.get(expv, self.ring.domain.zero)
- def coeff(self, element):
- """
- Returns the coefficient that stands next to the given monomial.
- Parameters
- ==========
- element : PolyElement (with ``is_monomial = True``) or 1
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y, z = ring("x,y,z", ZZ)
- >>> f = 3*x**2*y - x*y*z + 7*z**3 + 23
- >>> f.coeff(x**2*y)
- 3
- >>> f.coeff(x*y)
- 0
- >>> f.coeff(1)
- 23
- """
- if element == 1:
- return self._get_coeff(self.ring.zero_monom)
- elif isinstance(element, self.ring.dtype):
- terms = list(element.iterterms())
- if len(terms) == 1:
- monom, coeff = terms[0]
- if coeff == self.ring.domain.one:
- return self._get_coeff(monom)
- raise ValueError("expected a monomial, got %s" % element)
- def const(self):
- """Returns the constant coeffcient. """
- return self._get_coeff(self.ring.zero_monom)
- @property
- def LC(self):
- return self._get_coeff(self.leading_expv())
- @property
- def LM(self):
- expv = self.leading_expv()
- if expv is None:
- return self.ring.zero_monom
- else:
- return expv
- def leading_monom(self):
- """
- Leading monomial as a polynomial element.
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y = ring('x, y', ZZ)
- >>> (3*x*y + y**2).leading_monom()
- x*y
- """
- p = self.ring.zero
- expv = self.leading_expv()
- if expv:
- p[expv] = self.ring.domain.one
- return p
- @property
- def LT(self):
- expv = self.leading_expv()
- if expv is None:
- return (self.ring.zero_monom, self.ring.domain.zero)
- else:
- return (expv, self._get_coeff(expv))
- def leading_term(self):
- """Leading term as a polynomial element.
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y = ring('x, y', ZZ)
- >>> (3*x*y + y**2).leading_term()
- 3*x*y
- """
- p = self.ring.zero
- expv = self.leading_expv()
- if expv is not None:
- p[expv] = self[expv]
- return p
- def _sorted(self, seq, order):
- if order is None:
- order = self.ring.order
- else:
- order = OrderOpt.preprocess(order)
- if order is lex:
- return sorted(seq, key=lambda monom: monom[0], reverse=True)
- else:
- return sorted(seq, key=lambda monom: order(monom[0]), reverse=True)
- def coeffs(self, order=None):
- """Ordered list of polynomial coefficients.
- Parameters
- ==========
- order : :class:`~.MonomialOrder` or coercible, optional
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.orderings import lex, grlex
- >>> _, x, y = ring("x, y", ZZ, lex)
- >>> f = x*y**7 + 2*x**2*y**3
- >>> f.coeffs()
- [2, 1]
- >>> f.coeffs(grlex)
- [1, 2]
- """
- return [ coeff for _, coeff in self.terms(order) ]
- def monoms(self, order=None):
- """Ordered list of polynomial monomials.
- Parameters
- ==========
- order : :class:`~.MonomialOrder` or coercible, optional
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.orderings import lex, grlex
- >>> _, x, y = ring("x, y", ZZ, lex)
- >>> f = x*y**7 + 2*x**2*y**3
- >>> f.monoms()
- [(2, 3), (1, 7)]
- >>> f.monoms(grlex)
- [(1, 7), (2, 3)]
- """
- return [ monom for monom, _ in self.terms(order) ]
- def terms(self, order=None):
- """Ordered list of polynomial terms.
- Parameters
- ==========
- order : :class:`~.MonomialOrder` or coercible, optional
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> from sympy.polys.orderings import lex, grlex
- >>> _, x, y = ring("x, y", ZZ, lex)
- >>> f = x*y**7 + 2*x**2*y**3
- >>> f.terms()
- [((2, 3), 2), ((1, 7), 1)]
- >>> f.terms(grlex)
- [((1, 7), 1), ((2, 3), 2)]
- """
- return self._sorted(list(self.items()), order)
- def itercoeffs(self):
- """Iterator over coefficients of a polynomial. """
- return iter(self.values())
- def itermonoms(self):
- """Iterator over monomials of a polynomial. """
- return iter(self.keys())
- def iterterms(self):
- """Iterator over terms of a polynomial. """
- return iter(self.items())
- def listcoeffs(self):
- """Unordered list of polynomial coefficients. """
- return list(self.values())
- def listmonoms(self):
- """Unordered list of polynomial monomials. """
- return list(self.keys())
- def listterms(self):
- """Unordered list of polynomial terms. """
- return list(self.items())
- def imul_num(p, c):
- """multiply inplace the polynomial p by an element in the
- coefficient ring, provided p is not one of the generators;
- else multiply not inplace
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y = ring('x, y', ZZ)
- >>> p = x + y**2
- >>> p1 = p.imul_num(3)
- >>> p1
- 3*x + 3*y**2
- >>> p1 is p
- True
- >>> p = x
- >>> p1 = p.imul_num(3)
- >>> p1
- 3*x
- >>> p1 is p
- False
- """
- if p in p.ring._gens_set:
- return p*c
- if not c:
- p.clear()
- return
- for exp in p:
- p[exp] *= c
- return p
- def content(f):
- """Returns GCD of polynomial's coefficients. """
- domain = f.ring.domain
- cont = domain.zero
- gcd = domain.gcd
- for coeff in f.itercoeffs():
- cont = gcd(cont, coeff)
- return cont
- def primitive(f):
- """Returns content and a primitive polynomial. """
- cont = f.content()
- return cont, f.quo_ground(cont)
- def monic(f):
- """Divides all coefficients by the leading coefficient. """
- if not f:
- return f
- else:
- return f.quo_ground(f.LC)
- def mul_ground(f, x):
- if not x:
- return f.ring.zero
- terms = [ (monom, coeff*x) for monom, coeff in f.iterterms() ]
- return f.new(terms)
- def mul_monom(f, monom):
- monomial_mul = f.ring.monomial_mul
- terms = [ (monomial_mul(f_monom, monom), f_coeff) for f_monom, f_coeff in f.items() ]
- return f.new(terms)
- def mul_term(f, term):
- monom, coeff = term
- if not f or not coeff:
- return f.ring.zero
- elif monom == f.ring.zero_monom:
- return f.mul_ground(coeff)
- monomial_mul = f.ring.monomial_mul
- terms = [ (monomial_mul(f_monom, monom), f_coeff*coeff) for f_monom, f_coeff in f.items() ]
- return f.new(terms)
- def quo_ground(f, x):
- domain = f.ring.domain
- if not x:
- raise ZeroDivisionError('polynomial division')
- if not f or x == domain.one:
- return f
- if domain.is_Field:
- quo = domain.quo
- terms = [ (monom, quo(coeff, x)) for monom, coeff in f.iterterms() ]
- else:
- terms = [ (monom, coeff // x) for monom, coeff in f.iterterms() if not (coeff % x) ]
- return f.new(terms)
- def quo_term(f, term):
- monom, coeff = term
- if not coeff:
- raise ZeroDivisionError("polynomial division")
- elif not f:
- return f.ring.zero
- elif monom == f.ring.zero_monom:
- return f.quo_ground(coeff)
- term_div = f._term_div()
- terms = [ term_div(t, term) for t in f.iterterms() ]
- return f.new([ t for t in terms if t is not None ])
- def trunc_ground(f, p):
- if f.ring.domain.is_ZZ:
- terms = []
- for monom, coeff in f.iterterms():
- coeff = coeff % p
- if coeff > p // 2:
- coeff = coeff - p
- terms.append((monom, coeff))
- else:
- terms = [ (monom, coeff % p) for monom, coeff in f.iterterms() ]
- poly = f.new(terms)
- poly.strip_zero()
- return poly
- rem_ground = trunc_ground
- def extract_ground(self, g):
- f = self
- fc = f.content()
- gc = g.content()
- gcd = f.ring.domain.gcd(fc, gc)
- f = f.quo_ground(gcd)
- g = g.quo_ground(gcd)
- return gcd, f, g
- def _norm(f, norm_func):
- if not f:
- return f.ring.domain.zero
- else:
- ground_abs = f.ring.domain.abs
- return norm_func([ ground_abs(coeff) for coeff in f.itercoeffs() ])
- def max_norm(f):
- return f._norm(max)
- def l1_norm(f):
- return f._norm(sum)
- def deflate(f, *G):
- ring = f.ring
- polys = [f] + list(G)
- J = [0]*ring.ngens
- for p in polys:
- for monom in p.itermonoms():
- for i, m in enumerate(monom):
- J[i] = igcd(J[i], m)
- for i, b in enumerate(J):
- if not b:
- J[i] = 1
- J = tuple(J)
- if all(b == 1 for b in J):
- return J, polys
- H = []
- for p in polys:
- h = ring.zero
- for I, coeff in p.iterterms():
- N = [ i // j for i, j in zip(I, J) ]
- h[tuple(N)] = coeff
- H.append(h)
- return J, H
- def inflate(f, J):
- poly = f.ring.zero
- for I, coeff in f.iterterms():
- N = [ i*j for i, j in zip(I, J) ]
- poly[tuple(N)] = coeff
- return poly
- def lcm(self, g):
- f = self
- domain = f.ring.domain
- if not domain.is_Field:
- fc, f = f.primitive()
- gc, g = g.primitive()
- c = domain.lcm(fc, gc)
- h = (f*g).quo(f.gcd(g))
- if not domain.is_Field:
- return h.mul_ground(c)
- else:
- return h.monic()
- def gcd(f, g):
- return f.cofactors(g)[0]
- def cofactors(f, g):
- if not f and not g:
- zero = f.ring.zero
- return zero, zero, zero
- elif not f:
- h, cff, cfg = f._gcd_zero(g)
- return h, cff, cfg
- elif not g:
- h, cfg, cff = g._gcd_zero(f)
- return h, cff, cfg
- elif len(f) == 1:
- h, cff, cfg = f._gcd_monom(g)
- return h, cff, cfg
- elif len(g) == 1:
- h, cfg, cff = g._gcd_monom(f)
- return h, cff, cfg
- J, (f, g) = f.deflate(g)
- h, cff, cfg = f._gcd(g)
- return (h.inflate(J), cff.inflate(J), cfg.inflate(J))
- def _gcd_zero(f, g):
- one, zero = f.ring.one, f.ring.zero
- if g.is_nonnegative:
- return g, zero, one
- else:
- return -g, zero, -one
- def _gcd_monom(f, g):
- ring = f.ring
- ground_gcd = ring.domain.gcd
- ground_quo = ring.domain.quo
- monomial_gcd = ring.monomial_gcd
- monomial_ldiv = ring.monomial_ldiv
- mf, cf = list(f.iterterms())[0]
- _mgcd, _cgcd = mf, cf
- for mg, cg in g.iterterms():
- _mgcd = monomial_gcd(_mgcd, mg)
- _cgcd = ground_gcd(_cgcd, cg)
- h = f.new([(_mgcd, _cgcd)])
- cff = f.new([(monomial_ldiv(mf, _mgcd), ground_quo(cf, _cgcd))])
- cfg = f.new([(monomial_ldiv(mg, _mgcd), ground_quo(cg, _cgcd)) for mg, cg in g.iterterms()])
- return h, cff, cfg
- def _gcd(f, g):
- ring = f.ring
- if ring.domain.is_QQ:
- return f._gcd_QQ(g)
- elif ring.domain.is_ZZ:
- return f._gcd_ZZ(g)
- else: # TODO: don't use dense representation (port PRS algorithms)
- return ring.dmp_inner_gcd(f, g)
- def _gcd_ZZ(f, g):
- return heugcd(f, g)
- def _gcd_QQ(self, g):
- f = self
- ring = f.ring
- new_ring = ring.clone(domain=ring.domain.get_ring())
- cf, f = f.clear_denoms()
- cg, g = g.clear_denoms()
- f = f.set_ring(new_ring)
- g = g.set_ring(new_ring)
- h, cff, cfg = f._gcd_ZZ(g)
- h = h.set_ring(ring)
- c, h = h.LC, h.monic()
- cff = cff.set_ring(ring).mul_ground(ring.domain.quo(c, cf))
- cfg = cfg.set_ring(ring).mul_ground(ring.domain.quo(c, cg))
- return h, cff, cfg
- def cancel(self, g):
- """
- Cancel common factors in a rational function ``f/g``.
- Examples
- ========
- >>> from sympy.polys import ring, ZZ
- >>> R, x,y = ring("x,y", ZZ)
- >>> (2*x**2 - 2).cancel(x**2 - 2*x + 1)
- (2*x + 2, x - 1)
- """
- f = self
- ring = f.ring
- if not f:
- return f, ring.one
- domain = ring.domain
- if not (domain.is_Field and domain.has_assoc_Ring):
- _, p, q = f.cofactors(g)
- else:
- new_ring = ring.clone(domain=domain.get_ring())
- cq, f = f.clear_denoms()
- cp, g = g.clear_denoms()
- f = f.set_ring(new_ring)
- g = g.set_ring(new_ring)
- _, p, q = f.cofactors(g)
- _, cp, cq = new_ring.domain.cofactors(cp, cq)
- p = p.set_ring(ring)
- q = q.set_ring(ring)
- p = p.mul_ground(cp)
- q = q.mul_ground(cq)
- # Make canonical with respect to sign or quadrant in the case of ZZ_I
- # or QQ_I. This ensures that the LC of the denominator is canonical by
- # multiplying top and bottom by a unit of the ring.
- u = q.canonical_unit()
- if u == domain.one:
- p, q = p, q
- elif u == -domain.one:
- p, q = -p, -q
- else:
- p = p.mul_ground(u)
- q = q.mul_ground(u)
- return p, q
- def canonical_unit(f):
- domain = f.ring.domain
- return domain.canonical_unit(f.LC)
- def diff(f, x):
- """Computes partial derivative in ``x``.
- Examples
- ========
- >>> from sympy.polys.rings import ring
- >>> from sympy.polys.domains import ZZ
- >>> _, x, y = ring("x,y", ZZ)
- >>> p = x + x**2*y**3
- >>> p.diff(x)
- 2*x*y**3 + 1
- """
- ring = f.ring
- i = ring.index(x)
- m = ring.monomial_basis(i)
- g = ring.zero
- for expv, coeff in f.iterterms():
- if expv[i]:
- e = ring.monomial_ldiv(expv, m)
- g[e] = ring.domain_new(coeff*expv[i])
- return g
- def __call__(f, *values):
- if 0 < len(values) <= f.ring.ngens:
- return f.evaluate(list(zip(f.ring.gens, values)))
- else:
- raise ValueError("expected at least 1 and at most %s values, got %s" % (f.ring.ngens, len(values)))
- def evaluate(self, x, a=None):
- f = self
- if isinstance(x, list) and a is None:
- (X, a), x = x[0], x[1:]
- f = f.evaluate(X, a)
- if not x:
- return f
- else:
- x = [ (Y.drop(X), a) for (Y, a) in x ]
- return f.evaluate(x)
- ring = f.ring
- i = ring.index(x)
- a = ring.domain.convert(a)
- if ring.ngens == 1:
- result = ring.domain.zero
- for (n,), coeff in f.iterterms():
- result += coeff*a**n
- return result
- else:
- poly = ring.drop(x).zero
- for monom, coeff in f.iterterms():
- n, monom = monom[i], monom[:i] + monom[i+1:]
- coeff = coeff*a**n
- if monom in poly:
- coeff = coeff + poly[monom]
- if coeff:
- poly[monom] = coeff
- else:
- del poly[monom]
- else:
- if coeff:
- poly[monom] = coeff
- return poly
- def subs(self, x, a=None):
- f = self
- if isinstance(x, list) and a is None:
- for X, a in x:
- f = f.subs(X, a)
- return f
- ring = f.ring
- i = ring.index(x)
- a = ring.domain.convert(a)
- if ring.ngens == 1:
- result = ring.domain.zero
- for (n,), coeff in f.iterterms():
- result += coeff*a**n
- return ring.ground_new(result)
- else:
- poly = ring.zero
- for monom, coeff in f.iterterms():
- n, monom = monom[i], monom[:i] + (0,) + monom[i+1:]
- coeff = coeff*a**n
- if monom in poly:
- coeff = coeff + poly[monom]
- if coeff:
- poly[monom] = coeff
- else:
- del poly[monom]
- else:
- if coeff:
- poly[monom] = coeff
- return poly
- def compose(f, x, a=None):
- ring = f.ring
- poly = ring.zero
- gens_map = dict(list(zip(ring.gens, list(range(ring.ngens)))))
- if a is not None:
- replacements = [(x, a)]
- else:
- if isinstance(x, list):
- replacements = list(x)
- elif isinstance(x, dict):
- replacements = sorted(list(x.items()), key=lambda k: gens_map[k[0]])
- else:
- raise ValueError("expected a generator, value pair a sequence of such pairs")
- for k, (x, g) in enumerate(replacements):
- replacements[k] = (gens_map[x], ring.ring_new(g))
- for monom, coeff in f.iterterms():
- monom = list(monom)
- subpoly = ring.one
- for i, g in replacements:
- n, monom[i] = monom[i], 0
- if n:
- subpoly *= g**n
- subpoly = subpoly.mul_term((tuple(monom), coeff))
- poly += subpoly
- return poly
- # TODO: following methods should point to polynomial
- # representation independent algorithm implementations.
- def pdiv(f, g):
- return f.ring.dmp_pdiv(f, g)
- def prem(f, g):
- return f.ring.dmp_prem(f, g)
- def pquo(f, g):
- return f.ring.dmp_quo(f, g)
- def pexquo(f, g):
- return f.ring.dmp_exquo(f, g)
- def half_gcdex(f, g):
- return f.ring.dmp_half_gcdex(f, g)
- def gcdex(f, g):
- return f.ring.dmp_gcdex(f, g)
- def subresultants(f, g):
- return f.ring.dmp_subresultants(f, g)
- def resultant(f, g):
- return f.ring.dmp_resultant(f, g)
- def discriminant(f):
- return f.ring.dmp_discriminant(f)
- def decompose(f):
- if f.ring.is_univariate:
- return f.ring.dup_decompose(f)
- else:
- raise MultivariatePolynomialError("polynomial decomposition")
- def shift(f, a):
- if f.ring.is_univariate:
- return f.ring.dup_shift(f, a)
- else:
- raise MultivariatePolynomialError("polynomial shift")
- def sturm(f):
- if f.ring.is_univariate:
- return f.ring.dup_sturm(f)
- else:
- raise MultivariatePolynomialError("sturm sequence")
- def gff_list(f):
- return f.ring.dmp_gff_list(f)
- def sqf_norm(f):
- return f.ring.dmp_sqf_norm(f)
- def sqf_part(f):
- return f.ring.dmp_sqf_part(f)
- def sqf_list(f, all=False):
- return f.ring.dmp_sqf_list(f, all=all)
- def factor_list(f):
- return f.ring.dmp_factor_list(f)
|