Basic functionality of the module

Introduction

This tutorial tries to give an overview of the functionality concerning polynomials within SymPy. All code examples assume:

>>> from sympy import *
>>> x, y, z = symbols('x,y,z')
>>> init_printing(use_unicode=False)

Basic concepts

Polynomials

Given a family $(x_i)$ of symbols, or other suitable objects, including numbers, expressions derived from them by repeated addition, subtraction and multiplication are called polynomial expressions in the generators $x_i$.

By the distributive law it is possible to perform multiplications before additions and subtractions. The products of generators thus obtained are called monomials. They are usually written in the form $x_1^{{\nu }_1}{x}_2^{{\nu }_2}\cdots x_n^{{\nu }_n}$ where the exponents ${\nu }_i$ are nonnegative integers. It is often convenient to write this briefly as $x^{\nu }$ where $x=(x_1,x_2,\dots ,x_n)$ denotes the family of generators and $\nu =({\nu }_1,{\nu }_2,\dots ,{\nu }_n)$ is the family of exponents.

When all monomials having the same exponents are combined, the polynomial expression becomes a sum of products $c_{\nu }{x}^{\nu }$, called the terms of the polynomial, where the coefficients $c_{\nu }$ are integers. If some of the $x_i$ are manifest numbers, they are incorporated in the coefficients and not regarded as generators. Such coefficients are typically rational, real or complex numbers. Some symbolic numbers, e.g., pi, can be either coefficients or generators.

A polynomial expression that is a sum of terms with different monomials is uniquely determined by its family of coefficients $(c_{\nu })$. Such an expression is customarily called a polynomial, though, more properly, that name does stand for the coefficient family once the generators are given. SymPy implements polynomials by default as dictionaries with monomials as keys and coefficients as values. Another implementation consists of nested lists of coefficients.

The set of all polynomials with integer coefficients in the generators $x_i$ is a ring, i.e., the sums, differences and products of its elements are again polynomials in the same generators. This ring is denoted $\mathbb{Z}[x_1,x_2,\dots ,x_n]$, or $\mathbb{Z}[(x_i)]$, and called the ring of polynomials in the $x_i$ with integer coefficients.

More generally, the coefficients of a polynomial can be elements of any commutative ring $A$, and the corresponding polynomial ring is then denoted $A[x_1,x_2,\dots ,x_n]$. The ring $A$ can also be a polynomial ring. In SymPy, the coefficient ring is called the domain of the polynomial ring, and it can be given as a keyword parameter. By default, it is determined by the coefficients of the polynomial arguments.

Polynomial expressions can be transformed into polynomials by the method sympy.core.expr.Expr.as_poly:

>>> e = (x + y)*(y - 2*z)
>>> e.as_poly()
Poly(x*y - 2*x*z + y**2 - 2*y*z, x, y, z, domain='ZZ')

If a polynomial expression contains numbers that are not integers, they are regarded as coefficients and the coefficient ring is extended accordingly. In particular, division by integers leads to rational coefficients:

>>> e = (3*x/2 + y)*(z - 1)
>>> e.as_poly()
Poly(3/2*x*z - 3/2*x + y*z - y, x, y, z, domain='QQ')

Symbolic numbers are considered generators unless they are explicitly excluded, in which case they are adjoined to the coefficient ring:

>>> e = (x + 2*pi)*y
>>> e.as_poly()
Poly(x*y + 2*y*pi, x, y, pi, domain='ZZ')
>>> e.as_poly(x, y)
Poly(x*y + 2*pi*y, x, y, domain='ZZ[pi]')

Alternatively, the coefficient domain can be specified by means of a keyword argument:

>>> e = (x + 2*pi)*y
>>> e.as_poly(domain=ZZ[pi])
Poly(x*y + 2*pi*y, x, y, domain='ZZ[pi]')

Note that the ring $\mathbb{Z}[\pi ][x,y]$ of polynomials in $x$ and $y$ with coefficients in $\mathbb{Z}[\pi ]$ is mathematically equivalent to $\mathbb{Z}[\pi ,x,y]$, only their implementations differ.

If an expression contains functions of the generators, other than their positive integer powers, these are interpreted as new generators:

>>> e = x*sin(y) - y
>>> e.as_poly()
Poly(x*(sin(y)) - y, x, y, sin(y), domain='ZZ')

Since $y$ and $\sin (y)$ are algebraically independent they can both appear as generators in a polynomial. However, polynomial expressions must not contain negative powers of generators:

>>> e = x - 1/x
>>> e.as_poly()
Poly(x - (1/x), x, 1/x, domain='ZZ')

It is important to realize that the generators $x$ and $1/x=x^{-1}$ are treated as algebraically independent variables. In particular, their product is not equal to 1. Hence generators in denominators should be avoided even if they raise no error in the current implementation. This behavior is undesirable and may change in the future. Similar problems emerge with rational powers of generators. So, for example, $x$ and $\sqrt{x}=x^{1/2}$ are not recognized as algebraically dependent.

If there are algebraic numbers in an expression, it is possible to adjoin them to the coefficient ring by setting the keyword extension:

>>> e = x + sqrt(2)
>>> e.as_poly()
Poly(x + (sqrt(2)), x, sqrt(2), domain='ZZ')
>>> e.as_poly(extension=True)
Poly(x + sqrt(2), x, domain='QQ<sqrt(2)>')

With the default setting extension=False, both $x$ and $\sqrt{2}$ are incorrectly considered algebraically independent variables. With coefficients in the extension field $\mathbb{Q}(\sqrt{2})$ the square root is treated properly as an algebraic number. Setting extension=True whenever algebraic numbers are involved is definitely recommended even though it is not forced in the current implementation.

Divisibility

The fourth rational operation, division, or inverted multiplication, is not generally possible in rings. If $a$ and $b$ are two elements of a ring $A$, then there may exist a third element $q$ in $A$ such that $a=bq$. In fact, there may exist several such elements.

If also $a=b{q}^{\prime }$ for some $q^{\prime }$ in $A$, then $b(q-q^{\prime })=0$. Hence either $b$ or $q-q^{\prime }$ is zero, or they are both zero divisors, nonzero elements whose product is zero.

Integral domains

Commutative rings with no zero divisors are called integral domains. Most of the commonly encountered rings, the ring of integers, fields, and polynomial rings over integral domains are integral domains.

Assume now that $A$ is an integral domain, and consider the set $P$ of its nonzero elements, which is closed under multiplication. If $a$ and $b$ are in $P$, and there exists an element $q$ in $P$ such that $a=bq$, then $q$ is unique and called the quotient, $a/b$, of $a$ by $b$. Moreover, it is said that

• $a$ is divisible by $b$,

• $b$ is a divisor of $a$,

• $a$ is a multiple of $b$,

• $b$ is a factor of $a$.

An element $a$ of $P$ is a divisor of $1$ if and only if it is invertible in $A$, with the inverse $a^{-1}=1/a$. Such elements are called units. The units of the ring of integers are $1$ and $-1$. The invertible elements in a polynomial ring over a field are the nonzero constant polynomials.

If two elements of $P$, $a$ and $b$, are divisible by each other, then the quotient $a/b$ is invertible with inverse $b/a$, or equivalently, $b=ua$ where $u$ is a unit. Such elements are said to be associated with, or associates of, each other. The associates of an integer $n$ are $n$ and $-n$. In a polynomial ring over a field the associates of a polynomial are its constant multiples.

Each element of $P$ is divisible by its associates and the units. An element is irreducible if it has no other divisors and is not a unit. The irreducible elements in the ring of integers are the prime numbers $p$ and their opposites $-p$. In a field, every nonzero element is invertible and there are no irreducible elements.

Factorial domains

In the ring of integers, each nonzero element can be represented as a product of irreducible elements and optionally a unit $\pm 1$. Moreover, any two such products have the same number of irreducible factors which are associated with each other in a suitable order. Integral domains having this property are called factorial, or unique factorization domains. In addition to the ring of integers, all polynomial rings over a field are factorial, and so are more generally polynomial rings over any factorial domain. Fields are trivially factorial since there are only units. The irreducible elements of a factorial domain are usually called primes.

A family of integers has only a finite number of common divisors and the greatest of them is divisible by all of them. More generally, given a family of nonzero elements $(a_i)$ in an integral domain, a common divisor $d$ of the elements is called a greatest common divisor, abbreviated gcd, of the family if it is a multiple of all common divisors. A greatest common divisor, if it exists, is not unique in general; all of its associates have the same property. It is denoted by $d=gcd(a_1,\dots ,a_n)$ if there is no danger of confusion. A least common multiple, or lcm, of a family $(a_i)$ is defined analogously as a common multiple $m$ that divides all common multiples. It is denoted by $m=\mathrm{lcm}(a_1,\dots ,a_n)$.

In a factorial domain, greatest common divisors always exists. They can be found, at least in principle, by factoring each element of a family into a product of prime powers and an optional unit, and, for each prime, taking the least power that appears in the factorizations. The product of these prime powers is then a greatest common divisor. A least common multiple can be obtained from the same factorizations as the product of the greatest powers for each prime.

Euclidean domains

A practical algorithm for computing a greatest common divisor can be implemented in Euclidean domains. They are integral domains that can be endowed with a function $w$ assigning a nonnegative integer to each nonzero element of the domain and having the following property:

if $a$ and $b$ are nonzero, there are $q$ and $r$ that satisfy the division identity

$a=qb+r$

such that either $r=0$ or $w(r)<w(b)$.

The ring of integers and all univariate polynomial rings over fields are Euclidean domains with $w(a)=|a|$ resp. $w(a)=\mathrm{deg}(a)$.

The division identity for integers is implemented in Python as the built-in function divmod that can also be applied to SymPy Integers:

>>> divmod(Integer(53), Integer(7))
(7, 4)

For polynomials the division identity is given in SymPy by the function div():

>>> f = 5*x**2 + 10*x + 3
>>> g = 2*x + 2

>>> q, r = div(f, g, domain='QQ')
>>> q
5*x   5
--- + -
 2    2
>>> r
-2
>>> (q*g + r).expand()
   2
5*x  + 10*x + 3

The division identity can be used to determine the divisibility of elements in a Euclidean domain. If $r=0$ in the division identity, then $a$ is divisible by $b$. Conversely, if $a=cb$ for some element $c$, then $(c-q)b=r$. It follows that $c=q$ and $r=0$ if $w$ has the additional property:

if $a$ and $b$ are nonzero, then $w(ab)\ge w(b)$.

This is satisfied by the functions given above. (And it is always possible to redefine $w(a)$ by taking the minimum of the values $w(xa)$ for $x\ne 0$.)

The principal application of the division identity is the efficient computation of a greatest common divisor by means of the Euclidean algorithm. It applies to two elements of a Euclidean domain. A gcd of several elements can be obtained by iteration.

The function for computing the greatest common divisor of integers in SymPy is currently igcd():

>>> igcd(2, 4)
2
>>> igcd(5, 10, 15)
5

For univariate polynomials over a field the function has its common name gcd(), and the returned polynomial is monic:

>>> f = 4*x**2 - 1
>>> g = 8*x**3 + 1
>>> gcd(f, g, domain=QQ)
x + 1/2

Divisibility of polynomials

The ring $A=\mathbb{Z}[x]$ of univariate polynomials over the ring of integers is not Euclidean but it is still factorial. To see this, consider the divisibility in $A$.

Let $f$ and $g$ be two nonzero polynomials in $A$. If $f$ is divisible by $g$ in $A$, then it is also divisible in the ring $B=\mathbb{Q}[x]$ of polynomials with rational coefficients. Since $B$ is Euclidean, this can be determined by means of the division identity.

Assume, conversely, that $f=gh$ for some polynomial $h$ in $B$. Then $f$ is divisible by $g$ in $A$ if and only if the coefficients of $h$ are integers. To find out when this is true it is necessary to consider the divisibility of the coefficients.

For a polynomial $f$ in $A$, let $c$ be the greatest common divisor of its coefficients. Then $f$ is divisible by the constant polynomial $c$ in $A$, and the quotient $f/c=p$ is a polynomial whose coefficients are integers that have no common divisor apart from the units. Such polynomials are called primitive. A polynomial with rational coefficients can also be written as $f=cp$, where $c$ is a rational number and $p$ is a primitive polynomial. The constant $c$ is called the content of $f$, and $p$ is its primitive part. These components can be found by the method sympy.core.expr.Expr.as_content_primitive:

>>> f = 6*x**2 - 3*x + 9
>>> c, p = f.as_content_primitive()
>>> c, p
       2
(3, 2*x  - x + 3)
>>> f = x**2/3 - x/2 + 1
>>> c, p = f.as_content_primitive()
>>> c, p
         2
(1/6, 2*x  - 3*x + 6)

Let $f$, $f^{\prime }$ be polynomials with contents $c$, $c^{\prime }$ and primitive parts $p$, $p^{\prime }$. Then $f{f}^{\prime }=(c{c}^{\prime })(p{p}^{\prime })$ where the product $p{p}^{\prime }$ is primitive by Gauss’s lemma. It follows that

the content of a product of polynomials is the product of their contents and the primitive part of the product is the product of the primitive parts.

Returning to the divisibility in the ring $\mathbb{Z}[x]$, assume that $f$ and $g$ are two polynomials with integer coefficients such that the division identity in $\mathbb{Q}[x]$ yields the equality $f=gh$ for some polynomial $h$ with rational coefficients. Then the content of $f$ is equal to the content of $g$ multiplied by the content of $h$. As $h$ has integer coefficients if and only if its content is an integer, we get the following criterion:

$f$ is divisible by $g$ in the ring $\mathbb{Z}[x]$ if and only if

1. $f$ is divisible by $g$ in $\mathbb{Q}[x]$, and

2. the content of $f$ is divisible by the content of $g$ in $\mathbb{Z}$.

If $f=cp$ is irreducible in $\mathbb{Z}[x]$, then either $c$ or $p$ must be a unit. If $p$ is not a unit, it must be irreducible also in $\mathbb{Q}[x]$. For if it is a product of two polynomials, it is also the product of their primitive parts, and one of them must be a unit. Hence there are two kinds of irreducible elements in $\mathbb{Z}[x]$:

1. prime numbers of $\mathbb{Z}$, and

2. primitive polynomials that are irreducible in $\mathbb{Q}[x]$.

It follows that each polynomial in $\mathbb{Z}[x]$ is a product of irreducible elements. It suffices to factor its content and primitive part separately. These products are essentially unique; hence $\mathbb{Z}[x]$ is also factorial.

Another important consequence is that a greatest common divisor of two polynomials in $\mathbb{Z}[x]$ can be found efficiently by applying the Euclidean algorithm separately to their contents and primitive parts in the Euclidean domains $\mathbb{Z}$ and $\mathbb{Q}[x]$. This is also implemented in SymPy:

>>> f = 4*x**2 - 1
>>> g = 8*x**3 + 1
>>> gcd(f, g)
2*x + 1
>>> gcd(6*f, 3*g)
6*x + 3

Basic functionality

These functions provide different algorithms dealing with polynomials in the form of SymPy expression, like symbols, sums etc.

Division

The function div() provides division of polynomials with remainder. That is, for polynomials f and g, it computes q and r, such that $f=g\cdot q+r$ and $\mathrm{deg}(r)<\mathrm{deg}(q)$. For polynomials in one variables with coefficients in a field, say, the rational numbers, q and r are uniquely defined this way:

>>> f = 5*x**2 + 10*x + 3
>>> g = 2*x + 2

>>> q, r = div(f, g, domain='QQ')
>>> q
5*x   5
--- + -
 2    2
>>> r
-2
>>> (q*g + r).expand()
   2
5*x  + 10*x + 3

As you can see, q has a non-integer coefficient. If you want to do division only in the ring of polynomials with integer coefficients, you can specify an additional parameter:

>>> q, r = div(f, g, domain='ZZ')
>>> q
0
>>> r
   2
5*x  + 10*x + 3

But be warned, that this ring is no longer Euclidean and that the degree of the remainder doesn’t need to be smaller than that of f. Since 2 doesn’t divide 5, $2x$ doesn’t divide $5x^2$, even if the degree is smaller. But:

>>> g = 5*x + 1

>>> q, r = div(f, g, domain='ZZ')
>>> q
x
>>> r
9*x + 3
>>> (q*g + r).expand()
   2
5*x  + 10*x + 3

This also works for polynomials with multiple variables:

>>> f = x*y + y*z
>>> g = 3*x + 3*z

>>> q, r = div(f, g, domain='QQ')
>>> q
y
-
3
>>> r
0

In the last examples, all of the three variables x, y and z are assumed to be variables of the polynomials. But if you have some unrelated constant as coefficient, you can specify the variables explicitly:

>>> a, b, c = symbols('a,b,c')
>>> f = a*x**2 + b*x + c
>>> g = 3*x + 2
>>> q, r = div(f, g, domain='QQ')
>>> q
a*x   2*a   b
--- - --- + -
 3     9    3

>>> r
4*a   2*b
--- - --- + c
 9     3

GCD and LCM

With division, there is also the computation of the greatest common divisor and the least common multiple.

When the polynomials have integer coefficients, the contents’ gcd is also considered:

>>> f = (12*x + 12)*x
>>> g = 16*x**2
>>> gcd(f, g)
4*x

But if the polynomials have rational coefficients, then the returned polynomial is monic:

>>> f = 3*x**2/2
>>> g = 9*x/4
>>> gcd(f, g)
x

It also works with multiple variables. In this case, the variables are ordered alphabetically, be default, which has influence on the leading coefficient:

>>> f = x*y/2 + y**2
>>> g = 3*x + 6*y

>>> gcd(f, g)
x + 2*y

The lcm is connected with the gcd and one can be computed using the other:

>>> f = x*y**2 + x**2*y
>>> g = x**2*y**2
>>> gcd(f, g)
x*y
>>> lcm(f, g)
 3  2    2  3
x *y  + x *y
>>> (f*g).expand()
 4  3    3  4
x *y  + x *y
>>> (gcd(f, g, x, y)*lcm(f, g, x, y)).expand()
 4  3    3  4
x *y  + x *y

Square-free factorization

The square-free factorization of a univariate polynomial is the product of all factors (not necessarily irreducible) of degree 1, 2 etc.:

>>> f = 2*x**2 + 5*x**3 + 4*x**4 + x**5

>>> sqf_list(f)
                   2
(1, [(x + 2, 1), (x  + x, 2)])

>>> sqf(f)
                2
        / 2    \
(x + 2)*\x  + x/

Factorization

This function provides factorization of univariate and multivariate polynomials with rational coefficients:

>>> factor(x**4/2 + 5*x**3/12 - x**2/3)
 2
x *(2*x - 1)*(3*x + 4)
----------------------
          12

>>> factor(x**2 + 4*x*y + 4*y**2)
         2
(x + 2*y)

Groebner bases

Buchberger’s algorithm is implemented, supporting various monomial orders:

>>> groebner([x**2 + 1, y**4*x + x**3], x, y, order='lex')
             /[ 2       4    ]                            \
GroebnerBasis\[x  + 1, y  - 1], x, y, domain=ZZ, order=lex/


>>> groebner([x**2 + 1, y**4*x + x**3, x*y*z**3], x, y, z, order='grevlex')
             /[ 4       3   2    ]                                   \
GroebnerBasis\[y  - 1, z , x  + 1], x, y, z, domain=ZZ, order=grevlex/

Solving Equations

We have (incomplete) methods to find the complex or even symbolic roots of polynomials and to solve some systems of polynomial equations:

>>> from sympy import roots, solve_poly_system

>>> solve(x**3 + 2*x + 3, x)
           ____          ____
     1   \/ 11 *I  1   \/ 11 *I
[-1, - - --------, - + --------]
     2      2      2      2

>>> p = Symbol('p')
>>> q = Symbol('q')

>>> solve(x**2 + p*x + q, x)
          __________           __________
         /  2                 /  2
   p   \/  p  - 4*q     p   \/  p  - 4*q
[- - - -------------, - - + -------------]
   2         2          2         2

>>> solve_poly_system([y - x, x - 5], x, y)
[(5, 5)]

>>> solve_poly_system([y**2 - x**3 + 1, y*x], x, y)
                                   ___                 ___
                             1   \/ 3 *I         1   \/ 3 *I
[(0, -I), (0, I), (1, 0), (- - - -------, 0), (- - + -------, 0)]
                             2      2            2      2