Polynomials with numpy

In this example, we explore how to perform operations with polynomials. We will do this in two ways:

✍🏼 First, manually implementing the operations ourselves.
🔧 Then, using the functions provided by numpy.

Often, Python provides everything we need without requiring us to implement it ourselves. However, it is beneficial to understand how things work by trying to implement them ourselves, even if we later use an external package.

Let’s start by understanding what a polynomial is. A polynomial has the following form: 5x^3 + 3x^2 + 10x + 4. It has the following components:

❓ The independent variable. In our case, we use x as it is common.
🧮 Some coefficients. In our case, 5, 3, 10, 4. Each coefficient accompanies a power of x. A coefficient can be zero.
🌡️ The degree. It refers to the maximum exponent. In our case, the degree is 3.

We can express our polynomial in Python as follows. A simple list with the coefficients, ordered from the lowest to the highest degree.

p = [4, 10, 3, 5]

We assume that the first index p[0] corresponds to the minimum degree, that is, to the x^0 term, which is the 4 in 5x^3 + 3x^2 + 10x + 4.

Now that we know what a polynomial is and how to represent it in Python, we can define the following operations:

➕ Sum. Add two polynomials p and q.
➖ Subtract. Subtract two polynomials p and q.
✖️ Multiply. Multiply two polynomials p and q.
🧮 Evaluate at one point. Replace the x of p with a particular value.

Let’s start with addition and subtraction.

from itertools import zip_longest

def add(p, q):
    return [pp + qq for pp, qq in zip_longest(p, q, fillvalue=0)]

def subtract(p, q):
    return [pp - qq for pp, qq in zip_longest(p, q, fillvalue=0)]

It is important to note that we use zip_longest instead of zip since the latter assumes that both lists have the same length. In the case of polynomials, this is not the case. You can add polynomials with different degrees.

Multiplication is a bit more involved. It consists of taking each term of the first polynomial and multiplying it by all the terms of the second polynomial, adding the result.

def multiply(p, q):
    degree = len(p) + len(q) - 2
    result = [0] * (degree + 1)

    for i, pp in enumerate(p):
        for j, qq in enumerate(q):
            result[i + j] += pp * qq

    return result

Evaluating a polynomial at a point involves replacing x with a given value and performing the operations. There are different ways to do this. This is Horner’s method, which is quite efficient.

def evaluate(p, x):
    result = 0
    for coef in reversed(p):
        result = result * x + coef
    return result

Now let’s use our functions. We define two polynomials:

p: With the value x^3 + 3.
q: With the value 5x^2 + x + 2.

p = [3, 0, 1]
q = [2, 1, 5]

➕ We add them together. The result is interpreted as 6 + x + 5x^2.

print(add(p, q))
# [6, 1, 5]

➖ We subtract them.

print(subtract(p, q))
# [-4, -1, 1]

✖️ We multiply them.

print(multiply(p, q))
# [5, 1, 17, 3, 6]

🧮 And we evaluate p at the point 5.

print(evaluate(p, 5))
# 76

Although it is important to understand how things work, in practice, it may be better to use packages that offer these functions already implemented and tested. In many cases, they are also optimized.

For this, we can use numpy and its functions. There is no need to reinvent the wheel. As you can see, it offers functions to do everything. The result is the same.

import numpy as np

p = [1, 0, 3]
q = [5, 1, 2]

print(np.polynomial.polynomial.polyadd(p, q))
print(np.polynomial.polynomial.polysub(p, q))
print(np.polynomial.polynomial.polymul(p, q))
print(np.polynomial.polynomial.polyval(5, p))

✏️ Exercises:

Write a function that, given two polynomials p and q, performs their division, indicating their quotient and remainder. Verify with numpy.polydiv that you got it right.
Write a function that, given a polynomial p expressed as a list [3, 2, 1] returns it in the form 1x^2 + 2x + 3 as a string.