Source code for sympy.polys.rings

"""Sparse polynomial rings. """

from __future__ import print_function, division

from operator import add, mul, lt, le, gt, ge
from types import GeneratorType

from sympy.core.compatibility import is_sequence, reduce, string_types, range
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.magic import pollute

[docs]@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:, optional `Order` 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
[docs]@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:, optional `Order` 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)
[docs]@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:, optional `Order` 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 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
[docs]@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.core import symbols >>> from sympy.polys.rings import sring >>> from sympy.polys.domains import ZZ >>> from sympy.polys.orderings import lex >>> 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: # NOTE: this is inefficient because construct_domain() automatically # performs conversion to the target domain. It shouldn't do this. coeffs = sum([ list(rep.values()) for rep in reps ], []) opt.domain, _ = construct_domain(coeffs, opt=opt) _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, string_types): return _symbols(symbols, seq=True) if symbols else () elif isinstance(symbols, Expr): return (symbols,) elif is_sequence(symbols): if all(isinstance(s, string_types) 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 = {}
[docs]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 = lambda f: max(f) 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)
[docs] 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, string_types): 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): domain_new = self.domain_new poly = self.zero for monom, coeff in element.items(): coeff = domain_new(coeff) if coeff: poly[monom] = coeff return poly def from_terms(self, element): return self.from_dict(dict(element)) 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))) elif expr.is_Pow and expr.exp.is_Integer and expr.exp >= 0: return _rebuild(expr.base)**int(expr.exp) else: return 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)
[docs] 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, string_types): 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
[docs] 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
[docs] 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
[docs] 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
[docs] 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))
[docs] 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
[docs] 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))
[docs]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
[docs] 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) else: return new_ring.from_dict(self) 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
[docs] 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
[docs] 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 else: 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("can't 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("can't 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("can't 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 not mon 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(): positive = ring.domain.is_positive(coeff) sign = " + " if positive else " - " sexpvs.append(sign) if expv == zm: scoeff = printer._print(coeff) if scoeff.startswith("-"): scoeff = scoeff[1:] else: if not positive: coeff = -coeff if coeff != 1: 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 actuall 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 = list(multinomial_coefficients(len(self), n).items()) monomial_mulpow = self.ring.monomial_mulpow zero_monom = self.ring.zero_monom terms = list(self.iterterms()) 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 else: del poly[monom] return poly
[docs] 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__ = __div__ = __truediv__ __rfloordiv__ = __rdiv__ = __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
[docs] 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 any(not f for f in 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 any(not g for g in 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
[docs] 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() ])
[docs] 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()))))
[docs] 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() ])
[docs] 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()))))
[docs] 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)
[docs] 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)
[docs] 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
[docs] 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))
[docs] 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)
[docs] def coeffs(self, order=None): """Ordered list of polynomial coefficients. Parameters ========== order : :class:`Order` 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) ]
[docs] def monoms(self, order=None): """Ordered list of polynomial monomials. Parameters ========== order : :class:`Order` 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) ]
[docs] def terms(self, order=None): """Ordered list of polynomial terms. Parameters ========== order : :class:`Order` 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)
[docs] def itercoeffs(self): """Iterator over coefficients of a polynomial. """ return iter(self.values())
[docs] def itermonoms(self): """Iterator over monomials of a polynomial. """ return iter(self.keys())
[docs] def iterterms(self): """Iterator over terms of a polynomial. """ return iter(self.items())
[docs] def listcoeffs(self): """Unordered list of polynomial coefficients. """ return list(self.values())
[docs] def listmonoms(self): """Unordered list of polynomial monomials. """ return list(self.keys())
[docs] def listterms(self): """Unordered list of polynomial terms. """ return list(self.items())
[docs] 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
[docs] 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
[docs] def primitive(f): """Returns content and a primitive polynomial. """ cont = f.content() return cont, f.quo_ground(cont)
[docs] 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
[docs] 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) if q.is_negative: p, q = -p, -q 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_neg = p.is_negative q_neg = q.is_negative if p_neg and q_neg: p, q = -p, -q elif p_neg: cp, p = -cp, -p elif q_neg: cp, q = -cp, -q p = p.mul_ground(cp) q = q.mul_ground(cq) return p, q
[docs] 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)