Introduction

Welcome to the division of polynomials in Python using the Polynomial library. In this tutorial, we will be discussing the mathematical concept of polynomial division and how to implement it using the Polynomial library. We will also be providing a step-by-step guide on how to divide two polynomials and how to use the code in your own projects.

Before we begin, it's important to understand what a polynomial is and how it is represented in Python. A polynomial is a mathematical expression consisting of variables (also known as coefficients) and exponents. In Python, with the Polynomial library, a polynomial is represented as a list of coefficients, like [1, 2, 3, ...].

The Euclidean division of polynomials is not a simple algorithm even if it is completely inspired by the division of integers. In contrast, summing two polynomials or the product is a more intuitive and easy-to-understand algorithm. In this article, I present a code to create a polynomial as a list of numeric coefficients, to know the size of the list, to calculate the largest non-zero degree of the polynomial as well as addition and multiplication.

Background

To make this code, I was inspired by the way we declare a polynomial in Python and also I found many explanations on the Internet on how to divide two polynomials including a YouTube video explaining the calculations.

Even if it is well explained, it is nevertheless a difficult exercise to implement because there are differences between the way of dividing two polynomials on a paper compared to an algorithm to be carried out.

Using the Code

To test this code, you can go to talenha.boribar.com. An input window allows you to enter the code and test it. Even if difficulties in testing this code are possible, you can easily transcribe it into your favorite language.

Introduction

This allows us to perform operations on polynomials more easily, especially when it comes to calculating the derivative and integral of a polynomial.

In this article, we will see how to implement a division algorithm for polynomials using the Euclidean division. This method consists of finding the quotient and remainder of a polynomial division. It is a generalization of the division algorithm for integers and can be used to find the solutions of polynomial equations.

We will first introduce the mathematical concepts necessary to understand the algorithm and then implement it in Python using the Polynomial library. Finally, we will compare the results with those obtained using the built-in function for polynomial division in Python and discuss the advantages and limitations of our implementation.

Mathematical Concepts

The division of two polynomials is a way to express the result of the division of one polynomial by another as a quotient polynomial and a remainder polynomial. It is based on the long division algorithm that we learn in elementary school, but extended to polynomials.

The division of a polynomial f(x) by a polynomial g(x) can be written as follows:

f(x) = g(x) * q(x) + r(x)

where q(x) is the quotient polynomial and r(x) is the remainder polynomial.

If you want to get the divided formula, you just need to divide the equality side by side by g(x).

f(x) / g(x) =        q(x)          +         r(x) / g(x)

To perform the division, we start by dividing the leading coefficient of f(x) by the leading coefficient of g(x). This gives us the first coefficient of the quotient polynomial q(x). We then multiply g(x) by this coefficient and subtract it from f(x) to obtain a new polynomial, which we call the "reduced" polynomial. We repeat this process until the degree of the reduced polynomial is less than the degree of g(x). The final reduced polynomial is the remainder polynomial r(x).

For example, let's consider the division of the polynomial f(x) = x^3 + 2x^2 + 3x + 4 by the polynomial g(x) = x^2 + x + 1. We start by dividing the leading coefficient of f(x), which is 1, by the leading coefficient of g(x), which is also 1. This gives us the first coefficient of the quotient polynomial q(x), which is 1. We then multiply g(x) by this coefficient and subtract it from f(x) to obtain the reduced polynomial. It is important to note that the next iteration starts with the new function f(x) - C * g(x) and not f(x).

The Iterative Algorithm to Divide Two Polynomials

f(x)	g(x)	C	f(x) - C * g(x)
x^3	x^2	x	x^2 + 2x + 4
x^2+2x+4	x^2	1	x + 3

The result is finally x + 1 - (x+3) / (x^2+x+1).

When we expand this equation, we get an equation that is exactly equal to f(x).

Python Program

In Python, with the Polynomial library, a polynomial is represented as a list of coefficients, like [1, 2, 3, ...]. The order in which the coefficients are arranged is important, because they are arranged from left to right starting with the coefficient from degree n to degree 0 with respect to the polynomial object. If you defined a list of coefficients with:

np.poly1d([1, 2, 3])

they are arranged from left to right starting with the coefficient from degree 0 to degree n.

Now that we have a basic understanding of polynomials and how they are represented in Python, let's dive into the concept of polynomial division and how to implement it in Python.

Test in Python

To test in Python, we are going to use the web application replit.com.

from polynomial import poly
import numpy as np

# Print a polynomial object with the coefficients [1, 2, 3, 4]
print(poly([1, 2, 3, 4]))

# Print a polynomial object with the coefficients [1, 1, 1]
print(poly([1,1,1]))

# Divide the polynomial [ 1, 2, 3, 4 ] by the polynomial [ 1, 1, 1 ]
print(np.polydiv(poly([1,2,3,4]), poly([ 1, 1, 1 ])))

The output is exactly equal to the previous result.

[1, 2, 3, 4]
[1, 1, 1]
(array([1., 1.]), array([1., 3.]))

Implemented Functions in a Simple Programming Language

In this version of a polynomial, it's important to note that coefficients are arranged from left to right starting with the coefficient from degree n to degree 0.

For example, [ 1, 2, 3, 4 ] is printed as an equation by:

x^3 + 2x^2 + 3x + 4

There are functions that can be used:

poly@create(n): Creates a new polynomial of order n
poly@copy(p): Clones an existing polynomial
poly@size(p): Computes the size of the list of polynomial coefficients
poly@degree(p): Returns the highest degree whose coefficient is nonzero
poly@trim(p): Removes all leading 0s from the list
poly@isNull(p): Determines if it is a null polynomial
poly@nth(p, i): Returns the coefficient of the i^th degree of the polynomial
poly@add(p1, p2): Returns a polynomial which is the sum of p1 and p2
poly@sub(p1, p2): Returns a polynomial which is the subtraction of p1 and p2
poly@mul(p1,p2): Returns a polynomial which is the product of p1 and p2
poly@monome(i, c): Constructs a monome of degree i with the coefficient c
poly@div(p1, p2): Divides the polynomial p1 by p2 and return the quotient and the remainder of the division
poly@print(p): Prints the polynomial expression in algebraic form

[<<[
fun poly@create(n) {
     p = [0],
     while(n > 0) {
          p = p 0,
          n = n - 1
     },
     output(p)
},

fun poly@copy(p) {
     copy = [],
     count = 0,
     while(count < p.length) {
          copy = copy p(count),
          count = count + 1
     },
     output(copy)
},

fun poly@size(p) {
     output(p.length - 1)
},

fun poly@degree(p) {
     select {
         default {
            foreach(i from p) {
                if (p(i) == 0) {} else jump "break"
            },
            output(-1)
         },
         break { output(p.length - 1 - i) }
     }
},

fun poly@nth(p, i) {
     output(p(poly@size(p) - i))
},

fun poly@isNull(p) {
     output(poly@degree(p) == -1)
},

fun poly@add(p1, p2) {
     count1 = poly@size(p1),
     count2 = poly@size(p2),
     // init result
     result = poly@create(number@max(count1, count2)),
     while(count1 >= 0 && count2 >= 0) {
           count = number@max(count1,count2),
           result(count) = p1(count1) + p2(count2),
           count1 = count1 - 1,
           count2 = count2 - 1
     },
     while(count1 >= 0) {
           result(count1) = p1(count1),
           count1 = count1 - 1
     },
     while(count2 >= 0) {
           result(count2) = p2(count2),
           count2 = count2 - 1
     },
     output(result)
},

fun poly@sub(p1, p2) {
     count1 = poly@size(p1),
     count2 = poly@size(p2),
     // init result
     result = poly@create(number@max(count1, count2)),
     while(count1 >= 0 && count2 >= 0) {
           count = number@max(count1,count2),
           result(count) = p1(count1) - p2(count2),
           count1 = count1 - 1,
           count2 = count2 - 1
     },
     while(count1 >= 0) {
           result(count1) = p1(count1),
           count1 = count1 - 1
     },
     while(count2 >= 0) {
           result(count2) = -p2(count2),
           count2 = count2 - 1
     },
     output(result)
},

fun poly@mul(p1, p2) {
     result = poly@create(poly@size(p1) + poly@size(p2)),
     count2 = poly@size(p2),
     while(count2 >= 0) {
           count1 = poly@size(p1),
           while(count1 >= 0) {
                 pos = poly@size(result) - (poly@size(p1) - count1 + 
                       poly@size(p2) - count2),
                 result(pos) = result(pos) + p1(count1) * p2(count2),
                 count1 = count1 - 1
           },
           count2 = count2 - 1
     },
     output(result)
},

fun poly@monome(degree, coefficient) {
     count = degree,
     o = [0],
     while(count > 0) {
         o = o 0,
         count = count - 1
     },
     o(0) = coefficient,
     output(o)
},

fun poly@div(p1, p2) {
     if (poly@degree(p1) < poly@degree(p2)) {
          output({ quotient : [0], remainder : p1})
     } else {
          quotient = [0],
          remainder = poly@copy(p1),
          select {
              default {
                   pos = poly@size(p2),
                   loop = 0,
                   while(poly@degree(remainder) >= poly@degree(p2)) {
                          q = poly@monome(poly@degree(remainder) - poly@degree(p2),
                                          remainder(poly@size(remainder) - 
                                             poly@degree(remainder)) / p2(0)),
                          print("polynome quotient " poly@print(q)),
                          quotient = poly@add(quotient, q),
                          remainder = poly@sub(remainder, poly@mul(p2, q)),
                          print("polynome reste = " poly@print(remainder)),
                          if (poly@degree(remainder) == -1)
                                 jump "break",
                          pos = pos - 1,
                          if (loop > 5) jump "break",
                          loop = loop + 1
                   }
              },
              break {}
          },
          output({ quotient : quotient, remainder : remainder })
     }
},
fun poly@print(p) {
         if (p.length > 0) {
             count = poly@size(p),
             if (p(count) == 0) { o = "" } else o = p(count),
             count = count - 1,
             while(count >= 0) {
                 if (poly@size(p) - count == 1) {
                        degree = "x"
                 } else degree = "x^" (poly@size(p) - count),
                 if (o) { plus = " + " } else plus = "",
                 if (p(count) == 1)
                      { o = degree plus o }
                 else if (p(count) == 0) {}
                 else
                      o = p(count) degree plus o,
                 count = count - 1
             },
             if (o) {} else o = 0,
             output(o)
         }
}
]>>]

Explanations of the Function poly@div(p1, p2)

Start Condition

First, the degree of p1 must be greater than the degree of p2. If the degree of p1 is lesser than the degree of p2, then the quotient is 0, the remainder is p1 and the function ends.

While Loop Condition

The while condition loop poly@degree(remainder) >= poly@degree(p2) is false when the degree of the remainder is less than the degree of p2. For example:

remainder = 2x + 1          // <-- degree 1
p2 = x^2 + 2x + 1           // <--  degree 2
poly@degree(remainder) >= poly@degree(p2)      // <-- false, the while loop ends

The Calculation of the Quotient

q = poly@monome(poly@degree(remainder) - poly@degree(p2),
                remainder(poly@size(remainder) - poly@degree(remainder)) / p2(0))

The function poly@monome(n, c) takes a number n and a coefficient c then creates a new polynomial of degree n and coefficient c. The degree is calculated by the current degree of the remainder subtracted from the degree of the polynomial p2. This means that you calculate the degree of the quotient of the polynomial as if you were dividing x^3 by x.

Then you calculate the coefficient by dividing the n^th remainder coefficient by the first coefficient of p2. The coefficient is positioned at the first non-zero coefficient in the list (from left to right). Remember that p2 must be trimmed on the left: the first coefficient must be at position 0 in the list.

The Addition to the Final Quotient

The computed quotient is a monomial. It is just added to the final quotient by this operation:

quotient = poly@add(quotient, q)

The Computation of the Remainder

remainder = poly@sub(remainder, poly@mul(p2, q))

This is exactly the equation to compute the new remainder f(x) - C * g(x).

Note that to compute the Euclidean division, you have to implement the subtract and the product of two polynomials.

Points of Interest

I used ChatGPT to help me debug the source code shown here. Although he explained Euclidean division to me several times, it is not an easy algorithm. But, it is an important algorithm for completing functions of polynomials.

History

18^th April, 2023: First version

Euclidean Division of Polynomials