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.

1 files changed, 96 insertions, 38 deletions

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