Optimization: Polynomial division

[pventuzelo]

Easy improvement of division speed between polynomials

Executive Summary

Dividing two polynomials with significantly different degrees can take a long time (about 1 minute (to several hours) for degree ≃20.000 and 10.000).

Steps to Reproduce

Running this program may take around 60 seconds, depending on your system:

use std::time::Instant;

use lambdaworks_math::elliptic_curve::short_weierstrass::curves::bls12_381::curve::BLS12381FieldElement as F;
use lambdaworks_math::polynomial::Polynomial;

fn main() {
    let degree_poly_1 = 20_000; // Arbitrary value
    let mut coeffs1 = vec![F::zero()];
    for _ in 0..degree_poly_1 - 1 {
        coeffs1.push(F::zero());
        // coeffs1.push(F::one() + F::one()); // Or in place, this one lasts several hours
    }
    coeffs1.push(F::one());
    let p1 = Polynomial::new(&coeffs1);
    // Arbitrary polynomial, it could be computed with `new_monomial`, or in another way

    let degree_poly_2 = 10_000; // Arbitrary value
    let mut coeffs2 = vec![F::zero()];
    for _ in 0..degree_poly_2 - 1 {
        coeffs2.push(F::zero());
    }
    coeffs2.push(F::one());
    let p2 = Polynomial::new(&coeffs2);
    // Arbitrary polynomial, it could be computed with `new_monomial`, or in another way

    let start = Instant::now();
    let _q = p1 / p2; // We divide p2 by p2
    let duration = start.elapsed();

    // And display the time it took to compute it
    println!("Time to compute the quotient between a polynomial of degree {} and another of degree {}: {}s",
        degree_poly_1,
        degree_poly_2,
        duration.as_secs()
    );
}

Root Cause Analysis

Dividing a polynomial by another will call the mul_with_ref function, which computes the product of two polynomials in a natural way. But in long_division_with_remainder, this function is used with a monomial in argument, and so the mul_with_ref function wastes time computing useless products and sums (for each 0 coefficient of the monomial factor).

Recommendations

We can replace the current code by this one, which just check that both coefficients are not 0 to perform operations:

pub fn mul_with_ref(&self, factor: &Self) -> Self {
        let degree = self.degree() + factor.degree();
        let mut coefficients = vec![FieldElement::zero(); degree + 1];

        if self.coefficients.is_empty() || factor.coefficients.is_empty() {
            Polynomial::new(&[FieldElement::zero()])
        } else {
            for i in 0..=factor.degree() {
                if !factor.coefficients[i].is_zero() {
                    for j in 0..=self.degree() {
                        if !self.coefficients[j].is_zero() {
                            coefficients[i + j] += &factor.coefficients[i] * &self.coefficients[j];
                        }
                    }
                }
            }
            Polynomial::new(&coefficients)
        }
    }

I believe the overhead of checking for zero coefficients is negligible compared to the performance gains from avoiding unnecessary computations.

Or to keep a constant time property for classical multiplication, one can instead create the previous function only for use in the long_division_with_remainder function, removing the if !self.coefficients[j].is_zero() { condition:

pub fn mul_with_ref_for_long_division_with_remainer(&self, factor: &Self) -> Self {
        let degree = self.degree() + factor.degree();
        let mut coefficients = vec![FieldElement::zero(); degree + 1];

        if self.coefficients.is_empty() || factor.coefficients.is_empty() {
            Polynomial::new(&[FieldElement::zero()])
        } else {
            for i in 0..=factor.degree() {
                if !factor.coefficients[i].is_zero() {
                    for j in 0..=self.degree() {
                            coefficients[i + j] += &factor.coefficients[i] * &self.coefficients[j];
                    }
                }
            }
            Polynomial::new(&coefficients)
        }
    }

and keep the current mul_with_ref function as is.

Benchmarks

On my system, the first division operation (the one not commented) is reduced from about 1 minutes without the verifications to less than 1 second with verifications (the duration.as_secs() function returned 0 seconds in my tests, indicating that the operation now completes in less than 1 second), and the second one (the commented one) is reduced from several hour to about 1 minute.