"""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.expr import Expr
from sympy.core.symbol import Symbol, symbols as _symbols
from sympy.core.numbers import igcd, oo
from sympy.core.sympify import CantSympify, sympify
from sympy.core.compatibility import is_sequence, reduce, string_types, range
from sympy.ntheory.multinomial import multinomial_coefficients
from sympy.polys.monomials import MonomialOps
from sympy.polys.orderings import lex
from sympy.polys.heuristicgcd import heugcd
from sympy.polys.compatibility import IPolys
from sympy.polys.polyutils import expr_from_dict, _dict_reorder, _parallel_dict_from_expr
from sympy.polys.polyerrors import CoercionFailed, GeneratorsError, GeneratorsNeeded, ExactQuotientFailed, MultivariatePolynomialError
from sympy.polys.domains.domainelement import DomainElement
from sympy.polys.domains.polynomialring import PolynomialRing
from sympy.polys.polyoptions import Domain as DomainOpt, Order as OrderOpt, build_options
from sympy.polys.densebasic import dmp_to_dict, dmp_from_dict
from sympy.polys.constructor import construct_domain
from sympy.printing.defaults import DefaultPrinting
from sympy.utilities import public
from sympy.utilities.magic import pollute
@public
[docs]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:`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
@public
[docs]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:`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)
@public
[docs]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:`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
@public
[docs]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 not symbols:
raise GeneratorsNeeded("generators weren't specified")
if isinstance(symbols, string_types):
return _symbols(symbols, seq=True)
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 = hash((cls.__name__, symbols, ngens, domain, order))
obj = _ring_cache.get(_hash)
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 = _hash
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)]
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()
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] = 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 self is other
def __ne__(self, other):
return self is not 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:
i = 0
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.has_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, p1.ring.dtype):
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.__eq__(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_add = precedence["Add"]
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):
return f.ring.dmp_sqf_p(f)
@property
def is_irreducible(f):
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
p = ring.zero
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
p = ring.zero
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
p = ring.zero
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.has_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
domain = ring.domain
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
order = ring.order
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
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
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:
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.has_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.has_Field:
fc, f = f.primitive()
gc, g = g.primitive()
c = domain.lcm(fc, gc)
h = (f*g).quo(f.gcd(g))
if not domain.has_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.has_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] = 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)
