Polynomial Manipulation with SymPy

Featured image SymPy Polynomials Python

Certainly, one of the most abstract subjects in math classes: polynomials. After all, who has never been confused in a polynomial division, it is x everywhere. In this article, we will see how the SymPy library, for the Python language, helps us to deal with operations involving polynomials.

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

import sympy

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

Now, we are going to 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

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

expr.expand()

If you initially have an expanded polynomial, you can factor it with the factor method:

_.factor()

The factored form already points out 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

Checking equalities

Moving on, let’s create two new polynomial expressions:

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

Observe 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

Uh…?! Well, we have already written about this issue of equality in SymPy. The library analyzes equality in form and not mathematical equivalence. Therefore, as the forms of the two polynomials are distinct, the comparison returns False.

P - Q == 0
False

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

P - Q

See that SymPy effectively evaluates the expression without simplifying it. If we request the simplification with the simplify method, we see that the difference effectively equals zero:

sympy.simplify(P - Q) == 0
True

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

P.simplify() == Q
True
P == Q.factor()
True

Using the polynomial constructor

Until now, we have seen little news when compared to the article about expressions with SymPy. But let’s change that. Let’s see the type of the object Q:

Q
type(Q)
sympy.core.add.Add

Note that, internally, it is an object of the Add class of SymPy. That is, for SymPy, such object is simply the result of adding two other SymPy objects, which is correct. We are the ones giving meaning to the object as a polynomial. Let’s understand this last sentence.

The great advantage of using a high-level language like Python is that we can write code that is closer to our understanding using abstractions. The language interpreter is responsible for translating these abstractions into machine code. So far, we have used mere expressions as representations of polynomials, but it would be ideal for SymPy to effectively recognize such expressions as polynomials and not as additions. And this is possible and we will see the advantages of this approach.

Expressions have the as_poly method:

Q.as_poly()

Behold the return! Now SymPy recognizes it as a Poly type, recognizes that the variable is x, and that the domain is of integers (it can be changed, it is just the default when we do not explicitly pass the domain). Let’s associate this return 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 request what we want with methods closer to the semantics we use in mathematics. For example, we can request the coefficients of the polynomial:

Q_poly.coeffs()

The return is from the highest degree in x to the lowest. However, note that the return can be misleading. After all, we know that there is a term of degree 1 in x, only its coefficient is zero. The coeffs returns only non-zero coefficients. To avoid misunderstandings, I always recommend using the all_coeffs method:

Q_poly.all_coeffs()

Now, all coefficients.

Another concept that always appears in the study of polynomials: the degree of the polynomial. Now that we have an abstraction closer to our language, we can effectively request the degree with the degree method:

Q_poly.degree()

We can request the list of factors of the polynomial:

Q_poly.factor_list()

Also, we can request the roots:

Q_poly.all_roots()

Derivatives and integrals

In fact, we can perform more elaborate operations. Leaving high school mathematics for a moment, it is very common to need to derive or integrate a polynomial. Now, such operations become trivial. Let’s see the first derivative of our polynomial:

Q_poly.diff()

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

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

Or, simpler, with the eval method and passing the value of x:

Q_poly.diff().eval(2)

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

Q_poly.diff((x, 2))
Q_poly.diff((x, 3))

With the integrate method, we can obtain the indefinite integral of the polynomial:

Q_poly.integrate()

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

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

The Poly class and complex roots

We have already seen that there is the Poly class in SymPy. 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
type(A)
sympy.polys.polytools.Poly

Note that our polynomial A is of the third degree. Let’s see all its roots:

A.all_roots()

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

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

It makes sense, the first root is real, the others, complex.

This explains the existence of an all_roots method. After all, if there is a method that is concerned with returning all roots, it is because there must be some method that returns only some, right? Depending on the application in question, we may be interested only in real roots. Then we can use the real_roots method:

A.real_roots()

Polynomial division

Now, let’s check two Poly objects that we created throughout the article:

A
Q_poly

To facilitate understanding, let’s associate a variable B to Q_poly:

B = Q_poly
B

In the article call, we talked about polynomial division. The time has come, let’s see how to do such an operation with SymPy.

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

sympy.div(A, B)

The return 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
remainder

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

B * quotient + remainder
A == B * quotient + remainder
True

Simple as that.

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

If you want to know when new articles are available, follow Chemistry Programming on X. The complete list of articles about SymPy can be seen in the SymPy tag. And, if you would like to see all the tutorial series posts, see the SymPy Tutorials page.

See you next time.

Leave a Comment

Your email address will not be published. Required fields are marked *

Scroll to Top