diff --git a/examples/profile_factorization.rs b/examples/profile_factorization.rs index 56727cd..3947024 100644 --- a/examples/profile_factorization.rs +++ b/examples/profile_factorization.rs @@ -1,23 +1,29 @@ -use num_prime::factor::{pollard_rho, squfof, one_line}; +use num_prime::factor::{pollard_rho, squfof, one_line, SQUFOF_MULTIPLIERS}; use num_prime::RandPrime; use rand::random; fn main() { let mut rng = rand::thread_rng(); + // let p1: u64 = rng.gen_prime(40, None); + // let p2: u64 = rng.gen_prime(60, None); - let p1: u64 = rng.gen_prime(40, None); - let p2: u64 = rng.gen_prime(60, None); - - // let p1: u64 = 3486784447; - // let p2: u64 = 94143178889; + let p1: u64 = 3486784447; + let p2: u64 = 94143178889; // let p2: u64 = 3486784709; - let n = p1 as u128 * p2 as u128; println!("Semiprime: {} * {}", p1, p2); - println!("pollard rho result: {:?}", pollard_rho(&n, random(), random(), 400000)); - println!("squfof k=1 result: {:?}", squfof(&n, n, 40000)); - println!("squfof k=3*5*7*11 result: {:?}", squfof(&n, n * 3 * 5 * 7 * 11, 40000)); - println!("one_line k=1 result: {:?}", one_line(&n, n, 400000)); - println!("one_line k=480 result: {:?}", one_line(&n, n * 480, 400000)); + let n = p1 as u128 * p2 as u128; + + // let n: u128 = 133717415095455410877609739380293; + + const MAXITER: usize = 2 << 20; + for k in SQUFOF_MULTIPLIERS { + if let Some(kn) = n.checked_mul(k as u128) { + println!("squfof k={} result: {:?}", k, squfof(&n, kn, MAXITER)); + } + } + // println!("one_line k=1 result: {:?}", one_line(&n, n, MAXITER)); + // println!("one_line k=480 result: {:?}", one_line(&n, n * 480, MAXITER)); + // println!("pollard rho result: {:?}", pollard_rho(&n, random(), random(), MAXITER)); } diff --git a/src/factor.rs b/src/factor.rs index 90b81fd..efeacca 100644 --- a/src/factor.rs +++ b/src/factor.rs @@ -141,7 +141,7 @@ where } /// This function implements Shanks's square forms factorization (SQUFOF). It will assume that target -/// is not a perfect square and the multiplier is square-free. +/// is **not a perfect square**. /// /// The input is usually multiplied by a multiplier, and the multiplied integer should be put in /// the `mul_target` argument. The multiplier can be choosen from SQUFOF_MULTIPLIERS, or other square-free odd numbers. @@ -150,74 +150,63 @@ where /// The max iteration can be choosed as 2√(2√n), which is the theoretical upper limit for factorization. /// /// Reference: Gower, J., & Wagstaff Jr, S. (2008). Square form factorization. -/// In [Mathematics of Computation](https://homes.cerias.purdue.edu/~ssw/gowerthesis804/wthe.pdf) -/// or [thesis](https://homes.cerias.purdue.edu/~ssw/gowerthesis804/wthe.pdf) -pub fn squfof(target: &T, mul_target: T, max_iter: usize) -> (Option, usize) +/// In [1] [Mathematics of Computation](https://homes.cerias.purdue.edu/~ssw/gowerthesis804/wthe.pdf) +/// or [2] [his thesis](https://homes.cerias.purdue.edu/~ssw/gowerthesis804/wthe.pdf) +/// The code is from [3] [Rosetta code](https://rosettacode.org/wiki/Square_form_factorization) +pub fn squfof(target: &T, mul_target: T, max_iter: usize) -> (Option, usize) where for<'r> &'r T: RefNum, { assert!(&mul_target.is_multiple_of(&target), "mul_target should be multiples of target"); + let rd = Roots::sqrt(&mul_target); // root of k*N - // forward - let p0 = Roots::sqrt(&mul_target); - let mut pm1 = p0.clone(); - let mut p; // to be initialized in the first iteration - let mut qm1 = T::one(); - let mut q = &mul_target - &p0 * &p0; - let mut i = 1usize; - let qsqrt = loop { - let b = (&p0 + &pm1) / &q; - p = &b * &q - &pm1; - let qnext = if pm1 > p { - &qm1 + &b * (&pm1 - &p) + /// Reduction operator for binary quadratic forms. It's equivalent to + /// the one used in the `num-irrational` crate, in a little different form. + /// + /// This function reduces (a, b, c) = (qm1, p, q), updates qm1 and q, returns new p. + #[inline] + fn rho (rd: &T, p: &T, q: &mut T, qm1: &mut T) -> T where + for<'r> &'r T: RefNum, { + let b = (rd + p).div_floor(&*q); + let new_p = &b * &*q - p; + let new_q = if p > &new_p { + &*qm1 + b * (p - &new_p) } else { - &qm1 - &b * (&p - &pm1) + &*qm1 - b * (&new_p - p) }; + + *qm1 = std::mem::replace(q, new_q); + new_p + } + + // forward loop, search principal cycle + let (mut p, mut q, mut qm1) = (rd.clone(), &mul_target - &rd * &rd, T::one()); + for i in 1..max_iter { + p = rho(&rd, &p, &mut q, &mut qm1); if i.is_odd() { - if let Some(v) = qnext.sqrt_exact() { - break v; + if let Some(rq) = q.sqrt_exact() { + let b = (&rd - &p) / &rq; + let mut u = b * &rq + &p; + let (mut v, mut vm1) = ((&mul_target - &u * &u) / &rq, rq); + + // backward loop, search ambiguous cycle + loop { + let new_u = rho(&rd, &u, &mut v, &mut vm1); + if new_u == u { + break; + } else { + u = new_u + } + } + + let d = target.gcd(&u); + if d > T::one() && &d < target { + return (Some(d), i) + } } } - - pm1 = p; - qm1 = q; - q = qnext; - - i += 1; - if i == max_iter { - return (None, i); - } }; - - // backward - let b0 = (&p0 - &p) / &qsqrt; - pm1 = &b0 * &qsqrt + &p; - qm1 = qsqrt; - q = (&mul_target - &pm1 * &pm1) / &qm1; - - loop { - let b = (&p0 + &pm1) / &q; - p = &b * &q - &pm1; - if p == pm1 { - break; - } - - let qnext = if pm1 > p { - &qm1 + &b * (&pm1 - &p) - } else { - &qm1 - &b * (&p - &pm1) - }; - pm1 = p; - qm1 = q; - q = qnext; - } - - let d = target.gcd(&p); - if d > T::one() && &d < target { - (Some(d), i) - } else { - (None, i) - } + (None, max_iter) } // Square-free even numbers are suitable as SQUFOF multipliers @@ -323,6 +312,9 @@ mod tests { #[test] fn squfof_test() { assert_eq!(squfof(&11111u32, 11111u32, 100).0, Some(41)); + + // this case should success at step 276, from https://rosettacode.org/wiki/Talk:Square_form_factorization + assert!(matches!(squfof(&4558849u32, 4558849u32, 300).0, Some(_))); } #[test] diff --git a/src/nt_funcs.rs b/src/nt_funcs.rs index a08c335..a53dbbc 100644 --- a/src/nt_funcs.rs +++ b/src/nt_funcs.rs @@ -274,15 +274,15 @@ pub(crate) fn factorize64_advanced(cofactors: &[(u64, usize)]) -> Vec<(u64, usiz match i % NMETHODS { 0 => { - // Pollard's rho, allow 4x iterations since it's fast + // Pollard's rho let start = MontgomeryInt::new(random::(), target); let offset = start.convert(random::()); - if let (Some(p), _) = pollard_rho(&Mint::from(target), start.into(), offset.into(), max_iter * 4) { + if let (Some(p), _) = pollard_rho(&Mint::from(target), start.into(), offset.into(), max_iter) { break p.value(); } } 1 => { - // Hart's one-line, test 16 iterations + // Hart's one-line let mul_target = target.checked_mul(480).unwrap_or(target); if let (Some(p), _) = one_line(&target, mul_target, max_iter) { break p; @@ -290,17 +290,10 @@ pub(crate) fn factorize64_advanced(cofactors: &[(u64, usize)]) -> Vec<(u64, usiz } 2 => { // Shank's squfof - let n = i / NMETHODS; - if n >= SQUFOF_MULTIPLIERS.len() { - continue; - } - let mul_target = if let Some(kn) = target.checked_mul(SQUFOF_MULTIPLIERS[n] as u64) { - kn - } else { - continue; - }; - if let (Some(p), _) = squfof(&target, mul_target, max_iter) { - break p; + if let Some(mul_target) = target.checked_mul(SQUFOF_MULTIPLIERS[i % SQUFOF_MULTIPLIERS.len()] as u64) { + if let (Some(p), _) = squfof(&target, mul_target, max_iter) { + break p; + } } } _ => unreachable!(), @@ -431,12 +424,12 @@ pub(crate) fn factorize128_advanced(cofactors: &[(u128, usize)]) -> Vec<(u128, u // Pollard's rho let start = MontgomeryInt::new(random::(), target); let offset = start.convert(random::()); - if let (Some(p), _) = pollard_rho(&Mint::from(target), start.into(), offset.into(), max_iter * 4) { + if let (Some(p), _) = pollard_rho(&Mint::from(target), start.into(), offset.into(), max_iter) { break p.value(); } } 1 => { - // Hart's one-line, test 16 iterations + // Hart's one-line let mul_target = target.checked_mul(480).unwrap_or(target); if let (Some(p), _) = one_line(&target, mul_target, max_iter) { break p; @@ -444,17 +437,10 @@ pub(crate) fn factorize128_advanced(cofactors: &[(u128, usize)]) -> Vec<(u128, u } 2 => { // Shanks's squfof - let n = i / NMETHODS; - if n >= SQUFOF_MULTIPLIERS.len() { - continue; - } - let mul_target = if let Some(kn) = target.checked_mul(SQUFOF_MULTIPLIERS[n] as u128) { - kn - } else { - continue; - }; - if let (Some(p), _) = squfof(&target, mul_target, max_iter) { - break p; + if let Some(mul_target) = target.checked_mul(SQUFOF_MULTIPLIERS[i % NMETHODS] as u128) { + if let (Some(p), _) = squfof(&target, mul_target, max_iter) { + break p; + } } } _ => unreachable!(), @@ -464,7 +450,7 @@ pub(crate) fn factorize128_advanced(cofactors: &[(u128, usize)]) -> Vec<(u128, u // increase max iterations after trying all methods #[allow(unused_assignments)] if i % NMETHODS == 0 { - max_iter *= 2; + max_iter *= 4; } };