Implementing polynomial ring inverse

Purpose:

The purpose of this deliverable was to implement a component algorithm from [Hoffstein1998]. By working on this deliverable I gained an understanding of polynomial rings, modular polynomial arithmetic, the Extended Euclidean Algorithm, finding polynomial inverses, and the NTRU public key system.

Polynomial arithmetic implementation:

Cargo.toml

[package]
name = "polynomial"
version = "0.1.0"
authors = ["Michaela Molina <michaela.molina@sjsu.edu>"]
edition = "2018"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
gcd="2.0.1"
modinverse = "0.1.1"

main.rs

use std::process;
mod arithmetic;
use gcd::Gcd;
use modinverse::modinverse;

fn main() {
    
    // Test 1: N = 11, p = 7
    let mut p = 7;
    let mut ax = vec![1,0,0,0,0,0,0,0,0,0,0,-1];
    let mut test1 = vec![-1,1,0,0,1,0,-1,0,1,1,-1];
    let mut inverse1 = find_inverse(p, &ax, &test1);
    println!("inverse1 = {}", arithmetic::get_polynomial_string(&inverse1));

    // Test 2: N = 11, p = 97
    let mut p = 11;
    let mut ax = vec![1,0,0,0,0,0,0,0,0,0,0,-1];
    let mut test2 = vec![-1,1,0,0,1,0,-1,0,1,1,-1];
    let mut inverse2 = find_inverse(p, &ax, &test2);
    println!("inverse2 = {}", arithmetic::get_polynomial_string(&inverse2));
 
    // Test 3: N = 7, p = 37
    let mut p = 37;
    let mut ax = vec![1,0,0,0,0,0,0,-1];
    let mut test3 = vec![-1,1,0,0,1,0,-1,0,1,1,-1];
    let mut inverse3 = find_inverse(p, &ax, &test3);
    println!("inverse3 = {}", arithmetic::get_polynomial_string(&inverse3));

    // Test 4: N = 12, p = 73
    let mut p = 73;
    let mut ax = vec![1,0,0,0,0,0,0,0,0,0,0,-1];
    let mut test4 = vec![-1,1,-1,0,0,1,0,-1,1,0,-1];
    let mut inverse4 = find_inverse(p, &ax, &test4);
    println!("inverse4 = {}", arithmetic::get_polynomial_string(&inverse4));

    // Test 5: N = 15, p = 11
    let mut p = 11;
    let mut ax = vec![1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1];
    let mut test5 = vec![1,1,1,1,-1,-1,-1,-1,1,0,0,0,0];
    let mut inverse5 = find_inverse(p, &ax, &test5);
    println!("inverse5 = {}", arithmetic::get_polynomial_string(&inverse5));

    // Encryption case: N = 11, p = 7, q = 97
    let mut fx = vec![-1,1,0,0,1,0,-1,0,1,1,-1];
    let mut gx = vec![-1,0,-1,0,0,1,0,1,1,0,-1];
    let mut ax = vec![1,0,0,0,0,0,0,0,0,0,0,-1];
    let mut test1 = vec![-1,1,0,0,1,0,-1,0,1,1,-1];
    let mut bx = fx.clone(); 
    let mut p = 7;
    let mut q = 97;
    let mut f_p = find_inverse(p, &ax, &bx);
    let mut f_q = find_inverse(q, &ax, &bx);
    println!("inverse1 = {}", arithmetic::get_polynomial_string(&f_p));
    println!("inverse2 = {}", arithmetic::get_polynomial_string(&f_q));
    let mut h = arithmetic::multiply(&f_q, &gx, q as i64);
    let mut h = arithmetic::get_remainder(&mut h.clone(), &mut ax.clone(), q as i64);
    println!("h = {}", arithmetic::get_polynomial_string(&h));
    
}

// returns one factor 
fn find_inverse(p:u64, ax:&Vec<i64>, bx:&Vec<i64>) -> Vec<i64> {
    
    let (mut q, mut r1) = gcd(ax.clone(),bx.clone(),p as i64);
    // q is quotient list
    // r1 is remainder list
    let mut v   = compute_v(&q, &r1, p as i64);
    let mut w   = compute_w(&q, &r1, p as i64);
    // r2 is checked remainder list
    let mut r2 = check_remainder(ax, bx, &q, &r1, &v, &w, p as i64); 

   for i in 0..v.len() {
        println!("v_{} = {}", i, arithmetic::get_polynomial_string(&v[i])); 
        println!("w_{} = {}", i, arithmetic::get_polynomial_string(&w[i])); 
        println!("q_{} = {}", i, arithmetic::get_polynomial_string(&q[i])); 
        println!("r1_{} = {}", i, arithmetic::get_polynomial_string(&r1[i])); 
        println!("r2_{} = {}", i, arithmetic::get_polynomial_string(&r2[i])); 
        println!();
    }
    let mut pinverse = w[w.len()-1].clone();
    let mut first_constant_remainder = arithmetic::get_polynomial_string(&r1[r1.len()-1]);
    let mut constant_remainder = r1[r1.len()-1][0] as u64;
    if (constant_remainder != 1) {
        let mut gcd = p.gcd_euclid(constant_remainder);
        if gcd == 1 {
            let mut inverse = mod_inverse(constant_remainder, p);
            pinverse = arithmetic::polynomial_multiply(&pinverse.clone(), &vec![inverse as i64]);
            pinverse = arithmetic::get_remainder(&mut pinverse, &mut ax.clone(), p as i64);   
        }
    }
    pinverse
}

fn mod_inverse(p:u64,q:u64) -> u64 {
    let mut i = 0;
    while true {
        let mut product = p*i;
        if (product % q == 1) {
            break;
        }
        i += 1;
    }
    i
}

// returns (list of quotients, list of remainders)
fn gcd(ax:Vec<i64>, bx:Vec<i64>, modulus:i64) -> (Vec<Vec<i64>>, Vec<Vec<i64>>) {
    let mut quotient_array:Vec<Vec<i64>>  = Vec::new();
    let mut remainder_array:Vec<Vec<i64>> = Vec::new();
    quotient_array.push(Vec::new());  // array of polynomials, one quotient from each step
    remainder_array.push(Vec::new()); // array of polynomials, one remainder from each step
    let mut dividend = ax.clone();    // larger polynomial
    let mut divisor = bx.clone();     // smaller polynomial
    while true {
        let (mut quotient, mut remainder) = arithmetic::divide_polynomials(&mut dividend, 
                                                             &mut divisor, 
                                                             modulus);
        quotient_array.push(quotient.clone());
        remainder_array.push(remainder.clone());
        if (remainder.len() < 2) { // stop when remainder is a constant
            break;  
        } else {
            dividend = divisor.clone();
            divisor = remainder.clone();
        }
    }
    (quotient_array, remainder_array)
}

fn compute_v(quotient_array:&Vec<Vec<i64>>, 
             remainder_array:&Vec<Vec<i64>>, 
             modulus:i64) -> Vec<Vec<i64>> {
    let mut v:Vec<Vec<i64>> = Vec::new();
    v.push(vec![0]); // v_0 = 0
    v.push(vec![1]); // v_1 = 1
    for i in 2..quotient_array.len() {
        let mut product = arithmetic::polynomial_multiply(&quotient_array[i], &v[i-1]);
        arithmetic::make_polynomial_negative(&mut product);
        let mut value = arithmetic::polynomial_add(&v[i-2], &product);
        arithmetic::mod_polynomial_coefficients(&mut value, modulus);  
        v.push(value);
    }
    v
}

fn compute_w(quotient_array:&Vec<Vec<i64>>, 
             remainder_array:&Vec<Vec<i64>>,    
             modulus:i64) -> Vec<Vec<i64>> {
    let mut w:Vec<Vec<i64>> = Vec::new();
    w.push(vec![1]); // w_0 = 1
    let mut w_1 = quotient_array[1].clone();
    arithmetic::make_polynomial_negative(&mut w_1);
    w.push(w_1); // w_1 = - q_1
    for i in 2..quotient_array.len() {
        let mut product = arithmetic::polynomial_multiply(&quotient_array[i], &w[i-1]);
        arithmetic::make_polynomial_negative(&mut product);
        let mut value = arithmetic::polynomial_add(&w[i-2], &product);
        arithmetic::mod_polynomial_coefficients(&mut value, modulus);  
        w.push(value);
    }
    w
}

// remainder_i = ax*vx + bx*wx
fn check_remainder(ax:&Vec<i64>, 
                   bx:&Vec<i64>, 
                   quotient_array:&Vec<Vec<i64>>, 
                   remainder_array:&Vec<Vec<i64>>, 
                   v_array:&Vec<Vec<i64>>, 
                   w_array:&Vec<Vec<i64>>, 
                   modulus:i64) -> Vec<Vec<i64>> {
    let mut r:Vec<Vec<i64>> = Vec::new();
    r.push(vec![0]);
    for i in 1..remainder_array.len() {
        let mut p1 = arithmetic::polynomial_multiply(ax,&v_array[i]);
        let mut p2 = arithmetic::polynomial_multiply(bx,&w_array[i]);
        let mut sum = arithmetic::polynomial_add(&p1,&p2);
        arithmetic::mod_polynomial_coefficients(&mut sum,modulus);
        arithmetic::remove_leading_zeros(&mut sum);
        r.push(sum);
    }
    r
}

fn coefficients_are_zero(v:&Vec<i64>) -> bool {
    let mut result = true;
    for coeff in v {
        if *coeff != (0 as i64) {
            result = false;
            break;
        }     
    }
    result
}

fn print_quotient_array(polynomials:&Vec<Vec<i64>>) {
    for (index, p) in polynomials.iter().enumerate() {
        println!("quotient_{} = {}", index, arithmetic::get_polynomial_string(p));
    }
}
fn print_v_array(polynomials:&Vec<Vec<i64>>) {
    for (index, p) in polynomials.iter().enumerate() {
        println!("v_{} = {}", index, arithmetic::get_polynomial_string(p));
    }
}
fn print_w_array(polynomials:&Vec<Vec<i64>>) {
    for (index, p) in polynomials.iter().enumerate() {
        println!("w_{} = {}", index, arithmetic::get_polynomial_string(p));
    }
}

fn print_remainder_array(polynomials:&Vec<Vec<i64>>) {
    for (index, p) in polynomials.iter().enumerate() {
        println!("remainder_{} = {}", index, arithmetic::get_polynomial_string(p));
    }
}

arithmetic.rs

use std::process;

pub fn get_remainder(dividend:&mut Vec<i64>, divisor:&mut Vec<i64>, modulus:i64) -> Vec<i64> {
    let (mut quotient, mut remainder) = divide_polynomials(dividend, divisor, modulus);
    mod_polynomial_coefficients(&mut remainder, modulus); 
    remainder
}

// returns (quotient, reminder)
pub fn divide_polynomials(dividend:&mut Vec<i64>, 
                      divisor:&mut Vec<i64>, 
                      modulus:i64) -> (Vec<i64>,Vec<i64>) {
    mod_polynomial_coefficients(dividend, modulus);
    mod_polynomial_coefficients(divisor,modulus);
    let mut quotient_vector:Vec<Vec<i64>> = Vec::new();
    let (quotients, remainder) = long_division(dividend, divisor, modulus, &mut quotient_vector); 
    let quotient = polynomial_add_multiple(&quotients);
    //mod_polynomial_coefficients(
    (quotient, remainder)
}

// returns (array of quotients, remainder)
pub fn long_division(dividend:&mut Vec<i64>, 
                     divisor:&mut Vec<i64>, 
                 modulus:i64, 
                quotients:&mut Vec<Vec<i64>>) -> (Vec<Vec<i64>>,Vec<i64>) {
    let (mut quotient, mut remainder) =  divide_once(&dividend, &divisor, modulus);
    quotients.push(quotient);
    let mut remainder_degree:i64 = (remainder.len() as i64) - 1;
    let mut divisor_degree:i64 = (divisor.len() as i64) - 1;
    let mut pair = (Vec::new(), Vec::new());
    if (remainder.len() == 0) {
        pair = (quotients.clone(), remainder.clone());
    }
    else { 
        if (remainder_degree >= divisor_degree) {
            pair = long_division(&mut remainder,divisor,modulus, quotients); 
        } else {
            pair = (quotients.clone(),remainder.clone());
        }
    }
    pair
}
// returns (quotient, remainder)
pub fn divide_once(dividend:&Vec<i64>, 
               divisor:&Vec<i64>, 
               modulus:i64) -> (Vec<i64>,Vec<i64>) { 
    let mut modulus = modulus;
    let mut dividend_degree:i64 = dividend.len() as i64 -1;
    let mut divisor_degree:i64  = divisor.len() as i64 -1;
    let mut quotient_degree:i64 = 0;
    if (dividend_degree == 0 && divisor_degree == 0) {
        quotient_degree = 0;
    } else {
        quotient_degree = dividend_degree - divisor_degree;
    }
    let mut dividend_coefficient = dividend[0];
    let mut divisor_coefficient  = divisor[0];
    let mut quotient_coefficient = find_coefficient(divisor_coefficient, dividend_coefficient, modulus);
    let mut quotient:Vec<i64> = Vec::new();
    quotient.push(quotient_coefficient);
    while quotient_degree > 0 {
        quotient.push(0);
        quotient_degree -= 1;
    }
    let mut product = polynomial_multiply(&quotient,&divisor);
    make_polynomial_negative(&mut product);
    mod_polynomial_coefficients(&mut product, modulus);
    let mut remainder = polynomial_add(&product,&dividend);
    mod_polynomial_coefficients(&mut remainder, modulus);
    remove_leading_zeros(&mut remainder);
    (quotient,remainder)
}

pub fn remove_leading_zeros(v:&mut Vec<i64>) {
    while (v.len() != 0 && v[0] == 0) {
        v.remove(0);
    }
}

pub fn mod_polynomial_coefficients(v:&mut Vec<i64>, p:i64) {
    for coefficient in v {
        *coefficient = find_mod(*coefficient,p);
    } 
}

pub fn find_mod(x:i64,p:i64) -> i64 {
    let mut result:i64 = 0;
    if (x >=0) {
       result = x % p
    } else {
        let mut i:i64 = -1;
        while true {
            if (i*p) <= x {
                result = (i*p*-1)+x;
                break;
            } else {
                i -= 1;
            }
        }
    }
    result
}

pub fn make_polynomial_negative(v:&mut Vec<i64>) {
    for coefficient in v {
            *coefficient *= -1;
    } 
}

// divisor, dividend, modulus
pub fn find_coefficient(a:i64, b:i64, q:i64) -> i64 {
    let mut result: i64 = 0;
    let mut i = 0;
    while (true) {
        if (a * i % q) == b {
            result = i;
            break;
        }
         i += 1;
        if (i > 100) {
            println!("Error: no coefficient was found.");
            process::exit(1);        
        }
    }
    result
}

pub fn multiply(a:&Vec<i64>, b:&Vec<i64>, p:i64) -> Vec<i64> {
    let mut product = polynomial_multiply(a,b);
    mod_polynomial_coefficients(&mut product, p);
    product
}

// polynomials are arrays of coefficients from high order exponent to low order exponent
pub fn polynomial_multiply(a:&Vec<i64>, b:&Vec<i64>) -> Vec<i64> {
    let mut x:Vec<i64> = Vec::new();
    let mut y:Vec<i64> = Vec::new();
    if (a.len() >= b.len()) { // a is longer
        x = a.clone(); // make x the longer one
        y = b.clone();
    } else { // b is longer
        x = b.clone(); // make x the longer one
        y = a.clone();
    }
        x.reverse();
        y.reverse();
    let mut rows:Vec<Vec<i64>> = Vec::new();
    for i in 0..y.len() {
        let mut row:Vec<i64> = Vec::new();
        for k in 0..i {
            row.push(0);
        }
        for j in 0..x.len() {
            let mut coeff = y[i] * x[j];
            row.push(coeff);
        }
        rows.push(row);
    }
    let mut result = add_coefficients(rows);
    result.reverse();
    result
}

pub fn polynomial_add_multiple(polynomials:&Vec<Vec<i64>>) -> Vec<i64> {
    let mut sum:Vec<i64> = vec![0];
    for polynomial in polynomials {
        sum = polynomial_add(&sum,polynomial);
    }
    sum
}

// reverse them first and get a reversed result from add_coefficients
pub fn polynomial_add(a:&Vec<i64>, b:&Vec<i64>) -> Vec<i64> {
    let mut x = a.clone();
    let mut y = b.clone();
    x.reverse();
    y.reverse();
    let mut sum  = add_coefficients(vec![x,y]);
    sum.reverse();
    sum
}

// adds rows of coefficients for multiply function
pub fn add_coefficients(rows:Vec<Vec<i64>>) -> Vec<i64> {
    let mut max_row_length = 0;
    for row in &rows {
        if row.len() > max_row_length {
            max_row_length = row.len();
        }
    }
    let mut result:Vec<i64> = Vec::new();
    for i in 0..max_row_length {
        let mut column_sum = 0;
        for j in 0..rows.len() {
            column_sum += match rows[j].get(i) {
                Some(value) => *value,
                None        => 0
            }
        }
        result.push(column_sum);
    }
    result
}

// takes array from high exponent to low exponent
pub fn get_polynomial_string(v:&Vec<i64>) -> String {
    let mut result = String::new();
    let mut exponent:i64 = v.len() as i64 - 1;
    for coefficient in v {
        if (exponent != 0) {
            result.push_str(&format!("{}x^{} +", coefficient, exponent));
            exponent -= 1;
        } else {
            result.push_str(&format!("{}", coefficient));
        }
    }       
    result
}