Implementing the McEliece Cryptosystem

Purpose:

The purpose of this deliverable was to implement the public key system from [McEliece1978]. The system is a candidate for post-quantum cryptography.

Code:

main.rs

// main.rs 

mod text;
mod GF_2m_arithmetic;
mod GF_2mt_arithmetic;
mod matrix;
mod encrypt_decrypt;
mod test;
mod system;
use GF_2m_arithmetic as f1;
use GF_2mt_arithmetic as f2;
mod util;


fn main() {

    // Example use
    let mut message = String::from("!@#$%^&*()_+-=1234567890kjfkjhg");
    let (public_key, private_key) = system::create_system(32, 2);
    let iv                        = system::initialization_vector(&public_key); 
    let ciphertext                = system::encrypt(message.clone(), &public_key, &iv);
    let plaintext                 = system::decrypt(&ciphertext, &private_key, &iv);
    if message == plaintext {
        println!("True");
    } else {
        println!("False");
    }
    
}

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 irreducible 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 irreducible 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 {
    if a == 0 || b == 0 {
        return 0;
    }
    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 {
        if v == 0 {
            println!("Error in f1 gcd: b = {}", string(b));
            process::exit(1);
        }
        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
}

// 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;
}

// 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
}

// 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
}


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

// ========== helper methods ==========

// 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
}

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

// 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)   
}

// 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()
}

// ========== polynomial methods ==========

// Ben-Or Irreducibility Test
pub fn is_irreducible(f:u64) -> bool {
    let mut m = deg(f);
    let mut i = 0;
    let mut rhs = 2 as u64; // rhs = x
    rhs = multiply(rhs, rhs, f); // rhs = rhs^2 mod f
    i += 1;
    while i <= (m as f64/2 as f64).floor() as u64 {
        let mut rhs_plus_x = add(rhs,2);
        let mut gcd = gcd(f, rhs_plus_x);
        if gcd != 1 {
            return false;
        }
        rhs = multiply(rhs, rhs, f); // rhs = rhs^2 mod f
        i += 1;
    }
    return true; 
}

pub fn get_random_polynomial(m:u64) -> u64 {
    let mut f = 2_u64.pow(m as u32); // degree m 
    for index in 0..m {
        let mut coeff:u64 = rand::thread_rng().gen_range(0..2); // choose 0 or 1
        if coeff == 1 {
            f = add(f, 2_u64.pow(index as u32));
        }
    }
    f
}

pub fn get_irreducible_polynomial(m:u64) -> u64 {
    let mut f = get_random_polynomial(m);
    // test random polynomial
    let mut count = 0;
    while !is_irreducible(f) { 
        f = get_random_polynomial(m);
        count += 1;
    }
    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);
    }
    for element in L.clone() {
        if element != 0 {
            let mut inverse = inverse(element, f);
            println!("inverse of {} mod {} = {}\n", string(element), string(f), string(inverse));
        }
        
    }
    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 an irreducible 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 irreducible polynomial 
//      represented by g, and the irreducible 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 gcd(a:&Vec<u64>, b:&Vec<u64>, f:u64) -> Vec<u64> {
    if is_zero(&a) || is_zero(&b) {
        return vec![0];
    }
    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 > 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;
}

// 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);
    }
}

// Returns a^2 modulo g - do not reduce because 
// this operation is used in decoding and you don't 
// reduce it there 
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
}

// 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
}

// 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 {
        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
}

// ========== helper methods ==========

// 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 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 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;
    }
}

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

// 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
}

// ========== polynomial methods ==========

// Returns random polynomial of degree t over GF(2^m)
pub fn get_random_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);
    for index in 0..t {
        let mut constant = rand::thread_rng().gen_range(1..2_u64.pow(m as u32));
        g[index as usize] = constant;
    }
    g
}

// Returns irreducible polynomial of degree t over GF(2^m)
pub fn get_irreducible_polynomial(t:u64, f:u64) -> Vec<u64> {
    let mut g = get_random_polynomial(t, f);
    let mut count = 0;
    while !is_irreducible(&g, f) {
        g = get_random_polynomial(t, f);
        count += 1;
    }
    println!("{} tries to find irreducible g", count);
    g
}

// Ben-Or irreducibility test
pub fn is_irreducible(g:&Vec<u64>, f:u64) -> bool {
    let mut t = deg(&g);
    let mut m = f1::deg(f);
    let mut rhs = vec![0,1]; // rhs = x
    let t_over_2 = (t as f64/2 as f64).floor() as u64;
    let mut j = 0;
    for i in 0..m { // square rhs m times
        rhs = square(&rhs, &g, f);
        reduce(&mut rhs, &g, f);
    }
    j += 1;
    while j <= t_over_2 {
        let mut rhs_plus_x = add(&mut rhs, &vec![0 as u64, 1 as u64]);
        let mut gcd = gcd(&g, &rhs_plus_x, f);
        if !is_one(&gcd) {
            return false;
        } 
        for i in 0..m { // square rhs m times
            rhs = square(&rhs, &g, f);
            reduce(&mut rhs, &g, f);
        }
        j += 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;
use rand_chacha::ChaCha20Rng;
use rand_chacha::ChaCha8Rng;

// 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 the invertible matrix S
pub fn dense_nonsingular_matrix(k:u64) -> Matrix {
    let mut m = zeros(k as usize, k as usize);
    for i in 0..k {
        for j in 0..k {
            if i == j {
                m[(i as usize, j as usize)] = 1 as f64;
            }
        }
    }
    let mut density = k;
    let mut min_density = ((k as f64*k as f64)/2 as f64).floor() as u64;
    let mut no_density_increase = 0;
    while density < min_density {
        for row in 0..k {
            let mut first_density = density;
            let mut random_row = rand::thread_rng().gen_range(0..k);
            if row != random_row {
                let mut possible = true;
                for i in 0..k {
                    if m[(row as usize, i as usize)] == 1.0 && m[(random_row as usize, i as usize)] == 1.0 {
                        possible = false;
                        break;
                    }
                }
                if possible {
                    for i in 0..k {
                        let mut before = m[(row as usize, i as usize)];
                        m[(row as usize, i as usize)] = m[(row as usize, i as usize)] + m[(random_row as usize, i as usize)]; 
                        let mut after = m[(row as usize, i as usize)];
                        if after > before {
                            density += 1;
                        }
                    }
                }
            }
            let mut last_density = density;
            if first_density == last_density {
                no_density_increase += 1;
            } else {
                no_density_increase = 0;
            }
        }
        if no_density_increase > k*3 {
            break;
        }
    }
    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, n:u64, f:u64) -> Matrix {
    let mut L = L.clone();
    let mut y = zeros(t as usize, n as usize);
    for c in 0..n {
        y[(0,c as usize)] = 1 as f64;
    }
    for r in 1..t {
        for c in 0..n {
            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);
}

// 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 irreduicible polynomial for GF(2^mt) is represented by g
//    - the irreducible 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::*;
use itertools::Itertools;

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

// Returns encrypted message
// For use in testing
pub fn encrypt_with_error_vector(message:&Matrix, generator:&Matrix, length:usize, t:u64, 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);
    // low order to high order
    let mut syndrome_inv_plus_x = f2::add(&syndrome_inv, &vec![0 as u64, 1 as u64]); 
    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() {
        let mut denom = vec![L[i], 1]; // low order to high order 
        if yP_inv[(0,i)] == 1 as f64 {
            let mut inverse = f2::inverse(&denom, &g, f);
            sum = f2::add(&sum, &inverse);
        } 
    }
    sum
}

// Returns vector of all Choose(width, group_size) combinations of possible error vector
// Use only for debugging
pub 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 error vector with t errors 
// Must contain t errors or else it crashes
pub fn choose_error_vector(length:u64, t:u64) -> Matrix {
    let mut error_vector = zeros(1,length as usize);
    let mut indices = Vec::new();
    let mut count = 0;
    while count < t {
        let mut index = rand::thread_rng().gen_range(0..length);
        if indices.contains(&index) {
            continue;
        } else {
            indices.push(index);
            error_vector[(0,index as usize)] = 1 as f64;
            count += 1;
        }
    }
    error_vector
}

system.rs

// system.rs
// This file contains:
//    - method to create a cryptosystem with the specified parameters
//    - structs that hold the private and public keys
//    - implementation of cipher block chaining
//    - method to encrypt and decrypt a string
// Explanation: 
//    - we use this file for creating systems
//    - we also use it for testing the system on text

use std::process;
use peroxide::fuga::*;
use crate::GF_2m_arithmetic  as f1;
use crate::GF_2mt_arithmetic as f2;
use crate::matrix;
use crate::encrypt_decrypt;
use crate::test;
use crate::util;

// Struct used to hold the public key
pub struct public_key {
    pub generator: Matrix,
    pub length: u64,
    pub dimension: u64, 
    pub t: u64
}

// Struct used to hold the private key
pub struct private_key {
    pub S: Matrix,
    pub G: Matrix, 
    pub P: Matrix, 
    pub L: Vec<u64>,
    pub g: Vec<u64>,
    pub f: u64,
    pub length: u64,
    pub dimension: u64
}

// Returns tuple containing public_key and private_key
pub fn create_system(n:u64, t:u64) -> (public_key, private_key) {
    let m = (n as f64).log2();
    if m.fract() != 0.0 {
        println!("Error: n must be a power of 2");
        process::exit(1);
    } 
    let mut m = m as u64;
    let mut length = 2_u64.pow(m as u32); // set n = 2^m always
    test::check_precondition(length, t, m);
    let mut f = f1::get_irreducible_polynomial(m);
    let mut g = f2::get_irreducible_polynomial(t, f);
    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 mut public_key = public_key {
        generator: generator,
        length: length,
        dimension: dimension,
        t: t
    };
    let mut private_key = private_key {
        S: S,
        G: G,
        P: P,
        L: L,
        g: g, 
        f: f,
        length: length,
        dimension: dimension
    };
    (public_key, private_key)
}

// Returns initialization vector for cipher block chaining
pub fn initialization_vector(public_key:&public_key) -> Matrix {
    let mut dimension = public_key.dimension;
    let mut iv = zeros(1 as usize, dimension as usize);
    for i in 0..dimension {
        let mut value:u64 = rand::thread_rng().gen_range(0..2);
        iv[(0,i as usize)] = value as f64;
    }
    iv
}

// Encrypts a String using cipher block chaining
pub fn encrypt(message:String, public_key:&public_key, iv:&Matrix) -> Vec<u64> {
    let mut plaintext = util::string_to_vec_of_u64s(message);
    let     dimension = public_key.dimension;
    let        length = public_key.length;
    // make sure vector of 1 and 0 is divisible by dimension
    while plaintext.len() % dimension as usize != 0 { 
        plaintext.push(0);
    }
    let mut prev_ciphertext = iv.clone(); 
    let mut      ciphertext = Vec::new(); // vector of 1 x length matrices
    // split plaintext into 1 x dimension matrices
    let              blocks = vec_into_matrix_blocks(&plaintext, dimension); 
    let mut count = 1;
    let mut blocks_len = blocks.len();
    for block in blocks { // iterate through matrices
        let x = xor_matrix(&block, &prev_ciphertext, dimension); // xor with previous ciphertext
        let mut curr_ciphertext = encrypt_decrypt::encrypt(&x, 
                                                            &public_key.generator,  
                                                            public_key.length as usize, 
                                                            public_key.t);
        prev_ciphertext = curr_ciphertext.clone();
        ciphertext.push(curr_ciphertext.clone()); // ciphertext is a vector of encrypted matrices
        count += 1;
    }
    let mut ciphertext = matrix_blocks_into_vec(&ciphertext, length);
    ciphertext
}

// Decrypts a String using cipher block chaining
pub fn decrypt(ciphertext:&Vec<u64>, private_key:&private_key, iv:&Matrix) -> String {
    let mut ciphertext = vec_into_matrix_blocks(&ciphertext, private_key.length); 
    let mut plaintext = Vec::new();
    let mut prev_ciphertext = iv.clone();
    let mut count = 1;
    let mut ciphertext_len = ciphertext.len();
    for block in ciphertext {
        let mut plaintext_block = encrypt_decrypt::decrypt(&block, 
                                                                &private_key.S, 
                                                                &private_key.G, 
                                                                &private_key.P, 
                                                                &private_key.g, 
                                                                private_key.f, 
                                                                &private_key.L, 
                                                                private_key.length, 
                                                                private_key.dimension);
        let mut plaintext_block = xor_matrix(&plaintext_block, &prev_ciphertext, private_key.dimension);
        plaintext.push(plaintext_block);
        prev_ciphertext = block.clone();
        count += 1; 
    }
    let mut plaintext = matrix_blocks_into_vec(&plaintext, private_key.dimension);
    while plaintext.len() % 8 as usize != 0 {
        plaintext.pop();
    }
    let mut plaintext = util::vec_of_u64s_to_string(&plaintext);
    plaintext
}

// Converts matrix to vector
pub fn matrix(a:&Vec<u64>, dimension:u64) -> Matrix {
    let mut m = zeros(1 as usize, dimension as usize);
    for i in 0..dimension {
        m[(0,i as usize)] = a[i as usize] as f64;
    }
    m
}

// Converts vector of 1 x width matrices to one long vector 
pub fn matrix_blocks_into_vec(a:&Vec<Matrix>, width:u64) -> Vec<u64> {
    let mut v = Vec::new();
    for i in 0..a.len() {
        for j in 0..width as usize {
            v.push(a[i][(0,j)] as u64);
        }   
    }
    v
}

// Takes in vector of length that is divisible by dimension
// Returns vector of matrices of length equal to dimension
// Precondition: a.len()%dimension == 0
pub fn vec_into_matrix_blocks(a:&Vec<u64>, dimension:u64) -> Vec<Matrix> {
    let mut blocks = Vec::new();
    let mut i = 0;
    while i < a.len() {
        let mut block = Vec::new();
        while block.len() < dimension as usize {
            block.push(a[i]);
            i += 1;
        }
        blocks.push(matrix(&block, dimension));
    }
    blocks
}

// Takes in two 1 x dimension matrices
// Returns 1 x dimension matrix that represents their xor
pub fn xor_matrix(a:&Matrix, b:&Matrix, dimension:u64) -> Matrix {
    let mut c = zeros(1 as usize, dimension as usize);
    for i in 0..dimension as usize {
        if a[(0,i)] == 1.0 && b[(0,i)] == 1.0 {
            c[(0,i)] = 0 as f64;
        } else if a[(0,i)] == 1.0 || b[(0,i)] == 1.0 {
            c[(0,i)] = 1.0 as f64;
        } else {
            c[(0,i)] = 0.0 as f64;
        }
    }
    c
}

test.rs

// test.rs
// This file contains:
//    - method to test one instance of the system
//    - method to unit test all messages and all error vectors
//    - method to run unit tests on specified input parameters
// Explanation: 
//    - this file is used to unit test the system

use rand::{thread_rng, Rng};
use peroxide::fuga::*;
use std::process;
use itertools::Itertools;
use rand::seq::SliceRandom;
use crate::GF_2m_arithmetic  as f1;
use crate::GF_2mt_arithmetic as f2;
use crate::matrix;
use crate::encrypt_decrypt;
use std::time::Instant;
use std::time::Duration;
use get_size::GetSize;
use crate::system;

pub fn run_some_tests() {

    let (mut n, mut t)  = (16, 2);
    let (correct, incorrect, total) = test(n, t, 1);

    let (mut n, mut t)  = (16, 3);
    let (correct, incorrect, total) = test(n, t, 2);

    let (mut n, mut t)  = (32, 2);
    let (correct, incorrect, total) = test(n, t, 3);

    let (mut n, mut t)  = (32, 3);
    let (correct, incorrect, total) = test(n, t, 4);

    let (mut n, mut t)  = (64, 2);
    let (correct, incorrect, total) = test(n, t, 5); 

}

//  Runs unit tests on a system with the specified parameters
pub fn test(n:u64, t:u64, test_number:u64) -> (u64, u64, u64) {
    let (public_key, private_key) = system::create_system(n, t);
    let (correct, incorrect, total) = run_unit_tests(&public_key, &private_key); 
    println!("test {} results:\n  {}/{} correct, {}/{} incorrect\n", test_number, correct, total, incorrect, total);
    return (correct, incorrect, total);
}

// Runs encryption and decryption on a single instance
pub fn run_single_test() {
    let (public_key, private_key) = system::create_system(8, 2);
    let mut length = public_key.length;
    let mut dimension = public_key.dimension;
    let dimension = public_key.dimension;
    let length = public_key.length;
    let t = public_key.t;
    let generator = public_key.generator.clone();
    let S = private_key.S.clone();
    let G = private_key.G.clone();
    let P = private_key.P.clone();
    let g = private_key.g.clone();
    let f = private_key.f;
    let L = private_key.L.clone();
    let mut message = zeros(1,2);
    message[(0,0)] = 0.0;
    message[(0,1)] = 1.0;
    let mut y = encrypt_decrypt::encrypt(&message, &generator, length as usize, t);
    let mut decrypted = encrypt_decrypt::decrypt(&y, &S, &G, &P, &g, f, &L, length, dimension);
    println!("decrypted = {:?}", decrypted);
}

// Runs encryption and decryption on multiple messages and multiple error vectors
pub fn run_unit_tests(public_key:&system::public_key, private_key:&system::private_key) -> (u64, u64, u64) {
    let dimension = public_key.dimension;
    let length = public_key.length;
    let t = public_key.t;
    let generator = public_key.generator.clone();
    let S = private_key.S.clone();
    let G = private_key.G.clone();
    let P = private_key.P.clone();
    let g = private_key.g.clone();
    let f = private_key.f;
    let L = private_key.L.clone();
    let messages = get_random_messages(dimension, 100); 
    let mut number_of_error_vectors = 10; 
    let mut error_vectors = Vec::new();
    for i in 0..number_of_error_vectors {
        error_vectors.push(encrypt_decrypt::choose_error_vector(length, t));
    }
    let mut correct_count = 0;
    let mut incorrect_count = 0;
    let mut count = 0;
    for message in &messages {
        for error_vector in &error_vectors {
            let mut y = encrypt_decrypt::encrypt_with_error_vector(&message, &generator, length as usize, t, error_vector);
            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;
            }
            count += 1;
        }
    } 
    (correct_count, incorrect_count, messages.len() as u64 * number_of_error_vectors as u64)
}

// Checks the precondition that 2^m > mt
pub fn check_precondition(length:u64, t:u64, m:u64) {
    if length <= t*m {
        println!("Error: length must be > t*m; length = {}, t*m = {}", length, t*m);
        process::exit(1);
    }
}

// check if two vectors are equal
pub 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
pub 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 random messages with length elements
pub fn get_random_messages(length:u64, number:u64) -> Vec<Matrix> {
    let mut messages = Vec::new();
    for i in 0..number {
        let mut message = zeros(1, length as usize);
        for j in 0..length {
            message[(0,j as usize)] = rand::thread_rng().gen_range(0..2) as f64;
        }
        messages.push(message);
    }
    messages
}

// returns vector of all possible binary vectors of width width
pub 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
}

util.rs

// util.rs
// This file contains:
//    - method to convert a string to a binary vector 
//    - method to convert a binary vector to a string
// Explanation: 
//    - these methods are used in cipher block chaining

// convert string to vector of 0 and 1 of type u64
pub fn string_to_vec_of_u64s(s:String) -> Vec<u64> {
    let mut string_of_1_and_0 = String::from("");
    for character in s.clone().into_bytes() { // character is number in ascii
        let mut character_as_binary_string = format!("{:b}", character);
        while character_as_binary_string.len() < 8 {
            let mut new_str = String::from("0");
            new_str.push_str(&character_as_binary_string);
            character_as_binary_string = new_str;
        }
        string_of_1_and_0.push_str(&character_as_binary_string);
    }
    let mut vec_of_u64s = Vec::new();
    for ch in string_of_1_and_0.chars() {
        if ch == '1' {
            vec_of_u64s.push(1);
        } else if ch == '0' {
            vec_of_u64s.push(0);
        }
    }
    if s.len()*8 != vec_of_u64s.len() {
        use std::process;
        println!("Error in util.rs: translating string to vector did not work");
        process::exit(1);
    }
    vec_of_u64s
}

// Converts a vector of u64 ones and zeros to a String of ones and zeros
pub fn vec_of_u64s_to_string(v:&Vec<u64>) -> String {
    let mut vec_of_chars = Vec::new();
    let mut binary_string = String::from("");
    for i in 0..v.len() {
        if v[i] == 0 {
            binary_string += "0";
        } else if v[i] == 1 {
            binary_string += "1";
        }
        if binary_string.len() == 8 {
            vec_of_chars.push(binary_string.clone());
            binary_string = String::from("");
        }
    }
    let mut characters = Vec::new();
    for i in 0..vec_of_chars.len() {
        let bin = vec_of_chars[i].clone();
        characters.push(isize::from_str_radix(&bin, 2).unwrap() as u8);
    }
    while characters[characters.len()-1] == 0 { // get rid of nulls 
        characters.pop();
    }
    let mut characters = &characters[..];
    let st = std::str::from_utf8(&characters).unwrap();
    st.to_string()
}

Cargo.toml

[package]
name = "mceliece"
version = "0.1.0"
edition = "2021"

[dependencies]
rand = "0.8.5"
peroxide = "0.30"
itertools = "0.10.3"

Sources:

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

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