Coverage for sympy/polys/polytools.py : 25%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
"""User-friendly public interface to polynomial functions. """
from __future__ import print_function, division
from sympy.core import ( S, Basic, Expr, I, Integer, Add, Mul, Dummy, Tuple )
from sympy.core.mul import _keep_coeff from sympy.core.symbol import Symbol from sympy.core.basic import preorder_traversal from sympy.core.relational import Relational from sympy.core.sympify import sympify from sympy.core.decorators import _sympifyit from sympy.core.function import Derivative
from sympy.logic.boolalg import BooleanAtom
from sympy.polys.polyclasses import DMP
from sympy.polys.polyutils import ( basic_from_dict, _sort_gens, _unify_gens, _dict_reorder, _dict_from_expr, _parallel_dict_from_expr, )
from sympy.polys.rationaltools import together from sympy.polys.rootisolation import dup_isolate_real_roots_list from sympy.polys.groebnertools import groebner as _groebner from sympy.polys.fglmtools import matrix_fglm from sympy.polys.monomials import Monomial from sympy.polys.orderings import monomial_key
from sympy.polys.polyerrors import ( OperationNotSupported, DomainError, CoercionFailed, UnificationFailed, GeneratorsNeeded, PolynomialError, MultivariatePolynomialError, ExactQuotientFailed, PolificationFailed, ComputationFailed, GeneratorsError, )
from sympy.utilities import group, sift, public
import sympy.polys import mpmath from mpmath.libmp.libhyper import NoConvergence
from sympy.polys.domains import FF, QQ, ZZ from sympy.polys.constructor import construct_domain
from sympy.polys import polyoptions as options
from sympy.core.compatibility import iterable, range
@public class Poly(Expr): """Generic class for representing polynomial expressions. """
__slots__ = ['rep', 'gens']
is_commutative = True is_Poly = True
def __new__(cls, rep, *gens, **args): """Create a new polynomial instance out of something useful. """
raise NotImplementedError("'order' keyword is not implemented yet")
if isinstance(rep, dict): return cls._from_dict(rep, opt) else: return cls._from_list(list(rep), opt) else:
else:
@classmethod def new(cls, rep, *gens): """Construct :class:`Poly` instance from raw representation. """ raise PolynomialError( "invalid polynomial representation: %s" % rep) raise PolynomialError("invalid arguments: %s, %s" % (rep, gens))
@classmethod def from_dict(cls, rep, *gens, **args): """Construct a polynomial from a ``dict``. """
@classmethod def from_list(cls, rep, *gens, **args): """Construct a polynomial from a ``list``. """ opt = options.build_options(gens, args) return cls._from_list(rep, opt)
@classmethod def from_poly(cls, rep, *gens, **args): """Construct a polynomial from a polynomial. """ opt = options.build_options(gens, args) return cls._from_poly(rep, opt)
@classmethod def from_expr(cls, rep, *gens, **args): """Construct a polynomial from an expression. """ opt = options.build_options(gens, args) return cls._from_expr(rep, opt)
@classmethod def _from_dict(cls, rep, opt): """Construct a polynomial from a ``dict``. """
raise GeneratorsNeeded( "can't initialize from 'dict' without generators")
else:
@classmethod def _from_list(cls, rep, opt): """Construct a polynomial from a ``list``. """ gens = opt.gens
if not gens: raise GeneratorsNeeded( "can't initialize from 'list' without generators") elif len(gens) != 1: raise MultivariatePolynomialError( "'list' representation not supported")
level = len(gens) - 1 domain = opt.domain
if domain is None: domain, rep = construct_domain(rep, opt=opt) else: rep = list(map(domain.convert, rep))
return cls.new(DMP.from_list(rep, level, domain), *gens)
@classmethod def _from_poly(cls, rep, opt): """Construct a polynomial from a polynomial. """
else: rep = rep.reorder(*gens)
rep = rep.set_domain(domain)
@classmethod def _from_expr(cls, rep, opt): """Construct a polynomial from an expression. """
def _hashable_content(self): """Allow SymPy to hash Poly instances. """
def __hash__(self):
@property def free_symbols(self): """ Free symbols of a polynomial expression.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 1).free_symbols set([x]) >>> Poly(x**2 + y).free_symbols set([x, y]) >>> Poly(x**2 + y, x).free_symbols set([x, y])
""" symbols = set([])
for gen in self.gens: symbols |= gen.free_symbols
return symbols | self.free_symbols_in_domain
@property def free_symbols_in_domain(self): """ Free symbols of the domain of ``self``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 1).free_symbols_in_domain set() >>> Poly(x**2 + y).free_symbols_in_domain set() >>> Poly(x**2 + y, x).free_symbols_in_domain set([y])
""" domain, symbols = self.rep.dom, set()
if domain.is_Composite: for gen in domain.symbols: symbols |= gen.free_symbols elif domain.is_EX: for coeff in self.coeffs(): symbols |= coeff.free_symbols
return symbols
@property def args(self): """ Don't mess up with the core.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).args (x**2 + 1,)
""" return (self.as_expr(),)
@property def gen(self): """ Return the principal generator.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).gen x
"""
@property def domain(self): """Get the ground domain of ``self``. """
@property def zero(self): """Return zero polynomial with ``self``'s properties. """ return self.new(self.rep.zero(self.rep.lev, self.rep.dom), *self.gens)
@property def one(self): """Return one polynomial with ``self``'s properties. """ return self.new(self.rep.one(self.rep.lev, self.rep.dom), *self.gens)
@property def unit(self): """Return unit polynomial with ``self``'s properties. """ return self.new(self.rep.unit(self.rep.lev, self.rep.dom), *self.gens)
def unify(f, g): """ Make ``f`` and ``g`` belong to the same domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f, g = Poly(x/2 + 1), Poly(2*x + 1)
>>> f Poly(1/2*x + 1, x, domain='QQ') >>> g Poly(2*x + 1, x, domain='ZZ')
>>> F, G = f.unify(g)
>>> F Poly(1/2*x + 1, x, domain='QQ') >>> G Poly(2*x + 1, x, domain='QQ')
"""
def _unify(f, g):
try: return f.rep.dom, f.per, f.rep, f.rep.per(f.rep.dom.from_sympy(g)) except CoercionFailed: raise UnificationFailed("can't unify %s with %s" % (f, g))
f_monoms, f_coeffs = _dict_reorder( f.rep.to_dict(), f.gens, gens)
if f.rep.dom != dom: f_coeffs = [dom.convert(c, f.rep.dom) for c in f_coeffs]
F = DMP(dict(list(zip(f_monoms, f_coeffs))), dom, lev) else:
g_monoms, g_coeffs = _dict_reorder( g.rep.to_dict(), g.gens, gens)
if g.rep.dom != dom: g_coeffs = [dom.convert(c, g.rep.dom) for c in g_coeffs]
G = DMP(dict(list(zip(g_monoms, g_coeffs))), dom, lev) else: else: raise UnificationFailed("can't unify %s with %s" % (f, g))
def per(f, rep, gens=None, remove=None): """ Create a Poly out of the given representation.
Examples ========
>>> from sympy import Poly, ZZ >>> from sympy.abc import x, y
>>> from sympy.polys.polyclasses import DMP
>>> a = Poly(x**2 + 1)
>>> a.per(DMP([ZZ(1), ZZ(1)], ZZ), gens=[y]) Poly(y + 1, y, domain='ZZ')
"""
def set_domain(f, domain): """Set the ground domain of ``f``. """ opt = options.build_options(f.gens, {'domain': domain}) return f.per(f.rep.convert(opt.domain))
def get_domain(f): """Get the ground domain of ``f``. """
def set_modulus(f, modulus): """ Set the modulus of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(5*x**2 + 2*x - 1, x).set_modulus(2) Poly(x**2 + 1, x, modulus=2)
""" modulus = options.Modulus.preprocess(modulus) return f.set_domain(FF(modulus))
def get_modulus(f): """ Get the modulus of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, modulus=2).get_modulus() 2
""" domain = f.get_domain()
if domain.is_FiniteField: return Integer(domain.characteristic()) else: raise PolynomialError("not a polynomial over a Galois field")
def _eval_subs(f, old, new): """Internal implementation of :func:`subs`. """ if old in f.gens: if new.is_number: return f.eval(old, new) else: try: return f.replace(old, new) except PolynomialError: pass
return f.as_expr().subs(old, new)
def exclude(f): """ Remove unnecessary generators from ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import a, b, c, d, x
>>> Poly(a + x, a, b, c, d, x).exclude() Poly(a + x, a, x, domain='ZZ')
""" J, new = f.rep.exclude() gens = []
for j in range(len(f.gens)): if j not in J: gens.append(f.gens[j])
return f.per(new, gens=gens)
def replace(f, x, y=None): """ Replace ``x`` with ``y`` in generators list.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 1, x).replace(x, y) Poly(y**2 + 1, y, domain='ZZ')
""" if y is None: if f.is_univariate: x, y = f.gen, x else: raise PolynomialError( "syntax supported only in univariate case")
if x == y: return f
if x in f.gens and y not in f.gens: dom = f.get_domain()
if not dom.is_Composite or y not in dom.symbols: gens = list(f.gens) gens[gens.index(x)] = y return f.per(f.rep, gens=gens)
raise PolynomialError("can't replace %s with %s in %s" % (x, y, f))
def reorder(f, *gens, **args): """ Efficiently apply new order of generators.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + x*y**2, x, y).reorder(y, x) Poly(y**2*x + x**2, y, x, domain='ZZ')
"""
elif set(f.gens) != set(gens): raise PolynomialError( "generators list can differ only up to order of elements")
def ltrim(f, gen): """ Remove dummy generators from the "left" of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y, z
>>> Poly(y**2 + y*z**2, x, y, z).ltrim(y) Poly(y**2 + y*z**2, y, z, domain='ZZ')
"""
else: raise PolynomialError("can't left trim %s" % f)
def has_only_gens(f, *gens): """ Return ``True`` if ``Poly(f, *gens)`` retains ground domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y, z
>>> Poly(x*y + 1, x, y, z).has_only_gens(x, y) True >>> Poly(x*y + z, x, y, z).has_only_gens(x, y) False
""" indices = set([])
for gen in gens: try: index = f.gens.index(gen) except ValueError: raise GeneratorsError( "%s doesn't have %s as generator" % (f, gen)) else: indices.add(index)
for monom in f.monoms(): for i, elt in enumerate(monom): if i not in indices and elt: return False
return True
def to_ring(f): """ Make the ground domain a ring.
Examples ========
>>> from sympy import Poly, QQ >>> from sympy.abc import x
>>> Poly(x**2 + 1, domain=QQ).to_ring() Poly(x**2 + 1, x, domain='ZZ')
""" if hasattr(f.rep, 'to_ring'): result = f.rep.to_ring() else: # pragma: no cover raise OperationNotSupported(f, 'to_ring')
return f.per(result)
def to_field(f): """ Make the ground domain a field.
Examples ========
>>> from sympy import Poly, ZZ >>> from sympy.abc import x
>>> Poly(x**2 + 1, x, domain=ZZ).to_field() Poly(x**2 + 1, x, domain='QQ')
""" else: # pragma: no cover raise OperationNotSupported(f, 'to_field')
def to_exact(f): """ Make the ground domain exact.
Examples ========
>>> from sympy import Poly, RR >>> from sympy.abc import x
>>> Poly(x**2 + 1.0, x, domain=RR).to_exact() Poly(x**2 + 1, x, domain='QQ')
""" if hasattr(f.rep, 'to_exact'): result = f.rep.to_exact() else: # pragma: no cover raise OperationNotSupported(f, 'to_exact')
return f.per(result)
def retract(f, field=None): """ Recalculate the ground domain of a polynomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = Poly(x**2 + 1, x, domain='QQ[y]') >>> f Poly(x**2 + 1, x, domain='QQ[y]')
>>> f.retract() Poly(x**2 + 1, x, domain='ZZ') >>> f.retract(field=True) Poly(x**2 + 1, x, domain='QQ')
""" field=field, composite=f.domain.is_Composite or None)
def slice(f, x, m, n=None): """Take a continuous subsequence of terms of ``f``. """ if n is None: j, m, n = 0, x, m else: j = f._gen_to_level(x)
m, n = int(m), int(n)
if hasattr(f.rep, 'slice'): result = f.rep.slice(m, n, j) else: # pragma: no cover raise OperationNotSupported(f, 'slice')
return f.per(result)
def coeffs(f, order=None): """ Returns all non-zero coefficients from ``f`` in lex order.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 + 2*x + 3, x).coeffs() [1, 2, 3]
See Also ======== all_coeffs coeff_monomial nth
""" return [f.rep.dom.to_sympy(c) for c in f.rep.coeffs(order=order)]
def monoms(f, order=None): """ Returns all non-zero monomials from ``f`` in lex order.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).monoms() [(2, 0), (1, 2), (1, 1), (0, 1)]
See Also ======== all_monoms
"""
def terms(f, order=None): """ Returns all non-zero terms from ``f`` in lex order.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).terms() [((2, 0), 1), ((1, 2), 2), ((1, 1), 1), ((0, 1), 3)]
See Also ======== all_terms
"""
def all_coeffs(f): """ Returns all coefficients from a univariate polynomial ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 + 2*x - 1, x).all_coeffs() [1, 0, 2, -1]
"""
def all_monoms(f): """ Returns all monomials from a univariate polynomial ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 + 2*x - 1, x).all_monoms() [(3,), (2,), (1,), (0,)]
See Also ======== all_terms
""" return f.rep.all_monoms()
def all_terms(f): """ Returns all terms from a univariate polynomial ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 + 2*x - 1, x).all_terms() [((3,), 1), ((2,), 0), ((1,), 2), ((0,), -1)]
""" return [(m, f.rep.dom.to_sympy(c)) for m, c in f.rep.all_terms()]
def termwise(f, func, *gens, **args): """ Apply a function to all terms of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> def func(k, coeff): ... k = k[0] ... return coeff//10**(2-k)
>>> Poly(x**2 + 20*x + 400).termwise(func) Poly(x**2 + 2*x + 4, x, domain='ZZ')
"""
monom, coeff = result else:
else: raise PolynomialError( "%s monomial was generated twice" % monom)
def length(f): """ Returns the number of non-zero terms in ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 2*x - 1).length() 3
"""
def as_dict(f, native=False, zero=False): """ Switch to a ``dict`` representation.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 2*x*y**2 - y, x, y).as_dict() {(0, 1): -1, (1, 2): 2, (2, 0): 1}
""" else:
def as_list(f, native=False): """Switch to a ``list`` representation. """ if native: return f.rep.to_list() else: return f.rep.to_sympy_list()
def as_expr(f, *gens): """ Convert a Poly instance to an Expr instance.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> f = Poly(x**2 + 2*x*y**2 - y, x, y)
>>> f.as_expr() x**2 + 2*x*y**2 - y >>> f.as_expr({x: 5}) 10*y**2 - y + 25 >>> f.as_expr(5, 6) 379
"""
except ValueError: raise GeneratorsError( "%s doesn't have %s as generator" % (f, gen)) else:
def lift(f): """ Convert algebraic coefficients to rationals.
Examples ========
>>> from sympy import Poly, I >>> from sympy.abc import x
>>> Poly(x**2 + I*x + 1, x, extension=I).lift() Poly(x**4 + 3*x**2 + 1, x, domain='QQ')
""" if hasattr(f.rep, 'lift'): result = f.rep.lift() else: # pragma: no cover raise OperationNotSupported(f, 'lift')
return f.per(result)
def deflate(f): """ Reduce degree of ``f`` by mapping ``x_i**m`` to ``y_i``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**6*y**2 + x**3 + 1, x, y).deflate() ((3, 2), Poly(x**2*y + x + 1, x, y, domain='ZZ'))
""" if hasattr(f.rep, 'deflate'): J, result = f.rep.deflate() else: # pragma: no cover raise OperationNotSupported(f, 'deflate')
return J, f.per(result)
def inject(f, front=False): """ Inject ground domain generators into ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> f = Poly(x**2*y + x*y**3 + x*y + 1, x)
>>> f.inject() Poly(x**2*y + x*y**3 + x*y + 1, x, y, domain='ZZ') >>> f.inject(front=True) Poly(y**3*x + y*x**2 + y*x + 1, y, x, domain='ZZ')
"""
return f raise DomainError("can't inject generators over %s" % dom)
else: # pragma: no cover raise OperationNotSupported(f, 'inject')
gens = dom.symbols + f.gens else:
def eject(f, *gens): """ Eject selected generators into the ground domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> f = Poly(x**2*y + x*y**3 + x*y + 1, x, y)
>>> f.eject(x) Poly(x*y**3 + (x**2 + x)*y + 1, y, domain='ZZ[x]') >>> f.eject(y) Poly(y*x**2 + (y**3 + y)*x + 1, x, domain='ZZ[y]')
"""
raise DomainError("can't eject generators over %s" % dom)
_gens, front = f.gens[k:], True else: raise NotImplementedError( "can only eject front or back generators")
else: # pragma: no cover raise OperationNotSupported(f, 'eject')
def terms_gcd(f): """ Remove GCD of terms from the polynomial ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**6*y**2 + x**3*y, x, y).terms_gcd() ((3, 1), Poly(x**3*y + 1, x, y, domain='ZZ'))
""" else: # pragma: no cover raise OperationNotSupported(f, 'terms_gcd')
def add_ground(f, coeff): """ Add an element of the ground domain to ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x + 1).add_ground(2) Poly(x + 3, x, domain='ZZ')
""" if hasattr(f.rep, 'add_ground'): result = f.rep.add_ground(coeff) else: # pragma: no cover raise OperationNotSupported(f, 'add_ground')
return f.per(result)
def sub_ground(f, coeff): """ Subtract an element of the ground domain from ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x + 1).sub_ground(2) Poly(x - 1, x, domain='ZZ')
""" if hasattr(f.rep, 'sub_ground'): result = f.rep.sub_ground(coeff) else: # pragma: no cover raise OperationNotSupported(f, 'sub_ground')
return f.per(result)
def mul_ground(f, coeff): """ Multiply ``f`` by a an element of the ground domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x + 1).mul_ground(2) Poly(2*x + 2, x, domain='ZZ')
""" if hasattr(f.rep, 'mul_ground'): result = f.rep.mul_ground(coeff) else: # pragma: no cover raise OperationNotSupported(f, 'mul_ground')
return f.per(result)
def quo_ground(f, coeff): """ Quotient of ``f`` by a an element of the ground domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x + 4).quo_ground(2) Poly(x + 2, x, domain='ZZ')
>>> Poly(2*x + 3).quo_ground(2) Poly(x + 1, x, domain='ZZ')
""" if hasattr(f.rep, 'quo_ground'): result = f.rep.quo_ground(coeff) else: # pragma: no cover raise OperationNotSupported(f, 'quo_ground')
return f.per(result)
def exquo_ground(f, coeff): """ Exact quotient of ``f`` by a an element of the ground domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x + 4).exquo_ground(2) Poly(x + 2, x, domain='ZZ')
>>> Poly(2*x + 3).exquo_ground(2) Traceback (most recent call last): ... ExactQuotientFailed: 2 does not divide 3 in ZZ
""" if hasattr(f.rep, 'exquo_ground'): result = f.rep.exquo_ground(coeff) else: # pragma: no cover raise OperationNotSupported(f, 'exquo_ground')
return f.per(result)
def abs(f): """ Make all coefficients in ``f`` positive.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).abs() Poly(x**2 + 1, x, domain='ZZ')
""" if hasattr(f.rep, 'abs'): result = f.rep.abs() else: # pragma: no cover raise OperationNotSupported(f, 'abs')
return f.per(result)
def neg(f): """ Negate all coefficients in ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).neg() Poly(-x**2 + 1, x, domain='ZZ')
>>> -Poly(x**2 - 1, x) Poly(-x**2 + 1, x, domain='ZZ')
""" if hasattr(f.rep, 'neg'): result = f.rep.neg() else: # pragma: no cover raise OperationNotSupported(f, 'neg')
return f.per(result)
def add(f, g): """ Add two polynomials ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).add(Poly(x - 2, x)) Poly(x**2 + x - 1, x, domain='ZZ')
>>> Poly(x**2 + 1, x) + Poly(x - 2, x) Poly(x**2 + x - 1, x, domain='ZZ')
""" g = sympify(g)
if not g.is_Poly: return f.add_ground(g)
_, per, F, G = f._unify(g)
if hasattr(f.rep, 'add'): result = F.add(G) else: # pragma: no cover raise OperationNotSupported(f, 'add')
return per(result)
def sub(f, g): """ Subtract two polynomials ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).sub(Poly(x - 2, x)) Poly(x**2 - x + 3, x, domain='ZZ')
>>> Poly(x**2 + 1, x) - Poly(x - 2, x) Poly(x**2 - x + 3, x, domain='ZZ')
"""
return f.sub_ground(g)
else: # pragma: no cover raise OperationNotSupported(f, 'sub')
def mul(f, g): """ Multiply two polynomials ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).mul(Poly(x - 2, x)) Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ')
>>> Poly(x**2 + 1, x)*Poly(x - 2, x) Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ')
""" g = sympify(g)
if not g.is_Poly: return f.mul_ground(g)
_, per, F, G = f._unify(g)
if hasattr(f.rep, 'mul'): result = F.mul(G) else: # pragma: no cover raise OperationNotSupported(f, 'mul')
return per(result)
def sqr(f): """ Square a polynomial ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x - 2, x).sqr() Poly(x**2 - 4*x + 4, x, domain='ZZ')
>>> Poly(x - 2, x)**2 Poly(x**2 - 4*x + 4, x, domain='ZZ')
""" if hasattr(f.rep, 'sqr'): result = f.rep.sqr() else: # pragma: no cover raise OperationNotSupported(f, 'sqr')
return f.per(result)
def pow(f, n): """ Raise ``f`` to a non-negative power ``n``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x - 2, x).pow(3) Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ')
>>> Poly(x - 2, x)**3 Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ')
""" n = int(n)
if hasattr(f.rep, 'pow'): result = f.rep.pow(n) else: # pragma: no cover raise OperationNotSupported(f, 'pow')
return f.per(result)
def pdiv(f, g): """ Polynomial pseudo-division of ``f`` by ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).pdiv(Poly(2*x - 4, x)) (Poly(2*x + 4, x, domain='ZZ'), Poly(20, x, domain='ZZ'))
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'pdiv'): q, r = F.pdiv(G) else: # pragma: no cover raise OperationNotSupported(f, 'pdiv')
return per(q), per(r)
def prem(f, g): """ Polynomial pseudo-remainder of ``f`` by ``g``.
Caveat: The function prem(f, g, x) can be safely used to compute in Z[x] _only_ subresultant polynomial remainder sequences (prs's).
To safely compute Euclidean and Sturmian prs's in Z[x] employ anyone of the corresponding functions found in the module sympy.polys.subresultants_qq_zz. The functions in the module with suffix _pg compute prs's in Z[x] employing rem(f, g, x), whereas the functions with suffix _amv compute prs's in Z[x] employing rem_z(f, g, x).
The function rem_z(f, g, x) differs from prem(f, g, x) in that to compute the remainder polynomials in Z[x] it premultiplies the divident times the absolute value of the leading coefficient of the divisor raised to the power degree(f, x) - degree(g, x) + 1.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).prem(Poly(2*x - 4, x)) Poly(20, x, domain='ZZ')
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'prem'): result = F.prem(G) else: # pragma: no cover raise OperationNotSupported(f, 'prem')
return per(result)
def pquo(f, g): """ Polynomial pseudo-quotient of ``f`` by ``g``.
See the Caveat note in the function prem(f, g).
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).pquo(Poly(2*x - 4, x)) Poly(2*x + 4, x, domain='ZZ')
>>> Poly(x**2 - 1, x).pquo(Poly(2*x - 2, x)) Poly(2*x + 2, x, domain='ZZ')
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'pquo'): result = F.pquo(G) else: # pragma: no cover raise OperationNotSupported(f, 'pquo')
return per(result)
def pexquo(f, g): """ Polynomial exact pseudo-quotient of ``f`` by ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).pexquo(Poly(2*x - 2, x)) Poly(2*x + 2, x, domain='ZZ')
>>> Poly(x**2 + 1, x).pexquo(Poly(2*x - 4, x)) Traceback (most recent call last): ... ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'pexquo'): try: result = F.pexquo(G) except ExactQuotientFailed as exc: raise exc.new(f.as_expr(), g.as_expr()) else: # pragma: no cover raise OperationNotSupported(f, 'pexquo')
return per(result)
def div(f, g, auto=True): """ Polynomial division with remainder of ``f`` by ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x)) (Poly(1/2*x + 1, x, domain='QQ'), Poly(5, x, domain='QQ'))
>>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x), auto=False) (Poly(0, x, domain='ZZ'), Poly(x**2 + 1, x, domain='ZZ'))
""" dom, per, F, G = f._unify(g) retract = False
if auto and dom.has_Ring and not dom.has_Field: F, G = F.to_field(), G.to_field() retract = True
if hasattr(f.rep, 'div'): q, r = F.div(G) else: # pragma: no cover raise OperationNotSupported(f, 'div')
if retract: try: Q, R = q.to_ring(), r.to_ring() except CoercionFailed: pass else: q, r = Q, R
return per(q), per(r)
def rem(f, g, auto=True): """ Computes the polynomial remainder of ``f`` by ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x)) Poly(5, x, domain='ZZ')
>>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x), auto=False) Poly(x**2 + 1, x, domain='ZZ')
""" dom, per, F, G = f._unify(g) retract = False
if auto and dom.has_Ring and not dom.has_Field: F, G = F.to_field(), G.to_field() retract = True
if hasattr(f.rep, 'rem'): r = F.rem(G) else: # pragma: no cover raise OperationNotSupported(f, 'rem')
if retract: try: r = r.to_ring() except CoercionFailed: pass
return per(r)
def quo(f, g, auto=True): """ Computes polynomial quotient of ``f`` by ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).quo(Poly(2*x - 4, x)) Poly(1/2*x + 1, x, domain='QQ')
>>> Poly(x**2 - 1, x).quo(Poly(x - 1, x)) Poly(x + 1, x, domain='ZZ')
""" dom, per, F, G = f._unify(g) retract = False
if auto and dom.has_Ring and not dom.has_Field: F, G = F.to_field(), G.to_field() retract = True
if hasattr(f.rep, 'quo'): q = F.quo(G) else: # pragma: no cover raise OperationNotSupported(f, 'quo')
if retract: try: q = q.to_ring() except CoercionFailed: pass
return per(q)
def exquo(f, g, auto=True): """ Computes polynomial exact quotient of ``f`` by ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).exquo(Poly(x - 1, x)) Poly(x + 1, x, domain='ZZ')
>>> Poly(x**2 + 1, x).exquo(Poly(2*x - 4, x)) Traceback (most recent call last): ... ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
""" dom, per, F, G = f._unify(g) retract = False
if auto and dom.has_Ring and not dom.has_Field: F, G = F.to_field(), G.to_field() retract = True
if hasattr(f.rep, 'exquo'): try: q = F.exquo(G) except ExactQuotientFailed as exc: raise exc.new(f.as_expr(), g.as_expr()) else: # pragma: no cover raise OperationNotSupported(f, 'exquo')
if retract: try: q = q.to_ring() except CoercionFailed: pass
return per(q)
def _gen_to_level(f, gen): """Returns level associated with the given generator. """
else: else: raise PolynomialError("-%s <= gen < %s expected, got %s" % (length, length, gen)) else: except ValueError: raise PolynomialError( "a valid generator expected, got %s" % gen)
def degree(f, gen=0): """ Returns degree of ``f`` in ``x_j``.
The degree of 0 is negative infinity.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + y*x + 1, x, y).degree() 2 >>> Poly(x**2 + y*x + y, x, y).degree(y) 1 >>> Poly(0, x).degree() -oo
"""
else: # pragma: no cover raise OperationNotSupported(f, 'degree')
def degree_list(f): """ Returns a list of degrees of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + y*x + 1, x, y).degree_list() (2, 1)
""" else: # pragma: no cover raise OperationNotSupported(f, 'degree_list')
def total_degree(f): """ Returns the total degree of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + y*x + 1, x, y).total_degree() 2 >>> Poly(x + y**5, x, y).total_degree() 5
""" else: # pragma: no cover raise OperationNotSupported(f, 'total_degree')
def homogenize(f, s): """ Returns the homogeneous polynomial of ``f``.
A homogeneous polynomial is a polynomial whose all monomials with non-zero coefficients have the same total degree. If you only want to check if a polynomial is homogeneous, then use :func:`Poly.is_homogeneous`. If you want not only to check if a polynomial is homogeneous but also compute its homogeneous order, then use :func:`Poly.homogeneous_order`.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y, z
>>> f = Poly(x**5 + 2*x**2*y**2 + 9*x*y**3) >>> f.homogenize(z) Poly(x**5 + 2*x**2*y**2*z + 9*x*y**3*z, x, y, z, domain='ZZ')
""" if not isinstance(s, Symbol): raise TypeError("``Symbol`` expected, got %s" % type(s)) if s in f.gens: i = f.gens.index(s) gens = f.gens else: i = len(f.gens) gens = f.gens + (s,) if hasattr(f.rep, 'homogenize'): return f.per(f.rep.homogenize(i), gens=gens) raise OperationNotSupported(f, 'homogeneous_order')
def homogeneous_order(f): """ Returns the homogeneous order of ``f``.
A homogeneous polynomial is a polynomial whose all monomials with non-zero coefficients have the same total degree. This degree is the homogeneous order of ``f``. If you only want to check if a polynomial is homogeneous, then use :func:`Poly.is_homogeneous`.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> f = Poly(x**5 + 2*x**3*y**2 + 9*x*y**4) >>> f.homogeneous_order() 5
""" if hasattr(f.rep, 'homogeneous_order'): return f.rep.homogeneous_order() else: # pragma: no cover raise OperationNotSupported(f, 'homogeneous_order')
def LC(f, order=None): """ Returns the leading coefficient of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(4*x**3 + 2*x**2 + 3*x, x).LC() 4
""" if order is not None: return f.coeffs(order)[0]
if hasattr(f.rep, 'LC'): result = f.rep.LC() else: # pragma: no cover raise OperationNotSupported(f, 'LC')
return f.rep.dom.to_sympy(result)
def TC(f): """ Returns the trailing coefficient of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 + 2*x**2 + 3*x, x).TC() 0
""" if hasattr(f.rep, 'TC'): result = f.rep.TC() else: # pragma: no cover raise OperationNotSupported(f, 'TC')
return f.rep.dom.to_sympy(result)
def EC(f, order=None): """ Returns the last non-zero coefficient of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 + 2*x**2 + 3*x, x).EC() 3
""" if hasattr(f.rep, 'coeffs'): return f.coeffs(order)[-1] else: # pragma: no cover raise OperationNotSupported(f, 'EC')
def coeff_monomial(f, monom): """ Returns the coefficient of ``monom`` in ``f`` if there, else None.
Examples ========
>>> from sympy import Poly, exp >>> from sympy.abc import x, y
>>> p = Poly(24*x*y*exp(8) + 23*x, x, y)
>>> p.coeff_monomial(x) 23 >>> p.coeff_monomial(y) 0 >>> p.coeff_monomial(x*y) 24*exp(8)
Note that ``Expr.coeff()`` behaves differently, collecting terms if possible; the Poly must be converted to an Expr to use that method, however:
>>> p.as_expr().coeff(x) 24*y*exp(8) + 23 >>> p.as_expr().coeff(y) 24*x*exp(8) >>> p.as_expr().coeff(x*y) 24*exp(8)
See Also ======== nth: more efficient query using exponents of the monomial's generators
"""
def nth(f, *N): """ Returns the ``n``-th coefficient of ``f`` where ``N`` are the exponents of the generators in the term of interest.
Examples ========
>>> from sympy import Poly, sqrt >>> from sympy.abc import x, y
>>> Poly(x**3 + 2*x**2 + 3*x, x).nth(2) 2 >>> Poly(x**3 + 2*x*y**2 + y**2, x, y).nth(1, 2) 2 >>> Poly(4*sqrt(x)*y) Poly(4*y*(sqrt(x)), y, sqrt(x), domain='ZZ') >>> _.nth(1, 1) 4
See Also ======== coeff_monomial
""" raise ValueError('exponent of each generator must be specified') else: # pragma: no cover raise OperationNotSupported(f, 'nth')
def coeff(f, x, n=1, right=False): # the semantics of coeff_monomial and Expr.coeff are different; # if someone is working with a Poly, they should be aware of the # differences and chose the method best suited for the query. # Alternatively, a pure-polys method could be written here but # at this time the ``right`` keyword would be ignored because Poly # doesn't work with non-commutatives. raise NotImplementedError( 'Either convert to Expr with `as_expr` method ' 'to use Expr\'s coeff method or else use the ' '`coeff_monomial` method of Polys.')
def LM(f, order=None): """ Returns the leading monomial of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LM() x**2*y**0
"""
def EM(f, order=None): """ Returns the last non-zero monomial of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).EM() x**0*y**1
""" return Monomial(f.monoms(order)[-1], f.gens)
def LT(f, order=None): """ Returns the leading term of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LT() (x**2*y**0, 4)
""" monom, coeff = f.terms(order)[0] return Monomial(monom, f.gens), coeff
def ET(f, order=None): """ Returns the last non-zero term of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).ET() (x**0*y**1, 3)
""" monom, coeff = f.terms(order)[-1] return Monomial(monom, f.gens), coeff
def max_norm(f): """ Returns maximum norm of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(-x**2 + 2*x - 3, x).max_norm() 3
""" if hasattr(f.rep, 'max_norm'): result = f.rep.max_norm() else: # pragma: no cover raise OperationNotSupported(f, 'max_norm')
return f.rep.dom.to_sympy(result)
def l1_norm(f): """ Returns l1 norm of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(-x**2 + 2*x - 3, x).l1_norm() 6
""" if hasattr(f.rep, 'l1_norm'): result = f.rep.l1_norm() else: # pragma: no cover raise OperationNotSupported(f, 'l1_norm')
return f.rep.dom.to_sympy(result)
def clear_denoms(self, convert=False): """ Clear denominators, but keep the ground domain.
Examples ========
>>> from sympy import Poly, S, QQ >>> from sympy.abc import x
>>> f = Poly(x/2 + S(1)/3, x, domain=QQ)
>>> f.clear_denoms() (6, Poly(3*x + 2, x, domain='QQ')) >>> f.clear_denoms(convert=True) (6, Poly(3*x + 2, x, domain='ZZ'))
"""
dom = f.rep.dom.get_ring()
else: # pragma: no cover raise OperationNotSupported(f, 'clear_denoms')
else: return coeff, f.to_ring()
def rat_clear_denoms(self, g): """ Clear denominators in a rational function ``f/g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> f = Poly(x**2/y + 1, x) >>> g = Poly(x**3 + y, x)
>>> p, q = f.rat_clear_denoms(g)
>>> p Poly(x**2 + y, x, domain='ZZ[y]') >>> q Poly(y*x**3 + y**2, x, domain='ZZ[y]')
""" f = self
dom, per, f, g = f._unify(g)
f = per(f) g = per(g)
if not (dom.has_Field and dom.has_assoc_Ring): return f, g
a, f = f.clear_denoms(convert=True) b, g = g.clear_denoms(convert=True)
f = f.mul_ground(b) g = g.mul_ground(a)
return f, g
def integrate(self, *specs, **args): """ Computes indefinite integral of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 2*x + 1, x).integrate() Poly(1/3*x**3 + x**2 + x, x, domain='QQ')
>>> Poly(x*y**2 + x, x, y).integrate((0, 1), (1, 0)) Poly(1/2*x**2*y**2 + 1/2*x**2, x, y, domain='QQ')
""" f = self
if args.get('auto', True) and f.rep.dom.has_Ring: f = f.to_field()
if hasattr(f.rep, 'integrate'): if not specs: return f.per(f.rep.integrate(m=1))
rep = f.rep
for spec in specs: if type(spec) is tuple: gen, m = spec else: gen, m = spec, 1
rep = rep.integrate(int(m), f._gen_to_level(gen))
return f.per(rep) else: # pragma: no cover raise OperationNotSupported(f, 'integrate')
def diff(f, *specs, **kwargs): """ Computes partial derivative of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + 2*x + 1, x).diff() Poly(2*x + 2, x, domain='ZZ')
>>> Poly(x*y**2 + x, x, y).diff((0, 0), (1, 1)) Poly(2*x*y, x, y, domain='ZZ')
""" if not kwargs.get('evaluate', True): return Derivative(f, *specs, **kwargs)
if hasattr(f.rep, 'diff'): if not specs: return f.per(f.rep.diff(m=1))
rep = f.rep
for spec in specs: if type(spec) is tuple: gen, m = spec else: gen, m = spec, 1
rep = rep.diff(int(m), f._gen_to_level(gen))
return f.per(rep) else: # pragma: no cover raise OperationNotSupported(f, 'diff')
_eval_derivative = diff _eval_diff = diff
def eval(self, x, a=None, auto=True): """ Evaluate ``f`` at ``a`` in the given variable.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y, z
>>> Poly(x**2 + 2*x + 3, x).eval(2) 11
>>> Poly(2*x*y + 3*x + y + 2, x, y).eval(x, 2) Poly(5*y + 8, y, domain='ZZ')
>>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z)
>>> f.eval({x: 2}) Poly(5*y + 2*z + 6, y, z, domain='ZZ') >>> f.eval({x: 2, y: 5}) Poly(2*z + 31, z, domain='ZZ') >>> f.eval({x: 2, y: 5, z: 7}) 45
>>> f.eval((2, 5)) Poly(2*z + 31, z, domain='ZZ') >>> f(2, 5) Poly(2*z + 31, z, domain='ZZ')
"""
mapping = x
for gen, value in mapping.items(): f = f.eval(gen, value)
return f values = x
if len(values) > len(f.gens): raise ValueError("too many values provided")
for gen, value in zip(f.gens, values): f = f.eval(gen, value)
return f else: else:
if not hasattr(f.rep, 'eval'): # pragma: no cover raise OperationNotSupported(f, 'eval')
except CoercionFailed: if not auto: raise DomainError("can't evaluate at %s in %s" % (a, f.rep.dom)) else: a_domain, [a] = construct_domain([a]) new_domain = f.get_domain().unify_with_symbols(a_domain, f.gens)
f = f.set_domain(new_domain) a = new_domain.convert(a, a_domain)
result = f.rep.eval(a, j)
def __call__(f, *values): """ Evaluate ``f`` at the give values.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y, z
>>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z)
>>> f(2) Poly(5*y + 2*z + 6, y, z, domain='ZZ') >>> f(2, 5) Poly(2*z + 31, z, domain='ZZ') >>> f(2, 5, 7) 45
""" return f.eval(values)
def half_gcdex(f, g, auto=True): """ Half extended Euclidean algorithm of ``f`` and ``g``.
Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15 >>> g = x**3 + x**2 - 4*x - 4
>>> Poly(f).half_gcdex(Poly(g)) (Poly(-1/5*x + 3/5, x, domain='QQ'), Poly(x + 1, x, domain='QQ'))
""" dom, per, F, G = f._unify(g)
if auto and dom.has_Ring: F, G = F.to_field(), G.to_field()
if hasattr(f.rep, 'half_gcdex'): s, h = F.half_gcdex(G) else: # pragma: no cover raise OperationNotSupported(f, 'half_gcdex')
return per(s), per(h)
def gcdex(f, g, auto=True): """ Extended Euclidean algorithm of ``f`` and ``g``.
Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15 >>> g = x**3 + x**2 - 4*x - 4
>>> Poly(f).gcdex(Poly(g)) (Poly(-1/5*x + 3/5, x, domain='QQ'), Poly(1/5*x**2 - 6/5*x + 2, x, domain='QQ'), Poly(x + 1, x, domain='QQ'))
""" dom, per, F, G = f._unify(g)
if auto and dom.has_Ring: F, G = F.to_field(), G.to_field()
if hasattr(f.rep, 'gcdex'): s, t, h = F.gcdex(G) else: # pragma: no cover raise OperationNotSupported(f, 'gcdex')
return per(s), per(t), per(h)
def invert(f, g, auto=True): """ Invert ``f`` modulo ``g`` when possible.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).invert(Poly(2*x - 1, x)) Poly(-4/3, x, domain='QQ')
>>> Poly(x**2 - 1, x).invert(Poly(x - 1, x)) Traceback (most recent call last): ... NotInvertible: zero divisor
""" dom, per, F, G = f._unify(g)
if auto and dom.has_Ring: F, G = F.to_field(), G.to_field()
if hasattr(f.rep, 'invert'): result = F.invert(G) else: # pragma: no cover raise OperationNotSupported(f, 'invert')
return per(result)
def revert(f, n): """Compute ``f**(-1)`` mod ``x**n``. """ if hasattr(f.rep, 'revert'): result = f.rep.revert(int(n)) else: # pragma: no cover raise OperationNotSupported(f, 'revert')
return f.per(result)
def subresultants(f, g): """ Computes the subresultant PRS of ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 1, x).subresultants(Poly(x**2 - 1, x)) [Poly(x**2 + 1, x, domain='ZZ'), Poly(x**2 - 1, x, domain='ZZ'), Poly(-2, x, domain='ZZ')]
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'subresultants'): result = F.subresultants(G) else: # pragma: no cover raise OperationNotSupported(f, 'subresultants')
return list(map(per, result))
def resultant(f, g, includePRS=False): """ Computes the resultant of ``f`` and ``g`` via PRS.
If includePRS=True, it includes the subresultant PRS in the result. Because the PRS is used to calculate the resultant, this is more efficient than calling :func:`subresultants` separately.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = Poly(x**2 + 1, x)
>>> f.resultant(Poly(x**2 - 1, x)) 4 >>> f.resultant(Poly(x**2 - 1, x), includePRS=True) (4, [Poly(x**2 + 1, x, domain='ZZ'), Poly(x**2 - 1, x, domain='ZZ'), Poly(-2, x, domain='ZZ')])
"""
result, R = F.resultant(G, includePRS=includePRS) else: else: # pragma: no cover raise OperationNotSupported(f, 'resultant')
return (per(result, remove=0), list(map(per, R)))
def discriminant(f): """ Computes the discriminant of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + 2*x + 3, x).discriminant() -8
""" if hasattr(f.rep, 'discriminant'): result = f.rep.discriminant() else: # pragma: no cover raise OperationNotSupported(f, 'discriminant')
return f.per(result, remove=0)
def dispersionset(f, g=None): r"""Compute the *dispersion set* of two polynomials.
For two polynomials `f(x)` and `g(x)` with `\deg f > 0` and `\deg g > 0` the dispersion set `\operatorname{J}(f, g)` is defined as:
.. math:: \operatorname{J}(f, g) & := \{a \in \mathbb{N}_0 | \gcd(f(x), g(x+a)) \neq 1\} \\ & = \{a \in \mathbb{N}_0 | \deg \gcd(f(x), g(x+a)) \geq 1\}
For a single polynomial one defines `\operatorname{J}(f) := \operatorname{J}(f, f)`.
Examples ========
>>> from sympy import poly >>> from sympy.polys.dispersion import dispersion, dispersionset >>> from sympy.abc import x
Dispersion set and dispersion of a simple polynomial:
>>> fp = poly((x - 3)*(x + 3), x) >>> sorted(dispersionset(fp)) [0, 6] >>> dispersion(fp) 6
Note that the definition of the dispersion is not symmetric:
>>> fp = poly(x**4 - 3*x**2 + 1, x) >>> gp = fp.shift(-3) >>> sorted(dispersionset(fp, gp)) [2, 3, 4] >>> dispersion(fp, gp) 4 >>> sorted(dispersionset(gp, fp)) [] >>> dispersion(gp, fp) -oo
Computing the dispersion also works over field extensions:
>>> from sympy import sqrt >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>') >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>') >>> sorted(dispersionset(fp, gp)) [2] >>> sorted(dispersionset(gp, fp)) [1, 4]
We can even perform the computations for polynomials having symbolic coefficients:
>>> from sympy.abc import a >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x) >>> sorted(dispersionset(fp)) [0, 1]
See Also ========
dispersion
References ==========
1. [ManWright94]_ 2. [Koepf98]_ 3. [Abramov71]_ 4. [Man93]_ """ from sympy.polys.dispersion import dispersionset return dispersionset(f, g)
def dispersion(f, g=None): r"""Compute the *dispersion* of polynomials.
For two polynomials `f(x)` and `g(x)` with `\deg f > 0` and `\deg g > 0` the dispersion `\operatorname{dis}(f, g)` is defined as:
.. math:: \operatorname{dis}(f, g) & := \max\{ J(f,g) \cup \{0\} \} \\ & = \max\{ \{a \in \mathbb{N} | \gcd(f(x), g(x+a)) \neq 1\} \cup \{0\} \}
and for a single polynomial `\operatorname{dis}(f) := \operatorname{dis}(f, f)`.
Examples ========
>>> from sympy import poly >>> from sympy.polys.dispersion import dispersion, dispersionset >>> from sympy.abc import x
Dispersion set and dispersion of a simple polynomial:
>>> fp = poly((x - 3)*(x + 3), x) >>> sorted(dispersionset(fp)) [0, 6] >>> dispersion(fp) 6
Note that the definition of the dispersion is not symmetric:
>>> fp = poly(x**4 - 3*x**2 + 1, x) >>> gp = fp.shift(-3) >>> sorted(dispersionset(fp, gp)) [2, 3, 4] >>> dispersion(fp, gp) 4 >>> sorted(dispersionset(gp, fp)) [] >>> dispersion(gp, fp) -oo
Computing the dispersion also works over field extensions:
>>> from sympy import sqrt >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>') >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>') >>> sorted(dispersionset(fp, gp)) [2] >>> sorted(dispersionset(gp, fp)) [1, 4]
We can even perform the computations for polynomials having symbolic coefficients:
>>> from sympy.abc import a >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x) >>> sorted(dispersionset(fp)) [0, 1]
See Also ========
dispersionset
References ==========
1. [ManWright94]_ 2. [Koepf98]_ 3. [Abramov71]_ 4. [Man93]_ """ from sympy.polys.dispersion import dispersion return dispersion(f, g)
def cofactors(f, g): """ Returns the GCD of ``f`` and ``g`` and their cofactors.
Returns polynomials ``(h, cff, cfg)`` such that ``h = gcd(f, g)``, and ``cff = quo(f, h)`` and ``cfg = quo(g, h)`` are, so called, cofactors of ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).cofactors(Poly(x**2 - 3*x + 2, x)) (Poly(x - 1, x, domain='ZZ'), Poly(x + 1, x, domain='ZZ'), Poly(x - 2, x, domain='ZZ'))
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'cofactors'): h, cff, cfg = F.cofactors(G) else: # pragma: no cover raise OperationNotSupported(f, 'cofactors')
return per(h), per(cff), per(cfg)
def gcd(f, g): """ Returns the polynomial GCD of ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).gcd(Poly(x**2 - 3*x + 2, x)) Poly(x - 1, x, domain='ZZ')
"""
else: # pragma: no cover raise OperationNotSupported(f, 'gcd')
def lcm(f, g): """ Returns polynomial LCM of ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 1, x).lcm(Poly(x**2 - 3*x + 2, x)) Poly(x**3 - 2*x**2 - x + 2, x, domain='ZZ')
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'lcm'): result = F.lcm(G) else: # pragma: no cover raise OperationNotSupported(f, 'lcm')
return per(result)
def trunc(f, p): """ Reduce ``f`` modulo a constant ``p``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x**3 + 3*x**2 + 5*x + 7, x).trunc(3) Poly(-x**3 - x + 1, x, domain='ZZ')
""" p = f.rep.dom.convert(p)
if hasattr(f.rep, 'trunc'): result = f.rep.trunc(p) else: # pragma: no cover raise OperationNotSupported(f, 'trunc')
return f.per(result)
def monic(self, auto=True): """ Divides all coefficients by ``LC(f)``.
Examples ========
>>> from sympy import Poly, ZZ >>> from sympy.abc import x
>>> Poly(3*x**2 + 6*x + 9, x, domain=ZZ).monic() Poly(x**2 + 2*x + 3, x, domain='QQ')
>>> Poly(3*x**2 + 4*x + 2, x, domain=ZZ).monic() Poly(x**2 + 4/3*x + 2/3, x, domain='QQ')
"""
else: # pragma: no cover raise OperationNotSupported(f, 'monic')
def content(f): """ Returns the GCD of polynomial coefficients.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(6*x**2 + 8*x + 12, x).content() 2
""" if hasattr(f.rep, 'content'): result = f.rep.content() else: # pragma: no cover raise OperationNotSupported(f, 'content')
return f.rep.dom.to_sympy(result)
def primitive(f): """ Returns the content and a primitive form of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x**2 + 8*x + 12, x).primitive() (2, Poly(x**2 + 4*x + 6, x, domain='ZZ'))
""" else: # pragma: no cover raise OperationNotSupported(f, 'primitive')
def compose(f, g): """ Computes the functional composition of ``f`` and ``g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + x, x).compose(Poly(x - 1, x)) Poly(x**2 - x, x, domain='ZZ')
""" _, per, F, G = f._unify(g)
if hasattr(f.rep, 'compose'): result = F.compose(G) else: # pragma: no cover raise OperationNotSupported(f, 'compose')
return per(result)
def decompose(f): """ Computes a functional decomposition of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**4 + 2*x**3 - x - 1, x, domain='ZZ').decompose() [Poly(x**2 - x - 1, x, domain='ZZ'), Poly(x**2 + x, x, domain='ZZ')]
""" else: # pragma: no cover raise OperationNotSupported(f, 'decompose')
def shift(f, a): """ Efficiently compute Taylor shift ``f(x + a)``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 2*x + 1, x).shift(2) Poly(x**2 + 2*x + 1, x, domain='ZZ')
""" if hasattr(f.rep, 'shift'): result = f.rep.shift(a) else: # pragma: no cover raise OperationNotSupported(f, 'shift')
return f.per(result)
def sturm(self, auto=True): """ Computes the Sturm sequence of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 - 2*x**2 + x - 3, x).sturm() [Poly(x**3 - 2*x**2 + x - 3, x, domain='QQ'), Poly(3*x**2 - 4*x + 1, x, domain='QQ'), Poly(2/9*x + 25/9, x, domain='QQ'), Poly(-2079/4, x, domain='QQ')]
""" f = self
if auto and f.rep.dom.has_Ring: f = f.to_field()
if hasattr(f.rep, 'sturm'): result = f.rep.sturm() else: # pragma: no cover raise OperationNotSupported(f, 'sturm')
return list(map(f.per, result))
def gff_list(f): """ Computes greatest factorial factorization of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = x**5 + 2*x**4 - x**3 - 2*x**2
>>> Poly(f).gff_list() [(Poly(x, x, domain='ZZ'), 1), (Poly(x + 2, x, domain='ZZ'), 4)]
""" if hasattr(f.rep, 'gff_list'): result = f.rep.gff_list() else: # pragma: no cover raise OperationNotSupported(f, 'gff_list')
return [(f.per(g), k) for g, k in result]
def sqf_norm(f): """ Computes square-free norm of ``f``.
Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and ``r(x) = Norm(g(x))`` is a square-free polynomial over ``K``, where ``a`` is the algebraic extension of the ground domain.
Examples ========
>>> from sympy import Poly, sqrt >>> from sympy.abc import x
>>> s, f, r = Poly(x**2 + 1, x, extension=[sqrt(3)]).sqf_norm()
>>> s 1 >>> f Poly(x**2 - 2*sqrt(3)*x + 4, x, domain='QQ<sqrt(3)>') >>> r Poly(x**4 - 4*x**2 + 16, x, domain='QQ')
""" if hasattr(f.rep, 'sqf_norm'): s, g, r = f.rep.sqf_norm() else: # pragma: no cover raise OperationNotSupported(f, 'sqf_norm')
return s, f.per(g), f.per(r)
def sqf_part(f): """ Computes square-free part of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**3 - 3*x - 2, x).sqf_part() Poly(x**2 - x - 2, x, domain='ZZ')
""" if hasattr(f.rep, 'sqf_part'): result = f.rep.sqf_part() else: # pragma: no cover raise OperationNotSupported(f, 'sqf_part')
return f.per(result)
def sqf_list(f, all=False): """ Returns a list of square-free factors of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = 2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16
>>> Poly(f).sqf_list() (2, [(Poly(x + 1, x, domain='ZZ'), 2), (Poly(x + 2, x, domain='ZZ'), 3)])
>>> Poly(f).sqf_list(all=True) (2, [(Poly(1, x, domain='ZZ'), 1), (Poly(x + 1, x, domain='ZZ'), 2), (Poly(x + 2, x, domain='ZZ'), 3)])
""" if hasattr(f.rep, 'sqf_list'): coeff, factors = f.rep.sqf_list(all) else: # pragma: no cover raise OperationNotSupported(f, 'sqf_list')
return f.rep.dom.to_sympy(coeff), [(f.per(g), k) for g, k in factors]
def sqf_list_include(f, all=False): """ Returns a list of square-free factors of ``f``.
Examples ========
>>> from sympy import Poly, expand >>> from sympy.abc import x
>>> f = expand(2*(x + 1)**3*x**4) >>> f 2*x**7 + 6*x**6 + 6*x**5 + 2*x**4
>>> Poly(f).sqf_list_include() [(Poly(2, x, domain='ZZ'), 1), (Poly(x + 1, x, domain='ZZ'), 3), (Poly(x, x, domain='ZZ'), 4)]
>>> Poly(f).sqf_list_include(all=True) [(Poly(2, x, domain='ZZ'), 1), (Poly(1, x, domain='ZZ'), 2), (Poly(x + 1, x, domain='ZZ'), 3), (Poly(x, x, domain='ZZ'), 4)]
""" if hasattr(f.rep, 'sqf_list_include'): factors = f.rep.sqf_list_include(all) else: # pragma: no cover raise OperationNotSupported(f, 'sqf_list_include')
return [(f.per(g), k) for g, k in factors]
def factor_list(f): """ Returns a list of irreducible factors of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y
>>> Poly(f).factor_list() (2, [(Poly(x + y, x, y, domain='ZZ'), 1), (Poly(x**2 + 1, x, y, domain='ZZ'), 2)])
""" except DomainError: return S.One, [(f, 1)] else: # pragma: no cover raise OperationNotSupported(f, 'factor_list')
def factor_list_include(f): """ Returns a list of irreducible factors of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y
>>> Poly(f).factor_list_include() [(Poly(2*x + 2*y, x, y, domain='ZZ'), 1), (Poly(x**2 + 1, x, y, domain='ZZ'), 2)]
""" if hasattr(f.rep, 'factor_list_include'): try: factors = f.rep.factor_list_include() except DomainError: return [(f, 1)] else: # pragma: no cover raise OperationNotSupported(f, 'factor_list_include')
return [(f.per(g), k) for g, k in factors]
def intervals(f, all=False, eps=None, inf=None, sup=None, fast=False, sqf=False): """ Compute isolating intervals for roots of ``f``.
For real roots the Vincent-Akritas-Strzebonski (VAS) continued fractions method is used.
References: =========== 1. Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative Study of Two Real Root Isolation Methods . Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005. 2. Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. Vigklas: Improving the Performance of the Continued Fractions Method Using new Bounds of Positive Roots. Nonlinear Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 3, x).intervals() [((-2, -1), 1), ((1, 2), 1)] >>> Poly(x**2 - 3, x).intervals(eps=1e-2) [((-26/15, -19/11), 1), ((19/11, 26/15), 1)]
""" if eps is not None: eps = QQ.convert(eps)
if eps <= 0: raise ValueError("'eps' must be a positive rational")
if inf is not None: inf = QQ.convert(inf) if sup is not None: sup = QQ.convert(sup)
if hasattr(f.rep, 'intervals'): result = f.rep.intervals( all=all, eps=eps, inf=inf, sup=sup, fast=fast, sqf=sqf) else: # pragma: no cover raise OperationNotSupported(f, 'intervals')
if sqf: def _real(interval): s, t = interval return (QQ.to_sympy(s), QQ.to_sympy(t))
if not all: return list(map(_real, result))
def _complex(rectangle): (u, v), (s, t) = rectangle return (QQ.to_sympy(u) + I*QQ.to_sympy(v), QQ.to_sympy(s) + I*QQ.to_sympy(t))
real_part, complex_part = result
return list(map(_real, real_part)), list(map(_complex, complex_part)) else: def _real(interval): (s, t), k = interval return ((QQ.to_sympy(s), QQ.to_sympy(t)), k)
if not all: return list(map(_real, result))
def _complex(rectangle): ((u, v), (s, t)), k = rectangle return ((QQ.to_sympy(u) + I*QQ.to_sympy(v), QQ.to_sympy(s) + I*QQ.to_sympy(t)), k)
real_part, complex_part = result
return list(map(_real, real_part)), list(map(_complex, complex_part))
def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False): """ Refine an isolating interval of a root to the given precision.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 3, x).refine_root(1, 2, eps=1e-2) (19/11, 26/15)
""" if check_sqf and not f.is_sqf: raise PolynomialError("only square-free polynomials supported")
s, t = QQ.convert(s), QQ.convert(t)
if eps is not None: eps = QQ.convert(eps)
if eps <= 0: raise ValueError("'eps' must be a positive rational")
if steps is not None: steps = int(steps) elif eps is None: steps = 1
if hasattr(f.rep, 'refine_root'): S, T = f.rep.refine_root(s, t, eps=eps, steps=steps, fast=fast) else: # pragma: no cover raise OperationNotSupported(f, 'refine_root')
return QQ.to_sympy(S), QQ.to_sympy(T)
def count_roots(f, inf=None, sup=None): """ Return the number of roots of ``f`` in ``[inf, sup]`` interval.
Examples ========
>>> from sympy import Poly, I >>> from sympy.abc import x
>>> Poly(x**4 - 4, x).count_roots(-3, 3) 2 >>> Poly(x**4 - 4, x).count_roots(0, 1 + 3*I) 1
""" inf_real, sup_real = True, True
if inf is not None: inf = sympify(inf)
if inf is S.NegativeInfinity: inf = None else: re, im = inf.as_real_imag()
if not im: inf = QQ.convert(inf) else: inf, inf_real = list(map(QQ.convert, (re, im))), False
if sup is not None: sup = sympify(sup)
if sup is S.Infinity: sup = None else: re, im = sup.as_real_imag()
if not im: sup = QQ.convert(sup) else: sup, sup_real = list(map(QQ.convert, (re, im))), False
if inf_real and sup_real: if hasattr(f.rep, 'count_real_roots'): count = f.rep.count_real_roots(inf=inf, sup=sup) else: # pragma: no cover raise OperationNotSupported(f, 'count_real_roots') else: if inf_real and inf is not None: inf = (inf, QQ.zero)
if sup_real and sup is not None: sup = (sup, QQ.zero)
if hasattr(f.rep, 'count_complex_roots'): count = f.rep.count_complex_roots(inf=inf, sup=sup) else: # pragma: no cover raise OperationNotSupported(f, 'count_complex_roots')
return Integer(count)
def root(f, index, radicals=True): """ Get an indexed root of a polynomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = Poly(2*x**3 - 7*x**2 + 4*x + 4)
>>> f.root(0) -1/2 >>> f.root(1) 2 >>> f.root(2) 2 >>> f.root(3) Traceback (most recent call last): ... IndexError: root index out of [-3, 2] range, got 3
>>> Poly(x**5 + x + 1).root(0) CRootOf(x**3 - x**2 + 1, 0)
""" return sympy.polys.rootoftools.rootof(f, index, radicals=radicals)
def real_roots(f, multiple=True, radicals=True): """ Return a list of real roots with multiplicities.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x**3 - 7*x**2 + 4*x + 4).real_roots() [-1/2, 2, 2] >>> Poly(x**3 + x + 1).real_roots() [CRootOf(x**3 + x + 1, 0)]
"""
else: return group(reals, multiple=False)
def all_roots(f, multiple=True, radicals=True): """ Return a list of real and complex roots with multiplicities.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x**3 - 7*x**2 + 4*x + 4).all_roots() [-1/2, 2, 2] >>> Poly(x**3 + x + 1).all_roots() [CRootOf(x**3 + x + 1, 0), CRootOf(x**3 + x + 1, 1), CRootOf(x**3 + x + 1, 2)]
"""
else: return group(roots, multiple=False)
def nroots(f, n=15, maxsteps=50, cleanup=True): """ Compute numerical approximations of roots of ``f``.
Parameters ==========
n ... the number of digits to calculate maxsteps ... the maximum number of iterations to do
If the accuracy `n` cannot be reached in `maxsteps`, it will raise an exception. You need to rerun with higher maxsteps.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 3).nroots(n=15) [-1.73205080756888, 1.73205080756888] >>> Poly(x**2 - 3).nroots(n=30) [-1.73205080756887729352744634151, 1.73205080756887729352744634151]
""" if f.is_multivariate: raise MultivariatePolynomialError( "can't compute numerical roots of %s" % f)
if f.degree() <= 0: return []
# For integer and rational coefficients, convert them to integers only # (for accuracy). Otherwise just try to convert the coefficients to # mpmath.mpc and raise an exception if the conversion fails. if f.rep.dom is ZZ: coeffs = [int(coeff) for coeff in f.all_coeffs()] elif f.rep.dom is QQ: denoms = [coeff.q for coeff in f.all_coeffs()] from sympy.core.numbers import ilcm fac = ilcm(*denoms) coeffs = [int(coeff*fac) for coeff in f.all_coeffs()] else: coeffs = [coeff.evalf(n=n).as_real_imag() for coeff in f.all_coeffs()] try: coeffs = [mpmath.mpc(*coeff) for coeff in coeffs] except TypeError: raise DomainError("Numerical domain expected, got %s" % \ f.rep.dom)
dps = mpmath.mp.dps mpmath.mp.dps = n
try: # We need to add extra precision to guard against losing accuracy. # 10 times the degree of the polynomial seems to work well. roots = mpmath.polyroots(coeffs, maxsteps=maxsteps, cleanup=cleanup, error=False, extraprec=f.degree()*10)
# Mpmath puts real roots first, then complex ones (as does all_roots) # so we make sure this convention holds here, too. roots = list(map(sympify, sorted(roots, key=lambda r: (1 if r.imag else 0, r.real, r.imag)))) except NoConvergence: raise NoConvergence( 'convergence to root failed; try n < %s or maxsteps > %s' % ( n, maxsteps)) finally: mpmath.mp.dps = dps
return roots
def ground_roots(f): """ Compute roots of ``f`` by factorization in the ground domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**6 - 4*x**4 + 4*x**3 - x**2).ground_roots() {0: 2, 1: 2}
""" if f.is_multivariate: raise MultivariatePolynomialError( "can't compute ground roots of %s" % f)
roots = {}
for factor, k in f.factor_list()[1]: if factor.is_linear: a, b = factor.all_coeffs() roots[-b/a] = k
return roots
def nth_power_roots_poly(f, n): """ Construct a polynomial with n-th powers of roots of ``f``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = Poly(x**4 - x**2 + 1)
>>> f.nth_power_roots_poly(2) Poly(x**4 - 2*x**3 + 3*x**2 - 2*x + 1, x, domain='ZZ') >>> f.nth_power_roots_poly(3) Poly(x**4 + 2*x**2 + 1, x, domain='ZZ') >>> f.nth_power_roots_poly(4) Poly(x**4 + 2*x**3 + 3*x**2 + 2*x + 1, x, domain='ZZ') >>> f.nth_power_roots_poly(12) Poly(x**4 - 4*x**3 + 6*x**2 - 4*x + 1, x, domain='ZZ')
""" if f.is_multivariate: raise MultivariatePolynomialError( "must be a univariate polynomial")
N = sympify(n)
if N.is_Integer and N >= 1: n = int(N) else: raise ValueError("'n' must an integer and n >= 1, got %s" % n)
x = f.gen t = Dummy('t')
r = f.resultant(f.__class__.from_expr(x**n - t, x, t))
return r.replace(t, x)
def cancel(f, g, include=False): """ Cancel common factors in a rational function ``f/g``.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x)) (1, Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ'))
>>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x), include=True) (Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ'))
"""
else: # pragma: no cover raise OperationNotSupported(f, 'cancel')
else: return tuple(map(per, result))
@property def is_zero(f): """ Returns ``True`` if ``f`` is a zero polynomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(0, x).is_zero True >>> Poly(1, x).is_zero False
""" return f.rep.is_zero
@property def is_one(f): """ Returns ``True`` if ``f`` is a unit polynomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(0, x).is_one False >>> Poly(1, x).is_one True
""" return f.rep.is_one
@property def is_sqf(f): """ Returns ``True`` if ``f`` is a square-free polynomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 - 2*x + 1, x).is_sqf False >>> Poly(x**2 - 1, x).is_sqf True
""" return f.rep.is_sqf
@property def is_monic(f): """ Returns ``True`` if the leading coefficient of ``f`` is one.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x + 2, x).is_monic True >>> Poly(2*x + 2, x).is_monic False
""" return f.rep.is_monic
@property def is_primitive(f): """ Returns ``True`` if GCD of the coefficients of ``f`` is one.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(2*x**2 + 6*x + 12, x).is_primitive False >>> Poly(x**2 + 3*x + 6, x).is_primitive True
""" return f.rep.is_primitive
@property def is_ground(f): """ Returns ``True`` if ``f`` is an element of the ground domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x, x).is_ground False >>> Poly(2, x).is_ground True >>> Poly(y, x).is_ground True
"""
@property def is_linear(f): """ Returns ``True`` if ``f`` is linear in all its variables.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x + y + 2, x, y).is_linear True >>> Poly(x*y + 2, x, y).is_linear False
""" return f.rep.is_linear
@property def is_quadratic(f): """ Returns ``True`` if ``f`` is quadratic in all its variables.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x*y + 2, x, y).is_quadratic True >>> Poly(x*y**2 + 2, x, y).is_quadratic False
""" return f.rep.is_quadratic
@property def is_monomial(f): """ Returns ``True`` if ``f`` is zero or has only one term.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(3*x**2, x).is_monomial True >>> Poly(3*x**2 + 1, x).is_monomial False
"""
@property def is_homogeneous(f): """ Returns ``True`` if ``f`` is a homogeneous polynomial.
A homogeneous polynomial is a polynomial whose all monomials with non-zero coefficients have the same total degree. If you want not only to check if a polynomial is homogeneous but also compute its homogeneous order, then use :func:`Poly.homogeneous_order`.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + x*y, x, y).is_homogeneous True >>> Poly(x**3 + x*y, x, y).is_homogeneous False
""" return f.rep.is_homogeneous
@property def is_irreducible(f): """ Returns ``True`` if ``f`` has no factors over its domain.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> Poly(x**2 + x + 1, x, modulus=2).is_irreducible True >>> Poly(x**2 + 1, x, modulus=2).is_irreducible False
"""
@property def is_univariate(f): """ Returns ``True`` if ``f`` is a univariate polynomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + x + 1, x).is_univariate True >>> Poly(x*y**2 + x*y + 1, x, y).is_univariate False >>> Poly(x*y**2 + x*y + 1, x).is_univariate True >>> Poly(x**2 + x + 1, x, y).is_univariate False
"""
@property def is_multivariate(f): """ Returns ``True`` if ``f`` is a multivariate polynomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x, y
>>> Poly(x**2 + x + 1, x).is_multivariate False >>> Poly(x*y**2 + x*y + 1, x, y).is_multivariate True >>> Poly(x*y**2 + x*y + 1, x).is_multivariate False >>> Poly(x**2 + x + 1, x, y).is_multivariate True
"""
@property def is_cyclotomic(f): """ Returns ``True`` if ``f`` is a cyclotomic polnomial.
Examples ========
>>> from sympy import Poly >>> from sympy.abc import x
>>> f = x**16 + x**14 - x**10 + x**8 - x**6 + x**2 + 1
>>> Poly(f).is_cyclotomic False
>>> g = x**16 + x**14 - x**10 - x**8 - x**6 + x**2 + 1
>>> Poly(g).is_cyclotomic True
"""
def __abs__(f): return f.abs()
def __neg__(f): return f.neg()
@_sympifyit('g', NotImplemented) def __add__(f, g): if not g.is_Poly: try: g = f.__class__(g, *f.gens) except PolynomialError: return f.as_expr() + g
return f.add(g)
@_sympifyit('g', NotImplemented) def __radd__(f, g): if not g.is_Poly: try: g = f.__class__(g, *f.gens) except PolynomialError: return g + f.as_expr()
return g.add(f)
@_sympifyit('g', NotImplemented) def __sub__(f, g): try: g = f.__class__(g, *f.gens) except PolynomialError: return f.as_expr() - g
@_sympifyit('g', NotImplemented) def __rsub__(f, g): if not g.is_Poly: try: g = f.__class__(g, *f.gens) except PolynomialError: return g - f.as_expr()
return g.sub(f)
@_sympifyit('g', NotImplemented) def __mul__(f, g): if not g.is_Poly: try: g = f.__class__(g, *f.gens) except PolynomialError: return f.as_expr()*g
return f.mul(g)
@_sympifyit('g', NotImplemented) def __rmul__(f, g): if not g.is_Poly: try: g = f.__class__(g, *f.gens) except PolynomialError: return g*f.as_expr()
return g.mul(f)
@_sympifyit('n', NotImplemented) def __pow__(f, n): if n.is_Integer and n >= 0: return f.pow(n) else: return f.as_expr()**n
@_sympifyit('g', NotImplemented) def __divmod__(f, g): if not g.is_Poly: g = f.__class__(g, *f.gens)
return f.div(g)
@_sympifyit('g', NotImplemented) def __rdivmod__(f, g): if not g.is_Poly: g = f.__class__(g, *f.gens)
return g.div(f)
@_sympifyit('g', NotImplemented) def __mod__(f, g): if not g.is_Poly: g = f.__class__(g, *f.gens)
return f.rem(g)
@_sympifyit('g', NotImplemented) def __rmod__(f, g): if not g.is_Poly: g = f.__class__(g, *f.gens)
return g.rem(f)
@_sympifyit('g', NotImplemented) def __floordiv__(f, g): if not g.is_Poly: g = f.__class__(g, *f.gens)
return f.quo(g)
@_sympifyit('g', NotImplemented) def __rfloordiv__(f, g): if not g.is_Poly: g = f.__class__(g, *f.gens)
return g.quo(f)
@_sympifyit('g', NotImplemented) def __div__(f, g):
@_sympifyit('g', NotImplemented) def __rdiv__(f, g): return g.as_expr()/f.as_expr()
__truediv__ = __div__ __rtruediv__ = __rdiv__
@_sympifyit('other', NotImplemented) def __eq__(self, other):
try: g = f.__class__(g, f.gens, domain=f.get_domain()) except (PolynomialError, DomainError, CoercionFailed): return False
return False
try: dom = f.rep.dom.unify(g.rep.dom, f.gens) except UnificationFailed: return False
f = f.set_domain(dom) g = g.set_domain(dom)
@_sympifyit('g', NotImplemented) def __ne__(f, g): return not f.__eq__(g)
def __nonzero__(f): return not f.is_zero
__bool__ = __nonzero__
def eq(f, g, strict=False): if not strict: return f.__eq__(g) else: return f._strict_eq(sympify(g))
def ne(f, g, strict=False): return not f.eq(g, strict=strict)
def _strict_eq(f, g): return isinstance(g, f.__class__) and f.gens == g.gens and f.rep.eq(g.rep, strict=True)
@public class PurePoly(Poly): """Class for representing pure polynomials. """
def _hashable_content(self): """Allow SymPy to hash Poly instances. """
def __hash__(self):
@property def free_symbols(self): """ Free symbols of a polynomial.
Examples ========
>>> from sympy import PurePoly >>> from sympy.abc import x, y
>>> PurePoly(x**2 + 1).free_symbols set() >>> PurePoly(x**2 + y).free_symbols set() >>> PurePoly(x**2 + y, x).free_symbols set([y])
""" return self.free_symbols_in_domain
@_sympifyit('other', NotImplemented) def __eq__(self, other):
try: g = f.__class__(g, f.gens, domain=f.get_domain()) except (PolynomialError, DomainError, CoercionFailed): return False
return False
try: dom = f.rep.dom.unify(g.rep.dom, f.gens) except UnificationFailed: return False
f = f.set_domain(dom) g = g.set_domain(dom)
def _strict_eq(f, g): return isinstance(g, f.__class__) and f.rep.eq(g.rep, strict=True)
def _unify(f, g): g = sympify(g)
if not g.is_Poly: try: return f.rep.dom, f.per, f.rep, f.rep.per(f.rep.dom.from_sympy(g)) except CoercionFailed: raise UnificationFailed("can't unify %s with %s" % (f, g))
if len(f.gens) != len(g.gens): raise UnificationFailed("can't unify %s with %s" % (f, g))
if not (isinstance(f.rep, DMP) and isinstance(g.rep, DMP)): raise UnificationFailed("can't unify %s with %s" % (f, g))
cls = f.__class__ gens = f.gens
dom = f.rep.dom.unify(g.rep.dom, gens)
F = f.rep.convert(dom) G = g.rep.convert(dom)
def per(rep, dom=dom, gens=gens, remove=None): if remove is not None: gens = gens[:remove] + gens[remove + 1:]
if not gens: return dom.to_sympy(rep)
return cls.new(rep, *gens)
return dom, per, F, G
@public def poly_from_expr(expr, *gens, **args): """Construct a polynomial from an expression. """
def _poly_from_expr(expr, opt): """Construct a polynomial from an expression. """
raise PolificationFailed(opt, orig, expr)
except GeneratorsNeeded: raise PolificationFailed(opt, orig, expr)
else: coeffs = list(map(domain.from_sympy, coeffs))
@public def parallel_poly_from_expr(exprs, *gens, **args): """Construct polynomials from expressions. """
def _parallel_poly_from_expr(exprs, opt): """Construct polynomials from expressions. """
else:
else: failed = True
raise PolificationFailed(opt, origs, exprs, True)
# XXX: this is a temporary solution
raise PolynomialError("Piecewise generators do not make sense")
else: coeffs_list = list(map(domain.from_sympy, coeffs_list))
def _update_args(args, key, value): """Add a new ``(key, value)`` pair to arguments ``dict``. """ args = dict(args)
if key not in args: args[key] = value
return args
@public def degree(f, *gens, **args): """ Return the degree of ``f`` in the given variable.
The degree of 0 is negative infinity.
Examples ========
>>> from sympy import degree >>> from sympy.abc import x, y
>>> degree(x**2 + y*x + 1, gen=x) 2 >>> degree(x**2 + y*x + 1, gen=y) 1 >>> degree(0, x) -oo
"""
except PolificationFailed as exc: raise ComputationFailed('degree', 1, exc)
@public def degree_list(f, *gens, **args): """ Return a list of degrees of ``f`` in all variables.
Examples ========
>>> from sympy import degree_list >>> from sympy.abc import x, y
>>> degree_list(x**2 + y*x + 1) (2, 1)
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('degree_list', 1, exc)
degrees = F.degree_list()
return tuple(map(Integer, degrees))
@public def LC(f, *gens, **args): """ Return the leading coefficient of ``f``.
Examples ========
>>> from sympy import LC >>> from sympy.abc import x, y
>>> LC(4*x**2 + 2*x*y**2 + x*y + 3*y) 4
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('LC', 1, exc)
return F.LC(order=opt.order)
@public def LM(f, *gens, **args): """ Return the leading monomial of ``f``.
Examples ========
>>> from sympy import LM >>> from sympy.abc import x, y
>>> LM(4*x**2 + 2*x*y**2 + x*y + 3*y) x**2
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('LM', 1, exc)
monom = F.LM(order=opt.order) return monom.as_expr()
@public def LT(f, *gens, **args): """ Return the leading term of ``f``.
Examples ========
>>> from sympy import LT >>> from sympy.abc import x, y
>>> LT(4*x**2 + 2*x*y**2 + x*y + 3*y) 4*x**2
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('LT', 1, exc)
monom, coeff = F.LT(order=opt.order) return coeff*monom.as_expr()
@public def pdiv(f, g, *gens, **args): """ Compute polynomial pseudo-division of ``f`` and ``g``.
Examples ========
>>> from sympy import pdiv >>> from sympy.abc import x
>>> pdiv(x**2 + 1, 2*x - 4) (2*x + 4, 20)
""" options.allowed_flags(args, ['polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('pdiv', 2, exc)
q, r = F.pdiv(G)
if not opt.polys: return q.as_expr(), r.as_expr() else: return q, r
@public def prem(f, g, *gens, **args): """ Compute polynomial pseudo-remainder of ``f`` and ``g``.
Examples ========
>>> from sympy import prem >>> from sympy.abc import x
>>> prem(x**2 + 1, 2*x - 4) 20
""" options.allowed_flags(args, ['polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('prem', 2, exc)
r = F.prem(G)
if not opt.polys: return r.as_expr() else: return r
@public def pquo(f, g, *gens, **args): """ Compute polynomial pseudo-quotient of ``f`` and ``g``.
Examples ========
>>> from sympy import pquo >>> from sympy.abc import x
>>> pquo(x**2 + 1, 2*x - 4) 2*x + 4 >>> pquo(x**2 - 1, 2*x - 1) 2*x + 1
""" options.allowed_flags(args, ['polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('pquo', 2, exc)
try: q = F.pquo(G) except ExactQuotientFailed: raise ExactQuotientFailed(f, g)
if not opt.polys: return q.as_expr() else: return q
@public def pexquo(f, g, *gens, **args): """ Compute polynomial exact pseudo-quotient of ``f`` and ``g``.
Examples ========
>>> from sympy import pexquo >>> from sympy.abc import x
>>> pexquo(x**2 - 1, 2*x - 2) 2*x + 2
>>> pexquo(x**2 + 1, 2*x - 4) Traceback (most recent call last): ... ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
""" options.allowed_flags(args, ['polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('pexquo', 2, exc)
q = F.pexquo(G)
if not opt.polys: return q.as_expr() else: return q
@public def div(f, g, *gens, **args): """ Compute polynomial division of ``f`` and ``g``.
Examples ========
>>> from sympy import div, ZZ, QQ >>> from sympy.abc import x
>>> div(x**2 + 1, 2*x - 4, domain=ZZ) (0, x**2 + 1) >>> div(x**2 + 1, 2*x - 4, domain=QQ) (x/2 + 1, 5)
""" options.allowed_flags(args, ['auto', 'polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('div', 2, exc)
q, r = F.div(G, auto=opt.auto)
if not opt.polys: return q.as_expr(), r.as_expr() else: return q, r
@public def rem(f, g, *gens, **args): """ Compute polynomial remainder of ``f`` and ``g``.
Examples ========
>>> from sympy import rem, ZZ, QQ >>> from sympy.abc import x
>>> rem(x**2 + 1, 2*x - 4, domain=ZZ) x**2 + 1 >>> rem(x**2 + 1, 2*x - 4, domain=QQ) 5
""" options.allowed_flags(args, ['auto', 'polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('rem', 2, exc)
r = F.rem(G, auto=opt.auto)
if not opt.polys: return r.as_expr() else: return r
@public def quo(f, g, *gens, **args): """ Compute polynomial quotient of ``f`` and ``g``.
Examples ========
>>> from sympy import quo >>> from sympy.abc import x
>>> quo(x**2 + 1, 2*x - 4) x/2 + 1 >>> quo(x**2 - 1, x - 1) x + 1
""" options.allowed_flags(args, ['auto', 'polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('quo', 2, exc)
q = F.quo(G, auto=opt.auto)
if not opt.polys: return q.as_expr() else: return q
@public def exquo(f, g, *gens, **args): """ Compute polynomial exact quotient of ``f`` and ``g``.
Examples ========
>>> from sympy import exquo >>> from sympy.abc import x
>>> exquo(x**2 - 1, x - 1) x + 1
>>> exquo(x**2 + 1, 2*x - 4) Traceback (most recent call last): ... ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
""" options.allowed_flags(args, ['auto', 'polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('exquo', 2, exc)
q = F.exquo(G, auto=opt.auto)
if not opt.polys: return q.as_expr() else: return q
@public def half_gcdex(f, g, *gens, **args): """ Half extended Euclidean algorithm of ``f`` and ``g``.
Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``.
Examples ========
>>> from sympy import half_gcdex >>> from sympy.abc import x
>>> half_gcdex(x**4 - 2*x**3 - 6*x**2 + 12*x + 15, x**3 + x**2 - 4*x - 4) (-x/5 + 3/5, x + 1)
""" options.allowed_flags(args, ['auto', 'polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: domain, (a, b) = construct_domain(exc.exprs)
try: s, h = domain.half_gcdex(a, b) except NotImplementedError: raise ComputationFailed('half_gcdex', 2, exc) else: return domain.to_sympy(s), domain.to_sympy(h)
s, h = F.half_gcdex(G, auto=opt.auto)
if not opt.polys: return s.as_expr(), h.as_expr() else: return s, h
@public def gcdex(f, g, *gens, **args): """ Extended Euclidean algorithm of ``f`` and ``g``.
Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``.
Examples ========
>>> from sympy import gcdex >>> from sympy.abc import x
>>> gcdex(x**4 - 2*x**3 - 6*x**2 + 12*x + 15, x**3 + x**2 - 4*x - 4) (-x/5 + 3/5, x**2/5 - 6*x/5 + 2, x + 1)
""" options.allowed_flags(args, ['auto', 'polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: domain, (a, b) = construct_domain(exc.exprs)
try: s, t, h = domain.gcdex(a, b) except NotImplementedError: raise ComputationFailed('gcdex', 2, exc) else: return domain.to_sympy(s), domain.to_sympy(t), domain.to_sympy(h)
s, t, h = F.gcdex(G, auto=opt.auto)
if not opt.polys: return s.as_expr(), t.as_expr(), h.as_expr() else: return s, t, h
@public def invert(f, g, *gens, **args): """ Invert ``f`` modulo ``g`` when possible.
Examples ========
>>> from sympy import invert, S >>> from sympy.core.numbers import mod_inverse >>> from sympy.abc import x
>>> invert(x**2 - 1, 2*x - 1) -4/3
>>> invert(x**2 - 1, x - 1) Traceback (most recent call last): ... NotInvertible: zero divisor
For more efficient inversion of Rationals, use the ``mod_inverse`` function:
>>> mod_inverse(3, 5) 2 >>> (S(2)/5).invert(S(7)/3) 5/2
See Also ======== sympy.core.numbers.mod_inverse """ options.allowed_flags(args, ['auto', 'polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: domain, (a, b) = construct_domain(exc.exprs)
try: return domain.to_sympy(domain.invert(a, b)) except NotImplementedError: raise ComputationFailed('invert', 2, exc)
h = F.invert(G, auto=opt.auto)
if not opt.polys: return h.as_expr() else: return h
@public def subresultants(f, g, *gens, **args): """ Compute subresultant PRS of ``f`` and ``g``.
Examples ========
>>> from sympy import subresultants >>> from sympy.abc import x
>>> subresultants(x**2 + 1, x**2 - 1) [x**2 + 1, x**2 - 1, -2]
""" options.allowed_flags(args, ['polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('subresultants', 2, exc)
result = F.subresultants(G)
if not opt.polys: return [r.as_expr() for r in result] else: return result
@public def resultant(f, g, *gens, **args): """ Compute resultant of ``f`` and ``g``.
Examples ========
>>> from sympy import resultant >>> from sympy.abc import x
>>> resultant(x**2 + 1, x**2 - 1) 4
"""
except PolificationFailed as exc: raise ComputationFailed('resultant', 2, exc)
result, R = F.resultant(G, includePRS=includePRS) else:
return result.as_expr(), [r.as_expr() for r in R] else: if includePRS: return result, R return result
@public def discriminant(f, *gens, **args): """ Compute discriminant of ``f``.
Examples ========
>>> from sympy import discriminant >>> from sympy.abc import x
>>> discriminant(x**2 + 2*x + 3) -8
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('discriminant', 1, exc)
result = F.discriminant()
if not opt.polys: return result.as_expr() else: return result
@public def cofactors(f, g, *gens, **args): """ Compute GCD and cofactors of ``f`` and ``g``.
Returns polynomials ``(h, cff, cfg)`` such that ``h = gcd(f, g)``, and ``cff = quo(f, h)`` and ``cfg = quo(g, h)`` are, so called, cofactors of ``f`` and ``g``.
Examples ========
>>> from sympy import cofactors >>> from sympy.abc import x
>>> cofactors(x**2 - 1, x**2 - 3*x + 2) (x - 1, x + 1, x - 2)
""" options.allowed_flags(args, ['polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: domain, (a, b) = construct_domain(exc.exprs)
try: h, cff, cfg = domain.cofactors(a, b) except NotImplementedError: raise ComputationFailed('cofactors', 2, exc) else: return domain.to_sympy(h), domain.to_sympy(cff), domain.to_sympy(cfg)
h, cff, cfg = F.cofactors(G)
if not opt.polys: return h.as_expr(), cff.as_expr(), cfg.as_expr() else: return h, cff, cfg
@public def gcd_list(seq, *gens, **args): """ Compute GCD of a list of polynomials.
Examples ========
>>> from sympy import gcd_list >>> from sympy.abc import x
>>> gcd_list([x**3 - 1, x**2 - 1, x**2 - 3*x + 2]) x - 1
"""
return domain.zero
return None
options.allowed_flags(args, ['polys'])
try: polys, opt = parallel_poly_from_expr(seq, *gens, **args) except PolificationFailed as exc: result = try_non_polynomial_gcd(exc.exprs)
if result is not None: return result else: raise ComputationFailed('gcd_list', len(seq), exc)
if not polys: if not opt.polys: return S.Zero else: return Poly(0, opt=opt)
result, polys = polys[0], polys[1:]
for poly in polys: result = result.gcd(poly)
if result.is_one: break
if not opt.polys: return result.as_expr() else: return result
@public def gcd(f, g=None, *gens, **args): """ Compute GCD of ``f`` and ``g``.
Examples ========
>>> from sympy import gcd >>> from sympy.abc import x
>>> gcd(x**2 - 1, x**2 - 3*x + 2) x - 1
""" if g is not None: gens = (g,) + gens
return gcd_list(f, *gens, **args) raise TypeError("gcd() takes 2 arguments or a sequence of arguments")
except NotImplementedError: raise ComputationFailed('gcd', 2, exc)
else: return result
@public def lcm_list(seq, *gens, **args): """ Compute LCM of a list of polynomials.
Examples ========
>>> from sympy import lcm_list >>> from sympy.abc import x
>>> lcm_list([x**3 - 1, x**2 - 1, x**2 - 3*x + 2]) x**5 - x**4 - 2*x**3 - x**2 + x + 2
"""
return domain.one
return None
options.allowed_flags(args, ['polys'])
try: polys, opt = parallel_poly_from_expr(seq, *gens, **args) except PolificationFailed as exc: result = try_non_polynomial_lcm(exc.exprs)
if result is not None: return result else: raise ComputationFailed('lcm_list', len(seq), exc)
if not polys: if not opt.polys: return S.One else: return Poly(1, opt=opt)
result, polys = polys[0], polys[1:]
for poly in polys: result = result.lcm(poly)
if not opt.polys: return result.as_expr() else: return result
@public def lcm(f, g=None, *gens, **args): """ Compute LCM of ``f`` and ``g``.
Examples ========
>>> from sympy import lcm >>> from sympy.abc import x
>>> lcm(x**2 - 1, x**2 - 3*x + 2) x**3 - 2*x**2 - x + 2
""" gens = (g,) + gens
raise TypeError("lcm() takes 2 arguments or a sequence of arguments")
except NotImplementedError: raise ComputationFailed('lcm', 2, exc)
result = F.lcm(G)
if not opt.polys: return result.as_expr() else: return result
@public def terms_gcd(f, *gens, **args): """ Remove GCD of terms from ``f``.
If the ``deep`` flag is True, then the arguments of ``f`` will have terms_gcd applied to them.
If a fraction is factored out of ``f`` and ``f`` is an Add, then an unevaluated Mul will be returned so that automatic simplification does not redistribute it. The hint ``clear``, when set to False, can be used to prevent such factoring when all coefficients are not fractions.
Examples ========
>>> from sympy import terms_gcd, cos >>> from sympy.abc import x, y >>> terms_gcd(x**6*y**2 + x**3*y, x, y) x**3*y*(x**3*y + 1)
The default action of polys routines is to expand the expression given to them. terms_gcd follows this behavior:
>>> terms_gcd((3+3*x)*(x+x*y)) 3*x*(x*y + x + y + 1)
If this is not desired then the hint ``expand`` can be set to False. In this case the expression will be treated as though it were comprised of one or more terms:
>>> terms_gcd((3+3*x)*(x+x*y), expand=False) (3*x + 3)*(x*y + x)
In order to traverse factors of a Mul or the arguments of other functions, the ``deep`` hint can be used:
>>> terms_gcd((3 + 3*x)*(x + x*y), expand=False, deep=True) 3*x*(x + 1)*(y + 1) >>> terms_gcd(cos(x + x*y), deep=True) cos(x*(y + 1))
Rationals are factored out by default:
>>> terms_gcd(x + y/2) (2*x + y)/2
Only the y-term had a coefficient that was a fraction; if one does not want to factor out the 1/2 in cases like this, the flag ``clear`` can be set to False:
>>> terms_gcd(x + y/2, clear=False) x + y/2 >>> terms_gcd(x*y/2 + y**2, clear=False) y*(x/2 + y)
The ``clear`` flag is ignored if all coefficients are fractions:
>>> terms_gcd(x/3 + y/2, clear=False) (2*x + 3*y)/6
See Also ======== sympy.core.exprtools.gcd_terms, sympy.core.exprtools.factor_terms
""" from sympy.core.relational import Equality
orig = sympify(f) if not isinstance(f, Expr) or f.is_Atom: return orig
if args.get('deep', False): new = f.func(*[terms_gcd(a, *gens, **args) for a in f.args]) args.pop('deep') args['expand'] = False return terms_gcd(new, *gens, **args)
if isinstance(f, Equality): return f
clear = args.pop('clear', True) options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: return exc.expr
J, f = F.terms_gcd()
if opt.domain.has_Ring: if opt.domain.has_Field: denom, f = f.clear_denoms(convert=True)
coeff, f = f.primitive()
if opt.domain.has_Field: coeff /= denom else: coeff = S.One
term = Mul(*[x**j for x, j in zip(f.gens, J)]) if coeff == 1: coeff = S.One if term == 1: return orig
if clear: return _keep_coeff(coeff, term*f.as_expr()) # base the clearing on the form of the original expression, not # the (perhaps) Mul that we have now coeff, f = _keep_coeff(coeff, f.as_expr(), clear=False).as_coeff_Mul() return _keep_coeff(coeff, term*f, clear=False)
@public def trunc(f, p, *gens, **args): """ Reduce ``f`` modulo a constant ``p``.
Examples ========
>>> from sympy import trunc >>> from sympy.abc import x
>>> trunc(2*x**3 + 3*x**2 + 5*x + 7, 3) -x**3 - x + 1
""" options.allowed_flags(args, ['auto', 'polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('trunc', 1, exc)
result = F.trunc(sympify(p))
if not opt.polys: return result.as_expr() else: return result
@public def monic(f, *gens, **args): """ Divide all coefficients of ``f`` by ``LC(f)``.
Examples ========
>>> from sympy import monic >>> from sympy.abc import x
>>> monic(3*x**2 + 4*x + 2) x**2 + 4*x/3 + 2/3
""" options.allowed_flags(args, ['auto', 'polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('monic', 1, exc)
result = F.monic(auto=opt.auto)
if not opt.polys: return result.as_expr() else: return result
@public def content(f, *gens, **args): """ Compute GCD of coefficients of ``f``.
Examples ========
>>> from sympy import content >>> from sympy.abc import x
>>> content(6*x**2 + 8*x + 12) 2
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('content', 1, exc)
return F.content()
@public def primitive(f, *gens, **args): """ Compute content and the primitive form of ``f``.
Examples ========
>>> from sympy.polys.polytools import primitive >>> from sympy.abc import x
>>> primitive(6*x**2 + 8*x + 12) (2, 3*x**2 + 4*x + 6)
>>> eq = (2 + 2*x)*x + 2
Expansion is performed by default:
>>> primitive(eq) (2, x**2 + x + 1)
Set ``expand`` to False to shut this off. Note that the extraction will not be recursive; use the as_content_primitive method for recursive, non-destructive Rational extraction.
>>> primitive(eq, expand=False) (1, x*(2*x + 2) + 2)
>>> eq.as_content_primitive() (2, x*(x + 1) + 1)
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('primitive', 1, exc)
cont, result = F.primitive() if not opt.polys: return cont, result.as_expr() else: return cont, result
@public def compose(f, g, *gens, **args): """ Compute functional composition ``f(g)``.
Examples ========
>>> from sympy import compose >>> from sympy.abc import x
>>> compose(x**2 + x, x - 1) x**2 - x
""" options.allowed_flags(args, ['polys'])
try: (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('compose', 2, exc)
result = F.compose(G)
if not opt.polys: return result.as_expr() else: return result
@public def decompose(f, *gens, **args): """ Compute functional decomposition of ``f``.
Examples ========
>>> from sympy import decompose >>> from sympy.abc import x
>>> decompose(x**4 + 2*x**3 - x - 1) [x**2 - x - 1, x**2 + x]
"""
except PolificationFailed as exc: raise ComputationFailed('decompose', 1, exc)
else: return result
@public def sturm(f, *gens, **args): """ Compute Sturm sequence of ``f``.
Examples ========
>>> from sympy import sturm >>> from sympy.abc import x
>>> sturm(x**3 - 2*x**2 + x - 3) [x**3 - 2*x**2 + x - 3, 3*x**2 - 4*x + 1, 2*x/9 + 25/9, -2079/4]
""" options.allowed_flags(args, ['auto', 'polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('sturm', 1, exc)
result = F.sturm(auto=opt.auto)
if not opt.polys: return [r.as_expr() for r in result] else: return result
@public def gff_list(f, *gens, **args): """ Compute a list of greatest factorial factors of ``f``.
Examples ========
>>> from sympy import gff_list, ff >>> from sympy.abc import x
>>> f = x**5 + 2*x**4 - x**3 - 2*x**2
>>> gff_list(f) [(x, 1), (x + 2, 4)]
>>> (ff(x, 1)*ff(x + 2, 4)).expand() == f True
>>> f = x**12 + 6*x**11 - 11*x**10 - 56*x**9 + 220*x**8 + 208*x**7 - \ 1401*x**6 + 1090*x**5 + 2715*x**4 - 6720*x**3 - 1092*x**2 + 5040*x
>>> gff_list(f) [(x**3 + 7, 2), (x**2 + 5*x, 3)]
>>> ff(x**3 + 7, 2)*ff(x**2 + 5*x, 3) == f True
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('gff_list', 1, exc)
factors = F.gff_list()
if not opt.polys: return [(g.as_expr(), k) for g, k in factors] else: return factors
@public def gff(f, *gens, **args): """Compute greatest factorial factorization of ``f``. """ raise NotImplementedError('symbolic falling factorial')
@public def sqf_norm(f, *gens, **args): """ Compute square-free norm of ``f``.
Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and ``r(x) = Norm(g(x))`` is a square-free polynomial over ``K``, where ``a`` is the algebraic extension of the ground domain.
Examples ========
>>> from sympy import sqf_norm, sqrt >>> from sympy.abc import x
>>> sqf_norm(x**2 + 1, extension=[sqrt(3)]) (1, x**2 - 2*sqrt(3)*x + 4, x**4 - 4*x**2 + 16)
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('sqf_norm', 1, exc)
s, g, r = F.sqf_norm()
if not opt.polys: return Integer(s), g.as_expr(), r.as_expr() else: return Integer(s), g, r
@public def sqf_part(f, *gens, **args): """ Compute square-free part of ``f``.
Examples ========
>>> from sympy import sqf_part >>> from sympy.abc import x
>>> sqf_part(x**3 - 3*x - 2) x**2 - x - 2
""" options.allowed_flags(args, ['polys'])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('sqf_part', 1, exc)
result = F.sqf_part()
if not opt.polys: return result.as_expr() else: return result
def _sorted_factors(factors, method): """Sort a list of ``(expr, exp)`` pairs. """ def key(obj): poly, exp = obj rep = poly.rep.rep return (exp, len(rep), len(poly.gens), rep) else:
def _factors_product(factors): """Multiply a list of ``(expr, exp)`` pairs. """
def _symbolic_factor_list(expr, opt, method): """Helper function for :func:`_symbolic_factor`. """
for i in Mul.make_args(expr)] args.extend(arg.args) continue coeff *= arg continue factors.append((base, exp)) continue else:
except PolificationFailed as exc: factors.append((exc.expr, exp)) else:
factors.append((_coeff, exp)) else:
else:
else:
def _symbolic_factor(expr, opt, method): """Helper function for :func:`_factor`. """ return expr._eval_factor() elif hasattr(expr, 'args'): return expr.func(*[_symbolic_factor(arg, opt, method) for arg in expr.args]) elif hasattr(expr, '__iter__'): return expr.__class__([_symbolic_factor(arg, opt, method) for arg in expr]) else: return expr
def _generic_factor_list(expr, gens, args, method): """Helper function for :func:`sqf_list` and :func:`factor_list`. """
raise PolynomialError("a polynomial expected, got %s" % expr)
f, _ = _poly_from_expr(f, _opt) factors[i] = (f, k)
else: return coeff, fp, fq else: raise PolynomialError("a polynomial expected, got %s" % expr)
def _generic_factor(expr, gens, args, method): """Helper function for :func:`sqf` and :func:`factor`. """
def to_rational_coeffs(f): """ try to transform a polynomial to have rational coefficients
try to find a transformation ``x = alpha*y``
``f(x) = lc*alpha**n * g(y)`` where ``g`` is a polynomial with rational coefficients, ``lc`` the leading coefficient.
If this fails, try ``x = y + beta`` ``f(x) = g(y)``
Returns ``None`` if ``g`` not found; ``(lc, alpha, None, g)`` in case of rescaling ``(None, None, beta, g)`` in case of translation
Notes =====
Currently it transforms only polynomials without roots larger than 2.
Examples ========
>>> from sympy import sqrt, Poly, simplify >>> from sympy.polys.polytools import to_rational_coeffs >>> from sympy.abc import x >>> p = Poly(((x**2-1)*(x-2)).subs({x:x*(1 + sqrt(2))}), x, domain='EX') >>> lc, r, _, g = to_rational_coeffs(p) >>> lc, r (7 + 5*sqrt(2), -2*sqrt(2) + 2) >>> g Poly(x**3 + x**2 - 1/4*x - 1/4, x, domain='QQ') >>> r1 = simplify(1/r) >>> Poly(lc*r**3*(g.as_expr()).subs({x:x*r1}), x, domain='EX') == p True
""" from sympy.simplify.simplify import simplify
def _try_rescale(f, f1=None): """ try rescaling ``x -> alpha*x`` to convert f to a polynomial with rational coefficients. Returns ``alpha, f``; if the rescaling is successful, ``alpha`` is the rescaling factor, and ``f`` is the rescaled polynomial; else ``alpha`` is ``None``. """ from sympy.core.add import Add if not len(f.gens) == 1 or not (f.gens[0]).is_Atom: return None, f n = f.degree() lc = f.LC() f1 = f1 or f1.monic() coeffs = f1.all_coeffs()[1:] coeffs = [simplify(coeffx) for coeffx in coeffs] if coeffs[-2]: rescale1_x = simplify(coeffs[-2]/coeffs[-1]) coeffs1 = [] for i in range(len(coeffs)): coeffx = simplify(coeffs[i]*rescale1_x**(i + 1)) if not coeffx.is_rational: break coeffs1.append(coeffx) else: rescale_x = simplify(1/rescale1_x) x = f.gens[0] v = [x**n] for i in range(1, n + 1): v.append(coeffs1[i - 1]*x**(n - i)) f = Add(*v) f = Poly(f) return lc, rescale_x, f return None
def _try_translate(f, f1=None): """ try translating ``x -> x + alpha`` to convert f to a polynomial with rational coefficients. Returns ``alpha, f``; if the translating is successful, ``alpha`` is the translating factor, and ``f`` is the shifted polynomial; else ``alpha`` is ``None``. """ from sympy.core.add import Add if not len(f.gens) == 1 or not (f.gens[0]).is_Atom: return None, f n = f.degree() f1 = f1 or f1.monic() coeffs = f1.all_coeffs()[1:] c = simplify(coeffs[0]) if c and not c.is_rational: func = Add if c.is_Add: args = c.args func = c.func else: args = [c] sifted = sift(args, lambda z: z.is_rational) c1, c2 = sifted[True], sifted[False] alpha = -func(*c2)/n f2 = f1.shift(alpha) return alpha, f2 return None
def _has_square_roots(p): """ Return True if ``f`` is a sum with square roots but no other root """ from sympy.core.exprtools import Factors coeffs = p.coeffs() has_sq = False for y in coeffs: for x in Add.make_args(y): f = Factors(x).factors r = [wx.q for b, wx in f.items() if b.is_number and wx.is_Rational and wx.q >= 2] if not r: continue if min(r) == 2: has_sq = True if max(r) > 2: return False return has_sq
if f.get_domain().is_EX and _has_square_roots(f): f1 = f.monic() r = _try_rescale(f, f1) if r: return r[0], r[1], None, r[2] else: r = _try_translate(f, f1) if r: return None, None, r[0], r[1] return None
def _torational_factor_list(p, x): """ helper function to factor polynomial using to_rational_coeffs
Examples ========
>>> from sympy.polys.polytools import _torational_factor_list >>> from sympy.abc import x >>> from sympy import sqrt, expand, Mul >>> p = expand(((x**2-1)*(x-2)).subs({x:x*(1 + sqrt(2))})) >>> factors = _torational_factor_list(p, x); factors (-2, [(-x*(1 + sqrt(2))/2 + 1, 1), (-x*(1 + sqrt(2)) - 1, 1), (-x*(1 + sqrt(2)) + 1, 1)]) >>> expand(factors[0]*Mul(*[z[0] for z in factors[1]])) == p True >>> p = expand(((x**2-1)*(x-2)).subs({x:x + sqrt(2)})) >>> factors = _torational_factor_list(p, x); factors (1, [(x - 2 + sqrt(2), 1), (x - 1 + sqrt(2), 1), (x + 1 + sqrt(2), 1)]) >>> expand(factors[0]*Mul(*[z[0] for z in factors[1]])) == p True
""" from sympy.simplify.simplify import simplify p1 = Poly(p, x, domain='EX') n = p1.degree() res = to_rational_coeffs(p1) if not res: return None lc, r, t, g = res factors = factor_list(g.as_expr()) if lc: c = simplify(factors[0]*lc*r**n) r1 = simplify(1/r) a = [] for z in factors[1:][0]: a.append((simplify(z[0].subs({x: x*r1})), z[1])) else: c = factors[0] a = [] for z in factors[1:][0]: a.append((z[0].subs({x: x - t}), z[1])) return (c, a)
@public def sqf_list(f, *gens, **args): """ Compute a list of square-free factors of ``f``.
Examples ========
>>> from sympy import sqf_list >>> from sympy.abc import x
>>> sqf_list(2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16) (2, [(x + 1, 2), (x + 2, 3)])
""" return _generic_factor_list(f, gens, args, method='sqf')
@public def sqf(f, *gens, **args): """ Compute square-free factorization of ``f``.
Examples ========
>>> from sympy import sqf >>> from sympy.abc import x
>>> sqf(2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16) 2*(x + 1)**2*(x + 2)**3
""" return _generic_factor(f, gens, args, method='sqf')
@public def factor_list(f, *gens, **args): """ Compute a list of irreducible factors of ``f``.
Examples ========
>>> from sympy import factor_list >>> from sympy.abc import x, y
>>> factor_list(2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y) (2, [(x + y, 1), (x**2 + 1, 2)])
"""
@public def factor(f, *gens, **args): """ Compute the factorization of expression, ``f``, into irreducibles. (To factor an integer into primes, use ``factorint``.)
There two modes implemented: symbolic and formal. If ``f`` is not an instance of :class:`Poly` and generators are not specified, then the former mode is used. Otherwise, the formal mode is used.
In symbolic mode, :func:`factor` will traverse the expression tree and factor its components without any prior expansion, unless an instance of :class:`Add` is encountered (in this case formal factorization is used). This way :func:`factor` can handle large or symbolic exponents.
By default, the factorization is computed over the rationals. To factor over other domain, e.g. an algebraic or finite field, use appropriate options: ``extension``, ``modulus`` or ``domain``.
Examples ========
>>> from sympy import factor, sqrt >>> from sympy.abc import x, y
>>> factor(2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y) 2*(x + y)*(x**2 + 1)**2
>>> factor(x**2 + 1) x**2 + 1 >>> factor(x**2 + 1, modulus=2) (x + 1)**2 >>> factor(x**2 + 1, gaussian=True) (x - I)*(x + I)
>>> factor(x**2 - 2, extension=sqrt(2)) (x - sqrt(2))*(x + sqrt(2))
>>> factor((x**2 - 1)/(x**2 + 4*x + 4)) (x - 1)*(x + 1)/(x + 2)**2 >>> factor((x**2 + 4*x + 4)**10000000*(x**2 + 1)) (x + 2)**20000000*(x**2 + 1)
By default, factor deals with an expression as a whole:
>>> eq = 2**(x**2 + 2*x + 1) >>> factor(eq) 2**(x**2 + 2*x + 1)
If the ``deep`` flag is True then subexpressions will be factored:
>>> factor(eq, deep=True) 2**((x + 1)**2)
See Also ======== sympy.ntheory.factor_.factorint
"""
except PolynomialError as msg: if not f.is_commutative: from sympy.core.exprtools import factor_nc return factor_nc(f) else: raise PolynomialError(msg)
@public def intervals(F, all=False, eps=None, inf=None, sup=None, strict=False, fast=False, sqf=False): """ Compute isolating intervals for roots of ``f``.
Examples ========
>>> from sympy import intervals >>> from sympy.abc import x
>>> intervals(x**2 - 3) [((-2, -1), 1), ((1, 2), 1)] >>> intervals(x**2 - 3, eps=1e-2) [((-26/15, -19/11), 1), ((19/11, 26/15), 1)]
""" if not hasattr(F, '__iter__'): try: F = Poly(F) except GeneratorsNeeded: return []
return F.intervals(all=all, eps=eps, inf=inf, sup=sup, fast=fast, sqf=sqf) else: polys, opt = parallel_poly_from_expr(F, domain='QQ')
if len(opt.gens) > 1: raise MultivariatePolynomialError
for i, poly in enumerate(polys): polys[i] = poly.rep.rep
if eps is not None: eps = opt.domain.convert(eps)
if eps <= 0: raise ValueError("'eps' must be a positive rational")
if inf is not None: inf = opt.domain.convert(inf) if sup is not None: sup = opt.domain.convert(sup)
intervals = dup_isolate_real_roots_list(polys, opt.domain, eps=eps, inf=inf, sup=sup, strict=strict, fast=fast)
result = []
for (s, t), indices in intervals: s, t = opt.domain.to_sympy(s), opt.domain.to_sympy(t) result.append(((s, t), indices))
return result
@public def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False): """ Refine an isolating interval of a root to the given precision.
Examples ========
>>> from sympy import refine_root >>> from sympy.abc import x
>>> refine_root(x**2 - 3, 1, 2, eps=1e-2) (19/11, 26/15)
""" try: F = Poly(f) except GeneratorsNeeded: raise PolynomialError( "can't refine a root of %s, not a polynomial" % f)
return F.refine_root(s, t, eps=eps, steps=steps, fast=fast, check_sqf=check_sqf)
@public def count_roots(f, inf=None, sup=None): """ Return the number of roots of ``f`` in ``[inf, sup]`` interval.
If one of ``inf`` or ``sup`` is complex, it will return the number of roots in the complex rectangle with corners at ``inf`` and ``sup``.
Examples ========
>>> from sympy import count_roots, I >>> from sympy.abc import x
>>> count_roots(x**4 - 4, -3, 3) 2 >>> count_roots(x**4 - 4, 0, 1 + 3*I) 1
""" try: F = Poly(f, greedy=False) except GeneratorsNeeded: raise PolynomialError("can't count roots of %s, not a polynomial" % f)
return F.count_roots(inf=inf, sup=sup)
@public def real_roots(f, multiple=True): """ Return a list of real roots with multiplicities of ``f``.
Examples ========
>>> from sympy import real_roots >>> from sympy.abc import x
>>> real_roots(2*x**3 - 7*x**2 + 4*x + 4) [-1/2, 2, 2]
""" except GeneratorsNeeded: raise PolynomialError( "can't compute real roots of %s, not a polynomial" % f)
@public def nroots(f, n=15, maxsteps=50, cleanup=True): """ Compute numerical approximations of roots of ``f``.
Examples ========
>>> from sympy import nroots >>> from sympy.abc import x
>>> nroots(x**2 - 3, n=15) [-1.73205080756888, 1.73205080756888] >>> nroots(x**2 - 3, n=30) [-1.73205080756887729352744634151, 1.73205080756887729352744634151]
""" try: F = Poly(f, greedy=False) except GeneratorsNeeded: raise PolynomialError( "can't compute numerical roots of %s, not a polynomial" % f)
return F.nroots(n=n, maxsteps=maxsteps, cleanup=cleanup)
@public def ground_roots(f, *gens, **args): """ Compute roots of ``f`` by factorization in the ground domain.
Examples ========
>>> from sympy import ground_roots >>> from sympy.abc import x
>>> ground_roots(x**6 - 4*x**4 + 4*x**3 - x**2) {0: 2, 1: 2}
""" options.allowed_flags(args, [])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('ground_roots', 1, exc)
return F.ground_roots()
@public def nth_power_roots_poly(f, n, *gens, **args): """ Construct a polynomial with n-th powers of roots of ``f``.
Examples ========
>>> from sympy import nth_power_roots_poly, factor, roots >>> from sympy.abc import x
>>> f = x**4 - x**2 + 1 >>> g = factor(nth_power_roots_poly(f, 2))
>>> g (x**2 - x + 1)**2
>>> R_f = [ (r**2).expand() for r in roots(f) ] >>> R_g = roots(g).keys()
>>> set(R_f) == set(R_g) True
""" options.allowed_flags(args, [])
try: F, opt = poly_from_expr(f, *gens, **args) except PolificationFailed as exc: raise ComputationFailed('nth_power_roots_poly', 1, exc)
result = F.nth_power_roots_poly(n)
if not opt.polys: return result.as_expr() else: return result
@public def cancel(f, *gens, **args): """ Cancel common factors in a rational function ``f``.
Examples ========
>>> from sympy import cancel, sqrt, Symbol >>> from sympy.abc import x >>> A = Symbol('A', commutative=False)
>>> cancel((2*x**2 - 2)/(x**2 - 2*x + 1)) (2*x + 2)/(x - 1) >>> cancel((sqrt(3) + sqrt(15)*A)/(sqrt(2) + sqrt(10)*A)) sqrt(6)/2 """
elif len(f) == 2: p, q = f elif isinstance(f, Tuple): return factor_terms(f) else: raise ValueError('unexpected argument: %s' % f)
except PolificationFailed: if not isinstance(f, (tuple, Tuple)): return f else: return S.One, p, q except PolynomialError as msg: if f.is_commutative and not f.has(Piecewise): raise PolynomialError(msg) # Handling of noncommutative and/or piecewise expressions if f.is_Add or f.is_Mul: sifted = sift(f.args, lambda x: x.is_commutative is True and not x.has(Piecewise)) c, nc = sifted[True], sifted[False] nc = [cancel(i) for i in nc] return f.func(cancel(f.func._from_args(c)), *nc) else: reps = [] pot = preorder_traversal(f) next(pot) for e in pot: # XXX: This should really skip anything that's not Expr. if isinstance(e, (tuple, Tuple, BooleanAtom)): continue try: reps.append((e, cancel(e))) pot.skip() # this was handled successfully except NotImplementedError: pass return f.xreplace(dict(reps))
else: if not opt.polys: return c, P.as_expr(), Q.as_expr() else: return c, P, Q
@public def reduced(f, G, *gens, **args): """ Reduces a polynomial ``f`` modulo a set of polynomials ``G``.
Given a polynomial ``f`` and a set of polynomials ``G = (g_1, ..., g_n)``, computes a set of quotients ``q = (q_1, ..., q_n)`` and the remainder ``r`` such that ``f = q_1*g_1 + ... + q_n*g_n + r``, where ``r`` vanishes or ``r`` is a completely reduced polynomial with respect to ``G``.
Examples ========
>>> from sympy import reduced >>> from sympy.abc import x, y
>>> reduced(2*x**4 + y**2 - x**2 + y**3, [x**3 - x, y**3 - y]) ([2*x, 1], x**2 + y**2 + y)
""" options.allowed_flags(args, ['polys', 'auto'])
try: polys, opt = parallel_poly_from_expr([f] + list(G), *gens, **args) except PolificationFailed as exc: raise ComputationFailed('reduced', 0, exc)
domain = opt.domain retract = False
if opt.auto and domain.has_Ring and not domain.has_Field: opt = opt.clone(dict(domain=domain.get_field())) retract = True
from sympy.polys.rings import xring _ring, _ = xring(opt.gens, opt.domain, opt.order)
for i, poly in enumerate(polys): poly = poly.set_domain(opt.domain).rep.to_dict() polys[i] = _ring.from_dict(poly)
Q, r = polys[0].div(polys[1:])
Q = [Poly._from_dict(dict(q), opt) for q in Q] r = Poly._from_dict(dict(r), opt)
if retract: try: _Q, _r = [q.to_ring() for q in Q], r.to_ring() except CoercionFailed: pass else: Q, r = _Q, _r
if not opt.polys: return [q.as_expr() for q in Q], r.as_expr() else: return Q, r
@public def groebner(F, *gens, **args): """ Computes the reduced Groebner basis for a set of polynomials.
Use the ``order`` argument to set the monomial ordering that will be used to compute the basis. Allowed orders are ``lex``, ``grlex`` and ``grevlex``. If no order is specified, it defaults to ``lex``.
For more information on Groebner bases, see the references and the docstring of `solve_poly_system()`.
Examples ========
Example taken from [1].
>>> from sympy import groebner >>> from sympy.abc import x, y
>>> F = [x*y - 2*y, 2*y**2 - x**2]
>>> groebner(F, x, y, order='lex') GroebnerBasis([x**2 - 2*y**2, x*y - 2*y, y**3 - 2*y], x, y, domain='ZZ', order='lex') >>> groebner(F, x, y, order='grlex') GroebnerBasis([y**3 - 2*y, x**2 - 2*y**2, x*y - 2*y], x, y, domain='ZZ', order='grlex') >>> groebner(F, x, y, order='grevlex') GroebnerBasis([y**3 - 2*y, x**2 - 2*y**2, x*y - 2*y], x, y, domain='ZZ', order='grevlex')
By default, an improved implementation of the Buchberger algorithm is used. Optionally, an implementation of the F5B algorithm can be used. The algorithm can be set using ``method`` flag or with the :func:`setup` function from :mod:`sympy.polys.polyconfig`:
>>> F = [x**2 - x - 1, (2*x - 1) * y - (x**10 - (1 - x)**10)]
>>> groebner(F, x, y, method='buchberger') GroebnerBasis([x**2 - x - 1, y - 55], x, y, domain='ZZ', order='lex') >>> groebner(F, x, y, method='f5b') GroebnerBasis([x**2 - x - 1, y - 55], x, y, domain='ZZ', order='lex')
References ==========
1. [Buchberger01]_ 2. [Cox97]_
"""
@public def is_zero_dimensional(F, *gens, **args): """ Checks if the ideal generated by a Groebner basis is zero-dimensional.
The algorithm checks if the set of monomials not divisible by the leading monomial of any element of ``F`` is bounded.
References ==========
David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties and Algorithms, 3rd edition, p. 230
"""
@public class GroebnerBasis(Basic): """Represents a reduced Groebner basis. """
def __new__(cls, F, *gens, **args): """Compute a reduced Groebner basis for a system of polynomials. """
except PolificationFailed as exc: raise ComputationFailed('groebner', len(F), exc)
@classmethod def _new(cls, basis, options):
@property def args(self): return (Tuple(*self._basis), Tuple(*self._options.gens))
@property def exprs(self): return [poly.as_expr() for poly in self._basis]
@property def polys(self):
@property def gens(self):
@property def domain(self): return self._options.domain
@property def order(self): return self._options.order
def __len__(self):
def __iter__(self): else: return iter(self.exprs)
def __getitem__(self, item): else: basis = self.exprs
def __hash__(self): return hash((self._basis, tuple(self._options.items())))
def __eq__(self, other): if isinstance(other, self.__class__): return self._basis == other._basis and self._options == other._options elif iterable(other): return self.polys == list(other) or self.exprs == list(other) else: return False
def __ne__(self, other): return not self.__eq__(other)
@property def is_zero_dimensional(self): """ Checks if the ideal generated by a Groebner basis is zero-dimensional.
The algorithm checks if the set of monomials not divisible by the leading monomial of any element of ``F`` is bounded.
References ==========
David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties and Algorithms, 3rd edition, p. 230
"""
# If any element of the exponents vector is zero, then there's # a variable for which there's no degree bound and the ideal # generated by this Groebner basis isn't zero-dimensional.
def fglm(self, order): """ Convert a Groebner basis from one ordering to another.
The FGLM algorithm converts reduced Groebner bases of zero-dimensional ideals from one ordering to another. This method is often used when it is infeasible to compute a Groebner basis with respect to a particular ordering directly.
Examples ========
>>> from sympy.abc import x, y >>> from sympy import groebner
>>> F = [x**2 - 3*y - x + 1, y**2 - 2*x + y - 1] >>> G = groebner(F, x, y, order='grlex')
>>> list(G.fglm('lex')) [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7] >>> list(groebner(F, x, y, order='lex')) [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7]
References ==========
J.C. Faugere, P. Gianni, D. Lazard, T. Mora (1994). Efficient Computation of Zero-dimensional Groebner Bases by Change of Ordering
""" opt = self._options
src_order = opt.order dst_order = monomial_key(order)
if src_order == dst_order: return self
if not self.is_zero_dimensional: raise NotImplementedError("can't convert Groebner bases of ideals with positive dimension")
polys = list(self._basis) domain = opt.domain
opt = opt.clone(dict( domain=domain.get_field(), order=dst_order, ))
from sympy.polys.rings import xring _ring, _ = xring(opt.gens, opt.domain, src_order)
for i, poly in enumerate(polys): poly = poly.set_domain(opt.domain).rep.to_dict() polys[i] = _ring.from_dict(poly)
G = matrix_fglm(polys, _ring, dst_order) G = [Poly._from_dict(dict(g), opt) for g in G]
if not domain.has_Field: G = [g.clear_denoms(convert=True)[1] for g in G] opt.domain = domain
return self._new(G, opt)
def reduce(self, expr, auto=True): """ Reduces a polynomial modulo a Groebner basis.
Given a polynomial ``f`` and a set of polynomials ``G = (g_1, ..., g_n)``, computes a set of quotients ``q = (q_1, ..., q_n)`` and the remainder ``r`` such that ``f = q_1*f_1 + ... + q_n*f_n + r``, where ``r`` vanishes or ``r`` is a completely reduced polynomial with respect to ``G``.
Examples ========
>>> from sympy import groebner, expand >>> from sympy.abc import x, y
>>> f = 2*x**4 - x**2 + y**3 + y**2 >>> G = groebner([x**3 - x, y**3 - y])
>>> G.reduce(f) ([2*x, 1], x**2 + y**2 + y) >>> Q, r = _
>>> expand(sum(q*g for q, g in zip(Q, G)) + r) 2*x**4 - x**2 + y**3 + y**2 >>> _ == f True
""" poly = Poly._from_expr(expr, self._options) polys = [poly] + list(self._basis)
opt = self._options domain = opt.domain
retract = False
if auto and domain.has_Ring and not domain.has_Field: opt = opt.clone(dict(domain=domain.get_field())) retract = True
from sympy.polys.rings import xring _ring, _ = xring(opt.gens, opt.domain, opt.order)
for i, poly in enumerate(polys): poly = poly.set_domain(opt.domain).rep.to_dict() polys[i] = _ring.from_dict(poly)
Q, r = polys[0].div(polys[1:])
Q = [Poly._from_dict(dict(q), opt) for q in Q] r = Poly._from_dict(dict(r), opt)
if retract: try: _Q, _r = [q.to_ring() for q in Q], r.to_ring() except CoercionFailed: pass else: Q, r = _Q, _r
if not opt.polys: return [q.as_expr() for q in Q], r.as_expr() else: return Q, r
def contains(self, poly): """ Check if ``poly`` belongs the ideal generated by ``self``.
Examples ========
>>> from sympy import groebner >>> from sympy.abc import x, y
>>> f = 2*x**3 + y**3 + 3*y >>> G = groebner([x**2 + y**2 - 1, x*y - 2])
>>> G.contains(f) True >>> G.contains(f + 1) False
""" return self.reduce(poly)[1] == 0
@public def poly(expr, *gens, **args): """ Efficiently transform an expression into a polynomial.
Examples ========
>>> from sympy import poly >>> from sympy.abc import x
>>> poly(x*(x**2 + x - 1)**2) Poly(x**5 + 2*x**4 - x**3 - 2*x**2 + x, x, domain='ZZ')
"""
poly_factors.append(_poly(factor, opt)) poly_factors.append( _poly(factor.base, opt).pow(factor.exp)) else:
else: product = poly_factors[0]
for factor in poly_factors[1:]: product = product.mul(factor)
if factors: factor = Mul(*factors)
if factor.is_Number: product = product.mul(factor) else: product = product.mul(Poly._from_expr(factor, opt))
poly_terms.append(product)
else: result = poly_terms[0]
for term in poly_terms[1:]: result = result.add(term)
if terms: term = Add(*terms)
if term.is_Number: result = result.add(term) else: result = result.add(Poly._from_expr(term, opt))
return Poly(expr, *gens, **args)
|