Polynomials ring

by marco.p.v.vezzoli · February 5, 2022

Photo by Farhad Ibrahimzade on Unsplash Polynomials post list:

post 1: Polynomial Basics
post 2: Polynomials Ring

In the previous post I presented a representation of polynomials of a single variable (monovariate) with coefficients in a ring and their evaluation. Also polynomials, equipped with sum and multiplication can form a ring: let's recall briefly the ring properties

it is an abelian group respect to the sum, i.e.:
1. sum is associative
2. sum is commutative
3. there is a neutral element for the sum
4. all elements have an additive inverse
it is a monoid respect to the multiplication, i.e.:
1. multiplication is associative
2. there is a neutral element for the multiplication
multiplication is distributive respect to sum

An important consequence is that multiplicative inverse may not exist. This imply that we can easily implement three traits in rust: Add, Sub and Mult. If we restrict the coefficient to be a field (i.e. each element has a multiplicative inverse) we could perform division returning a quotient and a remainder. Today (Feb 2021) we are using some unstable features to implement the Fn* traits, the code in this post require a nightly channel and the following flags:

Name: poly-preamble.

#![feature(unboxed_closures)]
#![feature(fn_traits)]

recalling the implementation so far

The following implementation just adds the overload of == operator (PartialEq) and the definition of a degree method which returns the length of the coefficient vector supposing its last member (the highest power coefficient) is non zero.

Name: poly-basics.

use std::cmp::{min,max};
use std::ops::{Fn,FnOnce,FnMut};

// I need to guarantee that a type has a neutral element for sum
pub trait Zero{
    fn zero() -> Self;
}

#[derive(Debug, PartialEq, Clone)]
pub struct Poly <T> {
    coeff : Vec<T>
}

impl<T> Poly<T>{
    pub fn new(coeff : Vec<T>) -> Poly<T>{
        Poly{ coeff:coeff }
    }
    pub fn degree(self : &Poly<T>) -> usize {
        max(0, self.coeff.len() - 1)
    }
}

impl<T: fmt::Debug> fmt::Display for Poly<T> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let result : Vec<String> = self.coeff.iter()
            .enumerate()
            .map(|(i, c)| format!("({:?})^{}",c,i))
            .collect();
        write!(f, "{}", result.join(" + "))
    }
}

// here starts the function interface implementation
impl<T: Mul<Output = T> + Add<Output = T> + Copy + Zero> Fn<(T, )> for Poly<T> {
    extern "rust-call" fn call(&self, args: (T,)) -> T {
        self.coeff.iter()
            .rev()
            .fold(T::zero(), |acc, c| acc * args.0 + *c)
    }
}

impl<T: Mul<Output = T> + Add<Output = T> + Copy + Zero> FnMut<(T, )> for Poly<T> {
    extern "rust-call" fn call_mut(&mut self, args: (T,)) -> T {
        self.call(args)
    }
}

impl<T: Mul<Output = T> + Add<Output = T> + Copy + Zero> FnOnce<(T, )> for Poly<T> {
    type Output = T;

    extern "rust-call" fn call_once(self, args: (T,)) -> T {
        self.call(args)
    }
}
// here are a couple of implementation of general use
impl Zero for u32{
    fn zero() -> u32 {
        0
    }
}

impl Zero for i32{
    fn zero() -> i32 {
        0
    }
}

impl Zero for f32{
    fn zero() -> f32 {
        0.0
    }
}

A couple of simple unit test was also defined to check the implementation

Name: poly-basics-tests.

#[test]
fn test_evaluation() {
    let p : Poly<i32> = Poly::new(vec![1,2,3]);
    assert_eq!(p(10), 321);
}

#[test]
fn test_formatting() {
    let x : Poly<i32>= Poly::new(vec![1,2,3]);
    assert_eq!(format!("{}", x), "(1)^0 + (2)^1 + (3)^2");
    assert_eq!(x.degree(), 2);
}

Now it's time to define a few operations on them; I'll create an implementation for the reference type (which does not take ownership of the structure) and another for the plain type

sum of polynomials

The sum overload is obtained by implementing the std::ops::Add trait. Self::Output is the expected outcoming type; its default is the same type where the implementation is required all the coefficients are summed from both polynomials up to the smallest degree, the coefficients for higher degrees are then chained

Name: poly-add.

impl<T: Add<Output = T> + Copy> Add for Poly<T>{
    type Output = Poly<T>;
    fn add(self : Poly<T>, other : Poly<T>) -> Self::Output {
        &self + &other
    }
}

impl<T: Add<Output = T> + Copy> Add for &Poly<T>{
    type Output = Poly<T>;
    fn add(self, other : &Poly<T>) -> Self::Output {
        let common : Vec<T> = self.coeff.iter()
            .zip(other.coeff.iter())
            .map(|(c1, c2)| *c1 + *c2)
            .collect();
        let (_, rest) = if self.degree() > other.degree() {
            self.coeff.split_at(other.degree() + 1)
        } else {
            other.coeff.split_at(self.degree() + 1)
        };
        Poly{
            coeff: common.iter()
                .clone()
                .chain(
                    rest.iter()
                        .clone()
                ).map(|x| *x)
                .collect()
        }
    }
}

here are the unit tests: you can see where ownership has to be accounted.

Name: poly-sum-test.

#[test]
fn test_sum(){
    let x = Poly{coeff:vec![1,3,2]};
    let y = Poly{coeff:vec![8,2]};
    let expected = Poly{coeff:vec![9,5,2]};
    assert_eq!(x.clone() + y.clone(), expected.clone());
    assert_eq!(y.clone() + x.clone(), expected);
    let expected2 = Poly{coeff:vec![2,6,4]};
    assert_eq!(x.clone() + x, expected2);
}

#[test]
fn test_sum_deref(){
    let x = Poly{coeff:vec![1,3,2]};
    let y = Poly{coeff:vec![8,2]};
    let expected = Poly{coeff:vec![9,5,2]};
    assert_eq!(&x + &y, expected);
    assert_eq!(&y + &x, expected);
    let expected2 = Poly{coeff:vec![2,6,4]};
    assert_eq!(&x + &x, expected2);
}

difference of polynomials

The subtraction overload is obtained by implementing the std::ops::Sub trait. The code is almost identical: only the upper part of the coefficient vector must be transformed when the second operand has higher degree than the first. This code is not simplifying 0 coefficients at the highest degree: this require some zero check that depends on the context: e.g. floating points rounding

Name: poly-sub.

impl<T: Sub<Output = T> + Copy + Zero + fmt::Debug> Sub for Poly<T>{
    type Output = Poly<T>;
    fn sub(self : Poly<T>, other : Poly<T>) -> Self::Output {
        &self - &other
    }
}

impl<T: Sub<Output = T> + Copy + Zero + fmt::Debug> Sub for &Poly<T>{
    type Output = Poly<T>;
    fn sub(self, other : &Poly<T>) -> Self::Output {
        let common : Vec<T> = self.coeff.iter()
            .zip(other.coeff.iter())
            .map(|(c1, c2)| *c1 - *c2)
            .collect();
        let ((_, rest), invert) = if self.degree() > other.degree() {
            (self.coeff.split_at(other.degree() + 1), false)
        } else {
            (other.coeff.split_at(self.degree() + 1), true)
        };
        Poly{
            coeff: common.iter()
                .clone()
                .map(|x| *x )
                .chain(
                    rest.iter()
                        .map(|&x| if invert {T::zero() - x} else { x })
                        .clone()
                )
                .collect()
        }
    }
}

here are the unit tests: you can see where ownership has to be accounted.

Name: poly-sub-test.

#[test]
fn test_sub(){
    let x = Poly{coeff:vec![1,3,2]};
    let y = Poly{coeff:vec![8,2]};
    let expected1 = Poly{coeff:vec![-7,1,2]};
    let expected2 = Poly{coeff:vec![7,-1,-2]};
    assert_eq!(x.clone() - y.clone(), expected1);
    assert_eq!(y.clone() - x.clone(), expected2);
    let expected3 = Poly{coeff:vec![0,0,0]};
    assert_eq!(x.clone() - x, expected3);
}

#[test]
fn test_sub_deref(){
    let x = Poly{coeff:vec![1,3,2]};
    let y = Poly{coeff:vec![8,2]};
    let expected1 = Poly{coeff:vec![-7,1,2]};
    let expected2 = Poly{coeff:vec![7,-1,-2]};
    assert_eq!(&x - &y, expected1);
    assert_eq!(&y - &x, expected2);
    let expected3 = Poly{coeff:vec![0,0,0]};
    assert_eq!(&x - &x, expected3);
}

product of two polynomials

This is a little bit trickier than the sum. The i-th coefficient should be calculatead as

$A[x] = \sum_{j=0}^{n}a_{j}x^{j}$

$B[x] = \sum_{k=0}^{m}b_{k}x^{k}$

$C[x] = A[x]B[x]$

$c_i = \sum_{j+k=i}a_{j}b_{k}$

where:

n is the degree of polynomial A
m is the degree of polynomial B
i runs from 0 to n+m
j runs from 0 to n
k runs from 0 to m

the outer product of the coefficient vector looks like this table

	0	1	2	…	n
0	0	1	2	…	n
1	1	2	3	…	n+1
2	2	3	4	…	n+2
…	…	…	…	…	…
m	m	m+1	m+2	…	n+m

we could rewrite the equation as

$c_i = \sum_{j=max(0,i-m)}^{min(i,n)}a_{j}b_{n+m-j}$

Name: poly-mul.

impl<T: Mul<Output = T> + Add<Output = T> + fmt::Debug + Copy + PartialEq + Zero> Mul for Poly<T>{
    type Output = Poly<T>;
    fn mul(self : Poly<T>, other : Poly<T>) -> Self::Output {
        &self * &other
    }
}

impl<T: Mul<Output = T> + Add<Output = T> + fmt::Debug + Copy + PartialEq + Zero> Mul for &Poly<T>{
    type Output = Poly<T>;
    fn mul(self, other : &Poly<T>) -> Self::Output {
        let mydeg = self.degree() as i64;
        let otherdeg = other.degree() as i64;
        let resultdeg = mydeg + otherdeg;
        let resultcoeff : Vec<T>= (0..(resultdeg +1))
            .map(|i| ((max(0,i-otherdeg))..(min(i,mydeg) + 1))
                 .map(|j| self.coeff[j as usize] * other.coeff[(i-j) as usize])
                 .fold(T::zero(), |a, b| a + b))
            .collect();
        Poly{coeff:resultcoeff}
    }
}

also here let's add a few tests

Name: poly-mul-test.

#[test]
fn test_mul_deref(){
    let x = Poly{coeff:vec![1,3,2]};
    let y = Poly{coeff:vec![8,2]};
    println!("{:?}", &x * &y);
    let expected = Poly{coeff:vec![8,26,22,4]};
    assert_eq!(&x * &y, expected);
    assert_eq!(&y * &x, expected);
    let expected2 = Poly{coeff:vec![1,6,13,12,4]};
    assert_eq!(&x * &x, expected2);
}

scalar operations

My main objective is to write a division algorithm for polynomials. To do this I need a few utility functions. the following trait implementation is a convenient shortcut to multiply a polynomial by an element of the ring of its coefficients. I'm going to implement only the multiplication on the right side

Name: poly-scalar-mul.

impl<T: Mul<Output = T> + Copy> Mul<T> for Poly<T> where {
    type Output = Poly<T>;
    fn mul(self, rhs : T) -> Poly<T> {
        let coeff : Vec<T> = self.coeff.iter()
        .map(
            |c| rhs * (*c)
        ).collect();
        Poly{
            coeff : coeff
        }
    }
}

Polynomials ring

recalling the implementation so far

sum of polynomials

difference of polynomials

product of two polynomials

scalar operations

Like this:

Related

You may also like...

Polynomials ring

recalling the implementation so far

sum of polynomials

difference of polynomials

product of two polynomials

scalar operations

Share this:

Like this:

Related

You may also like...

Polynomial basics