Click here to Skip to main content
65,938 articles
CodeProject is changing. Read more.
Articles / programming / algorithm

Euclidean Division of Polynomials

1.00/5 (1 vote)
18 Apr 2023CPOL7 min read 10.7K  
This article describes how to divide two polynomials and shows the source code to calculate this division. There is also code to add, subtract and multiply two polynomials.
This article describes how to implement an algorithm to divide two polynomials. The code is written in a simple programming language that can be easily translated into other programming languages.

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:

Python
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.

Python
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:

Python
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 ith 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
ASM
[<<[
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:

Python
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

Python
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 nth 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:

Python
quotient = poly@add(quotient, q)

The Computation of the Remainder

Python
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

  • 18th April, 2023: First version

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)