Implementing McEliece Cryptosystem

Purpose:

The purpose of this deliverable was to implement the public key system from [McEliece1978]. By working on this deliverable I gained an understanding of finite field arithmetic.

Sources:

[1] Hankerson, Darrel, Alfred J. Menezes, and Scott Vanstone. Guide to elliptic curve cryptography. Springer Science & Business Media, 2006.

   For the following files, I used section 2.3 "Binary field arithmetic" which starts on page 47:

   In GF_2m_arithmetic.rs:
      add()         implemented pseudocode on page 47
      multiply()    implemented pseudocode on page 48
      reduce()      implemented modified pseudocode on page 53
      gcd()         implemented pseudocode on page 58
      inverse()     implemented pseudocode on page 58

   In GF_2mt_arithmetic.rs (all techniques had to be modified):
      add()         implemented modified pseudocode on page 47
      square()      implemented modified pseudocode on page 52
      reduce()      implemented modified pseudocode on page 53
      gcd()         implemented modified pseudocode on page 58
      get_bx_ax()   implemented modified pseudocode on page 58
      inverse()     implemented modified pseudocode on page 58

[2] Stallings, William. Cryptography and Network Security Principles and Practice. Sixth edition. Boston: Pearson, 2014.

    There is an explanation at the bottom of page 164 under "Multiplication" that helped me understand 
    the techniques used in [1].

[3] Wikipedia rref

Code:

main.rs

mod GF_2m_arithmetic;
mod GF_2mt_arithmetic;
mod matrix;
mod encrypt_decrypt;
use rand::{thread_rng, Rng};
use peroxide::fuga::*;
use std::process;
use itertools::Itertools;
use rand::seq::SliceRandom;
use GF_2m_arithmetic  as f1; // f1 for "field 1" 
use GF_2mt_arithmetic as f2; // f2 for "field 2"

fn main() {

    // test 1
    let n = 16;
    let m = (n as f64).log2() as u64;
    let t = 2;
    let length = 16;
    check_precondition(length, t, m);
    let mut f = f1::get_primitive_polynomial(m);
    let mut g = f2::get_primitive_polynomial(t, f);
    println!("test 1 polynomials:");
    println!("{}", f1::string(f));
    println!("{}", f2::string(&g));
    let mut L = f1::create_L(&g, f, length as usize);
    let (mut G, mut dimension, mut length) = matrix::generator(&g, f, &L);
    let mut S = matrix::dense_nonsingular_matrix(dimension as u64);
    let mut P = matrix::permutation_matrix(length as u64); 
    let mut generator = matrix::mod2_matrix(&(S.clone() * G.clone() * P.clone()), dimension as usize, length as usize);
    let (correct, incorrect, total) = run_tests(dimension, length, t, &generator, &S, &G, &P, &g, f, &L); 
    println!("test 1 results:");
    println!("   {}/{} messages decrypted correctly", correct, total);
    println!("   {}/{} messages decrypted incorrectly\n", incorrect, total);

    // test 2
    let n = 16;
    let m = (n as f64).log2() as u64;
    let t = 3;
    let length = 16;
    check_precondition(length, t, m);
    let mut f = f1::get_primitive_polynomial(m);
    let mut g = f2::get_primitive_polynomial(t, f);
    println!("test 2 polynomials:");
    println!("{}", f1::string(f));
    println!("{}", f2::string(&g));
    let mut L = f1::create_L(&g, f, length as usize);
    let (mut G, mut dimension, mut length) = matrix::generator(&g, f, &L);
    let mut S = matrix::dense_nonsingular_matrix(dimension as u64);
    let mut P = matrix::permutation_matrix(length as u64); 
    let mut generator = matrix::mod2_matrix(&(S.clone() * G.clone() * P.clone()), dimension as usize, length as usize);
    let (correct, incorrect, total) = run_tests(dimension, length, t, &generator, &S, &G, &P, &g, f, &L); 
    println!("test 2 results:");
    println!("   {}/{} messages decrypted correctly", correct, total);
    println!("   {}/{} messages decrypted incorrectly\n", incorrect, total);

    // test 3
    let n = 32;
    let m = (n as f64).log2() as u64;
    let t = 2;
    let length = 20;
    check_precondition(length, t, m);
    let mut f = f1::get_primitive_polynomial(m);
    let mut g = f2::get_primitive_polynomial(t, f);
    println!("test 3 polynomials:");
    println!("{}", f1::string(f));
    println!("{}", f2::string(&g));
    let mut L = f1::create_L(&g, f, length as usize);
    let (mut G, mut dimension, mut length) = matrix::generator(&g, f, &L);
    let mut S = matrix::dense_nonsingular_matrix(dimension as u64);
    let mut P = matrix::permutation_matrix(length as u64); 
    let mut generator = matrix::mod2_matrix(&(S.clone() * G.clone() * P.clone()), dimension as usize, length as usize);
    let (correct, incorrect, total) = run_tests(dimension, length, t, &generator, &S, &G, &P, &g, f, &L); 
    println!("test 3 results:");
    println!("   {}/{} messages decrypted correctly", correct, total);
    println!("   {}/{} messages decrypted incorrectly\n", incorrect, total);

    // test 4
    let n = 32;
    let m = (n as f64).log2() as u64;
    let t = 3;
    let length = 20;  
    check_precondition(length, t, m);
    let mut f = f1::get_primitive_polynomial(m);
    let mut g = f2::get_primitive_polynomial(t, f);
    println!("test 4 polynomials:");
    println!("{}", f1::string(f));
    println!("{}", f2::string(&g));
    let mut L = f1::create_L(&g, f, length as usize);
    let (mut G, mut dimension, mut length) = matrix::generator(&g, f, &L);
    let mut S = matrix::dense_nonsingular_matrix(dimension as u64);
    let mut P = matrix::permutation_matrix(length as u64); 
    let mut generator = matrix::mod2_matrix(&(S.clone() * G.clone() * P.clone()), dimension as usize, length as usize);
    let (correct, incorrect, total) = run_tests(dimension, length, t, &generator, &S, &G, &P, &g, f, &L); 
    println!("test 4 results:");
    println!("   {}/{} messages decrypted correctly", correct, total);
    println!("   {}/{} messages decrypted incorrectly\n", incorrect, total);

    // test 5
    let n = 64;
    let m = (n as f64).log2() as u64;
    let t = 2;
    let length = 20;  
    check_precondition(length, t, m);
    let mut f = f1::get_primitive_polynomial(m);
    let mut g = f2::get_primitive_polynomial(t, f);
    println!("test 5 polynomials:");
    println!("{}", f1::string(f));
    println!("{}", f2::string(&g));
    let mut L = f1::create_L(&g, f, length as usize);
    let (mut G, mut dimension, mut length) = matrix::generator(&g, f, &L);
    let mut S = matrix::dense_nonsingular_matrix(dimension as u64);
    let mut P = matrix::permutation_matrix(length as u64); 
    let mut generator = matrix::mod2_matrix(&(S.clone() * G.clone() * P.clone()), dimension as usize, length as usize);
    let (correct, incorrect, total) = run_tests(dimension, length, t, &generator, &S, &G, &P, &g, f, &L); 
    println!("test 5 results:");
    println!("   {}/{} messages decrypted correctly", correct, total);
    println!("   {}/{} messages decrypted incorrectly\n", incorrect, total);
}

fn check_precondition(length:u64, t:u64, m:u64) {
    if length <= t*m {
        println!("Error: length must be > t*m");
        process::exit(1);
    }
}

fn run_tests(dimension:u64, length:u64, t:u64, generator:&Matrix, S:&Matrix, G:&Matrix, 
             P:&Matrix, g:&Vec<u64>, f:u64, L:&Vec<u64>) -> (u64, u64, u64) {
    let mut messages      = all_possible_messages(dimension);
    let mut error_vectors = all_possible_error_vectors(length, t); // takes too long to do all of them
    let mut error_indices: Vec<u32> = (0..error_vectors.len() as u32).collect();
    let mut number_of_error_vectors = 3;
    let mut chosen_error_indices:Vec<u32> = error_indices.choose_multiple(&mut rand::thread_rng(), 
                                                  number_of_error_vectors).cloned().collect();
    println!("chosen_error_indices = {:?}", chosen_error_indices);
    let mut correct_count = 0;
    let mut incorrect_count = 0;
    for message in &messages {
        for error_index in &chosen_error_indices {
            let mut y = encrypt_decrypt::encrypt(message, &generator, length as usize, 
                                                 &error_vectors[(*error_index) as usize]);
            let mut decrypted = encrypt_decrypt::decrypt(&y, &S, &G, &P, &g, f, &L, length, dimension);
            if is_equal(&message, &decrypted, dimension as u64) {
                correct_count += 1;
            } else {
                incorrect_count += 1;
            }
        }
    } 
    (correct_count, incorrect_count, messages.len() as u64 * number_of_error_vectors as u64)
}    

// check if two vectors are equal
fn is_equal(a:&Matrix, b:&Matrix, width:u64) -> bool {
    for i in 0..width as usize {
        if a[(0,i)] != b[(0,i)] {
            return false;
        } 
    }
    return true;
}

// returns vector of all Choose(width, group_size) combinations of possible error vector
fn all_possible_error_vectors(width:u64, group_size:u64) -> Vec<Matrix> {
    let it = (0..width).combinations(group_size as usize);
    let mut vectors = Vec::new();
    for element in it {
        let mut vec = zeros(1 as usize, width as usize);
        for index in element {
            vec[(0 as usize,index as usize)] = 1 as f64;
        }
        vectors.push(vec);
    }
    vectors
}

// returns vector of all possible binary vectors of width width
fn all_possible_messages(width:u64) -> Vec<Matrix> {
    let mut vectors = Vec::new();
    for i in 0..2_u64.pow(width as u32) {
        let mut vec = f1::vec(i, 2_u64.pow(width as u32));
        let mut matrix = zeros(1 as usize, width as usize);
        for j in 0..vec.len() {
            matrix[(0 as usize, j as usize)] = vec[j as usize];
        } 
        vectors.push(matrix);
    }
    vectors
}

GF_2m_arithmetic.rs

// GF_2m_arithmetic.rs 
// Used as "f1" for "field 1"
// This file contains:
//    - arithmetic operations add, multiply, reduce, inverse, power, and gcd for elements of GF(2^m)
//    - method to generate a primitive polynomial of degree m over GF(2)
// Explanation: 
//    - elements of GF(2^m) are polynomials of degree m-1 over GF(2).
//    - since these polynomials can be represented by binary strings, 
//      they can be represented by by binary u64s, with higher order 
//      coefficients corresponding to higher order bits.
//    - elements of GF(2^m) are represented by the names a, b, or u, 
//      and the primitive polynomial is represented by f

use crate::GF_2mt_arithmetic as f2;
use std::process;
use rand::{thread_rng, Rng};
use peroxide::fuga::*;
use itertools::Itertools;

// Returns gcd of two polynomials a and b 
pub fn gcd(a:u64, b:u64) -> u64 {
    let mut u = a.clone();
    let mut v = b.clone();
    let mut g1 = 1 as u64;
    let mut g2 = 0 as u64;
    let mut h1 = 0 as u64;
    let mut h2 = 1 as u64;
    while u != 0 {
        let mut j:i64 = deg(u) as i64 - deg(v) as i64;
        if j < 0 {
            swap(&mut u, &mut v);
            swap(&mut g1, &mut g2);
            swap(&mut h1, &mut h2);
            j *= -1;
        }
        let mut p1 = v.clone();
        let mut p2 = g2.clone();
        let mut p3 = h2.clone();
        shift_left_i(&mut p1, j as u64);
        shift_left_i(&mut p2, j as u64);
        shift_left_i(&mut p3, j as u64);
        u = add(u, p1);
        g1 = add(g1, p2);
        h1 = add(h1, p3);
    }
    let mut d = v.clone();
    d
}
        

// Reduces polynomial u modulo f
pub fn reduce(u:&mut u64, f:u64) {
    let m = deg(f);
    let n = deg(*u);
    let mut rz = rz(f);
    shift_left_i(&mut rz, 63-m);
    for i in (m..64).rev() {
        if 2_u64.pow(i as u32) & (*u) != 0 {
            *u = add(*u, rz);
            *u = add(*u, 2_u64.pow(i as u32));
        } 
        shift_right(&mut rz);
    }
}

// Return u^i modulo f for u a polynomial in GF(2^m)
pub fn pow(u:u64, i:u64, f:u64) -> u64 {
    let mut p = 1;
    for i in 0..i {
        p = multiply(p,u,f);
    }
    p
}

// Returns the vector representation of polynomial u
// The low order bits go into high order vector index
pub fn vec(u:u64, f:u64) -> Vec<f64> {
    let mut vec = vec![0 as f64; deg(f) as usize];
    let mut u = u.clone();
    for i in (0..vec.len()).rev() {
        if u%2 == 0 {
            vec[i] = 0 as f64;
        } else {
            vec[i] = 1 as f64;
        }
        shift_right(&mut u);
    }
    vec
}

// Returns the inverse of a modulo f
pub fn inverse(a:u64, f:u64) -> u64 {
    if a == 0 as u64 {
        println!("Error in f1::inverse - a was 0");
        process::exit(1);
    }
    let mut u = a.clone();
    let mut v = f.clone();
    let mut g1 = 1;
    let mut g2 = 0;
    while u != 1 {
        if u == 0 {
            println!("Error in f1: polynomial was not irreducible");
            println!("a = {:#066b}, f = {:#066b}", a, f);
            process::exit(1);
        } 
        let mut j:i64 = deg(u) as i64 - deg(v) as i64;
        if j < 0 as i64 {
            swap(&mut u, &mut v);            
            swap(&mut g1, &mut g2);            
            j *= -1 as i64
        }
        let mut p1 = v.clone();
        let mut p2 = g2.clone();
        shift_left_i(&mut p1, j as u64);
        shift_left_i(&mut p2, j as u64);
        u = add(u, p1);
        g1 = add(g1, p2);
    }
    return g1;
}

// Swaps two polynomial references
fn swap(a:&mut u64, b:&mut u64) {
    let mut temp = *a;
    *a = *b;
    *b = temp;
}

// Right-to-left shift-and-add field multiplication 
// Multiplies two polynomials a, b while keeping result less than degree of f
// Precondition: a, b in GF(2^m)
pub fn multiply(a:u64, b:u64, f:u64) -> u64 {
    let mut m = deg(f);
    let mut b = b.clone();
    let mut c = 0;
    if a%2 == 1 {
        c = b;
    }
    for i in 1..m {
        shift_left(&mut b);
        if deg(b) == m {
            b = add(b, rz(f));
            b = add(b, 2_u64.pow(m as u32) as u64);
        }
        if bit(a, i) == 1 {
            c = add(c, b);
        }
    }
    c
}

// Returns true if coefficient at position bit is 1 in polynomial a, false otherwise
fn bit(a:u64, bit:u64) -> u64 {
    if a & (1 << bit) != 0 {
        return 1;
    } else {
        return 0;
    }
}

// Returns poynomial f, except with highest order coefficient removed
// Used with the following facts in multiply, reduce:
// f(z) = z^m + r(z) 
// z^m mod f(z) = r(z)
pub fn rz(f:u64) -> u64 {
    let mut x = 1;
    shift_left_i(&mut x, deg(f));
    add(f, x)   
}

// Adds two polynomials a and b over GF(2) - same as bitwise XOR
pub fn add(a:u64, b:u64) -> u64 {
    a ^ b
}

// Returns degree of polynomial a
pub fn deg(a:u64) -> u64 {
    if a == 0 {
        return 0;
    }
    return ((a as f64).log2().floor()+1 as f64) as u64 - 1;
}

// Returns polynomial a as binary string
pub fn string(a:u64) -> String {
    return format!("{:#066b}", a);
}

// Returns polynomial a as a text string: {0,1}x^m-1 + ... + {0,1}x + {0,1}
pub fn poly_string(a:u64) -> String {
    let mut result = String::from("");
    let mut deg = deg(a) as i64;
    while deg >= 0 {
        if a & 2_u64.pow(deg as u32) != 0 {
            if deg != 0 {
                result.push_str(&format!("1x^{} + ", deg));
            } else {
                result.push_str(&format!("1x^{}", deg));
            }
        } else {
            if deg != 0 {
                result.push_str(&format!("0x^{} + ", deg));
            } else {
                result.push_str(&format!("0x^{}", deg));
            }
        }
        deg -= 1;
    }
    result
}

// Divide polynomial a by x
pub fn shift_right(a:&mut u64) -> u64 {
    *a = *a >> 1;
    a.clone()
}

// Multiply polynomial a by x
pub fn shift_left(a:&mut u64) -> u64 {
    *a = *a << 1;
    a.clone()
}

// Divide polynomial a by x^i
pub fn shift_right_i(a:&mut u64, i:u64) -> u64 {
    for i in 0..i {
        shift_right(a); 
    }
    a.clone()
}

// Multiply polynomial a by x^i 
pub fn shift_left_i(a:&mut u64, i:u64) -> u64 {
    for i in 0..i {
        shift_left(a); 
    }
    a.clone()
}

// Returns true if f is irreducible over GF(2), false otherwise
// Precondition: f has degree at most 12 
fn is_irreducible(f:u64) -> bool {
    let mut n = deg(f);
    let mut q = 2 as u64;
    for i in 1..(n as f64/2 as f64).floor() as u64 + 1 {
        let mut exp = q.pow(i as u32);
        let mut rhs = pow(2 as u64, exp, f);
        let mut rhs = add(rhs, 2 as u64);
        let mut gcd = gcd(f, rhs);
        if gcd != 1 {
            return false;
        }
    }
    return true;
}

// Returns true if f is a primitive polynomial, false otherwise
// Precondition: f is an irreducible polynomial of degree at most 6 
fn is_primitive(f:u64) -> bool {
    let m = deg(f);
    let mut n = 2_u64.pow(m as u32)-1; // n = p^m - 1
    let mut remainder = add(2_u64.pow(n as u32), 1 as u64); // x^n + 1
    reduce(&mut remainder, f);
    if remainder != 0 {
        return false;
    }
    n -= 1;
    while n > 0 {
        let mut remainder = add(2_u64.pow(n as u32), 1 as u64); // x^n + 1
        reduce(&mut remainder, f);
        if remainder == 0 {
            return false;
        }
        n -= 1;
    }
    return true;
}   

// Returns a primitive (irreducible) polynomial of degree m
pub fn get_primitive_polynomial(m:u64) -> u64 {
    let mut f = 2_u64.pow(m as u32); // degree m 
    f = add(f, 1 as u64); // constant is 1
    let mut index = rand::thread_rng().gen_range(1..m); // choose 3rd index
    f = add(f, 2_u64.pow(index as u32));
    let mut count = 0;
    while !is_irreducible(f) || !is_primitive(f) {
        f = add(f, 2_u64.pow(index as u32)); // undo third index 
        let mut index = rand::thread_rng().gen_range(1..m); // choose new 3rd index
        f = add(f, 2_u64.pow(index as u32)); // add 3rd index
        count += 1;
        if count > 100 {
            println!("Could not find primitive polynomial over GF(2) of degree {}", m);
            process::exit(1);
        }
    }   
    return f
}

// Returns length random elements of GF(2^m) in a shuffled vector
// If one of them is a root, polynomial g was not irreducible
pub fn create_L(g:&Vec<u64>, f:u64, length:usize) -> Vec<u64> {
    let mut L = Vec::new();
    for i in 0..2_u64.pow(deg(f) as u32) {
        if f2::compute(g, f, i) == 0 {
            println!("Error in choosing L: g(x) is not irreducible over GF(2^{})", deg(f));
            process::exit(1);
        }
        L.push(i);
    }
    L.shuffle(&mut thread_rng());   
    return (&L[0..length]).to_vec();
}

GF_2mt_arithmetic.rs

// GF_2mt_arithmetic.rs
// Used as "f2" for "field 2"
// This file contains:
//    - arithmetic operations add, multiply, reduce, inverse, square root, and gcd for elements of GF(2^mt)
//    - method to generate a primitive polynomial of degree t over GF(2^m)
// Explanation: 
//    - elements of GF(2^mt) are polynomials of degree t-1 over GF(2^m).
//    - they are represented as vectors of u64s, with the u64 being the coefficient in GF(2^m),
//      and low order indices corresponding to low order coefficients
//    - elements of GF(2^mt) are represented by the names a or b, with the primitive polynomial 
//      represented by g, and the primitive polynomial for the field GF(2^m) represented by f
//    - for visualization, the higher order coefficients are on the left so that a multiply by x
//      is a shift left 

use crate::GF_2m_arithmetic as f1;
use std::process;
use std::mem;
use rand::{thread_rng, Rng};
use peroxide::fuga::*;
use itertools::Itertools;

// Returns gcd of a and b
pub fn regular_gcd(a:&Vec<u64>, b:&Vec<u64>, f:u64) -> Vec<u64> {
    let mut u = a.clone();
    let mut v = b.clone();
    let mut g1 = vec![1 as u64];
    let mut g2 = vec![0 as u64];
    let mut h1 = vec![0 as u64];
    let mut h2 = vec![1 as u64];
    let mut t = deg(&b);
    let mut count = 0;
    while !is_zero(&u) {
        /*if count%100 == 0 {
            println!("f2::gcd iteration {}", count);
        }*/
        if count > 5000 {
            println!("Error: max iterations in f2 gcd");
            process::exit(1);
        }
        if is_zero(&v) {
            println!("Error in f2 gcd: g = {}", string(&b));
            process::exit(1);
        }            
        let mut j:i64 = deg(&u) as i64 - deg(&v) as i64;
        if j<0 as i64 {
            mem::swap(&mut u, &mut v);
            mem::swap(&mut g1, &mut g2);
            mem::swap(&mut h1, &mut h2);
            j *= -1 as i64;
        }
        let mut p1 = v.clone();
        let mut p2 = g2.clone();
        let mut p3 = h2.clone();
        shift_left_i(&mut p1, j as usize);
        shift_left_i(&mut p2, j as usize);
        shift_left_i(&mut p3, j as usize);
        let mut u_coeff = u[u.len()-1];
        let mut v_coeff = v[v.len()-1];
        let mut coeff = f1::multiply(f1::inverse(v_coeff, f), u_coeff, f); 
        p1 = multiply_by_constant(&p1, coeff, f);
        p2 = multiply_by_constant(&p2, coeff, f);
        p3 = multiply_by_constant(&p3, coeff, f);
        let mut len1 = u.len();
        u  = add(&u, &p1);
        let mut len2 = u.len();
        g1 = add(&g1, &p2);
        h1 = add(&h1, &p3);
        count += 1;
    }
    let mut d = v;
    d    
}

// Returns (b(x), a(x)) such that v(x)*b(x) = a(x) mod g(x)
// Pass in (v(x), g(x)) and use extended euclidean algorithm
// Used in decode 
pub fn get_bx_ax(vx:&Vec<u64>, gx:&Vec<u64>, f:u64) -> (Vec<u64>, Vec<u64>) {
    let mut u = vx.clone();
    let mut v = gx.clone();
    let mut g1 = vec![1 as u64];
    let mut g2 = vec![0 as u64];
    let mut h1 = vec![0 as u64];
    let mut h2 = vec![1 as u64];
    let mut t = deg(&gx);
    let mut u_bound = ((t as f64)/(2 as f64)).floor() as u64;
    let mut g1_bound = ((t as f64 - 1 as f64)/(2 as f64)).floor() as u64;
    let mut count = 0;
    while !is_zero(&u) {
        if count > 100 {
            println!("Error: max iterations in f2 get_bx_ax");
            process::exit(1);
        }
        if is_zero(&v) {
            println!("Error in f2 get_bx_ax: g = {}", string(&gx));
            process::exit(1);
        }            
        let mut j:i64 = deg(&u) as i64 - deg(&v) as i64;
        if j<0 as i64 {
            mem::swap(&mut u, &mut v);
            mem::swap(&mut g1, &mut g2);
            mem::swap(&mut h1, &mut h2);
            j *= -1 as i64;
        }
        let mut p1 = v.clone();
        let mut p2 = g2.clone();
        let mut p3 = h2.clone();
        shift_left_i(&mut p1, j as usize);
        shift_left_i(&mut p2, j as usize);
        shift_left_i(&mut p3, j as usize);
        let mut u_coeff = u[u.len()-1];
        let mut v_coeff = v[v.len()-1];
        let mut coeff = f1::multiply(f1::inverse(v_coeff, f), u_coeff, f); 
        p1 = multiply_by_constant(&p1, coeff, f);
        p2 = multiply_by_constant(&p2, coeff, f);
        p3 = multiply_by_constant(&p3, coeff, f);
        u  = add(&u, &p1);
        g1 = add(&g1, &p2);
        h1 = add(&h1, &p3);
        if deg(&u) > u_bound || deg(&g1) > g1_bound {
            if deg(&v) <= u_bound && deg(&g2) <= g1_bound { // check prev 
                break;
            }
        }
        count += 1;
    }
    return (g2, v);
}

// Returns the inverse of a modulo g
// Coefficient arithmetic done in GF(2^m) using f
pub fn inverse(a:&Vec<u64>, g:&Vec<u64>, f:u64) -> Vec<u64> {
    let mut u = a.clone();
    let mut v = g.clone();
    let mut g1 = vec![1 as u64];
    let mut g2 = vec![0 as u64];
    let mut count = 0;
    while !is_constant(&u) {
        if count > 100 {
            println!("Error: max iterations in f2 inverse");
            process::exit(1);
        }
        if is_zero(&u) || is_zero(&v) {
            println!("Error in f2: g = {}", string(&g));
            process::exit(1);
        }            
        let mut j:i64 = deg(&u) as i64 - deg(&v) as i64;
        if j<0 as i64 {
            mem::swap(&mut u, &mut v);
            mem::swap(&mut g1, &mut g2);
            j *= -1 as i64;
        }
        let mut p1 = v.clone();
        let mut p2 = g2.clone();
        shift_left_i(&mut p1, j as usize);
        shift_left_i(&mut p2, j as usize);
        let mut u_coeff = u[u.len()-1];
        let mut v_coeff = v[v.len()-1];
        let mut p1_coeff = f1::multiply(f1::inverse(v_coeff, f), u_coeff, f); 
        p1 = multiply_by_constant(&p1, p1_coeff, f);
        p2 = multiply_by_constant(&p2, p1_coeff, f);
        u = add(&u, &p1);
        g1 = add(&g1, &p2);
        count += 1;
    }
    let mut inv = f1::inverse(u[0], f);
    g1 = multiply_by_constant(&g1, inv, f);
    return g1;
}

// Returns true if polynomial is the constant 1, false otherwise
pub fn is_one(a:&Vec<u64>) -> bool {
    return a.len() == 1 && a[a.len()-1] == 1;
}

// Returns a^2 modulo g, coefficients reduced modulo f 
pub fn square(a:&Vec<u64>, g:&Vec<u64>, f:u64) -> Vec<u64> {
    let mut b = Vec::new();
    let mut i = 0;
    while i < a.len() {
        b.push(f1::multiply(a[i],a[i],f));
        b.push(0 as u64);
        i += 1;
    }
    b.pop(); 
    b
}

// Returns square root of a modulo g, coefficients reduced modulo f
pub fn sqrt(a:&Vec<u64>, g:&Vec<u64>, f:u64) -> Vec<u64> {  
    let mut t = deg(&g);
    let mut m = f1::deg(f);
    let mut n = t*m-1;
    let mut result = a.clone();
    for i in 0..n {
        result = square(&result, &g, f);
        reduce(&mut result, &g, f);
    }
    result
}

// Reduces polynomial a modulo polynomial g, coefficients reduced modulo f 
pub fn reduce(a:&mut Vec<u64>, g:&Vec<u64>, f:u64) { 
    let mut deg_a = deg(&a) as usize;
    let mut deg_g = deg(&g) as usize;
    if deg_a < deg_g {
        return;
    }
    let mut rz = g.clone();
    rz.pop();
    shift_left_i(&mut rz, deg_a - deg_g);
    for i in (deg_g..deg_a+1).rev() {
        if a[i] != 0 {
            let mut to_add = multiply_by_constant(&rz, a[i], f);
            *a = add(&a, &to_add);
            a.pop();
        } else {
            a.pop();
        }
        shift_right(&mut rz);
    }
}

// Multiplies polynomial a by constant c and reduces coefficients modulo f
pub fn multiply_by_constant(a:&Vec<u64>, c:u64, f:u64) -> Vec<u64> {
    let mut b = a.clone();
    for i in 0..a.len() {
        b[i] = f1::multiply(b[i], c, f);
    }
    b
}

// Returns true is polynomial a is only a constant value
pub fn is_constant(a:&Vec<u64>) -> bool {
    if a.len() == 1 && a[0] != 0 {
        return true;
    }
    return false;
}

// Returns true is polynomial a is zero, false otherwise
pub fn is_zero(a:&Vec<u64>) -> bool {
    for i in 0..a.len() {
        if a[i] != 0 {
            return false;
        } 
    }
    return true;
}

// Multiplies polynomial a by x
pub fn shift_left(a:&mut Vec<u64>) {
    a.insert(0,0 as u64);
}      

// Divides polynomial a by x
pub fn shift_right(a:&mut Vec<u64>) {
    a.remove(0);
}

// Multiplies polynomial a by x^i
pub fn shift_left_i(a:&mut Vec<u64>, i:usize) {
    let mut i = i;
    while i > 0 {
        shift_left(a);
        i -= 1;
    }
}

// Adds two polynomials a and b in GF(2^tm)        
pub fn add(a:&Vec<u64>, b:&Vec<u64>) -> Vec<u64> {
    let mut c = Vec::new();
    let mut i = 0;
    while i < a.len() && i < b.len() { 
        c.push(f1::add(a[i], b[i]));
        i += 1;
    }
    while i < a.len() {
        c.push(a[i]);
        i += 1;
    }
    while i < b.len() {
        c.push(b[i]);
        i += 1;
    }
    while c[c.len()-1] == 0 && c.len() > 1 {
        c.pop();
    }
    c
}

// Returns degree of polynomial g
pub fn deg(g:&Vec<u64>) -> u64 {
    return g.len() as u64-1;
}


// Plugs x, an element of GF(2^m), into g(x)
// Coefficients are reduced modulo f
pub fn compute(g:&Vec<u64>, f:u64, x:u64) -> u64 {
    let mut sum = 0;
    let mut power = deg(g) as usize;
    while power > 0 {
        let mut variable = f1::pow(x, power as u64, f);
        let mut coeff = g[power];
        let mut prod = f1::multiply(coeff, variable, f);
        sum = f1::add(f1::multiply(g[power], f1::pow(x, power as u64, f), f), sum);
        power -= 1;
    }
    sum = f1::add(sum, g[power]); // add constant term
    sum
}

// Returns polynomial a as a string: {GF(2^m)}x^{t-1} + ... + {GF(2^m)}x + {GF(2^m)}
pub fn string(a:&Vec<u64>) -> String {
    let mut result = String::from("");
    let mut deg = deg(&a) as i64;
    for u in a.iter().rev() {
        result.push_str(&format!("x^{}: {:#066b}", deg, u));
        if deg != 0 {
            result.push_str("\n");
        }
        deg -= 1;
    }
    result
}

// Returns a primitive (irreducible) polynomial of degree t over GF(2^m)
pub fn get_primitive_polynomial(t:u64, f:u64) -> Vec<u64> {
    let mut m = f1::deg(f);
    let mut g = vec![1 as u64];
    shift_left_i(&mut g, t as usize);
    let mut constant = rand::thread_rng().gen_range(1..2_u64.pow(m as u32));
    g[0] = constant;
    let mut index = rand::thread_rng().gen_range(1..t) as usize; // choose third index
    g[index] = rand::thread_rng().gen_range(1..2_u64.pow(m as u32));
    let mut count = 0;
    while !is_irreducible(&g, f) { // || !is_primitive(&g, f) { // takes too long for deg(g) > 2 
        g[0] = 0; // clear constant
        g[index] = 0; // clear index
        let mut constant = rand::thread_rng().gen_range(1..2_u64.pow(m as u32));
        g[0] = constant;
        let mut index = rand::thread_rng().gen_range(1..t) as usize; // choose third index
        g[index] = rand::thread_rng().gen_range(1..2_u64.pow(m as u32));
        count += 1;
        if count > 500 {
            println!("Could not find primitive polynomial over GF(2^{}) of degree {}", m, t);
            process::exit(1);
        }   
    }
    g
}


// Returns true if g is irreducible over GF(2^m), false otherwise
fn is_irreducible(g:&Vec<u64>, f:u64) -> bool {
    let mut n = deg(&g);
    let mut q = 2_u64.pow(f1::deg(f) as u32); // q = 2^m
    for i in 1..(n as f64/2 as f64).floor() as u64 + 1 {
        let mut exp = q.pow(i as u32);
        let mut rhs = vec![1 as u64]; 
        shift_left_i(&mut rhs, exp as usize); // x^{q^i}
        let mut rhs = add(&mut rhs, &vec![0 as u64, 1 as u64]); // x^{q^i}-x
        let mut gcd = regular_gcd(&g, &rhs, f);
        if !is_one(&gcd) {
            return false;
        } 
    }
    return true;
}

// Returns true if g is a primitive polynomial, false otherwise
// Precondition: g is an irreducible polynomial of degree at most 2
// Degree t > 2 takes too long 
fn is_primitive(g:&Vec<u64>, f:u64) -> bool {
    let mut t = deg(&g);
    let mut exp = t * f1::deg(f); // exp = mt
    let mut n = 2_u64.pow(exp as u32)-1; // 2^{mt}-1
    let mut remainder = vec![1 as u64];
    shift_left_i(&mut remainder, n as usize);
    remainder = add(&remainder, &vec![1 as u64]); 
    reduce(&mut remainder, &g, f);
    if !is_zero(&remainder) {
        return false;
    }
    n -= 1;
    while n > 0 {
        let mut remainder = vec![1 as u64];
        shift_left_i(&mut remainder, n as usize);
        remainder = add(&remainder, &vec![1 as u64]); 
        reduce(&mut remainder, &g, f);
        if is_zero(&remainder) {
            return false;
        }
        n -= 1;
    }
    return true;
} 

matrix.rs

// matrix.rs
// This file contains:
//    - method to create generator matrix for a binary irreducible Goppa code
//      using an irreducible polynomial g 
//    - method to generate the matrices S and P for the private key, 
//      used to scramble the generator matrix
//    - method used to convert matrix of elements of GF(2^m) to binary matrix
//    - method to multiply two matrices using field multiplication from GF(2^m)
// Explanation: 
//    - the generator matrix for the code is found by taking the transpose of the nullspace 
//      of the parity check matrix for the code
//    - the parity check matrix is H = xyz where x,y,and z are matrices formed 
//      using g, f, and L

use crate::GF_2m_arithmetic as f1;
use crate::GF_2mt_arithmetic as f2;
use peroxide::fuga::*;
use std::process;

// Returns the generator matrix for the binary irreducible Goppa code corresponding
// to the irreducible polynomial g 
// Returns (generator, height, width)
pub fn generator(g:&Vec<u64>, f:u64, L:&Vec<u64>) -> (Matrix, u64, u64) {
    let t  = f2::deg(g);
    let n  = L.len();
    let m  = f1::deg(f); 
    let x  = x(g, t);
    let y  = y(L, t, n as u64, f);
    let z  = z(&g, &L, n as u64, f);
    let xy = multiply(&x, &y, t as usize, t as usize, t as usize, n, f);
    let h  = multiply(&xy, &z, t as usize, n, n, n, f);
    let mut h  = binary_matrix(&h, t as usize, n, m as usize, f);
    let (generator, dimension, length) = nullspace_t(&mut h, t*m as u64, n as u64); 
    return (generator, dimension, length);
}   

// Append two matrices together
pub fn append(a:&Matrix, b:&Matrix, h1:u64, w1:u64, h2:u64, w2:u64) -> Matrix {
    let mut m = zeros(h1 as usize, (w1+w2) as usize);
    let mut r1 = 0 as usize;
    let mut c1 = 0 as usize;
    while r1 < h1 as usize {
        c1 = 0;
        while c1 < w1 as usize {
            m[(r1,c1)] = a[(r1,c1)];
            c1 += 1;
        }
        r1 += 1;
    }
    let mut r2 = 0 as usize;
    while r2 < h2 as usize {
        let mut c2 = 0 as usize;
        while c2 < w2 as usize {
            m[(r2, c1+c2)] = b[(r2,c2)];
            c2 += 1;
        }
        r2 += 1;
    }
    m
}


// Returns dense nonsingular matrix S for use in private key
// Make sure determinant is 1 or -1 so inverse contains whole numbers
// Also, if matrix is too large with high density, library cannot find the matrix
// Also Rust library sometimes crashes while finding determinant
pub fn dense_nonsingular_matrix(k:u64) -> Matrix {
    let size:u64 = k*k;
    let min_density = (size as f64/2 as f64).ceil() as u64 + 1; // 1/2
    let density = min_density; // crashes for higher density
    let mut det:f64 = 0.0;
    let mut m = zeros(k as usize, k as usize);
    let mut iteration = 0;
    loop {
        m = zeros(k as usize, k as usize);
        let mut count = 0;
        while count < density {
            let c = rand::thread_rng().gen_range(0..k);
            let r = rand::thread_rng().gen_range(0..k);
            if m[(c as usize, r as usize)] != 1 as f64 {
                m[(c as usize, r as usize)] = 1 as f64;
                count += 1;
            }
        }
        if m.det() == 1 as f64 || m.det() == -1 as f64 { // Rust library sometimes crashes here
                break;
        }
        iteration += 1;
    }
    m
}

// Returns random permutation matrix P for use in private key
pub fn permutation_matrix(n:u64) -> Matrix  {
    let mut ones: Vec<u64> = (0..n).collect();
    ones.shuffle(&mut thread_rng());
    let mut m = zeros(n as usize,n as usize);
    let mut i:usize = 0;
    for i in 0..n {
        m[(i as usize,ones[i as usize] as usize)] = 1 as f64;
    }
    m
}

// Convert a matrix of elements of GF(2^m) to a binary matrix
// by converting each element of GF(2^m) to a vertical binary string where
// the top number is the high order bit
pub fn binary_matrix(h:&Matrix, t:usize, d:usize, m:usize, f:u64) -> Matrix {
    let mut b = zeros(t*m, d);
    for c in 0..d {
        for r in 0..t {
            let mut vec = f1::vec(h[(r,c)] as u64, f); // high order to low order
            for i in 0..m {
                b[(r*m+i, c)] = vec[i];
            }
        }
    }
    b
}

// Multiplies two matrices together
// Since the elements of the matrix are elements of GF(2^m), regular
// multiplication cannot be used and field multiplication in GF(2^m) must be used 
pub fn multiply(x:&Matrix, y:&Matrix, xrows:usize, xcols:usize, yrows:usize, ycols:usize, f:u64) -> Matrix {
    if xcols != yrows {
        println!("matrix multiplication error");
        process::exit(1);
    }
    let mut z = zeros(xrows, ycols);
    for r in 0..xrows {
        for c in 0..ycols {
            for t in 0..xcols {
                z[(r,c)] = f1::add(z[(r,c)] as u64, f1::multiply(x[(r,t)] as u64, y[(t,c)] as u64, f)) as f64;
            }
        }
    }
    z
}             

// Return z in (xyz = parity check matrix)
// It is a diagonal matrix where the elements are 
// the inverse of g(l) modulo g(x) for l in L
pub fn z(g:&Vec<u64>, L:&Vec<u64>, d:u64, f:u64) -> Matrix {
    let mut z = zeros(d as usize, d as usize);
    for r in 0..d {
        for c in 0..d {
            if r == c {
                z[(r as usize, c as usize)] = f1::inverse(f2::compute(&g, f, L[r as usize]), f) as f64;
            }
        }
    }
    z
}

// Returns y in (xyz = parity check matrix)
// It contains the elements of L raised to powers and reduced modulo f
pub fn y(L:&Vec<u64>, t:u64, d:u64, f:u64) -> Matrix {
    let mut L = L.clone();
    let mut y = zeros(t as usize, d as usize);
    for c in 0..d {
        y[(0,c as usize)] = 1 as f64;
    }
    for r in 1..t {
        for c in 0..d {
            y[(r as usize,c as usize)] = f1::pow(L[c as usize],r,f) as f64;
        }
    }
    y
}

// Returns x in (xyz = parity check matrix)
// It is a lower triangular matrix made 
// using the coefficients of irreducible polynomial g(x)
pub fn x(g:&Vec<u64>, t:u64) -> Matrix {
    let mut g = g.clone();
    g.remove(0);
    let mut x = zeros(t as usize, t as usize);
    for r in (0..t).rev() {
        for c in 0..t {
            x[(r as usize, c as usize)] = g[c as usize] as f64;
        }
        g.remove(0);
        g.push(0);
    }
    x
}

// Modified rref, used in last decryption step
pub fn rref_for_decrypt(M:&mut Matrix, m:u64, n:u64, rows_to_include:u64, columns_to_include:u64) {
    let mut lead:usize = 0;
    let mut rowCount:usize = rows_to_include as usize;
    let mut columnCount:usize = columns_to_include as usize;
    for r in 0..rowCount {
        if columnCount <= lead {
            return 
        }
        let mut i = r;
        while M[(i, lead)] == 0 as f64{
            i = i+1;
            if rowCount == i {
                i = r;
                lead = lead + 1;
                if columnCount == lead {
                    return 
                }
            }
        }
        if i != r {
            unsafe {
                M.swap(i, r, Row);
            }
        }
        for i in 0..rowCount {
            if i != r && M[(i, lead)] != 0 as f64 {
                for c in 0..(n as  usize) {
                    M[(i,c)] = f1::add(M[(i,c)] as u64, M[(r,c)] as u64) as f64;    
                }
            }
        }
        lead = lead + 1;
    } 
}

// Puts a matrix into reduced row echelon form modulo 2
// Modified from wikipedia rref
pub fn rref_mod_2(M:&mut Matrix, m:u64, n:u64) {
    let mut lead:usize = 0;
    let mut rowCount:usize = m as usize;
    let mut columnCount:usize = n as usize;
    for r in 0..rowCount {
        if columnCount <= lead {
            return 
        }
        let mut i = r;
        while M[(i, lead)] == 0 as f64{
            i = i+1;
            if rowCount == i {
                i = r;
                lead = lead + 1;
                if columnCount == lead {
                    return 
                }
            }
        }
        if i != r {
            unsafe {
                M.swap(i, r, Row);
            }
        }
        for i in 0..rowCount {
            if i != r && M[(i, lead)] != 0 as f64 {
                // add row r to row i
                for c in 0..columnCount {
                    M[(i,c)] = f1::add(M[(i,c)] as u64, M[(r,c)] as u64) as f64;    
                }
            }
        }
        lead = lead + 1;
    } 
}

// Reduces an element of a matrix modulo 2
pub fn mod2(a:f64) -> f64 {
    let result = (a.round() as i64).rem_euclid(2) as f64;
    result
}

// Reduces all the elements of a matrix modulo 2
pub fn mod2_matrix(m:&Matrix, rows:usize, cols:usize) -> Matrix {
    let mut result = m.clone();
    for row in 0..rows {
        for col in 0..cols {
            result[(row, col)] = mod2(result[(row,col)]);
        }
    }
    result
}

// Returns true if row is all zeros, false otherwise
pub fn row_is_zero(M:&Matrix, i:u64) -> bool {
    let mut row = M.row(i as usize);
    for e in row {
        if e != 0 as f64 {
            return false;
        }
    }
    return true;
}
// Check if a column is all zeros except a 1 at position col
pub fn is_col(vec:&Vec<f64>, col:usize) -> bool {
    for i in 0..vec.len() {
        if i == col {
            if vec[i] != 1  as f64 {
                return false;
            }
        } else {
            if vec[i] != 0 as f64 {
                return false;
            }
        }
    }
    return true;
}

// Find column to swap with in swap_columns()
pub fn find_col(M:&Matrix, rows:usize, cols:usize, i:usize) -> (bool, usize) {
    for c in 0..cols {
        let mut vec = M.col(c);
        if is_col(&vec, i) {
            return (true,c );
        }
    }
    return (false,0);
}

// Function to swap columns in rref 
pub fn swap_columns(M:&mut Matrix, rows:u64, cols:u64) {
    for row in 0..rows {
        for col in 0..cols {
            if row == col && M[(row as usize, col as usize)] != 1 as f64 {
                let (mut is_col, index) = find_col(&M, rows as usize, cols as usize, col as usize);
                if is_col {
                    unsafe {
                        M.swap(row as usize, index as usize, Col);
                    }
                } else {
                }
            }
        }
    }
}

// Returns the transpose of the nullspace of matrix M along with its height, width
// Arithmetic done modulo 2
pub fn nullspace_t(M:&mut Matrix, rows:u64, cols:u64) -> (Matrix, u64, u64) {  
    let mut swap = Vec::new();
    rref_mod_2(M, rows, cols);
    // swap columns
    for row in 0..rows {
        for col in 0..cols {
            if row == col && M[(row as usize, col as usize)] != 1 as f64 {
                let (mut is_col, index) = find_col(&M, rows as usize, cols as usize, col as usize);
                if is_col {
                    unsafe {
                        M.swap(row as usize, index as usize, Col);
                        swap.push((row as usize, index as usize));
                    }
                } else {
                }
            }
        }
    }
    // end swap columns
    // find where identity matrix ends
    let mut last_nonzero_row = rows-1;
    for row in (0..rows).rev() {
        if !row_is_zero(&M, row) {
            last_nonzero_row = row;
            break;
        }
    }
    let mut first_column = last_nonzero_row+1;
    // if rref form is an identity matrix, there is no nullspace
    if (cols - first_column) < 1 {
        println!("Error: Nullspace is zero or too small");
        process::exit(1);
    }   
    // nullspace has one row for each column and one column for each column after the identity matrix
    let mut nullspace = zeros(cols as usize, cols as usize-first_column as usize);
    for row in 0..rows {
        for col in 0..(cols as usize - first_column as usize) {
            nullspace[(row as usize, col as usize)] = M[(row as usize, col as usize +first_column as usize)];
        }
    }
    // fill in identity matrix
    let mut row = first_column as usize;
    let mut col = 0 as usize;
    while col < (cols as usize - first_column as usize) {
        nullspace[(row, col)] = 1 as f64;
        col += 1;
        row += 1;
    }
    // swap rows back in reverse order
    if swap.len() != 0 {
        for &(r1, r2) in swap.iter().rev() {
            unsafe {
                nullspace.swap(r1,r2,Row);
            }
        }
    }
    // end swap rows back
    let mut nullspace = nullspace.transpose();
    return (nullspace, cols as u64 - first_column as u64, cols as u64);
}

encrypt_decrypt.rs

// encrypt_decrypt.rs
// This file contains:
//    - methods to encrypt and decrypt a message 
//    - includes Patterson's error correcting algorithm, method to 
//      compute the syndrome, and method to compute error vector 
// Explanation: 
//    - the primitive polynomial for GF(2^mt) is represented by g
//    - the primitive polynomial for GF(2^m) is represented by f
//    - the public key is the scambled generator matrix
//    - the private key is S, G, P, L, g, and f

use crate::GF_2m_arithmetic as f1;
use crate::GF_2mt_arithmetic as f2;
use crate::matrix;
use peroxide::fuga::*;

// Returns encrypted message
// generator is the public key
pub fn encrypt(message:&Matrix, generator:&Matrix, length:usize, error_vector:&Matrix) -> Matrix {
    let mut mG = matrix::mod2_matrix(&(message.clone() * generator.clone()), 1 as usize, length as usize);
    let mut y = matrix::mod2_matrix(&(mG.clone() + error_vector.clone()), 1 as usize, length as usize);
    y
}

// Returns decrypted message
// y is the received, encrypted message
pub fn decrypt(y:&Matrix, S:&Matrix, G:&Matrix, P:&Matrix, g:&Vec<u64>, f:u64, L:&Vec<u64>, 
                                                                    length:u64, dimension:u64) -> Matrix {
    let mut yP_inv = y.clone()*P.inv();
    let mut error_vector = error_vector(&yP_inv, &g, f, &L);
    let mut mSG = matrix::mod2_matrix(&(yP_inv + error_vector), 1 as usize, length as usize);
    let mut find_mS = matrix::append(&G.transpose(), &mSG.transpose(), length, dimension, length, 1 as u64);
    matrix::rref_for_decrypt(&mut find_mS, length, dimension+1, length, dimension);
    let mut mS = zeros(1 as usize, dimension as usize);
    let mut col = dimension as usize;
    for row in 0..dimension {
        mS[(0 as usize,row as usize)] = find_mS[(row as usize, col as usize)];
    }
    let mut S_inv = matrix::mod2_matrix(&(S.inv()), dimension as usize, dimension as usize);
    let mut decoded = matrix::mod2_matrix(&(mS.clone() * S_inv), 1 as usize, dimension as usize);
    decoded
}

// Returns error vector which contains moved errors
pub fn error_vector(yP_inv:&Matrix, g:&Vec<u64>, f:u64, L:&Vec<u64>) -> Matrix {
    let mut syndrome = syndrome(yP_inv, g, f, L); 
    let mut sigma = sigma(&syndrome, g, f);
    let mut error_vector = zeros(1 as usize, L.len() as usize);
    for i in 0..L.len() {
        let mut result = f2::compute(&sigma, f, L[i]);
        if result == 0 as u64 {
            error_vector[(0,i)] = 1 as f64;
        } else {
            error_vector[(0,i)] = 0 as f64;
        }
    }
    error_vector
}

// Returns polynomial to correct errors
// Uses Patterson's error correcting algorithm
pub fn sigma(syndrome:&Vec<u64>, g:&Vec<u64>, f:u64) -> Vec<u64> {
    let syndrome_inv = f2::inverse(syndrome, g, f);
    let mut syndrome_inv_plus_x = f2::add(&syndrome_inv, &vec![0 as u64, 1 as u64]); // low order to high order
    let mut sqrt = f2::sqrt(&syndrome_inv_plus_x, &g, f);
    let (mut bx, mut ax) = f2::get_bx_ax(&sqrt, &g, f);
    let mut bx2 = f2::square(&bx, &g, f);
    let mut ax2 = f2::square(&ax, &g, f);
    f2::shift_left(&mut bx2); // multiply b(x) by x
    let mut sigma = f2::add(&ax2, &bx2);
    sigma
}
    
// Returns syndrome which is a polynomial in GF(2^tm)
pub fn syndrome(yP_inv:&Matrix, g:&Vec<u64>, f:u64, L:&Vec<u64>) -> Vec<u64> {
    let mut sum = vec![0 as u64];
    for i in 0..L.len() {
        if yP_inv[(0,i)] == 1 as f64 {
            let mut denom = vec![L[i], 1]; // low order to high order 
            let mut inverse = f2::inverse(&denom, &g, f);
            sum = f2::add(&sum, &inverse);
        } 
    }
    sum
}