rings.py 67 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467
  1. """Sparse polynomial rings. """
  2. from typing import Any, Dict as tDict
  3. from operator import add, mul, lt, le, gt, ge
  4. from functools import reduce
  5. from types import GeneratorType
  6. from sympy.core.expr import Expr
  7. from sympy.core.numbers import igcd, oo
  8. from sympy.core.symbol import Symbol, symbols as _symbols
  9. from sympy.core.sympify import CantSympify, sympify
  10. from sympy.ntheory.multinomial import multinomial_coefficients
  11. from sympy.polys.compatibility import IPolys
  12. from sympy.polys.constructor import construct_domain
  13. from sympy.polys.densebasic import dmp_to_dict, dmp_from_dict
  14. from sympy.polys.domains.domainelement import DomainElement
  15. from sympy.polys.domains.polynomialring import PolynomialRing
  16. from sympy.polys.heuristicgcd import heugcd
  17. from sympy.polys.monomials import MonomialOps
  18. from sympy.polys.orderings import lex
  19. from sympy.polys.polyerrors import (
  20. CoercionFailed, GeneratorsError,
  21. ExactQuotientFailed, MultivariatePolynomialError)
  22. from sympy.polys.polyoptions import (Domain as DomainOpt,
  23. Order as OrderOpt, build_options)
  24. from sympy.polys.polyutils import (expr_from_dict, _dict_reorder,
  25. _parallel_dict_from_expr)
  26. from sympy.printing.defaults import DefaultPrinting
  27. from sympy.utilities import public
  28. from sympy.utilities.iterables import is_sequence
  29. from sympy.utilities.magic import pollute
  30. @public
  31. def ring(symbols, domain, order=lex):
  32. """Construct a polynomial ring returning ``(ring, x_1, ..., x_n)``.
  33. Parameters
  34. ==========
  35. symbols : str
  36. Symbol/Expr or sequence of str, Symbol/Expr (non-empty)
  37. domain : :class:`~.Domain` or coercible
  38. order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``
  39. Examples
  40. ========
  41. >>> from sympy.polys.rings import ring
  42. >>> from sympy.polys.domains import ZZ
  43. >>> from sympy.polys.orderings import lex
  44. >>> R, x, y, z = ring("x,y,z", ZZ, lex)
  45. >>> R
  46. Polynomial ring in x, y, z over ZZ with lex order
  47. >>> x + y + z
  48. x + y + z
  49. >>> type(_)
  50. <class 'sympy.polys.rings.PolyElement'>
  51. """
  52. _ring = PolyRing(symbols, domain, order)
  53. return (_ring,) + _ring.gens
  54. @public
  55. def xring(symbols, domain, order=lex):
  56. """Construct a polynomial ring returning ``(ring, (x_1, ..., x_n))``.
  57. Parameters
  58. ==========
  59. symbols : str
  60. Symbol/Expr or sequence of str, Symbol/Expr (non-empty)
  61. domain : :class:`~.Domain` or coercible
  62. order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``
  63. Examples
  64. ========
  65. >>> from sympy.polys.rings import xring
  66. >>> from sympy.polys.domains import ZZ
  67. >>> from sympy.polys.orderings import lex
  68. >>> R, (x, y, z) = xring("x,y,z", ZZ, lex)
  69. >>> R
  70. Polynomial ring in x, y, z over ZZ with lex order
  71. >>> x + y + z
  72. x + y + z
  73. >>> type(_)
  74. <class 'sympy.polys.rings.PolyElement'>
  75. """
  76. _ring = PolyRing(symbols, domain, order)
  77. return (_ring, _ring.gens)
  78. @public
  79. def vring(symbols, domain, order=lex):
  80. """Construct a polynomial ring and inject ``x_1, ..., x_n`` into the global namespace.
  81. Parameters
  82. ==========
  83. symbols : str
  84. Symbol/Expr or sequence of str, Symbol/Expr (non-empty)
  85. domain : :class:`~.Domain` or coercible
  86. order : :class:`~.MonomialOrder` or coercible, optional, defaults to ``lex``
  87. Examples
  88. ========
  89. >>> from sympy.polys.rings import vring
  90. >>> from sympy.polys.domains import ZZ
  91. >>> from sympy.polys.orderings import lex
  92. >>> vring("x,y,z", ZZ, lex)
  93. Polynomial ring in x, y, z over ZZ with lex order
  94. >>> x + y + z # noqa:
  95. x + y + z
  96. >>> type(_)
  97. <class 'sympy.polys.rings.PolyElement'>
  98. """
  99. _ring = PolyRing(symbols, domain, order)
  100. pollute([ sym.name for sym in _ring.symbols ], _ring.gens)
  101. return _ring
  102. @public
  103. def sring(exprs, *symbols, **options):
  104. """Construct a ring deriving generators and domain from options and input expressions.
  105. Parameters
  106. ==========
  107. exprs : :class:`~.Expr` or sequence of :class:`~.Expr` (sympifiable)
  108. symbols : sequence of :class:`~.Symbol`/:class:`~.Expr`
  109. options : keyword arguments understood by :class:`~.Options`
  110. Examples
  111. ========
  112. >>> from sympy import sring, symbols
  113. >>> x, y, z = symbols("x,y,z")
  114. >>> R, f = sring(x + 2*y + 3*z)
  115. >>> R
  116. Polynomial ring in x, y, z over ZZ with lex order
  117. >>> f
  118. x + 2*y + 3*z
  119. >>> type(_)
  120. <class 'sympy.polys.rings.PolyElement'>
  121. """
  122. single = False
  123. if not is_sequence(exprs):
  124. exprs, single = [exprs], True
  125. exprs = list(map(sympify, exprs))
  126. opt = build_options(symbols, options)
  127. # TODO: rewrite this so that it doesn't use expand() (see poly()).
  128. reps, opt = _parallel_dict_from_expr(exprs, opt)
  129. if opt.domain is None:
  130. coeffs = sum([ list(rep.values()) for rep in reps ], [])
  131. opt.domain, coeffs_dom = construct_domain(coeffs, opt=opt)
  132. coeff_map = dict(zip(coeffs, coeffs_dom))
  133. reps = [{m: coeff_map[c] for m, c in rep.items()} for rep in reps]
  134. _ring = PolyRing(opt.gens, opt.domain, opt.order)
  135. polys = list(map(_ring.from_dict, reps))
  136. if single:
  137. return (_ring, polys[0])
  138. else:
  139. return (_ring, polys)
  140. def _parse_symbols(symbols):
  141. if isinstance(symbols, str):
  142. return _symbols(symbols, seq=True) if symbols else ()
  143. elif isinstance(symbols, Expr):
  144. return (symbols,)
  145. elif is_sequence(symbols):
  146. if all(isinstance(s, str) for s in symbols):
  147. return _symbols(symbols)
  148. elif all(isinstance(s, Expr) for s in symbols):
  149. return symbols
  150. raise GeneratorsError("expected a string, Symbol or expression or a non-empty sequence of strings, Symbols or expressions")
  151. _ring_cache = {} # type: tDict[Any, Any]
  152. class PolyRing(DefaultPrinting, IPolys):
  153. """Multivariate distributed polynomial ring. """
  154. def __new__(cls, symbols, domain, order=lex):
  155. symbols = tuple(_parse_symbols(symbols))
  156. ngens = len(symbols)
  157. domain = DomainOpt.preprocess(domain)
  158. order = OrderOpt.preprocess(order)
  159. _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
  160. obj = _ring_cache.get(_hash_tuple)
  161. if obj is None:
  162. if domain.is_Composite and set(symbols) & set(domain.symbols):
  163. raise GeneratorsError("polynomial ring and it's ground domain share generators")
  164. obj = object.__new__(cls)
  165. obj._hash_tuple = _hash_tuple
  166. obj._hash = hash(_hash_tuple)
  167. obj.dtype = type("PolyElement", (PolyElement,), {"ring": obj})
  168. obj.symbols = symbols
  169. obj.ngens = ngens
  170. obj.domain = domain
  171. obj.order = order
  172. obj.zero_monom = (0,)*ngens
  173. obj.gens = obj._gens()
  174. obj._gens_set = set(obj.gens)
  175. obj._one = [(obj.zero_monom, domain.one)]
  176. if ngens:
  177. # These expect monomials in at least one variable
  178. codegen = MonomialOps(ngens)
  179. obj.monomial_mul = codegen.mul()
  180. obj.monomial_pow = codegen.pow()
  181. obj.monomial_mulpow = codegen.mulpow()
  182. obj.monomial_ldiv = codegen.ldiv()
  183. obj.monomial_div = codegen.div()
  184. obj.monomial_lcm = codegen.lcm()
  185. obj.monomial_gcd = codegen.gcd()
  186. else:
  187. monunit = lambda a, b: ()
  188. obj.monomial_mul = monunit
  189. obj.monomial_pow = monunit
  190. obj.monomial_mulpow = lambda a, b, c: ()
  191. obj.monomial_ldiv = monunit
  192. obj.monomial_div = monunit
  193. obj.monomial_lcm = monunit
  194. obj.monomial_gcd = monunit
  195. if order is lex:
  196. obj.leading_expv = max
  197. else:
  198. obj.leading_expv = lambda f: max(f, key=order)
  199. for symbol, generator in zip(obj.symbols, obj.gens):
  200. if isinstance(symbol, Symbol):
  201. name = symbol.name
  202. if not hasattr(obj, name):
  203. setattr(obj, name, generator)
  204. _ring_cache[_hash_tuple] = obj
  205. return obj
  206. def _gens(self):
  207. """Return a list of polynomial generators. """
  208. one = self.domain.one
  209. _gens = []
  210. for i in range(self.ngens):
  211. expv = self.monomial_basis(i)
  212. poly = self.zero
  213. poly[expv] = one
  214. _gens.append(poly)
  215. return tuple(_gens)
  216. def __getnewargs__(self):
  217. return (self.symbols, self.domain, self.order)
  218. def __getstate__(self):
  219. state = self.__dict__.copy()
  220. del state["leading_expv"]
  221. for key, value in state.items():
  222. if key.startswith("monomial_"):
  223. del state[key]
  224. return state
  225. def __hash__(self):
  226. return self._hash
  227. def __eq__(self, other):
  228. return isinstance(other, PolyRing) and \
  229. (self.symbols, self.domain, self.ngens, self.order) == \
  230. (other.symbols, other.domain, other.ngens, other.order)
  231. def __ne__(self, other):
  232. return not self == other
  233. def clone(self, symbols=None, domain=None, order=None):
  234. return self.__class__(symbols or self.symbols, domain or self.domain, order or self.order)
  235. def monomial_basis(self, i):
  236. """Return the ith-basis element. """
  237. basis = [0]*self.ngens
  238. basis[i] = 1
  239. return tuple(basis)
  240. @property
  241. def zero(self):
  242. return self.dtype()
  243. @property
  244. def one(self):
  245. return self.dtype(self._one)
  246. def domain_new(self, element, orig_domain=None):
  247. return self.domain.convert(element, orig_domain)
  248. def ground_new(self, coeff):
  249. return self.term_new(self.zero_monom, coeff)
  250. def term_new(self, monom, coeff):
  251. coeff = self.domain_new(coeff)
  252. poly = self.zero
  253. if coeff:
  254. poly[monom] = coeff
  255. return poly
  256. def ring_new(self, element):
  257. if isinstance(element, PolyElement):
  258. if self == element.ring:
  259. return element
  260. elif isinstance(self.domain, PolynomialRing) and self.domain.ring == element.ring:
  261. return self.ground_new(element)
  262. else:
  263. raise NotImplementedError("conversion")
  264. elif isinstance(element, str):
  265. raise NotImplementedError("parsing")
  266. elif isinstance(element, dict):
  267. return self.from_dict(element)
  268. elif isinstance(element, list):
  269. try:
  270. return self.from_terms(element)
  271. except ValueError:
  272. return self.from_list(element)
  273. elif isinstance(element, Expr):
  274. return self.from_expr(element)
  275. else:
  276. return self.ground_new(element)
  277. __call__ = ring_new
  278. def from_dict(self, element, orig_domain=None):
  279. domain_new = self.domain_new
  280. poly = self.zero
  281. for monom, coeff in element.items():
  282. coeff = domain_new(coeff, orig_domain)
  283. if coeff:
  284. poly[monom] = coeff
  285. return poly
  286. def from_terms(self, element, orig_domain=None):
  287. return self.from_dict(dict(element), orig_domain)
  288. def from_list(self, element):
  289. return self.from_dict(dmp_to_dict(element, self.ngens-1, self.domain))
  290. def _rebuild_expr(self, expr, mapping):
  291. domain = self.domain
  292. def _rebuild(expr):
  293. generator = mapping.get(expr)
  294. if generator is not None:
  295. return generator
  296. elif expr.is_Add:
  297. return reduce(add, list(map(_rebuild, expr.args)))
  298. elif expr.is_Mul:
  299. return reduce(mul, list(map(_rebuild, expr.args)))
  300. else:
  301. # XXX: Use as_base_exp() to handle Pow(x, n) and also exp(n)
  302. # XXX: E can be a generator e.g. sring([exp(2)]) -> ZZ[E]
  303. base, exp = expr.as_base_exp()
  304. if exp.is_Integer and exp > 1:
  305. return _rebuild(base)**int(exp)
  306. else:
  307. return self.ground_new(domain.convert(expr))
  308. return _rebuild(sympify(expr))
  309. def from_expr(self, expr):
  310. mapping = dict(list(zip(self.symbols, self.gens)))
  311. try:
  312. poly = self._rebuild_expr(expr, mapping)
  313. except CoercionFailed:
  314. raise ValueError("expected an expression convertible to a polynomial in %s, got %s" % (self, expr))
  315. else:
  316. return self.ring_new(poly)
  317. def index(self, gen):
  318. """Compute index of ``gen`` in ``self.gens``. """
  319. if gen is None:
  320. if self.ngens:
  321. i = 0
  322. else:
  323. i = -1 # indicate impossible choice
  324. elif isinstance(gen, int):
  325. i = gen
  326. if 0 <= i and i < self.ngens:
  327. pass
  328. elif -self.ngens <= i and i <= -1:
  329. i = -i - 1
  330. else:
  331. raise ValueError("invalid generator index: %s" % gen)
  332. elif isinstance(gen, self.dtype):
  333. try:
  334. i = self.gens.index(gen)
  335. except ValueError:
  336. raise ValueError("invalid generator: %s" % gen)
  337. elif isinstance(gen, str):
  338. try:
  339. i = self.symbols.index(gen)
  340. except ValueError:
  341. raise ValueError("invalid generator: %s" % gen)
  342. else:
  343. raise ValueError("expected a polynomial generator, an integer, a string or None, got %s" % gen)
  344. return i
  345. def drop(self, *gens):
  346. """Remove specified generators from this ring. """
  347. indices = set(map(self.index, gens))
  348. symbols = [ s for i, s in enumerate(self.symbols) if i not in indices ]
  349. if not symbols:
  350. return self.domain
  351. else:
  352. return self.clone(symbols=symbols)
  353. def __getitem__(self, key):
  354. symbols = self.symbols[key]
  355. if not symbols:
  356. return self.domain
  357. else:
  358. return self.clone(symbols=symbols)
  359. def to_ground(self):
  360. # TODO: should AlgebraicField be a Composite domain?
  361. if self.domain.is_Composite or hasattr(self.domain, 'domain'):
  362. return self.clone(domain=self.domain.domain)
  363. else:
  364. raise ValueError("%s is not a composite domain" % self.domain)
  365. def to_domain(self):
  366. return PolynomialRing(self)
  367. def to_field(self):
  368. from sympy.polys.fields import FracField
  369. return FracField(self.symbols, self.domain, self.order)
  370. @property
  371. def is_univariate(self):
  372. return len(self.gens) == 1
  373. @property
  374. def is_multivariate(self):
  375. return len(self.gens) > 1
  376. def add(self, *objs):
  377. """
  378. Add a sequence of polynomials or containers of polynomials.
  379. Examples
  380. ========
  381. >>> from sympy.polys.rings import ring
  382. >>> from sympy.polys.domains import ZZ
  383. >>> R, x = ring("x", ZZ)
  384. >>> R.add([ x**2 + 2*i + 3 for i in range(4) ])
  385. 4*x**2 + 24
  386. >>> _.factor_list()
  387. (4, [(x**2 + 6, 1)])
  388. """
  389. p = self.zero
  390. for obj in objs:
  391. if is_sequence(obj, include=GeneratorType):
  392. p += self.add(*obj)
  393. else:
  394. p += obj
  395. return p
  396. def mul(self, *objs):
  397. """
  398. Multiply a sequence of polynomials or containers of polynomials.
  399. Examples
  400. ========
  401. >>> from sympy.polys.rings import ring
  402. >>> from sympy.polys.domains import ZZ
  403. >>> R, x = ring("x", ZZ)
  404. >>> R.mul([ x**2 + 2*i + 3 for i in range(4) ])
  405. x**8 + 24*x**6 + 206*x**4 + 744*x**2 + 945
  406. >>> _.factor_list()
  407. (1, [(x**2 + 3, 1), (x**2 + 5, 1), (x**2 + 7, 1), (x**2 + 9, 1)])
  408. """
  409. p = self.one
  410. for obj in objs:
  411. if is_sequence(obj, include=GeneratorType):
  412. p *= self.mul(*obj)
  413. else:
  414. p *= obj
  415. return p
  416. def drop_to_ground(self, *gens):
  417. r"""
  418. Remove specified generators from the ring and inject them into
  419. its domain.
  420. """
  421. indices = set(map(self.index, gens))
  422. symbols = [s for i, s in enumerate(self.symbols) if i not in indices]
  423. gens = [gen for i, gen in enumerate(self.gens) if i not in indices]
  424. if not symbols:
  425. return self
  426. else:
  427. return self.clone(symbols=symbols, domain=self.drop(*gens))
  428. def compose(self, other):
  429. """Add the generators of ``other`` to ``self``"""
  430. if self != other:
  431. syms = set(self.symbols).union(set(other.symbols))
  432. return self.clone(symbols=list(syms))
  433. else:
  434. return self
  435. def add_gens(self, symbols):
  436. """Add the elements of ``symbols`` as generators to ``self``"""
  437. syms = set(self.symbols).union(set(symbols))
  438. return self.clone(symbols=list(syms))
  439. class PolyElement(DomainElement, DefaultPrinting, CantSympify, dict):
  440. """Element of multivariate distributed polynomial ring. """
  441. def new(self, init):
  442. return self.__class__(init)
  443. def parent(self):
  444. return self.ring.to_domain()
  445. def __getnewargs__(self):
  446. return (self.ring, list(self.iterterms()))
  447. _hash = None
  448. def __hash__(self):
  449. # XXX: This computes a hash of a dictionary, but currently we don't
  450. # protect dictionary from being changed so any use site modifications
  451. # will make hashing go wrong. Use this feature with caution until we
  452. # figure out how to make a safe API without compromising speed of this
  453. # low-level class.
  454. _hash = self._hash
  455. if _hash is None:
  456. self._hash = _hash = hash((self.ring, frozenset(self.items())))
  457. return _hash
  458. def copy(self):
  459. """Return a copy of polynomial self.
  460. Polynomials are mutable; if one is interested in preserving
  461. a polynomial, and one plans to use inplace operations, one
  462. can copy the polynomial. This method makes a shallow copy.
  463. Examples
  464. ========
  465. >>> from sympy.polys.domains import ZZ
  466. >>> from sympy.polys.rings import ring
  467. >>> R, x, y = ring('x, y', ZZ)
  468. >>> p = (x + y)**2
  469. >>> p1 = p.copy()
  470. >>> p2 = p
  471. >>> p[R.zero_monom] = 3
  472. >>> p
  473. x**2 + 2*x*y + y**2 + 3
  474. >>> p1
  475. x**2 + 2*x*y + y**2
  476. >>> p2
  477. x**2 + 2*x*y + y**2 + 3
  478. """
  479. return self.new(self)
  480. def set_ring(self, new_ring):
  481. if self.ring == new_ring:
  482. return self
  483. elif self.ring.symbols != new_ring.symbols:
  484. terms = list(zip(*_dict_reorder(self, self.ring.symbols, new_ring.symbols)))
  485. return new_ring.from_terms(terms, self.ring.domain)
  486. else:
  487. return new_ring.from_dict(self, self.ring.domain)
  488. def as_expr(self, *symbols):
  489. if symbols and len(symbols) != self.ring.ngens:
  490. raise ValueError("not enough symbols, expected %s got %s" % (self.ring.ngens, len(symbols)))
  491. else:
  492. symbols = self.ring.symbols
  493. return expr_from_dict(self.as_expr_dict(), *symbols)
  494. def as_expr_dict(self):
  495. to_sympy = self.ring.domain.to_sympy
  496. return {monom: to_sympy(coeff) for monom, coeff in self.iterterms()}
  497. def clear_denoms(self):
  498. domain = self.ring.domain
  499. if not domain.is_Field or not domain.has_assoc_Ring:
  500. return domain.one, self
  501. ground_ring = domain.get_ring()
  502. common = ground_ring.one
  503. lcm = ground_ring.lcm
  504. denom = domain.denom
  505. for coeff in self.values():
  506. common = lcm(common, denom(coeff))
  507. poly = self.new([ (k, v*common) for k, v in self.items() ])
  508. return common, poly
  509. def strip_zero(self):
  510. """Eliminate monomials with zero coefficient. """
  511. for k, v in list(self.items()):
  512. if not v:
  513. del self[k]
  514. def __eq__(p1, p2):
  515. """Equality test for polynomials.
  516. Examples
  517. ========
  518. >>> from sympy.polys.domains import ZZ
  519. >>> from sympy.polys.rings import ring
  520. >>> _, x, y = ring('x, y', ZZ)
  521. >>> p1 = (x + y)**2 + (x - y)**2
  522. >>> p1 == 4*x*y
  523. False
  524. >>> p1 == 2*(x**2 + y**2)
  525. True
  526. """
  527. if not p2:
  528. return not p1
  529. elif isinstance(p2, PolyElement) and p2.ring == p1.ring:
  530. return dict.__eq__(p1, p2)
  531. elif len(p1) > 1:
  532. return False
  533. else:
  534. return p1.get(p1.ring.zero_monom) == p2
  535. def __ne__(p1, p2):
  536. return not p1 == p2
  537. def almosteq(p1, p2, tolerance=None):
  538. """Approximate equality test for polynomials. """
  539. ring = p1.ring
  540. if isinstance(p2, ring.dtype):
  541. if set(p1.keys()) != set(p2.keys()):
  542. return False
  543. almosteq = ring.domain.almosteq
  544. for k in p1.keys():
  545. if not almosteq(p1[k], p2[k], tolerance):
  546. return False
  547. return True
  548. elif len(p1) > 1:
  549. return False
  550. else:
  551. try:
  552. p2 = ring.domain.convert(p2)
  553. except CoercionFailed:
  554. return False
  555. else:
  556. return ring.domain.almosteq(p1.const(), p2, tolerance)
  557. def sort_key(self):
  558. return (len(self), self.terms())
  559. def _cmp(p1, p2, op):
  560. if isinstance(p2, p1.ring.dtype):
  561. return op(p1.sort_key(), p2.sort_key())
  562. else:
  563. return NotImplemented
  564. def __lt__(p1, p2):
  565. return p1._cmp(p2, lt)
  566. def __le__(p1, p2):
  567. return p1._cmp(p2, le)
  568. def __gt__(p1, p2):
  569. return p1._cmp(p2, gt)
  570. def __ge__(p1, p2):
  571. return p1._cmp(p2, ge)
  572. def _drop(self, gen):
  573. ring = self.ring
  574. i = ring.index(gen)
  575. if ring.ngens == 1:
  576. return i, ring.domain
  577. else:
  578. symbols = list(ring.symbols)
  579. del symbols[i]
  580. return i, ring.clone(symbols=symbols)
  581. def drop(self, gen):
  582. i, ring = self._drop(gen)
  583. if self.ring.ngens == 1:
  584. if self.is_ground:
  585. return self.coeff(1)
  586. else:
  587. raise ValueError("Cannot drop %s" % gen)
  588. else:
  589. poly = ring.zero
  590. for k, v in self.items():
  591. if k[i] == 0:
  592. K = list(k)
  593. del K[i]
  594. poly[tuple(K)] = v
  595. else:
  596. raise ValueError("Cannot drop %s" % gen)
  597. return poly
  598. def _drop_to_ground(self, gen):
  599. ring = self.ring
  600. i = ring.index(gen)
  601. symbols = list(ring.symbols)
  602. del symbols[i]
  603. return i, ring.clone(symbols=symbols, domain=ring[i])
  604. def drop_to_ground(self, gen):
  605. if self.ring.ngens == 1:
  606. raise ValueError("Cannot drop only generator to ground")
  607. i, ring = self._drop_to_ground(gen)
  608. poly = ring.zero
  609. gen = ring.domain.gens[0]
  610. for monom, coeff in self.iterterms():
  611. mon = monom[:i] + monom[i+1:]
  612. if mon not in poly:
  613. poly[mon] = (gen**monom[i]).mul_ground(coeff)
  614. else:
  615. poly[mon] += (gen**monom[i]).mul_ground(coeff)
  616. return poly
  617. def to_dense(self):
  618. return dmp_from_dict(self, self.ring.ngens-1, self.ring.domain)
  619. def to_dict(self):
  620. return dict(self)
  621. def str(self, printer, precedence, exp_pattern, mul_symbol):
  622. if not self:
  623. return printer._print(self.ring.domain.zero)
  624. prec_mul = precedence["Mul"]
  625. prec_atom = precedence["Atom"]
  626. ring = self.ring
  627. symbols = ring.symbols
  628. ngens = ring.ngens
  629. zm = ring.zero_monom
  630. sexpvs = []
  631. for expv, coeff in self.terms():
  632. negative = ring.domain.is_negative(coeff)
  633. sign = " - " if negative else " + "
  634. sexpvs.append(sign)
  635. if expv == zm:
  636. scoeff = printer._print(coeff)
  637. if negative and scoeff.startswith("-"):
  638. scoeff = scoeff[1:]
  639. else:
  640. if negative:
  641. coeff = -coeff
  642. if coeff != self.ring.one:
  643. scoeff = printer.parenthesize(coeff, prec_mul, strict=True)
  644. else:
  645. scoeff = ''
  646. sexpv = []
  647. for i in range(ngens):
  648. exp = expv[i]
  649. if not exp:
  650. continue
  651. symbol = printer.parenthesize(symbols[i], prec_atom, strict=True)
  652. if exp != 1:
  653. if exp != int(exp) or exp < 0:
  654. sexp = printer.parenthesize(exp, prec_atom, strict=False)
  655. else:
  656. sexp = exp
  657. sexpv.append(exp_pattern % (symbol, sexp))
  658. else:
  659. sexpv.append('%s' % symbol)
  660. if scoeff:
  661. sexpv = [scoeff] + sexpv
  662. sexpvs.append(mul_symbol.join(sexpv))
  663. if sexpvs[0] in [" + ", " - "]:
  664. head = sexpvs.pop(0)
  665. if head == " - ":
  666. sexpvs.insert(0, "-")
  667. return "".join(sexpvs)
  668. @property
  669. def is_generator(self):
  670. return self in self.ring._gens_set
  671. @property
  672. def is_ground(self):
  673. return not self or (len(self) == 1 and self.ring.zero_monom in self)
  674. @property
  675. def is_monomial(self):
  676. return not self or (len(self) == 1 and self.LC == 1)
  677. @property
  678. def is_term(self):
  679. return len(self) <= 1
  680. @property
  681. def is_negative(self):
  682. return self.ring.domain.is_negative(self.LC)
  683. @property
  684. def is_positive(self):
  685. return self.ring.domain.is_positive(self.LC)
  686. @property
  687. def is_nonnegative(self):
  688. return self.ring.domain.is_nonnegative(self.LC)
  689. @property
  690. def is_nonpositive(self):
  691. return self.ring.domain.is_nonpositive(self.LC)
  692. @property
  693. def is_zero(f):
  694. return not f
  695. @property
  696. def is_one(f):
  697. return f == f.ring.one
  698. @property
  699. def is_monic(f):
  700. return f.ring.domain.is_one(f.LC)
  701. @property
  702. def is_primitive(f):
  703. return f.ring.domain.is_one(f.content())
  704. @property
  705. def is_linear(f):
  706. return all(sum(monom) <= 1 for monom in f.itermonoms())
  707. @property
  708. def is_quadratic(f):
  709. return all(sum(monom) <= 2 for monom in f.itermonoms())
  710. @property
  711. def is_squarefree(f):
  712. if not f.ring.ngens:
  713. return True
  714. return f.ring.dmp_sqf_p(f)
  715. @property
  716. def is_irreducible(f):
  717. if not f.ring.ngens:
  718. return True
  719. return f.ring.dmp_irreducible_p(f)
  720. @property
  721. def is_cyclotomic(f):
  722. if f.ring.is_univariate:
  723. return f.ring.dup_cyclotomic_p(f)
  724. else:
  725. raise MultivariatePolynomialError("cyclotomic polynomial")
  726. def __neg__(self):
  727. return self.new([ (monom, -coeff) for monom, coeff in self.iterterms() ])
  728. def __pos__(self):
  729. return self
  730. def __add__(p1, p2):
  731. """Add two polynomials.
  732. Examples
  733. ========
  734. >>> from sympy.polys.domains import ZZ
  735. >>> from sympy.polys.rings import ring
  736. >>> _, x, y = ring('x, y', ZZ)
  737. >>> (x + y)**2 + (x - y)**2
  738. 2*x**2 + 2*y**2
  739. """
  740. if not p2:
  741. return p1.copy()
  742. ring = p1.ring
  743. if isinstance(p2, ring.dtype):
  744. p = p1.copy()
  745. get = p.get
  746. zero = ring.domain.zero
  747. for k, v in p2.items():
  748. v = get(k, zero) + v
  749. if v:
  750. p[k] = v
  751. else:
  752. del p[k]
  753. return p
  754. elif isinstance(p2, PolyElement):
  755. if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
  756. pass
  757. elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
  758. return p2.__radd__(p1)
  759. else:
  760. return NotImplemented
  761. try:
  762. cp2 = ring.domain_new(p2)
  763. except CoercionFailed:
  764. return NotImplemented
  765. else:
  766. p = p1.copy()
  767. if not cp2:
  768. return p
  769. zm = ring.zero_monom
  770. if zm not in p1.keys():
  771. p[zm] = cp2
  772. else:
  773. if p2 == -p[zm]:
  774. del p[zm]
  775. else:
  776. p[zm] += cp2
  777. return p
  778. def __radd__(p1, n):
  779. p = p1.copy()
  780. if not n:
  781. return p
  782. ring = p1.ring
  783. try:
  784. n = ring.domain_new(n)
  785. except CoercionFailed:
  786. return NotImplemented
  787. else:
  788. zm = ring.zero_monom
  789. if zm not in p1.keys():
  790. p[zm] = n
  791. else:
  792. if n == -p[zm]:
  793. del p[zm]
  794. else:
  795. p[zm] += n
  796. return p
  797. def __sub__(p1, p2):
  798. """Subtract polynomial p2 from p1.
  799. Examples
  800. ========
  801. >>> from sympy.polys.domains import ZZ
  802. >>> from sympy.polys.rings import ring
  803. >>> _, x, y = ring('x, y', ZZ)
  804. >>> p1 = x + y**2
  805. >>> p2 = x*y + y**2
  806. >>> p1 - p2
  807. -x*y + x
  808. """
  809. if not p2:
  810. return p1.copy()
  811. ring = p1.ring
  812. if isinstance(p2, ring.dtype):
  813. p = p1.copy()
  814. get = p.get
  815. zero = ring.domain.zero
  816. for k, v in p2.items():
  817. v = get(k, zero) - v
  818. if v:
  819. p[k] = v
  820. else:
  821. del p[k]
  822. return p
  823. elif isinstance(p2, PolyElement):
  824. if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
  825. pass
  826. elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
  827. return p2.__rsub__(p1)
  828. else:
  829. return NotImplemented
  830. try:
  831. p2 = ring.domain_new(p2)
  832. except CoercionFailed:
  833. return NotImplemented
  834. else:
  835. p = p1.copy()
  836. zm = ring.zero_monom
  837. if zm not in p1.keys():
  838. p[zm] = -p2
  839. else:
  840. if p2 == p[zm]:
  841. del p[zm]
  842. else:
  843. p[zm] -= p2
  844. return p
  845. def __rsub__(p1, n):
  846. """n - p1 with n convertible to the coefficient domain.
  847. Examples
  848. ========
  849. >>> from sympy.polys.domains import ZZ
  850. >>> from sympy.polys.rings import ring
  851. >>> _, x, y = ring('x, y', ZZ)
  852. >>> p = x + y
  853. >>> 4 - p
  854. -x - y + 4
  855. """
  856. ring = p1.ring
  857. try:
  858. n = ring.domain_new(n)
  859. except CoercionFailed:
  860. return NotImplemented
  861. else:
  862. p = ring.zero
  863. for expv in p1:
  864. p[expv] = -p1[expv]
  865. p += n
  866. return p
  867. def __mul__(p1, p2):
  868. """Multiply two polynomials.
  869. Examples
  870. ========
  871. >>> from sympy.polys.domains import QQ
  872. >>> from sympy.polys.rings import ring
  873. >>> _, x, y = ring('x, y', QQ)
  874. >>> p1 = x + y
  875. >>> p2 = x - y
  876. >>> p1*p2
  877. x**2 - y**2
  878. """
  879. ring = p1.ring
  880. p = ring.zero
  881. if not p1 or not p2:
  882. return p
  883. elif isinstance(p2, ring.dtype):
  884. get = p.get
  885. zero = ring.domain.zero
  886. monomial_mul = ring.monomial_mul
  887. p2it = list(p2.items())
  888. for exp1, v1 in p1.items():
  889. for exp2, v2 in p2it:
  890. exp = monomial_mul(exp1, exp2)
  891. p[exp] = get(exp, zero) + v1*v2
  892. p.strip_zero()
  893. return p
  894. elif isinstance(p2, PolyElement):
  895. if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
  896. pass
  897. elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
  898. return p2.__rmul__(p1)
  899. else:
  900. return NotImplemented
  901. try:
  902. p2 = ring.domain_new(p2)
  903. except CoercionFailed:
  904. return NotImplemented
  905. else:
  906. for exp1, v1 in p1.items():
  907. v = v1*p2
  908. if v:
  909. p[exp1] = v
  910. return p
  911. def __rmul__(p1, p2):
  912. """p2 * p1 with p2 in the coefficient domain of p1.
  913. Examples
  914. ========
  915. >>> from sympy.polys.domains import ZZ
  916. >>> from sympy.polys.rings import ring
  917. >>> _, x, y = ring('x, y', ZZ)
  918. >>> p = x + y
  919. >>> 4 * p
  920. 4*x + 4*y
  921. """
  922. p = p1.ring.zero
  923. if not p2:
  924. return p
  925. try:
  926. p2 = p.ring.domain_new(p2)
  927. except CoercionFailed:
  928. return NotImplemented
  929. else:
  930. for exp1, v1 in p1.items():
  931. v = p2*v1
  932. if v:
  933. p[exp1] = v
  934. return p
  935. def __pow__(self, n):
  936. """raise polynomial to power `n`
  937. Examples
  938. ========
  939. >>> from sympy.polys.domains import ZZ
  940. >>> from sympy.polys.rings import ring
  941. >>> _, x, y = ring('x, y', ZZ)
  942. >>> p = x + y**2
  943. >>> p**3
  944. x**3 + 3*x**2*y**2 + 3*x*y**4 + y**6
  945. """
  946. ring = self.ring
  947. if not n:
  948. if self:
  949. return ring.one
  950. else:
  951. raise ValueError("0**0")
  952. elif len(self) == 1:
  953. monom, coeff = list(self.items())[0]
  954. p = ring.zero
  955. if coeff == 1:
  956. p[ring.monomial_pow(monom, n)] = coeff
  957. else:
  958. p[ring.monomial_pow(monom, n)] = coeff**n
  959. return p
  960. # For ring series, we need negative and rational exponent support only
  961. # with monomials.
  962. n = int(n)
  963. if n < 0:
  964. raise ValueError("Negative exponent")
  965. elif n == 1:
  966. return self.copy()
  967. elif n == 2:
  968. return self.square()
  969. elif n == 3:
  970. return self*self.square()
  971. elif len(self) <= 5: # TODO: use an actual density measure
  972. return self._pow_multinomial(n)
  973. else:
  974. return self._pow_generic(n)
  975. def _pow_generic(self, n):
  976. p = self.ring.one
  977. c = self
  978. while True:
  979. if n & 1:
  980. p = p*c
  981. n -= 1
  982. if not n:
  983. break
  984. c = c.square()
  985. n = n // 2
  986. return p
  987. def _pow_multinomial(self, n):
  988. multinomials = multinomial_coefficients(len(self), n).items()
  989. monomial_mulpow = self.ring.monomial_mulpow
  990. zero_monom = self.ring.zero_monom
  991. terms = self.items()
  992. zero = self.ring.domain.zero
  993. poly = self.ring.zero
  994. for multinomial, multinomial_coeff in multinomials:
  995. product_monom = zero_monom
  996. product_coeff = multinomial_coeff
  997. for exp, (monom, coeff) in zip(multinomial, terms):
  998. if exp:
  999. product_monom = monomial_mulpow(product_monom, monom, exp)
  1000. product_coeff *= coeff**exp
  1001. monom = tuple(product_monom)
  1002. coeff = product_coeff
  1003. coeff = poly.get(monom, zero) + coeff
  1004. if coeff:
  1005. poly[monom] = coeff
  1006. elif monom in poly:
  1007. del poly[monom]
  1008. return poly
  1009. def square(self):
  1010. """square of a polynomial
  1011. Examples
  1012. ========
  1013. >>> from sympy.polys.rings import ring
  1014. >>> from sympy.polys.domains import ZZ
  1015. >>> _, x, y = ring('x, y', ZZ)
  1016. >>> p = x + y**2
  1017. >>> p.square()
  1018. x**2 + 2*x*y**2 + y**4
  1019. """
  1020. ring = self.ring
  1021. p = ring.zero
  1022. get = p.get
  1023. keys = list(self.keys())
  1024. zero = ring.domain.zero
  1025. monomial_mul = ring.monomial_mul
  1026. for i in range(len(keys)):
  1027. k1 = keys[i]
  1028. pk = self[k1]
  1029. for j in range(i):
  1030. k2 = keys[j]
  1031. exp = monomial_mul(k1, k2)
  1032. p[exp] = get(exp, zero) + pk*self[k2]
  1033. p = p.imul_num(2)
  1034. get = p.get
  1035. for k, v in self.items():
  1036. k2 = monomial_mul(k, k)
  1037. p[k2] = get(k2, zero) + v**2
  1038. p.strip_zero()
  1039. return p
  1040. def __divmod__(p1, p2):
  1041. ring = p1.ring
  1042. if not p2:
  1043. raise ZeroDivisionError("polynomial division")
  1044. elif isinstance(p2, ring.dtype):
  1045. return p1.div(p2)
  1046. elif isinstance(p2, PolyElement):
  1047. if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
  1048. pass
  1049. elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
  1050. return p2.__rdivmod__(p1)
  1051. else:
  1052. return NotImplemented
  1053. try:
  1054. p2 = ring.domain_new(p2)
  1055. except CoercionFailed:
  1056. return NotImplemented
  1057. else:
  1058. return (p1.quo_ground(p2), p1.rem_ground(p2))
  1059. def __rdivmod__(p1, p2):
  1060. return NotImplemented
  1061. def __mod__(p1, p2):
  1062. ring = p1.ring
  1063. if not p2:
  1064. raise ZeroDivisionError("polynomial division")
  1065. elif isinstance(p2, ring.dtype):
  1066. return p1.rem(p2)
  1067. elif isinstance(p2, PolyElement):
  1068. if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
  1069. pass
  1070. elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
  1071. return p2.__rmod__(p1)
  1072. else:
  1073. return NotImplemented
  1074. try:
  1075. p2 = ring.domain_new(p2)
  1076. except CoercionFailed:
  1077. return NotImplemented
  1078. else:
  1079. return p1.rem_ground(p2)
  1080. def __rmod__(p1, p2):
  1081. return NotImplemented
  1082. def __truediv__(p1, p2):
  1083. ring = p1.ring
  1084. if not p2:
  1085. raise ZeroDivisionError("polynomial division")
  1086. elif isinstance(p2, ring.dtype):
  1087. if p2.is_monomial:
  1088. return p1*(p2**(-1))
  1089. else:
  1090. return p1.quo(p2)
  1091. elif isinstance(p2, PolyElement):
  1092. if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
  1093. pass
  1094. elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
  1095. return p2.__rtruediv__(p1)
  1096. else:
  1097. return NotImplemented
  1098. try:
  1099. p2 = ring.domain_new(p2)
  1100. except CoercionFailed:
  1101. return NotImplemented
  1102. else:
  1103. return p1.quo_ground(p2)
  1104. def __rtruediv__(p1, p2):
  1105. return NotImplemented
  1106. __floordiv__ = __truediv__
  1107. __rfloordiv__ = __rtruediv__
  1108. # TODO: use // (__floordiv__) for exquo()?
  1109. def _term_div(self):
  1110. zm = self.ring.zero_monom
  1111. domain = self.ring.domain
  1112. domain_quo = domain.quo
  1113. monomial_div = self.ring.monomial_div
  1114. if domain.is_Field:
  1115. def term_div(a_lm_a_lc, b_lm_b_lc):
  1116. a_lm, a_lc = a_lm_a_lc
  1117. b_lm, b_lc = b_lm_b_lc
  1118. if b_lm == zm: # apparently this is a very common case
  1119. monom = a_lm
  1120. else:
  1121. monom = monomial_div(a_lm, b_lm)
  1122. if monom is not None:
  1123. return monom, domain_quo(a_lc, b_lc)
  1124. else:
  1125. return None
  1126. else:
  1127. def term_div(a_lm_a_lc, b_lm_b_lc):
  1128. a_lm, a_lc = a_lm_a_lc
  1129. b_lm, b_lc = b_lm_b_lc
  1130. if b_lm == zm: # apparently this is a very common case
  1131. monom = a_lm
  1132. else:
  1133. monom = monomial_div(a_lm, b_lm)
  1134. if not (monom is None or a_lc % b_lc):
  1135. return monom, domain_quo(a_lc, b_lc)
  1136. else:
  1137. return None
  1138. return term_div
  1139. def div(self, fv):
  1140. """Division algorithm, see [CLO] p64.
  1141. fv array of polynomials
  1142. return qv, r such that
  1143. self = sum(fv[i]*qv[i]) + r
  1144. All polynomials are required not to be Laurent polynomials.
  1145. Examples
  1146. ========
  1147. >>> from sympy.polys.rings import ring
  1148. >>> from sympy.polys.domains import ZZ
  1149. >>> _, x, y = ring('x, y', ZZ)
  1150. >>> f = x**3
  1151. >>> f0 = x - y**2
  1152. >>> f1 = x - y
  1153. >>> qv, r = f.div((f0, f1))
  1154. >>> qv[0]
  1155. x**2 + x*y**2 + y**4
  1156. >>> qv[1]
  1157. 0
  1158. >>> r
  1159. y**6
  1160. """
  1161. ring = self.ring
  1162. ret_single = False
  1163. if isinstance(fv, PolyElement):
  1164. ret_single = True
  1165. fv = [fv]
  1166. if not all(fv):
  1167. raise ZeroDivisionError("polynomial division")
  1168. if not self:
  1169. if ret_single:
  1170. return ring.zero, ring.zero
  1171. else:
  1172. return [], ring.zero
  1173. for f in fv:
  1174. if f.ring != ring:
  1175. raise ValueError('self and f must have the same ring')
  1176. s = len(fv)
  1177. qv = [ring.zero for i in range(s)]
  1178. p = self.copy()
  1179. r = ring.zero
  1180. term_div = self._term_div()
  1181. expvs = [fx.leading_expv() for fx in fv]
  1182. while p:
  1183. i = 0
  1184. divoccurred = 0
  1185. while i < s and divoccurred == 0:
  1186. expv = p.leading_expv()
  1187. term = term_div((expv, p[expv]), (expvs[i], fv[i][expvs[i]]))
  1188. if term is not None:
  1189. expv1, c = term
  1190. qv[i] = qv[i]._iadd_monom((expv1, c))
  1191. p = p._iadd_poly_monom(fv[i], (expv1, -c))
  1192. divoccurred = 1
  1193. else:
  1194. i += 1
  1195. if not divoccurred:
  1196. expv = p.leading_expv()
  1197. r = r._iadd_monom((expv, p[expv]))
  1198. del p[expv]
  1199. if expv == ring.zero_monom:
  1200. r += p
  1201. if ret_single:
  1202. if not qv:
  1203. return ring.zero, r
  1204. else:
  1205. return qv[0], r
  1206. else:
  1207. return qv, r
  1208. def rem(self, G):
  1209. f = self
  1210. if isinstance(G, PolyElement):
  1211. G = [G]
  1212. if not all(G):
  1213. raise ZeroDivisionError("polynomial division")
  1214. ring = f.ring
  1215. domain = ring.domain
  1216. zero = domain.zero
  1217. monomial_mul = ring.monomial_mul
  1218. r = ring.zero
  1219. term_div = f._term_div()
  1220. ltf = f.LT
  1221. f = f.copy()
  1222. get = f.get
  1223. while f:
  1224. for g in G:
  1225. tq = term_div(ltf, g.LT)
  1226. if tq is not None:
  1227. m, c = tq
  1228. for mg, cg in g.iterterms():
  1229. m1 = monomial_mul(mg, m)
  1230. c1 = get(m1, zero) - c*cg
  1231. if not c1:
  1232. del f[m1]
  1233. else:
  1234. f[m1] = c1
  1235. ltm = f.leading_expv()
  1236. if ltm is not None:
  1237. ltf = ltm, f[ltm]
  1238. break
  1239. else:
  1240. ltm, ltc = ltf
  1241. if ltm in r:
  1242. r[ltm] += ltc
  1243. else:
  1244. r[ltm] = ltc
  1245. del f[ltm]
  1246. ltm = f.leading_expv()
  1247. if ltm is not None:
  1248. ltf = ltm, f[ltm]
  1249. return r
  1250. def quo(f, G):
  1251. return f.div(G)[0]
  1252. def exquo(f, G):
  1253. q, r = f.div(G)
  1254. if not r:
  1255. return q
  1256. else:
  1257. raise ExactQuotientFailed(f, G)
  1258. def _iadd_monom(self, mc):
  1259. """add to self the monomial coeff*x0**i0*x1**i1*...
  1260. unless self is a generator -- then just return the sum of the two.
  1261. mc is a tuple, (monom, coeff), where monomial is (i0, i1, ...)
  1262. Examples
  1263. ========
  1264. >>> from sympy.polys.rings import ring
  1265. >>> from sympy.polys.domains import ZZ
  1266. >>> _, x, y = ring('x, y', ZZ)
  1267. >>> p = x**4 + 2*y
  1268. >>> m = (1, 2)
  1269. >>> p1 = p._iadd_monom((m, 5))
  1270. >>> p1
  1271. x**4 + 5*x*y**2 + 2*y
  1272. >>> p1 is p
  1273. True
  1274. >>> p = x
  1275. >>> p1 = p._iadd_monom((m, 5))
  1276. >>> p1
  1277. 5*x*y**2 + x
  1278. >>> p1 is p
  1279. False
  1280. """
  1281. if self in self.ring._gens_set:
  1282. cpself = self.copy()
  1283. else:
  1284. cpself = self
  1285. expv, coeff = mc
  1286. c = cpself.get(expv)
  1287. if c is None:
  1288. cpself[expv] = coeff
  1289. else:
  1290. c += coeff
  1291. if c:
  1292. cpself[expv] = c
  1293. else:
  1294. del cpself[expv]
  1295. return cpself
  1296. def _iadd_poly_monom(self, p2, mc):
  1297. """add to self the product of (p)*(coeff*x0**i0*x1**i1*...)
  1298. unless self is a generator -- then just return the sum of the two.
  1299. mc is a tuple, (monom, coeff), where monomial is (i0, i1, ...)
  1300. Examples
  1301. ========
  1302. >>> from sympy.polys.rings import ring
  1303. >>> from sympy.polys.domains import ZZ
  1304. >>> _, x, y, z = ring('x, y, z', ZZ)
  1305. >>> p1 = x**4 + 2*y
  1306. >>> p2 = y + z
  1307. >>> m = (1, 2, 3)
  1308. >>> p1 = p1._iadd_poly_monom(p2, (m, 3))
  1309. >>> p1
  1310. x**4 + 3*x*y**3*z**3 + 3*x*y**2*z**4 + 2*y
  1311. """
  1312. p1 = self
  1313. if p1 in p1.ring._gens_set:
  1314. p1 = p1.copy()
  1315. (m, c) = mc
  1316. get = p1.get
  1317. zero = p1.ring.domain.zero
  1318. monomial_mul = p1.ring.monomial_mul
  1319. for k, v in p2.items():
  1320. ka = monomial_mul(k, m)
  1321. coeff = get(ka, zero) + v*c
  1322. if coeff:
  1323. p1[ka] = coeff
  1324. else:
  1325. del p1[ka]
  1326. return p1
  1327. def degree(f, x=None):
  1328. """
  1329. The leading degree in ``x`` or the main variable.
  1330. Note that the degree of 0 is negative infinity (the SymPy object -oo).
  1331. """
  1332. i = f.ring.index(x)
  1333. if not f:
  1334. return -oo
  1335. elif i < 0:
  1336. return 0
  1337. else:
  1338. return max([ monom[i] for monom in f.itermonoms() ])
  1339. def degrees(f):
  1340. """
  1341. A tuple containing leading degrees in all variables.
  1342. Note that the degree of 0 is negative infinity (the SymPy object -oo)
  1343. """
  1344. if not f:
  1345. return (-oo,)*f.ring.ngens
  1346. else:
  1347. return tuple(map(max, list(zip(*f.itermonoms()))))
  1348. def tail_degree(f, x=None):
  1349. """
  1350. The tail degree in ``x`` or the main variable.
  1351. Note that the degree of 0 is negative infinity (the SymPy object -oo)
  1352. """
  1353. i = f.ring.index(x)
  1354. if not f:
  1355. return -oo
  1356. elif i < 0:
  1357. return 0
  1358. else:
  1359. return min([ monom[i] for monom in f.itermonoms() ])
  1360. def tail_degrees(f):
  1361. """
  1362. A tuple containing tail degrees in all variables.
  1363. Note that the degree of 0 is negative infinity (the SymPy object -oo)
  1364. """
  1365. if not f:
  1366. return (-oo,)*f.ring.ngens
  1367. else:
  1368. return tuple(map(min, list(zip(*f.itermonoms()))))
  1369. def leading_expv(self):
  1370. """Leading monomial tuple according to the monomial ordering.
  1371. Examples
  1372. ========
  1373. >>> from sympy.polys.rings import ring
  1374. >>> from sympy.polys.domains import ZZ
  1375. >>> _, x, y, z = ring('x, y, z', ZZ)
  1376. >>> p = x**4 + x**3*y + x**2*z**2 + z**7
  1377. >>> p.leading_expv()
  1378. (4, 0, 0)
  1379. """
  1380. if self:
  1381. return self.ring.leading_expv(self)
  1382. else:
  1383. return None
  1384. def _get_coeff(self, expv):
  1385. return self.get(expv, self.ring.domain.zero)
  1386. def coeff(self, element):
  1387. """
  1388. Returns the coefficient that stands next to the given monomial.
  1389. Parameters
  1390. ==========
  1391. element : PolyElement (with ``is_monomial = True``) or 1
  1392. Examples
  1393. ========
  1394. >>> from sympy.polys.rings import ring
  1395. >>> from sympy.polys.domains import ZZ
  1396. >>> _, x, y, z = ring("x,y,z", ZZ)
  1397. >>> f = 3*x**2*y - x*y*z + 7*z**3 + 23
  1398. >>> f.coeff(x**2*y)
  1399. 3
  1400. >>> f.coeff(x*y)
  1401. 0
  1402. >>> f.coeff(1)
  1403. 23
  1404. """
  1405. if element == 1:
  1406. return self._get_coeff(self.ring.zero_monom)
  1407. elif isinstance(element, self.ring.dtype):
  1408. terms = list(element.iterterms())
  1409. if len(terms) == 1:
  1410. monom, coeff = terms[0]
  1411. if coeff == self.ring.domain.one:
  1412. return self._get_coeff(monom)
  1413. raise ValueError("expected a monomial, got %s" % element)
  1414. def const(self):
  1415. """Returns the constant coeffcient. """
  1416. return self._get_coeff(self.ring.zero_monom)
  1417. @property
  1418. def LC(self):
  1419. return self._get_coeff(self.leading_expv())
  1420. @property
  1421. def LM(self):
  1422. expv = self.leading_expv()
  1423. if expv is None:
  1424. return self.ring.zero_monom
  1425. else:
  1426. return expv
  1427. def leading_monom(self):
  1428. """
  1429. Leading monomial as a polynomial element.
  1430. Examples
  1431. ========
  1432. >>> from sympy.polys.rings import ring
  1433. >>> from sympy.polys.domains import ZZ
  1434. >>> _, x, y = ring('x, y', ZZ)
  1435. >>> (3*x*y + y**2).leading_monom()
  1436. x*y
  1437. """
  1438. p = self.ring.zero
  1439. expv = self.leading_expv()
  1440. if expv:
  1441. p[expv] = self.ring.domain.one
  1442. return p
  1443. @property
  1444. def LT(self):
  1445. expv = self.leading_expv()
  1446. if expv is None:
  1447. return (self.ring.zero_monom, self.ring.domain.zero)
  1448. else:
  1449. return (expv, self._get_coeff(expv))
  1450. def leading_term(self):
  1451. """Leading term as a polynomial element.
  1452. Examples
  1453. ========
  1454. >>> from sympy.polys.rings import ring
  1455. >>> from sympy.polys.domains import ZZ
  1456. >>> _, x, y = ring('x, y', ZZ)
  1457. >>> (3*x*y + y**2).leading_term()
  1458. 3*x*y
  1459. """
  1460. p = self.ring.zero
  1461. expv = self.leading_expv()
  1462. if expv is not None:
  1463. p[expv] = self[expv]
  1464. return p
  1465. def _sorted(self, seq, order):
  1466. if order is None:
  1467. order = self.ring.order
  1468. else:
  1469. order = OrderOpt.preprocess(order)
  1470. if order is lex:
  1471. return sorted(seq, key=lambda monom: monom[0], reverse=True)
  1472. else:
  1473. return sorted(seq, key=lambda monom: order(monom[0]), reverse=True)
  1474. def coeffs(self, order=None):
  1475. """Ordered list of polynomial coefficients.
  1476. Parameters
  1477. ==========
  1478. order : :class:`~.MonomialOrder` or coercible, optional
  1479. Examples
  1480. ========
  1481. >>> from sympy.polys.rings import ring
  1482. >>> from sympy.polys.domains import ZZ
  1483. >>> from sympy.polys.orderings import lex, grlex
  1484. >>> _, x, y = ring("x, y", ZZ, lex)
  1485. >>> f = x*y**7 + 2*x**2*y**3
  1486. >>> f.coeffs()
  1487. [2, 1]
  1488. >>> f.coeffs(grlex)
  1489. [1, 2]
  1490. """
  1491. return [ coeff for _, coeff in self.terms(order) ]
  1492. def monoms(self, order=None):
  1493. """Ordered list of polynomial monomials.
  1494. Parameters
  1495. ==========
  1496. order : :class:`~.MonomialOrder` or coercible, optional
  1497. Examples
  1498. ========
  1499. >>> from sympy.polys.rings import ring
  1500. >>> from sympy.polys.domains import ZZ
  1501. >>> from sympy.polys.orderings import lex, grlex
  1502. >>> _, x, y = ring("x, y", ZZ, lex)
  1503. >>> f = x*y**7 + 2*x**2*y**3
  1504. >>> f.monoms()
  1505. [(2, 3), (1, 7)]
  1506. >>> f.monoms(grlex)
  1507. [(1, 7), (2, 3)]
  1508. """
  1509. return [ monom for monom, _ in self.terms(order) ]
  1510. def terms(self, order=None):
  1511. """Ordered list of polynomial terms.
  1512. Parameters
  1513. ==========
  1514. order : :class:`~.MonomialOrder` or coercible, optional
  1515. Examples
  1516. ========
  1517. >>> from sympy.polys.rings import ring
  1518. >>> from sympy.polys.domains import ZZ
  1519. >>> from sympy.polys.orderings import lex, grlex
  1520. >>> _, x, y = ring("x, y", ZZ, lex)
  1521. >>> f = x*y**7 + 2*x**2*y**3
  1522. >>> f.terms()
  1523. [((2, 3), 2), ((1, 7), 1)]
  1524. >>> f.terms(grlex)
  1525. [((1, 7), 1), ((2, 3), 2)]
  1526. """
  1527. return self._sorted(list(self.items()), order)
  1528. def itercoeffs(self):
  1529. """Iterator over coefficients of a polynomial. """
  1530. return iter(self.values())
  1531. def itermonoms(self):
  1532. """Iterator over monomials of a polynomial. """
  1533. return iter(self.keys())
  1534. def iterterms(self):
  1535. """Iterator over terms of a polynomial. """
  1536. return iter(self.items())
  1537. def listcoeffs(self):
  1538. """Unordered list of polynomial coefficients. """
  1539. return list(self.values())
  1540. def listmonoms(self):
  1541. """Unordered list of polynomial monomials. """
  1542. return list(self.keys())
  1543. def listterms(self):
  1544. """Unordered list of polynomial terms. """
  1545. return list(self.items())
  1546. def imul_num(p, c):
  1547. """multiply inplace the polynomial p by an element in the
  1548. coefficient ring, provided p is not one of the generators;
  1549. else multiply not inplace
  1550. Examples
  1551. ========
  1552. >>> from sympy.polys.rings import ring
  1553. >>> from sympy.polys.domains import ZZ
  1554. >>> _, x, y = ring('x, y', ZZ)
  1555. >>> p = x + y**2
  1556. >>> p1 = p.imul_num(3)
  1557. >>> p1
  1558. 3*x + 3*y**2
  1559. >>> p1 is p
  1560. True
  1561. >>> p = x
  1562. >>> p1 = p.imul_num(3)
  1563. >>> p1
  1564. 3*x
  1565. >>> p1 is p
  1566. False
  1567. """
  1568. if p in p.ring._gens_set:
  1569. return p*c
  1570. if not c:
  1571. p.clear()
  1572. return
  1573. for exp in p:
  1574. p[exp] *= c
  1575. return p
  1576. def content(f):
  1577. """Returns GCD of polynomial's coefficients. """
  1578. domain = f.ring.domain
  1579. cont = domain.zero
  1580. gcd = domain.gcd
  1581. for coeff in f.itercoeffs():
  1582. cont = gcd(cont, coeff)
  1583. return cont
  1584. def primitive(f):
  1585. """Returns content and a primitive polynomial. """
  1586. cont = f.content()
  1587. return cont, f.quo_ground(cont)
  1588. def monic(f):
  1589. """Divides all coefficients by the leading coefficient. """
  1590. if not f:
  1591. return f
  1592. else:
  1593. return f.quo_ground(f.LC)
  1594. def mul_ground(f, x):
  1595. if not x:
  1596. return f.ring.zero
  1597. terms = [ (monom, coeff*x) for monom, coeff in f.iterterms() ]
  1598. return f.new(terms)
  1599. def mul_monom(f, monom):
  1600. monomial_mul = f.ring.monomial_mul
  1601. terms = [ (monomial_mul(f_monom, monom), f_coeff) for f_monom, f_coeff in f.items() ]
  1602. return f.new(terms)
  1603. def mul_term(f, term):
  1604. monom, coeff = term
  1605. if not f or not coeff:
  1606. return f.ring.zero
  1607. elif monom == f.ring.zero_monom:
  1608. return f.mul_ground(coeff)
  1609. monomial_mul = f.ring.monomial_mul
  1610. terms = [ (monomial_mul(f_monom, monom), f_coeff*coeff) for f_monom, f_coeff in f.items() ]
  1611. return f.new(terms)
  1612. def quo_ground(f, x):
  1613. domain = f.ring.domain
  1614. if not x:
  1615. raise ZeroDivisionError('polynomial division')
  1616. if not f or x == domain.one:
  1617. return f
  1618. if domain.is_Field:
  1619. quo = domain.quo
  1620. terms = [ (monom, quo(coeff, x)) for monom, coeff in f.iterterms() ]
  1621. else:
  1622. terms = [ (monom, coeff // x) for monom, coeff in f.iterterms() if not (coeff % x) ]
  1623. return f.new(terms)
  1624. def quo_term(f, term):
  1625. monom, coeff = term
  1626. if not coeff:
  1627. raise ZeroDivisionError("polynomial division")
  1628. elif not f:
  1629. return f.ring.zero
  1630. elif monom == f.ring.zero_monom:
  1631. return f.quo_ground(coeff)
  1632. term_div = f._term_div()
  1633. terms = [ term_div(t, term) for t in f.iterterms() ]
  1634. return f.new([ t for t in terms if t is not None ])
  1635. def trunc_ground(f, p):
  1636. if f.ring.domain.is_ZZ:
  1637. terms = []
  1638. for monom, coeff in f.iterterms():
  1639. coeff = coeff % p
  1640. if coeff > p // 2:
  1641. coeff = coeff - p
  1642. terms.append((monom, coeff))
  1643. else:
  1644. terms = [ (monom, coeff % p) for monom, coeff in f.iterterms() ]
  1645. poly = f.new(terms)
  1646. poly.strip_zero()
  1647. return poly
  1648. rem_ground = trunc_ground
  1649. def extract_ground(self, g):
  1650. f = self
  1651. fc = f.content()
  1652. gc = g.content()
  1653. gcd = f.ring.domain.gcd(fc, gc)
  1654. f = f.quo_ground(gcd)
  1655. g = g.quo_ground(gcd)
  1656. return gcd, f, g
  1657. def _norm(f, norm_func):
  1658. if not f:
  1659. return f.ring.domain.zero
  1660. else:
  1661. ground_abs = f.ring.domain.abs
  1662. return norm_func([ ground_abs(coeff) for coeff in f.itercoeffs() ])
  1663. def max_norm(f):
  1664. return f._norm(max)
  1665. def l1_norm(f):
  1666. return f._norm(sum)
  1667. def deflate(f, *G):
  1668. ring = f.ring
  1669. polys = [f] + list(G)
  1670. J = [0]*ring.ngens
  1671. for p in polys:
  1672. for monom in p.itermonoms():
  1673. for i, m in enumerate(monom):
  1674. J[i] = igcd(J[i], m)
  1675. for i, b in enumerate(J):
  1676. if not b:
  1677. J[i] = 1
  1678. J = tuple(J)
  1679. if all(b == 1 for b in J):
  1680. return J, polys
  1681. H = []
  1682. for p in polys:
  1683. h = ring.zero
  1684. for I, coeff in p.iterterms():
  1685. N = [ i // j for i, j in zip(I, J) ]
  1686. h[tuple(N)] = coeff
  1687. H.append(h)
  1688. return J, H
  1689. def inflate(f, J):
  1690. poly = f.ring.zero
  1691. for I, coeff in f.iterterms():
  1692. N = [ i*j for i, j in zip(I, J) ]
  1693. poly[tuple(N)] = coeff
  1694. return poly
  1695. def lcm(self, g):
  1696. f = self
  1697. domain = f.ring.domain
  1698. if not domain.is_Field:
  1699. fc, f = f.primitive()
  1700. gc, g = g.primitive()
  1701. c = domain.lcm(fc, gc)
  1702. h = (f*g).quo(f.gcd(g))
  1703. if not domain.is_Field:
  1704. return h.mul_ground(c)
  1705. else:
  1706. return h.monic()
  1707. def gcd(f, g):
  1708. return f.cofactors(g)[0]
  1709. def cofactors(f, g):
  1710. if not f and not g:
  1711. zero = f.ring.zero
  1712. return zero, zero, zero
  1713. elif not f:
  1714. h, cff, cfg = f._gcd_zero(g)
  1715. return h, cff, cfg
  1716. elif not g:
  1717. h, cfg, cff = g._gcd_zero(f)
  1718. return h, cff, cfg
  1719. elif len(f) == 1:
  1720. h, cff, cfg = f._gcd_monom(g)
  1721. return h, cff, cfg
  1722. elif len(g) == 1:
  1723. h, cfg, cff = g._gcd_monom(f)
  1724. return h, cff, cfg
  1725. J, (f, g) = f.deflate(g)
  1726. h, cff, cfg = f._gcd(g)
  1727. return (h.inflate(J), cff.inflate(J), cfg.inflate(J))
  1728. def _gcd_zero(f, g):
  1729. one, zero = f.ring.one, f.ring.zero
  1730. if g.is_nonnegative:
  1731. return g, zero, one
  1732. else:
  1733. return -g, zero, -one
  1734. def _gcd_monom(f, g):
  1735. ring = f.ring
  1736. ground_gcd = ring.domain.gcd
  1737. ground_quo = ring.domain.quo
  1738. monomial_gcd = ring.monomial_gcd
  1739. monomial_ldiv = ring.monomial_ldiv
  1740. mf, cf = list(f.iterterms())[0]
  1741. _mgcd, _cgcd = mf, cf
  1742. for mg, cg in g.iterterms():
  1743. _mgcd = monomial_gcd(_mgcd, mg)
  1744. _cgcd = ground_gcd(_cgcd, cg)
  1745. h = f.new([(_mgcd, _cgcd)])
  1746. cff = f.new([(monomial_ldiv(mf, _mgcd), ground_quo(cf, _cgcd))])
  1747. cfg = f.new([(monomial_ldiv(mg, _mgcd), ground_quo(cg, _cgcd)) for mg, cg in g.iterterms()])
  1748. return h, cff, cfg
  1749. def _gcd(f, g):
  1750. ring = f.ring
  1751. if ring.domain.is_QQ:
  1752. return f._gcd_QQ(g)
  1753. elif ring.domain.is_ZZ:
  1754. return f._gcd_ZZ(g)
  1755. else: # TODO: don't use dense representation (port PRS algorithms)
  1756. return ring.dmp_inner_gcd(f, g)
  1757. def _gcd_ZZ(f, g):
  1758. return heugcd(f, g)
  1759. def _gcd_QQ(self, g):
  1760. f = self
  1761. ring = f.ring
  1762. new_ring = ring.clone(domain=ring.domain.get_ring())
  1763. cf, f = f.clear_denoms()
  1764. cg, g = g.clear_denoms()
  1765. f = f.set_ring(new_ring)
  1766. g = g.set_ring(new_ring)
  1767. h, cff, cfg = f._gcd_ZZ(g)
  1768. h = h.set_ring(ring)
  1769. c, h = h.LC, h.monic()
  1770. cff = cff.set_ring(ring).mul_ground(ring.domain.quo(c, cf))
  1771. cfg = cfg.set_ring(ring).mul_ground(ring.domain.quo(c, cg))
  1772. return h, cff, cfg
  1773. def cancel(self, g):
  1774. """
  1775. Cancel common factors in a rational function ``f/g``.
  1776. Examples
  1777. ========
  1778. >>> from sympy.polys import ring, ZZ
  1779. >>> R, x,y = ring("x,y", ZZ)
  1780. >>> (2*x**2 - 2).cancel(x**2 - 2*x + 1)
  1781. (2*x + 2, x - 1)
  1782. """
  1783. f = self
  1784. ring = f.ring
  1785. if not f:
  1786. return f, ring.one
  1787. domain = ring.domain
  1788. if not (domain.is_Field and domain.has_assoc_Ring):
  1789. _, p, q = f.cofactors(g)
  1790. else:
  1791. new_ring = ring.clone(domain=domain.get_ring())
  1792. cq, f = f.clear_denoms()
  1793. cp, g = g.clear_denoms()
  1794. f = f.set_ring(new_ring)
  1795. g = g.set_ring(new_ring)
  1796. _, p, q = f.cofactors(g)
  1797. _, cp, cq = new_ring.domain.cofactors(cp, cq)
  1798. p = p.set_ring(ring)
  1799. q = q.set_ring(ring)
  1800. p = p.mul_ground(cp)
  1801. q = q.mul_ground(cq)
  1802. # Make canonical with respect to sign or quadrant in the case of ZZ_I
  1803. # or QQ_I. This ensures that the LC of the denominator is canonical by
  1804. # multiplying top and bottom by a unit of the ring.
  1805. u = q.canonical_unit()
  1806. if u == domain.one:
  1807. p, q = p, q
  1808. elif u == -domain.one:
  1809. p, q = -p, -q
  1810. else:
  1811. p = p.mul_ground(u)
  1812. q = q.mul_ground(u)
  1813. return p, q
  1814. def canonical_unit(f):
  1815. domain = f.ring.domain
  1816. return domain.canonical_unit(f.LC)
  1817. def diff(f, x):
  1818. """Computes partial derivative in ``x``.
  1819. Examples
  1820. ========
  1821. >>> from sympy.polys.rings import ring
  1822. >>> from sympy.polys.domains import ZZ
  1823. >>> _, x, y = ring("x,y", ZZ)
  1824. >>> p = x + x**2*y**3
  1825. >>> p.diff(x)
  1826. 2*x*y**3 + 1
  1827. """
  1828. ring = f.ring
  1829. i = ring.index(x)
  1830. m = ring.monomial_basis(i)
  1831. g = ring.zero
  1832. for expv, coeff in f.iterterms():
  1833. if expv[i]:
  1834. e = ring.monomial_ldiv(expv, m)
  1835. g[e] = ring.domain_new(coeff*expv[i])
  1836. return g
  1837. def __call__(f, *values):
  1838. if 0 < len(values) <= f.ring.ngens:
  1839. return f.evaluate(list(zip(f.ring.gens, values)))
  1840. else:
  1841. raise ValueError("expected at least 1 and at most %s values, got %s" % (f.ring.ngens, len(values)))
  1842. def evaluate(self, x, a=None):
  1843. f = self
  1844. if isinstance(x, list) and a is None:
  1845. (X, a), x = x[0], x[1:]
  1846. f = f.evaluate(X, a)
  1847. if not x:
  1848. return f
  1849. else:
  1850. x = [ (Y.drop(X), a) for (Y, a) in x ]
  1851. return f.evaluate(x)
  1852. ring = f.ring
  1853. i = ring.index(x)
  1854. a = ring.domain.convert(a)
  1855. if ring.ngens == 1:
  1856. result = ring.domain.zero
  1857. for (n,), coeff in f.iterterms():
  1858. result += coeff*a**n
  1859. return result
  1860. else:
  1861. poly = ring.drop(x).zero
  1862. for monom, coeff in f.iterterms():
  1863. n, monom = monom[i], monom[:i] + monom[i+1:]
  1864. coeff = coeff*a**n
  1865. if monom in poly:
  1866. coeff = coeff + poly[monom]
  1867. if coeff:
  1868. poly[monom] = coeff
  1869. else:
  1870. del poly[monom]
  1871. else:
  1872. if coeff:
  1873. poly[monom] = coeff
  1874. return poly
  1875. def subs(self, x, a=None):
  1876. f = self
  1877. if isinstance(x, list) and a is None:
  1878. for X, a in x:
  1879. f = f.subs(X, a)
  1880. return f
  1881. ring = f.ring
  1882. i = ring.index(x)
  1883. a = ring.domain.convert(a)
  1884. if ring.ngens == 1:
  1885. result = ring.domain.zero
  1886. for (n,), coeff in f.iterterms():
  1887. result += coeff*a**n
  1888. return ring.ground_new(result)
  1889. else:
  1890. poly = ring.zero
  1891. for monom, coeff in f.iterterms():
  1892. n, monom = monom[i], monom[:i] + (0,) + monom[i+1:]
  1893. coeff = coeff*a**n
  1894. if monom in poly:
  1895. coeff = coeff + poly[monom]
  1896. if coeff:
  1897. poly[monom] = coeff
  1898. else:
  1899. del poly[monom]
  1900. else:
  1901. if coeff:
  1902. poly[monom] = coeff
  1903. return poly
  1904. def compose(f, x, a=None):
  1905. ring = f.ring
  1906. poly = ring.zero
  1907. gens_map = dict(list(zip(ring.gens, list(range(ring.ngens)))))
  1908. if a is not None:
  1909. replacements = [(x, a)]
  1910. else:
  1911. if isinstance(x, list):
  1912. replacements = list(x)
  1913. elif isinstance(x, dict):
  1914. replacements = sorted(list(x.items()), key=lambda k: gens_map[k[0]])
  1915. else:
  1916. raise ValueError("expected a generator, value pair a sequence of such pairs")
  1917. for k, (x, g) in enumerate(replacements):
  1918. replacements[k] = (gens_map[x], ring.ring_new(g))
  1919. for monom, coeff in f.iterterms():
  1920. monom = list(monom)
  1921. subpoly = ring.one
  1922. for i, g in replacements:
  1923. n, monom[i] = monom[i], 0
  1924. if n:
  1925. subpoly *= g**n
  1926. subpoly = subpoly.mul_term((tuple(monom), coeff))
  1927. poly += subpoly
  1928. return poly
  1929. # TODO: following methods should point to polynomial
  1930. # representation independent algorithm implementations.
  1931. def pdiv(f, g):
  1932. return f.ring.dmp_pdiv(f, g)
  1933. def prem(f, g):
  1934. return f.ring.dmp_prem(f, g)
  1935. def pquo(f, g):
  1936. return f.ring.dmp_quo(f, g)
  1937. def pexquo(f, g):
  1938. return f.ring.dmp_exquo(f, g)
  1939. def half_gcdex(f, g):
  1940. return f.ring.dmp_half_gcdex(f, g)
  1941. def gcdex(f, g):
  1942. return f.ring.dmp_gcdex(f, g)
  1943. def subresultants(f, g):
  1944. return f.ring.dmp_subresultants(f, g)
  1945. def resultant(f, g):
  1946. return f.ring.dmp_resultant(f, g)
  1947. def discriminant(f):
  1948. return f.ring.dmp_discriminant(f)
  1949. def decompose(f):
  1950. if f.ring.is_univariate:
  1951. return f.ring.dup_decompose(f)
  1952. else:
  1953. raise MultivariatePolynomialError("polynomial decomposition")
  1954. def shift(f, a):
  1955. if f.ring.is_univariate:
  1956. return f.ring.dup_shift(f, a)
  1957. else:
  1958. raise MultivariatePolynomialError("polynomial shift")
  1959. def sturm(f):
  1960. if f.ring.is_univariate:
  1961. return f.ring.dup_sturm(f)
  1962. else:
  1963. raise MultivariatePolynomialError("sturm sequence")
  1964. def gff_list(f):
  1965. return f.ring.dmp_gff_list(f)
  1966. def sqf_norm(f):
  1967. return f.ring.dmp_sqf_norm(f)
  1968. def sqf_part(f):
  1969. return f.ring.dmp_sqf_part(f)
  1970. def sqf_list(f, all=False):
  1971. return f.ring.dmp_sqf_list(f, all=all)
  1972. def factor_list(f):
  1973. return f.ring.dmp_factor_list(f)