summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/lib.rs103
-rw-r--r--src/main.rs134
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
+}