use num::{ bigint::{BigInt, ToBigInt}, Signed, }; use berlekamp::Poly; /// 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() } #[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!(); } #[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(); if std::env::args().len() >= 4 { // Search for primes modulo which the polynomial is irreducible instead. let primes = generate_primes(p); println!("Searching {} primes up to {p}...", primes.len()); for p in primes { let factors = berlekamp::factors(&poly, p); let total = factors.iter().fold(0usize, |sum, (deg, _)| sum + deg); if total == 1 { println!("The polynomial {} is irreducible mod {p}", Poly(&poly)); return; } } println!("The polynomial is reducible modulo primes up to {p}."); return; } // println!("poly is:"); // print_poly(&poly); // println!("prime is {p}"); // let poly1 = [-20, -21, 0, 1].map(conversion); // let poly2 = [1, 3].map(conversion); // let composition = berlekamp::composition(&poly1, &poly2); // 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); print!("{} ≡ ", Poly(&poly)); let mut total = 0usize; for (index, (mul, f)) in factors.into_iter().enumerate() { total += mul; match mul { 0 => { dbg!("zero multiplicity?"); } 1 => { print!("{}({})", if index > 0 { " * " } else { "" }, Poly(&f)); } _ => { print!( "{}({})^{}", if index > 0 { " * " } else { "" }, Poly(&f), mul ); } } } println!(" (mod {p})"); if total == 1 { println!("the polynomial is irreducible"); } } #[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!("]"); } #[allow(unused)] /// Return a vector of primes <= BOUND. /// /// This is naïve and slow, so BOUND should not be too large. fn generate_primes(bound: usize) -> Vec { if bound <= 1 { return Vec::new(); } let sqrt = (bound as f64).sqrt() as usize; let mut result: Vec = Vec::with_capacity(bound); let mut records: Vec = std::iter::repeat(true).take(bound).collect(); records[0] = false; records[1] = false; for i in 2..=(sqrt + 1) { if !matches!(records.get(i), Some(true)) { continue; } let mut multiple = i + i; while let Some(rec) = records.get_mut(multiple) { *rec = false; multiple += i; } } result.extend( records .iter() .enumerate() .filter_map(|(n, rec)| rec.then_some(n)), ); result }