RefinedHEADmaster

* src/lib.rs: Now the `factors` function directly finds all
  irreducible factors of the polynomial without resoting to manual
  inspections when the polynomial is not square-free.  Also some
  quality-of-life improvements were made.

* src/main.rs: Now has the option of searching for primes modulo which
  the polynomial is irreducible, up till a given bound.  This is
  useful if I am curious to know whether such small primes exist.

2 files changed, 179 insertions, 58 deletions

diff --git a/src/lib.rs b/src/lib.rs
index db3186a..8d8c3a0 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -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(&quotient)
-            );
-
-            let factors =
-                berlekamp::factors(&quotient, 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
+}