From 5d7ae639739e1b5389bc94c227aeb6661cb8eab1 Mon Sep 17 00:00:00 2001 From: JSDurand Date: Wed, 6 Dec 2023 03:44:28 +0800 Subject: initial commit A manual implementation of the Berlekamp algorithm for finding irreducible factors of polynomials over finite fields. --- src/lib.rs | 842 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/main.rs | 317 +++++++++++++++++++++++ 2 files changed, 1159 insertions(+) create mode 100644 src/lib.rs create mode 100644 src/main.rs (limited to 'src') diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..db3186a --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,842 @@ +use num::{ + bigint::{BigInt, ToBigInt}, + rational::BigRational, + traits::Euclid, + Signed, +}; + +/// A simple conversion function from anything that can be converted +/// to a big integer to a big integer. +fn conversion(n: impl ToBigInt) -> BigInt { + n.to_bigint().unwrap() +} + +/// A simple conversion function from anything that can be converted +/// to a big integer to a big rational number. +#[allow(unused)] +fn ratconversion(n: impl ToBigInt) -> BigRational { + BigRational::from_integer(n.to_bigint().unwrap()) +} + +/// Return the inverse of `a` modulo `p`. Here `p` is assumed to be +/// relatively prime with `a`. +/// +/// This is based on Euclidean algorithm. In fact, it returns a +/// coefficient `x` such that ax + py = gcd(a, p), for some y. +pub fn modulo_inverse(a: BigInt, p: BigInt) -> BigInt { + // perform Euclidean algorithm of a by p + let mut quotients: Vec = Vec::new(); + let mut lhss: Vec = Vec::new(); + + let mut divider: BigInt = p.clone(); + let mut divided: BigInt = a.clone(); + + if divider <= conversion(1) { + panic!("1 or 0 is not a prime"); + } + + while divider > conversion(1) { + // The mighty compiler should simplify this to a single + // instruction on a modern "enough" CPU. + let div: BigInt = divided.div_euclid(÷r); + let rem: BigInt = divided.rem_euclid(÷r); + + // for debugging + + // #[cfg(debug_assertions)] + // println!("{divided} = {div} * {divider} + {rem}"); + + quotients.push(div); + lhss.push(divided); + + divided = divider; + divider = rem; + } + + if divider != conversion(1) { + dbg!(÷r, &a, &p); + panic!("they are not relatively prime"); + } + + assert!(!quotients.is_empty()); + assert!(!lhss.is_empty()); + + let mut rhs = quotients.pop().unwrap(); + rhs *= conversion(-1); + + let mut lhs = conversion(1); + + // #[cfg(debug_assertions)] + // println!("lhs = {lhs}, rhs = {rhs}"); + + while let Some(q) = quotients.pop() { + // already: a * lhs + b * rhs = 1 + + // now: b = x - q * a + + // so 1 = a * lhs + rhs (x - q * a) = rhs * x + (lhs - rhs * q) * a + + let temp = rhs.clone(); + + rhs = lhs - rhs * q; + + lhs = temp; + + // #[cfg(debug_assertions)] + // println!("lhs = {lhs}, rhs = {rhs}"); + } + + // finally lhs * a + rhs * p = 1 + + lhs.rem_euclid(&p) +} + +/// Return the remainder of the Euclidean division of a by b, over the +/// finite field of p elements, where p is assumed to be a natural +/// prime number. +pub fn poly_rem(a: &[BigInt], b: &[BigInt], p: BigInt) -> Vec { + if a.len() < b.len() { + return a.to_vec(); + } + + if b.is_empty() { + panic!("cannot divide by 0"); + } + + let b_leading_coefficient = b[b.len() - 1].clone(); + + let b_leading_inverse = modulo_inverse(b_leading_coefficient, p.clone()); + + let mut candidate: Vec = a.iter().map(|c| c.rem_euclid(&p)).collect(); + + let mut diff: usize = a.len() - b.len(); + + loop { + let cand_leading = candidate[b.len() + diff - 1].clone(); + + let leading_coefficient = b_leading_inverse.clone() * cand_leading; + + if leading_coefficient != conversion(0) { + for (i, cand_ref) in candidate.iter_mut().skip(diff).take(b.len()).enumerate() { + // println!( + // "bi = {}, leading_coefficient = {leading_coefficient}, cand = {cand_ref}", + // b[i] + // ); + *cand_ref -= b[i].clone() * leading_coefficient.clone(); + // println!("cand = {cand_ref}"); + *cand_ref = cand_ref.rem_euclid(&p); + } + } + + if diff > 0 { + diff -= 1; + } else { + break; + } + } + + let mut max_non_zero = 0usize; + let mut non_zero = false; + + // dbg!(); + // print_poly(&candidate); + + for (index, cand) in candidate.iter().enumerate() { + if *cand != conversion(0) { + max_non_zero = index; + non_zero = true; + } + } + + if !non_zero { + return Vec::new(); + } + + candidate.into_iter().take(max_non_zero + 1).collect() +} + +/// Return both the quotient and the remainder of the Euclidean +/// division of a by b, over the finite field of p elements, where p +/// is assumed to be a natural prime number. +pub fn poly_division(a: &[BigInt], b: &[BigInt], p: BigInt) -> (Vec, Vec) { + if a.len() < b.len() { + return (Vec::new(), a.to_vec()); + } + + if b.is_empty() { + panic!("cannot divide by 0"); + } + + let b_leading_coefficient = b[b.len() - 1].clone(); + + let b_leading_inverse = modulo_inverse(b_leading_coefficient, p.clone()); + + let mut candidate: Vec = a.iter().map(|c| c.rem_euclid(&p)).collect(); + + let mut diff: usize = a.len() - b.len(); + + let mut quotient: Vec = Vec::with_capacity(diff); + + loop { + let cand_leading = candidate[b.len() + diff - 1].clone(); + + let leading_coefficient = b_leading_inverse.clone() * cand_leading; + + quotient.push(leading_coefficient.clone()); + + if leading_coefficient != conversion(0) { + for (i, cand_ref) in candidate.iter_mut().skip(diff).take(b.len()).enumerate() { + // println!( + // "bi = {}, leading_coefficient = {leading_coefficient}, cand = {cand_ref}", + // b[i] + // ); + *cand_ref -= b[i].clone() * leading_coefficient.clone(); + // println!("cand = {cand_ref}"); + *cand_ref = cand_ref.rem_euclid(&p); + } + } + + if diff > 0 { + diff -= 1; + } else { + break; + } + } + + quotient = quotient.into_iter().rev().collect(); + + let mut max_non_zero = 0usize; + let mut non_zero = false; + + // dbg!(); + // print_poly(&candidate); + + for (index, cand) in candidate.iter().enumerate() { + if *cand != conversion(0) { + max_non_zero = index; + non_zero = true; + } + } + + if !non_zero { + (quotient, Vec::new()) + } else { + ( + quotient, + candidate.into_iter().take(max_non_zero + 1).collect(), + ) + } +} + +pub fn poly_gcd(a: &[BigInt], b: &[BigInt], p: BigInt) -> Vec { + let mut a = a.to_vec(); + let mut b = b.to_vec(); + + let mut rem = poly_rem(&a, &b, p.clone()); + + while !rem.is_empty() { + a = b; + b = rem; + + rem = poly_rem(&a, &b, p.clone()); + } + + b +} + +/// Return a vector of left-null-vectors over the finite field of p +/// elements, where p is assumed to be a natural prime number. +/// +/// It is assumed that T is the square of X. +pub fn find_left_null_space(x: usize, t: usize, matrix: &[BigInt], p: usize) -> Vec> { + assert_eq!(t, x * x); + assert_eq!(matrix.len(), t); + + let mut matrix = matrix.to_vec(); + let mut result: Vec> = Vec::new(); + + let mut used_pivots_rows: Vec> = + std::iter::repeat(None).take(x).collect::>(); + + // also perform row-operations on the identity matrix + + let mut identity: Vec = std::iter::repeat(conversion(0)).take(t).collect::>(); + + for i in 0..x { + *identity.get_mut(i * x + i).unwrap() = conversion(1); + } + + for i in 0..x { + for j in 0..x { + if used_pivots_rows.get(j).unwrap().is_some() { + continue; + } + + let mat_ij = matrix.get(i * x + j).unwrap().clone(); + + if mat_ij == conversion(0) { + continue; + } + + *used_pivots_rows.get_mut(j).unwrap() = Some(i); + + let mat_ij_inverse = modulo_inverse(mat_ij.clone(), conversion(p)); + + for k in 0..x { + if k == i || matrix.get(k * x + j).unwrap() == &conversion(0) { + continue; + } + + let mat_kj = matrix.get(k * x + j).unwrap().clone(); + + for ell in 0..x { + let ident_iell = identity.get(i * x + ell).unwrap().clone(); + let mat_iell = matrix.get(i * x + ell).unwrap().clone(); + + *identity.get_mut(k * x + ell).unwrap() -= + mat_ij_inverse.clone() * mat_kj.clone() * ident_iell; + + *identity.get_mut(k * x + ell).unwrap() = identity + .get(k * x + ell) + .unwrap() + .rem_euclid(&conversion(p)); + + *matrix.get_mut(k * x + ell).unwrap() -= + mat_ij_inverse.clone() * mat_kj.clone() * mat_iell; + + *matrix.get_mut(k * x + ell).unwrap() = + matrix.get(k * x + ell).unwrap().rem_euclid(&conversion(p)); + } + } + + break; + } + } + + // println!("pivot: {used_pivots_rows:?}"); + + let mut used_rows: Vec = std::iter::repeat(false) + .take(x) + .collect::>() + .as_slice() + .try_into() + .unwrap(); + + for i in used_pivots_rows.iter().filter_map(|x| *x) { + *used_rows.get_mut(i).unwrap() = true; + } + + for (index, used) in used_rows.iter().copied().enumerate() { + if !used { + let mut newrow: Vec = Vec::with_capacity(x); + + for j in 0..x { + newrow.push(identity.get(index * x + j).unwrap().clone()); + } + + result.push(newrow); + } + } + + result +} + +fn derivative(poly: &[BigInt], p: usize) -> Vec { + let big_p = conversion(p); + + let untrimmed = poly + .iter() + .enumerate() + .skip(1) + .map(|(n, a)| (conversion(n) * a.clone()).rem_euclid(&big_p)); + + let mut max_non_zero = 0usize; + let mut non_zero_p = false; + + for (i, c) in untrimmed.clone().enumerate() { + if c != conversion(0) { + non_zero_p = true; + max_non_zero = i; + } + } + + if !non_zero_p { + Vec::new() + } else { + untrimmed.take(max_non_zero + 1).collect() + } +} + +/// Return true if and only if the polynomial has no repeated factors +/// modulo the prime p. +pub fn is_square_free(poly: &[BigInt], p: usize) -> Result<(), Vec> { + let derivative = derivative(poly, p); + + if derivative.is_empty() { + return Err(poly.to_vec()); + } + + let gcd = poly_gcd(poly, &derivative, conversion(p)); + + if gcd.len() <= 1 { + Ok(()) + } else { + Err(gcd) + } +} + +#[derive(Debug, Clone, Default)] +pub enum BeRe { + Prod(Vec, Vec), + #[default] + Irreducible, + NonSquareFree, +} + +impl std::fmt::Display for BeRe { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + BeRe::Prod(a, b) => write!(f, "product of {} and {}", Poly(a), Poly(b)), + BeRe::Irreducible => write!(f, "irreducible"), + BeRe::NonSquareFree => write!(f, "not square-free"), + } + } +} + +/// Return an irreducible polynomial that divides the polynomial, over +/// the finite field specified by p, which is assumed to be a positive +/// prime number. +/// +/// Also the polynomial should be square-free. +pub fn berlekamp(poly: &[BigInt], p: usize) -> BeRe { + // reject non-square-free polynomials + if is_square_free(poly, p).is_err() { + return BeRe::NonSquareFree; + } + + // compute the basis matrix + let n = poly.len(); + + if n == 0 { + return BeRe::Prod(Vec::default(), Vec::default()); + } + + // print!("poly: "); + // print_poly(poly); + // println!(); + + let mut xqi_mat: Vec<_> = std::iter::repeat_with(|| conversion(0)) + .take(n * n) + .collect(); + + let mut xq = vec![conversion(1)]; + + let rem = poly_rem(&xq, poly, conversion(p)); + + for (co, result_mut) in rem.iter().zip(xqi_mat.iter_mut()) { + *result_mut = co.clone(); + } + + xq = rem; + + for i in 1..=(n - 1) { + let mut temp: Vec<_> = std::iter::repeat_with(|| conversion(0)).take(p).collect(); + + temp.append(&mut xq); + + xq = temp; + + let rem = poly_rem(&xq, poly, conversion(p)); + + // print!("i = {i}, rem = "); + // print_poly(&rem); + + for (co, result_mut) in rem.iter().zip(xqi_mat.iter_mut().skip(n * i)) { + *result_mut = co.clone(); + } + } + + // now subtract the identity matrix + + for i in 0..n { + *xqi_mat.get_mut(i * n + i).unwrap() -= conversion(1); + } + + // println!("haha"); + // print_matrix::<3, 3, 9>(xqi_mat.as_slice().try_into().unwrap()); + + // now find a basis for the left-null-space + let basis = find_left_null_space(n, n * n, &xqi_mat, p); + + assert!(!basis.is_empty(), "the kernel should not be zero"); + + // for (i, v) in basis.clone().into_iter().enumerate() { + // print!("{i}: "); + // print_matrix::<1, 3, 3>(v.as_slice().try_into().unwrap()); + // } + + if basis.len() == 1 { + return BeRe::default(); + } + + // now find the corresponding polynomials + + // The first element correponds always to the constant polynomial + // and serves no purpose to us, so we omit it here. + + let _irr_factors_num = basis.len(); + + // dbg!(irr_factors_num); + + for v in basis.clone().into_iter().skip(1) { + let mut max_non_zero = 0usize; + let mut non_zero = false; + + for (index, cand) in v.iter().enumerate() { + if *cand != conversion(0) { + max_non_zero = index; + non_zero = true; + } + } + + if !non_zero { + panic!("this should not happen"); + } + + let mut v: Vec<_> = v.iter().take(max_non_zero + 1).cloned().collect(); + + for _ in 0..p { + *v.get_mut(0).unwrap() -= 1; + + // print_poly(&v); + + // Now take the gcd of the polynomial with v + + let gcd = poly_gcd(poly, &v, conversion(p)); + + if (2..poly.len()).contains(&gcd.len()) { + let (quotient, _) = poly_division(poly, &gcd, conversion(p)); + + return BeRe::Prod(gcd, quotient); + } + } + } + + println!("weird poly = {}", Poly(poly)); + + for (n, v) in basis.clone().into_iter().enumerate() { + println!("{n}-th basis: {}", Poly(&v)); + } + + panic!("this should not happen"); +} + +/// Return all irreducible factors of the polynomial modulo the prime +/// p. +/// +/// If the polynomial has repeated factors, return an error. +#[allow(clippy::result_unit_err)] +pub fn factors(poly: &[BigInt], p: usize) -> Result>, Vec> { + if let Err(gcd) = is_square_free(poly, p) { + return Err(gcd); + } + + let mut result = Vec::new(); + let mut candidates: Vec> = vec![poly.to_vec()]; + + while let Some(poly) = candidates.pop() { + let bere = berlekamp(&poly, p); + + match bere { + BeRe::Prod(x, y) => { + candidates.push(x); + candidates.push(y); + } + BeRe::Irreducible => result.push(poly), + BeRe::NonSquareFree => return Err(Vec::new()), + } + } + + Ok(result) +} + +/// Return a sub-slice that trims the leading zeroes. +pub fn trim(p: &[BigInt]) -> &[BigInt] { + let mut max_non_zero = 0usize; + let mut non_zero_p = false; + + for (n, c) in p.iter().enumerate() { + if c != &conversion(0) { + non_zero_p = true; + max_non_zero = n; + } + } + + if !non_zero_p { + &p[0..0] + } else { + &p[0..=max_non_zero] + } +} + +/// Return the product of two polynomials. +pub fn poly_prod(p1: &[BigInt], p2: &[BigInt]) -> Vec { + let p1 = trim(p1); + + let p2 = trim(p2); + + // degrees sum up, so the lengths have to subtract one + let mut result: Vec = std::iter::repeat_with(|| conversion(0)) + .take(p1.len() + p2.len() - 1) + .collect(); + + let mut degree_map: std::collections::HashMap = + std::collections::HashMap::with_capacity(p1.len() + p2.len()); + + for (n1, c1) in p1.iter().enumerate() { + for (n2, c2) in p2.iter().enumerate() { + let new_degree = n1 + n2; + + if let Some(orig) = degree_map.get(&new_degree) { + degree_map.insert(new_degree, orig + c1 * c2); + } else { + degree_map.insert(new_degree, c1 * c2); + } + } + } + + for (n, c) in degree_map { + *result.get_mut(n).unwrap() = c; + } + + result +} + +pub fn power(p: &[BigInt], n: usize) -> Vec { + if n == 0 { + return Vec::new(); + } + + let mut result = p.to_vec(); + + for _ in 0..(n - 1) { + result = poly_prod(&result, p); + } + + result +} + +/// Return c * p + q. +pub fn poly_sum(p: &[BigInt], q: &[BigInt], cp: BigInt) -> Vec { + let p = trim(p); + let q = trim(q); + + if p.is_empty() { + return q.to_vec(); + } + + if q.is_empty() { + return p.iter().map(|c| c * &cp).collect(); + } + + if p.len() >= q.len() { + let mut result: Vec<_> = p.iter().map(|c| c * &cp).collect(); + + for (n, c) in q.iter().enumerate() { + *result.get_mut(n).unwrap() += c.clone(); + } + + result + } else { + let mut result = q.to_vec(); + + for (n, c) in p.iter().enumerate() { + *result.get_mut(n).unwrap() += c.clone() * cp.clone(); + } + + result + } +} + +/// Return the composition of p and q, i.e. p ( q ( x ) ) as a +/// polynomial in x. +pub fn composition(p: &[BigInt], q: &[BigInt]) -> Vec { + let p = trim(p); + let q = trim(q); + + let mut power: Vec<_> = vec![conversion(1)]; + + let mut result: Vec = Vec::new(); + + for c in p.iter() { + result = poly_sum(&power, &result, c.clone()); + + power = poly_prod(&power, q); + } + + result +} + +/// Print an array as a square matrix nicely. +/// +/// This has the same requirement as the function `power`. Similarly, +/// one can use the macro `print_square_matrix` to correctly fill in +/// the dimensions automatically. +#[allow(unused)] +fn print_matrix(matrix: &[BigInt; T]) { + if X * Y != T { + panic!("dimensions do not match: {X} * {Y} is not {T}"); + } + + println!("["); + + let mut max_lens: [usize; Y] = [0; Y]; + + for j in 0..Y { + for i in 0..X { + let entry_str = format!("{}", matrix.get(Y * i + j).unwrap()); + + *max_lens.get_mut(j).unwrap() = + std::cmp::max(*max_lens.get(j).unwrap(), entry_str.len()); + } + } + + for i in 0..X { + print!(" "); + + for j in 0..Y { + let entry_str = format!("{}", matrix.get(Y * i + j).unwrap()); + + let fill_in_space: String = std::iter::repeat(32 as char) + .take(*max_lens.get(j).unwrap() - entry_str.len()) + .collect(); + + print!("{}{} ", fill_in_space, entry_str); + } + + println!(); + } + println!("]"); +} + +pub struct Poly<'a>(pub &'a [BigInt]); + +impl<'a> std::fmt::Display for Poly<'a> { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let poly = self.0; + + struct Monomial(usize); + + impl std::fmt::Display for Monomial { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self.0 { + 0 => Ok(()), + 1 => { + write!(f, "x") + } + _ => { + write!(f, "x^{}", self.0) + } + } + } + } + + struct Coefficient(BigInt); + + impl std::fmt::Display for Coefficient { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match &self.0 { + n if n == &conversion(0) || n == &conversion(1) || n == &conversion(-1) => { + Ok(()) + } + _ => { + write!(f, "{}", self.0.abs()) + } + } + } + } + + for (n, coeff) in poly.iter().enumerate().rev() { + if coeff == &conversion(0) { + continue; + } + + write!( + f, + "{}{}{}", + if n + 1 >= poly.len() { + "" + } else if coeff >= &conversion(0) { + " + " + } else { + " - " + }, + Coefficient(coeff.clone()), + Monomial(n), + )?; + + if coeff.abs() == conversion(1) && n == 0 { + write!(f, "1")?; + } + } + + Ok(()) + } +} + +#[allow(unused)] +fn print_poly(poly: &[BigInt]) { + struct Monomial(usize); + + impl std::fmt::Display for Monomial { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self.0 { + 0 => Ok(()), + 1 => { + write!(f, "x") + } + _ => { + write!(f, "x^{}", self.0) + } + } + } + } + + struct Coefficient(BigInt); + + impl std::fmt::Display for Coefficient { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match &self.0 { + n if n == &conversion(0) || n == &conversion(1) || n == &conversion(-1) => Ok(()), + _ => { + write!(f, "{}", self.0.abs()) + } + } + } + } + + for (n, coeff) in poly.iter().enumerate().rev() { + if coeff == &conversion(0) { + continue; + } + + print!( + "{}{}{}", + if n + 1 >= poly.len() { + "" + } else if coeff >= &conversion(0) { + " + " + } else { + " - " + }, + Coefficient(coeff.clone()), + Monomial(n), + ); + + if coeff.abs() == conversion(1) && n == 0 { + print!("1"); + } + } + + println!(); +} diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..b6ce57c --- /dev/null +++ b/src/main.rs @@ -0,0 +1,317 @@ +use num::{ + bigint::{BigInt, ToBigInt}, + Signed, +}; + +/// A simple conversion function from anything that can be converted +/// to a big integer to a big rational number. +fn conversion(n: impl ToBigInt) -> BigInt { + n.to_bigint().unwrap() +} + +fn print_poly(poly: &[BigInt]) { + struct Monomial(usize); + + impl std::fmt::Display for Monomial { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self.0 { + 0 => Ok(()), + 1 => { + write!(f, "x") + } + _ => { + write!(f, "x^{}", self.0) + } + } + } + } + + struct Coefficient(BigInt); + + impl std::fmt::Display for Coefficient { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match &self.0 { + n if n == &conversion(0) || n == &conversion(1) || n == &conversion(-1) => Ok(()), + _ => { + write!(f, "{}", self.0.abs()) + } + } + } + } + + for (n, coeff) in poly.iter().enumerate().rev() { + if coeff == &conversion(0) { + continue; + } + + print!( + "{}{}{}", + if n + 1 == poly.len() { + "" + } else if coeff >= &conversion(0) { + " + " + } else { + " - " + }, + Coefficient(coeff.clone()), + Monomial(n), + ); + + if coeff.abs() == conversion(1) && n == 0 { + print!("1"); + } + } + + println!(); +} + +#[allow(unused)] +macro_rules! print_square_matrix { + ($s:literal, $g: expr) => { + print_matrix::<$s, $s, { $s * $s }>($g); + }; +} + +fn main() { + let poly = read_poly(&std::env::args().nth(1).expect("Enter a polynomial")); + let p = std::env::args() + .nth(2) + .expect("Enter a prime") + .parse::() + .unwrap(); + + println!("poly is:"); + print_poly(&poly); + + println!("prime is {p}"); + + // let poly1 = [12, 12, 0, 1].map(conversion); + // let poly2 = [-4, 1].map(conversion); + + // let composition = berlekamp::composition(&poly1, &poly2); + + use berlekamp::Poly; + + // println!( + // "the composition of {} and {} is {}", + // Poly(&poly1), + // Poly(&poly2), + // Poly(&composition) + // ); + + // println!("poly2 is {factor}"); + + // let poly = [-25, 0, 15, 0, -3, 0, 1].map(conversion); + + // let p = 5; + + let factors = berlekamp::factors(&poly, p); + + match factors { + Ok(factors) => { + println!("the irreducible factors are:"); + + for factor in factors { + print_poly(&factor); + } + } + Err(gcd) => { + println!( + "the gcd of the polynomial with its derivative is {}", + Poly(&gcd) + ); + + if gcd.len() == poly.len() { + // the derivative is zero + println!("the polynomial is a {p}-th power of some polynomial"); + + return; + } + + let (quotient, _) = berlekamp::poly_division(&poly, &gcd, conversion(p)); + + println!( + "the quotient of the polynomial by the above gcd is {}", + Poly("ient) + ); + + let factors = + berlekamp::factors("ient, p).expect("this quotient should be square-free"); + + println!("the irreducible factors of the quotient, and hence of the polynomial, are:"); + + for factor in factors { + print_poly(&factor); + } + } + } +} + +#[allow(unused)] +fn read_poly(s: &str) -> Vec { + let mut degree: usize; + let mut coefficient = 1isize; + let mut pending_coefficient = false; + + let mut iter = s.chars().peekable(); + + let mut degrees_and_cos: Vec<(usize, BigInt)> = Vec::new(); + + while let Some(c) = iter.next() { + match c { + '\n' | '\t' | ' ' => {} + '+' => { + if pending_coefficient { + degrees_and_cos.push((0usize, conversion(coefficient))); + } + + pending_coefficient = false; + coefficient = 1isize; + } + 'x' => { + pending_coefficient = false; + + if !matches!(iter.peek(), Some(c) if *c == '^') { + degrees_and_cos.push((1usize, conversion(coefficient))); + continue; + } + + let _ = iter.next(); + + degree = 0usize; + + while matches!(iter.peek(), Some(c) if c.is_digit(10)) { + degree *= 10; + degree += iter.next().unwrap() as usize - '0' as usize; + } + + degrees_and_cos.push((degree, conversion(coefficient))); + } + c if c.is_digit(10) || c == '-' => { + pending_coefficient = true; + + let mut negative = false; + + coefficient = if c.is_digit(10) { + c as isize - '0' as isize + } else { + negative = true; + -1isize + }; + + let mut first_negative = negative; + + while matches!(iter.peek(), Some(c) if c.is_digit(10)) { + if first_negative { + coefficient = '0' as isize - iter.next().unwrap() as isize; + + first_negative = false; + + continue; + } + + let factor = if negative { + '0' as isize - iter.next().unwrap() as isize + } else { + iter.next().unwrap() as isize - '0' as isize + }; + + coefficient *= 10isize; + coefficient += factor; + } + } + _ => { + panic!("invalid: {c}"); + } + } + } + + if pending_coefficient { + degrees_and_cos.push((0usize, conversion(coefficient))); + } + + let mut degree_co_map: std::collections::HashMap = Default::default(); + + for (degree, coefficient) in degrees_and_cos { + if let Some(orig) = degree_co_map.get(°ree) { + degree_co_map.insert(degree, orig.clone() + coefficient); + } else { + degree_co_map.insert(degree, coefficient); + } + } + + // degree_co_map.extend(degrees_and_cos); + + let mut max_degree = 0usize; + let mut non_zero_p = false; + + for (d, c) in degree_co_map.iter() { + if c == &conversion(0) { + continue; + } + + non_zero_p = true; + + max_degree = std::cmp::max(*d, max_degree); + } + + if !non_zero_p { + Vec::new() + } else { + let mut result: Vec<_> = std::iter::repeat_with(|| conversion(0)) + .take(max_degree + 1) + .collect(); + + for (d, c) in degree_co_map { + if c == conversion(0) { + continue; + } + + *result.get_mut(d).unwrap() = c; + } + + result + } +} + +/// Print an array as a square matrix nicely. +/// +/// This has the same requirement as the function `power`. Similarly, +/// one can use the macro `print_square_matrix` to correctly fill in +/// the dimensions automatically. +#[allow(unused)] +fn print_matrix(matrix: &[BigInt; T]) { + if X * Y != T { + panic!("dimensions do not match: {X} * {Y} is not {T}"); + } + + println!("["); + + let mut max_lens: [usize; Y] = [0; Y]; + + for j in 0..Y { + for i in 0..X { + let entry_str = format!("{}", matrix.get(Y * i + j).unwrap()); + + *max_lens.get_mut(j).unwrap() = + std::cmp::max(*max_lens.get(j).unwrap(), entry_str.len()); + } + } + + for i in 0..X { + print!(" "); + + for j in 0..Y { + let entry_str = format!("{}", matrix.get(Y * i + j).unwrap()); + + let fill_in_space: String = std::iter::repeat(32 as char) + .take(*max_lens.get(j).unwrap() - entry_str.len()) + .collect(); + + print!("{}{} ", fill_in_space, entry_str); + } + + println!(); + } + println!("]"); +} -- cgit v1.2.3-18-g5258