diff options
| -rw-r--r-- | src/lib.rs | 103 | ||||
| -rw-r--r-- | src/main.rs | 134 |
2 files changed, 179 insertions, 58 deletions
@@ -182,7 +182,7 @@ pub fn poly_division(a: &[BigInt], b: &[BigInt], p: BigInt) -> (Vec<BigInt>, Vec let leading_coefficient = b_leading_inverse.clone() * cand_leading; - quotient.push(leading_coefficient.clone()); + quotient.push(leading_coefficient.clone().rem_euclid(&p)); if leading_coefficient != conversion(0) { for (i, cand_ref) in candidate.iter_mut().skip(diff).take(b.len()).enumerate() { @@ -367,8 +367,9 @@ fn derivative(poly: &[BigInt], p: usize) -> Vec<BigInt> { } } -/// Return true if and only if the polynomial has no repeated factors -/// modulo the prime p. +/// Return Ok() if and only if the polynomial has no repeated factors +/// modulo the prime p, otherwise return an error containing the gcd +/// of the polynomial with its derivative. pub fn is_square_free(poly: &[BigInt], p: usize) -> Result<(), Vec<BigInt>> { let derivative = derivative(poly, p); @@ -535,29 +536,91 @@ pub fn berlekamp(poly: &[BigInt], p: usize) -> BeRe { /// p. /// /// If the polynomial has repeated factors, return an error. -#[allow(clippy::result_unit_err)] -pub fn factors(poly: &[BigInt], p: usize) -> Result<Vec<Vec<BigInt>>, Vec<BigInt>> { - if let Err(gcd) = is_square_free(poly, p) { - return Err(gcd); - } +pub fn factors(poly: &[BigInt], p: usize) -> Vec<(usize, Vec<BigInt>)> { + let poly = trim(poly); + + let mut result: Vec<(usize, Vec<BigInt>)> = Vec::new(); + let mut candidates: Vec<(usize, Vec<BigInt>)> = vec![(1, poly.to_vec())]; + + 'candidate_loop: while let Some((multiplicity, poly)) = candidates.pop() { + if let Err(gcd) = is_square_free(&poly, p) { + if gcd.len() == poly.len() { + // the derivative is zero + let pth_root = fake_p_th_root(&poly, p); + + if pth_root.is_empty() { + return Vec::new(); + } + + candidates.push((multiplicity * p, pth_root)); + + continue 'candidate_loop; + } + + let (quotient, _) = poly_division(&poly, &gcd, conversion(p)); + + candidates.push((multiplicity, quotient)); + candidates.push((multiplicity, gcd)); - let mut result = Vec::new(); - let mut candidates: Vec<Vec<BigInt>> = vec![poly.to_vec()]; + continue 'candidate_loop; + } - while let Some(poly) = candidates.pop() { let bere = berlekamp(&poly, p); match bere { BeRe::Prod(x, y) => { - candidates.push(x); - candidates.push(y); + candidates.push((multiplicity, x)); + candidates.push((multiplicity, y)); + } + BeRe::Irreducible => { + result.push((multiplicity, poly)); } - BeRe::Irreducible => result.push(poly), - BeRe::NonSquareFree => return Err(Vec::new()), + BeRe::NonSquareFree => panic!("this should be impossible"), } } - Ok(result) + result +} + +/// Return the p-th root of the polynomial, if it is a p-th power +/// indeed. +/// +/// Precisely speaking, a monomial cx^{pn+k} in the polynomial +/// contributes to the term cx^n in the new polynomial, where k is any +/// integer such that 0 <= k < p. +fn fake_p_th_root(poly: &[BigInt], p: usize) -> Vec<BigInt> { + let big_p = conversion(p); + let big_zero = conversion(0); + + let iter = poly.iter().enumerate().filter_map(|(n, c)| { + let rem = c.rem_euclid(&big_p); + + (&rem != &big_zero).then(|| (n.div_euclid(p), rem)) + }); + + let mut degrees_map: std::collections::HashMap<usize, BigInt> = Default::default(); + + degrees_map.extend(iter); + + if degrees_map.is_empty() { + return Vec::new(); + } + + let mut max = 0usize; + + for (n, _) in degrees_map.iter() { + max = std::cmp::max(max, *n); + } + + let mut result: Vec<BigInt> = std::iter::repeat_with(|| conversion(0)) + .take(max + 1) + .collect(); + + for (n, c) in degrees_map { + *result.get_mut(n).unwrap() = c; + } + + result } /// Return a sub-slice that trims the leading zeroes. @@ -764,12 +827,12 @@ impl<'a> std::fmt::Display for Poly<'a> { write!( f, "{}{}{}", - if n + 1 >= poly.len() { + if coeff < &conversion(0) { + " - " + } else if n + 1 >= poly.len() { "" - } else if coeff >= &conversion(0) { - " + " } else { - " - " + " + " }, Coefficient(coeff.clone()), Monomial(n), diff --git a/src/main.rs b/src/main.rs index b6ce57c..988a6bd 100644 --- a/src/main.rs +++ b/src/main.rs @@ -3,12 +3,15 @@ use num::{ 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); @@ -80,17 +83,39 @@ fn main() { .parse::<usize>() .unwrap(); - println!("poly is:"); - print_poly(&poly); + if std::env::args().len() >= 4 { + // Search for primes modulo which the polynomial is irreducible instead. - println!("prime is {p}"); + let primes = generate_primes(p); - // let poly1 = [12, 12, 0, 1].map(conversion); - // let poly2 = [-4, 1].map(conversion); + println!("Searching {} primes up to {p}...", primes.len()); - // let composition = berlekamp::composition(&poly1, &poly2); + 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; + } - use berlekamp::Poly; + // 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 {}", @@ -107,44 +132,36 @@ fn main() { let factors = berlekamp::factors(&poly, p); - match factors { - Ok(factors) => { - println!("the irreducible factors are:"); + print!("{} ≡ ", Poly(&poly)); - for factor in factors { - print_poly(&factor); - } - } - Err(gcd) => { - println!( - "the gcd of the polynomial with its derivative is {}", - Poly(&gcd) - ); + let mut total = 0usize; - if gcd.len() == poly.len() { - // the derivative is zero - println!("the polynomial is a {p}-th power of some polynomial"); + for (index, (mul, f)) in factors.into_iter().enumerate() { + total += mul; - return; + match mul { + 0 => { + dbg!("zero multiplicity?"); } - - 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); + 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)] @@ -315,3 +332,44 @@ fn print_matrix<const X: usize, const Y: usize, const T: usize>(matrix: &[BigInt } 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<usize> { + if bound <= 1 { + return Vec::new(); + } + + let sqrt = (bound as f64).sqrt() as usize; + + let mut result: Vec<usize> = Vec::with_capacity(bound); + + let mut records: Vec<bool> = 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 +} |