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, wrap_line=False)
Basic concepts#
Polynomials#
Given a family 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
.
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 where the exponents
are nonnegative integers. It is often
convenient to write this briefly as
where
denotes the family of generators and
is
the family of exponents.
When all monomials having the same exponents are combined, the polynomial
expression becomes a sum of products , called the terms of the
polynomial, where the coefficients
are integers. If some of the
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 . 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 is
a ring, i.e., the sums, differences and products of its elements are again
polynomials in the same generators. This ring is denoted
, or
, and called the ring of polynomials in
the
with integer coefficients.
More generally, the coefficients of a polynomial can be elements of any
commutative ring , and the corresponding polynomial ring is then denoted
. The ring
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 of polynomials in
and
with
coefficients in
is mathematically equivalent to
, 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 and
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 and
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,
and
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 and
are
incorrectly considered algebraically independent variables. With coefficients in
the extension field
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 and
are two elements of a ring
, then
there may exist a third element
in
such that
. In fact, there
may exist several such elements.
If also for some
in
, then
. Hence either
or
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 is an integral domain, and consider the set
of its
nonzero elements, which is closed under multiplication. If
and
are in
, and there exists an element
in
such that
, then
is
unique and called the quotient,
, of
by
. Moreover, it is said
that
is divisible by
,
is a divisor of
,
is a multiple of
,
is a factor of
.
An element of
is a divisor of
if and only if it is invertible in
, with the inverse
. Such elements are called units. The
units of the ring of integers are
and
. The invertible elements in a
polynomial ring over a field are the nonzero constant polynomials.
If two elements of ,
and
, are divisible by each other, then the
quotient
is invertible with inverse
, or equivalently,
where
is a unit. Such elements are said to be associated with, or associates
of, each other. The associates of an integer
are
and
. In a
polynomial ring over a field the associates of a polynomial are its constant
multiples.
Each element of 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
and their opposites
. 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 . 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 in an integral domain, a common divisor
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
if there is no danger of confusion.
A least common multiple, or lcm, of a family
is defined analogously
as a common multiple
that divides all common multiples. It is denoted by
.
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 assigning a nonnegative integer to each nonzero element of the
domain and having the following property:
if
and
are nonzero, there are
and
that satisfy the division identity
such that either
or
.
The ring of integers and all univariate polynomial rings over fields are
Euclidean domains with resp.
.
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 in the division identity, then
is divisible by
. Conversely, if
for some element
, then
. It
follows that
and
if
has the additional property:
if
and
are nonzero, then
.
This is satisfied by the functions given above. (And it is always possible to
redefine by taking the minimum of the values
for
.)
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 of univariate polynomials over the ring of integers
is not Euclidean but it is still factorial. To see this, consider the
divisibility in
.
Let and
be two nonzero polynomials in
. If
is divisible by
in
, then it is also divisible in the ring
of polynomials
with rational coefficients. Since
is Euclidean, this can be determined by
means of the division identity.
Assume, conversely, that for some polynomial
in
. Then
is
divisible by
in
if and only if the coefficients of
are integers. To
find out when this is true it is necessary to consider the divisibility of the
coefficients.
For a polynomial in
, let
be the greatest common divisor of its
coefficients. Then
is divisible by the constant polynomial
in
, and
the quotient
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
, where
is a rational number and
is a primitive polynomial. The constant
is
called the content of
, and
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 ,
be polynomials with contents
,
and primitive parts
,
. Then
where the product
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 , assume that
and
are two polynomials with integer coefficients such that the division
identity in
yields the equality
for some polynomial
with rational coefficients. Then the content of
is equal to the content of
multiplied by the content of
. As
has integer coefficients if and
only if its content is an integer, we get the following criterion:
is divisible by
in the ring
if and only if
is divisible by
in
, and
the content of
is divisible by the content of
in
.
If is irreducible in
, then either
or
must be a
unit. If
is not a unit, it must be irreducible also in
. 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
:
prime numbers of
, and
primitive polynomials that are irreducible in
.
It follows that each polynomial in is a product of irreducible
elements. It suffices to factor its content and primitive part separately. These
products are essentially unique; hence
is also factorial.
Another important consequence is that a greatest common divisor of two
polynomials in can be found efficiently by applying the
Euclidean algorithm separately to their contents and primitive parts in the
Euclidean domains
and
. 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 and
. 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, doesn’t divide
, 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