Polynomials with SymPy - Calculations with Python

Certainly one of the most abstract subjects in math class: polynomials. After all, let whoever has never been confused by a polynomial division cast the first variable: there is x everywhere you look. In this article, we will see how the SymPy library for Python helps us handle operations involving polynomials.

Creating polynomials manually

#

Let’s start by importing the library we will use throughout the article:

import sympy

# setup for better outputs in the article, can be ignored
sympy.init_printing(use_latex='png', scale=1.0, order='grlex',
                    forecolor='Black', backcolor='White',)

To begin, let’s use the most basic knowledge we acquired in the first article of the SymPy series. Let’s create a symbol for our variable and an expression that will represent a polynomial:

x = sympy.Symbol('x')

expr = (x - 1) * (x - 2) * (x - 3)

expr

\(\displaystyle \left(x - 3\right) \left(x - 2\right) \left(x - 1\right)\)

Note that we created our polynomial expression representing a polynomial in factored form. To obtain the expanded form, that is, the form where all the distributive multiplications are carried out, we use the expand method:

expr.expand()

\(\displaystyle x^{3} - 6 x^{2} + 11 x - 6\)

If we start with an expanded polynomial, we can factor it with the factor method:

_.factor()

\(\displaystyle \left(x - 3\right) \left(x - 2\right) \left(x - 1\right)\)

The factored form already makes it evident that the roots of the polynomial are 1, 2, and 3. But we can obtain the roots programmatically using SymPy’s solve:

roots = sympy.solve(expr, x)

roots

\(\displaystyle \left[ 1, \ 2, \ 3\right]\)

Checking equalities

#

Let’s create two new polynomial expressions:

P = (x - 5) * (x + 5)
Q = x**2 - 25

Note that P is simply the factored form of Q. Therefore, they are equivalent, and if we check the equality between them the result will be True, right?

P == Q

False

Huh…?! Well, we have already written about this equality issue in SymPy. The library checks equality of form, not mathematical equivalence. Therefore, since the forms of the two polynomials are different, the comparison returns False.

Likewise, we know that mathematically the difference between the polynomials is zero, since they are equivalent. But the result of the following comparison is also False:

P - Q == 0

False

Why? Because SymPy is not proactive about simplifications, as we have already seen in other articles. Let’s look at the representation of the difference:

P - Q

\(\displaystyle - x^{2} + \left(x - 5\right) \left(x + 5\right) + 25\)

Note that SymPy indeed stores the expression without simplifying it. If we request the simplification with the simplify method, we see that the difference does equal zero:

sympy.simplify(P - Q) == 0

True

And we can see that one polynomial is indeed the factored version of the other with the following comparisons:

P.simplify() == Q

True

P == Q.factor()

True

Using the polynomial constructor

#

So far, we have seen little that is new compared to the article about expressions with SymPy. But let’s change that. Let’s look at the type of the object Q:

Q

\(\displaystyle x^{2} - 25\)

type(Q)

sympy.core.add.Add

Note that, internally, it is an object of SymPy’s Add class. That is, for SymPy such an object is simply the result of adding two other SymPy objects, which is correct. The one giving the object the meaning of a polynomial is us. Let’s unpack that last sentence.

The great advantage of using a high-level language such as Python is that we can write code that is closer to how we think, by using abstractions. The language interpreter takes care of translating such abstractions for the machine. So far we have used mere expressions to represent polynomials, but ideally SymPy would recognize such expressions as polynomials and not as additions. And that is possible, and we will see the advantages of this approach.

Expressions have the as_poly method:

Q.as_poly()

\(\displaystyle \operatorname{Poly}{\left( x^{2} - 25, x, domain=\mathbb{Z} \right)}\)

Look at the result! Now SymPy recognizes it as a Poly type, recognizes that the variable is x and that the domain is the integers (it can be changed; it is just the default when we do not pass the domain explicitly). Let’s assign this result to a variable and confirm the type:

Q_poly = Q.as_poly()

type(Q_poly)

sympy.polys.polytools.Poly

And what is the advantage? Now we can ask for what we want with methods closer to the semantics we use in mathematics. For example, we can ask for the coefficients of the polynomial:

Q_poly.coeffs()

\(\displaystyle \left[ 1, \ -25\right]\)

The result goes from the highest degree in x to the lowest. However, notice that the result can be misleading. After all, we know there is a term of degree 1 in x; it is just that its coefficient is zero. coeffs returns only nonzero coefficients. To avoid confusion, I recommend always using the all_coeffs method:

Q_poly.all_coeffs()

\(\displaystyle \left[ 1, \ 0, \ -25\right]\)

Now we have it, all the coefficients.

Another concept that always shows up in the study of polynomials: the degree of the polynomial. Now that we have an abstraction closer to our language, we can simply ask for the degree with the degree method:

Q_poly.degree()

\(\displaystyle 2\)

We can ask for the list of factors of the polynomial:

Q_poly.factor_list()

\(\displaystyle \left( 1, \ \left[ \left( \operatorname{Poly}{\left( x - 5, x, domain=\mathbb{Z} \right)}, \ 1\right), \ \left( \operatorname{Poly}{\left( x + 5, x, domain=\mathbb{Z} \right)}, \ 1\right)\right]\right)\)

Ask for the roots:

Q_poly.all_roots()

\(\displaystyle \left[ -5, \ 5\right]\)

Derivatives and integrals

#

We can even do more elaborate operations. Moving a bit beyond high school math, it is very common to need to differentiate or integrate a polynomial. Now, such operations become trivial. Let’s look at the first derivative of our polynomial:

Q_poly.diff()

\(\displaystyle \operatorname{Poly}{\left( 2 x, x, domain=\mathbb{Z} \right)}\)

We can obtain the value of the derivative at a given value of x using the subs method already seen in other articles:

Q_poly.diff().subs({x: 2})

\(\displaystyle 4\)

Or, more simply, with the eval method, passing the value of x:

Q_poly.diff().eval(2)

\(\displaystyle 4\)

For higher-order derivatives, just pass a tuple with the variable and the desired order:

Q_poly.diff((x, 2))

\(\displaystyle \operatorname{Poly}{\left( 2, x, domain=\mathbb{Z} \right)}\)

Q_poly.diff((x, 3))

\(\displaystyle \operatorname{Poly}{\left( 0, x, domain=\mathbb{Z} \right)}\)

To obtain the integral, we use the integrate method:

Q_poly.integrate()

\(\displaystyle \operatorname{Poly}{\left( \frac{1}{3} x^{3} - 25 x, x, domain=\mathbb{Q} \right)}\)

Likewise, we can obtain the value of the integral at a given value of x with the subs method or the eval method:

Q_poly.integrate().subs({x: 1})

\(\displaystyle - \frac{74}{3}\)

Q_poly.integrate().eval(1)

\(\displaystyle - \frac{74}{3}\)

The Poly class and complex roots

#

We have already seen that SymPy has the Poly class. We can pass an expression directly to the class instead of using the as_poly method:

expr = x**3 + 2*x + 3

A = sympy.Poly(expr)
A

\(\displaystyle \operatorname{Poly}{\left( x^{3} + 2 x + 3, x, domain=\mathbb{Z} \right)}\)

type(A)

sympy.polys.polytools.Poly

Notice that our polynomial A has degree three. Let’s look at all of its roots:

A.all_roots()

\(\displaystyle \left[ -1, \ \frac{1}{2} - \frac{\sqrt{11} i}{2}, \ \frac{1}{2} + \frac{\sqrt{11} i}{2}\right]\)

We have complex roots now. Let’s look at the result of the is_real attribute for each root:

for root in A.all_roots():
    print(root.is_real)

True
False
False

Makes sense: the first root is real, the others are complex.

This explains the existence of an all_roots method. After all, if there is a method that takes care to return all the roots, there must be some method that returns only some of them, right? Depending on the application, we may be interested only in the real roots. Then we can use the real_roots method:

A.real_roots()

\(\displaystyle \left[ -1\right]\)

Polynomial division

#

Let’s check two Poly objects that we created along the article:

A

\(\displaystyle \operatorname{Poly}{\left( x^{3} + 2 x + 3, x, domain=\mathbb{Z} \right)}\)

Q_poly

\(\displaystyle \operatorname{Poly}{\left( x^{2} - 25, x, domain=\mathbb{Z} \right)}\)

To make things easier to follow, let’s assign a variable B to Q_poly:

B = Q_poly
B

\(\displaystyle \operatorname{Poly}{\left( x^{2} - 25, x, domain=\mathbb{Z} \right)}\)

In the intro of the article, we talked about polynomial division. The time has come: let’s see how to do that operation with SymPy.

It is quite simple, we call the div method passing the dividend and the divisor:

sympy.div(A, B)

\(\displaystyle \left( \operatorname{Poly}{\left( x, x, domain=\mathbb{Z} \right)}, \ \operatorname{Poly}{\left( 27 x + 3, x, domain=\mathbb{Z} \right)}\right)\)

The result is a tuple with two polynomials, the quotient and the remainder. We can obtain each one separately by unpacking the resulting tuple:

quotient, remainder = sympy.div(A, B)

quotient

\(\displaystyle \operatorname{Poly}{\left( x, x, domain=\mathbb{Z} \right)}\)

remainder

\(\displaystyle \operatorname{Poly}{\left( 27 x + 3, x, domain=\mathbb{Z} \right)}\)

Since this is a division, we expect that divisor × quotient + remainder = dividend. Let’s see:

B * quotient + remainder

\(\displaystyle \operatorname{Poly}{\left( x^{3} + 2 x + 3, x, domain=\mathbb{Z} \right)}\)

A == B * quotient + remainder

True

Simple.

Conclusion and more articles about SymPy

#

In this article we saw the basics of how to handle polynomials in SymPy, and I hope it was possible to see the power of this library.

See you next time.