diff --git a/src/buffer.rs b/src/buffer.rs index 1563d36..2189234 100644 --- a/src/buffer.rs +++ b/src/buffer.rs @@ -113,6 +113,8 @@ pub trait PrimeBufferExt: for<'a> PrimeBuffer<'a> { /// factorization failed, then a list of found factors (not necessarily primes) will be returned. A prime /// factor will repeat if its exponent is larget than one, and it's ensured that the product of the list of /// factors is equal to the original target. + /// + /// TODO(v0.next): Return two lists when failed, one for prime factors, another one for remaining cofactors fn factors( &self, target: T, @@ -475,6 +477,7 @@ impl NaiveBuffer { #[cfg(test)] mod tests { use super::*; + use crate::detail::Mint; #[cfg(feature = "num-bigint")] use core::str::FromStr; #[cfg(feature = "num-bigint")] @@ -536,23 +539,27 @@ mod tests { } #[test] - fn pb_is_prime_test() { + fn is_prime_test() { // test for is_prime let pb = NaiveBuffer::new(); // some mersenne numbers - assert!(matches!( + assert_eq!( pb.is_prime(&(2u32.pow(19) - 1), None), Primality::Yes - )); - assert!(matches!( + ); + assert_eq!( pb.is_prime(&(2u32.pow(23) - 1), None), Primality::No - )); + ); assert!(matches!( pb.is_prime(&(2u128.pow(89) - 1), None), Primality::Probable(_) )); + assert!(matches!( + pb.is_prime(&Mint::from(2u128.pow(89) - 1), None), + Primality::Probable(_) + )); // test against small prime assertion for _ in 0..100 { @@ -567,12 +574,27 @@ mod tests { } // test large numbers + const P: u128 = 18699199384836356663; // https://golang.org/issue/638 + assert!(matches!( + pb.is_prime(&P, None), + Primality::Probable(_) + )); + assert!(matches!( + pb.is_prime(&P, Some(PrimalityTestConfig::bpsw())), + Primality::Probable(_) + )); + assert!(matches!( + pb.is_prime(&Mint::from(P), None), + Primality::Probable(_) + )); + assert!(matches!( + pb.is_prime(&Mint::from(P), Some(PrimalityTestConfig::bpsw())), + Primality::Probable(_) + )); + #[cfg(feature = "num-bigint")] { let large_primes = [ - // https://golang.org/issue/638 - "18699199384836356663", - "98920366548084643601728869055592650835572950932266967461790948584315647051443", "94560208308847015747498523884063394671606671904944666360068158221458669711639", @@ -587,6 +609,9 @@ mod tests { "9850501549098619803069760025035903451269934817616361666987073351061430442874302652853566563721228910201656997576599", // E-382: 2^382-105 "42307582002575910332922579714097346549017899709713998034217522897561970639123926132812109468141778230245837569601494931472367", // Curve41417: 2^414-17 "6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151", // E-521: 2^521-1 + + // https://github.com/AtropineTears/num-primes/issues/1#issuecomment-934629597 + "169511182982703321453314585423962898651587669459838234386506572286328885534468792292646838949809616446341407457141008401355628947670484184607678853094537849610289912805960069455687743151708433319901176932959509872662610091644590437761688516626993416011399330087939042347256922771590903190536793274742859624657" ]; for pstr in large_primes { assert!( diff --git a/src/factor.rs b/src/factor.rs index 17aeb8f..7fce830 100644 --- a/src/factor.rs +++ b/src/factor.rs @@ -103,6 +103,9 @@ 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. +/// +/// Note that the multiplier is usually selected from the following sequence: +/// `[1, 3, 5, 7, 11, 3*5, 3*7, 3*11, .., 3*5*7*11]` pub fn squfof(target: &T, multiplier: T) -> Option where for<'r> &'r T: RefNum, @@ -176,8 +179,12 @@ where } } -// TODO: ECM, Quadratic sieve / Prime field sieve, Fermat(https://en.wikipedia.org/wiki/Fermat%27s_factorization_method) +// TODO: ECM, One line, Quadratic sieve / Prime field sieve, Fermat(https://en.wikipedia.org/wiki/Fermat%27s_factorization_method) // REF: https://pypi.org/project/primefac/ +// http://flintlib.org/doc/ulong_extras.html#factorisation +// https://github.com/zademn/facto-rs/ +// https://github.com/elmomoilanen/prime-factorization +// https://cseweb.ucsd.edu/~ethome/teaching/2022-cse-291-14/ fn pollard_pp1() {} fn williams_pp1() {} diff --git a/src/lib.rs b/src/lib.rs index ef14019..366dbe6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -71,9 +71,10 @@ pub use traits::*; pub mod detail { //! Implementation details for this crate. //! - //! The structs and traits in this module are exposed - //! for public use, but they are not designed for such usage. User-friendly is not a goal here. - //! Some traits in this module can be used to extend `num-prime` with new backends. + //! The structs and traits in this module are exposed for public use, although they are no + //! designed for such usage. User-friendly is not a goal and backward-compatilibity is not + //! strictly maintained here. Some traits in this module can be used to extend `num-prime` + //! with new backends. pub use super::primality::{LucasUtils, PrimalityBase, PrimalityRefBase}; pub use super::tables::SMALL_PRIMES; pub use super::mint::Mint; diff --git a/src/mint.rs b/src/mint.rs index 1973e1a..1ded1bc 100644 --- a/src/mint.rs +++ b/src/mint.rs @@ -733,18 +733,27 @@ where T::Inv: Clone, { #[inline] fn checked_jacobi(&self, n: &Self) -> Option { - let (a, n) = left_ref_only(self, n); - a.checked_jacobi(n) + match (&self.0, &n.0) { + (Left(a), Left(n)) => a.checked_jacobi(n), + (Right(a), Left(n)) => a.residue().checked_jacobi(n), + (_, Right(_)) => unreachable!() + } } #[inline] fn checked_legendre(&self, n: &Self) -> Option { - let (a, n) = left_ref_only(self, n); - a.checked_legendre(n) + match (&self.0, &n.0) { + (Left(a), Left(n)) => a.checked_legendre(n), + (Right(a), Left(n)) => a.residue().checked_legendre(n), + (_, Right(_)) => unreachable!() + } } #[inline] fn kronecker(&self, n: &Self) -> i8 { - let (a, n) = left_ref_only(self, n); - a.kronecker(n) + match (&self.0, &n.0) { + (Left(a), Left(n)) => a.kronecker(n), + (Right(a), Left(n)) => a.residue().kronecker(n), + (_, Right(_)) => unreachable!() + } } } @@ -762,7 +771,6 @@ where })) } } -// TODO: implement ModularRefOps #[cfg(test)] mod tests { diff --git a/src/nt_funcs.rs b/src/nt_funcs.rs index b96534d..ed5b284 100644 --- a/src/nt_funcs.rs +++ b/src/nt_funcs.rs @@ -56,6 +56,13 @@ pub fn is_prime64(target: u64) -> bool { } // Then do a deterministic Miller-rabin test + is_prime64_miller(target) +} + +// Primality test for u64 with only miller-rabin tests, used during factorization. +// It assumes the target is odd, not too small and cannot be divided small primes +#[cfg(not(feature = "big-table"))] +fn is_prime64_miller(target: u64) -> bool { // The collection of witnesses are from http://miller-rabin.appspot.com/ if let Ok(u) = u16::try_from(target) { // 2, 3 for u16 range @@ -86,6 +93,8 @@ pub fn is_prime64(target: u64) -> bool { if target & 1 == 0 { return target == 2; } + + // trial division if target < SMALL_PRIMES_NEXT { // find in the prime list if the target is small enough return SMALL_PRIMES.binary_search(&(target as u16)).is_ok(); @@ -97,6 +106,13 @@ pub fn is_prime64(target: u64) -> bool { } } + is_prime64_miller(target) +} + +// Primality test for u64 with only miller-rabin tests, used during factorization. +// It assumes the target is odd, not too small and cannot be divided small primes +#[cfg(feature = "big-table")] +fn is_prime64_miller(target: u64) -> bool { // 32bit test const MAGIC: u32 = 0xAD625B89; if let Ok(u) = u32::try_from(target) { @@ -159,16 +175,15 @@ pub fn factorize64(target: u64) -> BTreeMap { let mut factored = false; #[cfg(not(feature = "big-table"))] - for &p in SMALL_PRIMES.iter().skip(1) { - let p64 = p as u64; - if p64 > tsqrt { + for p in SMALL_PRIMES.iter().skip(1).map(|&v| v as u64) { + if p > tsqrt { factored = true; break; } - while residual % p64 == 0 { - residual = residual / p64; - *result.entry(p64).or_insert(0) += 1; + while residual % p == 0 { + residual = residual / p; + *result.entry(p).or_insert(0) += 1; } if residual == 1 { factored = true; @@ -239,9 +254,7 @@ pub fn factorize64(target: u64) -> BTreeMap { 3 * 5 * 7 * 11, ]; while let Some(target) = todo.pop() { - // TODO: add a separate method (is_prime64_mint?) which skips the trial division part, as it's - // already performed by division above - if is_prime64(target) { + if is_prime64_miller(target) { *result.entry(target).or_insert(0) += 1; } else { let mut i = 1usize; @@ -268,7 +281,54 @@ pub fn factorize64(target: u64) -> BTreeMap { result } -// TODO: support factorize128, as we have efficient modular arithmetic for u128 +pub fn factorize128(target: u128) -> BTreeMap { + // shortcut for u64 + if target < (1u128 << 64) { + return factorize64(target as u64) + .into_iter() + .map(|(k, v)| (k as u128, v)) + .collect(); + } + + let mut result = BTreeMap::new(); + let f2 = target.trailing_zeros(); // quick check on factors of 2 + if f2 != 0 { + result.insert(2, f2 as usize); + } + + // trial division using primes in the table + // TODO(v0.3.2): speed up this by precompute tables + let mut residual = target >> f2; + for p in SMALL_PRIMES.iter().skip(1).map(|&v| v as u128) { + while residual % p == 0 { + residual = residual / p; + *result.entry(p).or_insert(0) += 1; + } + } + if residual == 1 { + return result; + } + + // then try pollard's rho and SQUFOF methods util fully factored + let mut todo = vec![residual]; + while let Some(target) = todo.pop() { + if is_prime(&Mint::from(target), Some(PrimalityTestConfig::bpsw())).probably() { + *result.entry(target).or_insert(0) += 1; + } else { + let divisor = loop { + // TODO: only pollard rho is used by now, select better methods + let start = MontgomeryInt::new(random::(), target); + let offset = start.convert(random::()); + if let Some(p) = pollard_rho(&Mint::from(target), start.into(), offset.into()) { + break p.value(); + } + }; + todo.push(divisor); + todo.push(target / divisor); + } + } + result +} /// This function re-exports [PrimeBufferExt::is_prime()][crate::buffer::PrimeBufferExt::is_prime()] with a default buffer distance pub fn is_prime(target: &T, config: Option) -> Primality @@ -569,8 +629,8 @@ where { let buf = NaiveBuffer::new(); let config = Some(PrimalityTestConfig::strict()); - // TODO: use miller-rabin for large numbers (more than 256 bits?), as BPSW could be too slow - // the NIST recommends 5 rounds for 512 and 1024 bits. For 1536 bits, the recommendation is 4 rounds. + // XXX: use miller-rabin for large numbers (more than 256 bits?), as BPSW could be too slow (need check) + // the NIST recommends 5 rounds for 512 and 1024 bits. For 1536 bits, the recommendation is 4 rounds. // test (n-1)/2 first since its smaller let sophie_p = buf.is_prime(&(target >> 1), config); @@ -986,14 +1046,30 @@ mod tests { } #[test] - #[cfg(feature = "big-int")] - fn is_prime_test() { - #[cfg(feature = "num-bigint")] - { - use num_bigint::BigUint; - // https://github.com/AtropineTears/num-primes/issues/1#issuecomment-934629597 - let p = BigUint::parse_bytes(b"169511182982703321453314585423962898651587669459838234386506572286328885534468792292646838949809616446341407457141008401355628947670484184607678853094537849610289912805960069455687743151708433319901176932959509872662610091644590437761688516626993416011399330087939042347256922771590903190536793274742859624657", 10).unwrap(); - assert!(is_prime(&p, None).probably()); + fn factorize128_test() { + // some known cases + let fac_primorial19 = BTreeMap::from_iter(SMALL_PRIMES.iter().take(19).map(|&p| (p as u128, 1))); + let fac = factorize128(7858321551080267055879090); + assert_eq!(fac, fac_primorial19); + + let fac_smallbig = BTreeMap::from_iter([(167, 1), (2417851639229258349412369, 1)]); + let fac = factorize128(403781223751286144351865623); + assert_eq!(fac, fac_smallbig); + + // random factorization tests + for _ in 0..1 { // TODO: run more tests when other factorization methods are implemented + let x = random(); + let fac = factorize128(x); + let mut prod = 1; + for (p, exp) in fac { + assert!( + is_prime(&p, None).probably(), + "factorization result should have prime factors! (get {})", + p + ); + prod *= p.pow(exp as u32); + } + assert_eq!(x, prod, "factorization check failed! ({} != {})", x, prod); } } diff --git a/src/primality.rs b/src/primality.rs index 07ccc34..f872077 100644 --- a/src/primality.rs +++ b/src/primality.rs @@ -44,8 +44,8 @@ where T::from_isize(-q).unwrap().negm(&m) }; - let mut uk = T::zero(); // U(k) - let mut uk1 = T::one(); // U(k+1) + let mut uk = T::zero() % &m; // U(k), mod m for montgomery form + let mut uk1 = T::one() % &m; // U(k+1) for i in (0..n.bits()).rev() { if n.bit(i) { diff --git a/src/traits.rs b/src/traits.rs index 546f8c0..2ae6e18 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -18,7 +18,7 @@ pub trait BitTest { } /// This enum describes the result of primality checks -#[derive(Debug, Clone, Copy)] +#[derive(Debug, Clone, Copy, PartialEq)] pub enum Primality { /// The number passes deterministic primality check. Yes, @@ -82,6 +82,7 @@ impl BitOr for Primality { pub struct PrimalityTestConfig { // TODO: add option to divides small primes in the table // and this option should be enabled if the probabilistic test is used for strict config + /// Number of strong probable prime test, starting from base 2 pub sprp_trials: usize, @@ -209,6 +210,7 @@ pub trait ExactRoots: Roots + Pow + Clone { // factor.c `divexact_21` // TODO: implement quick is_x_power (specifically is_235_power) +// This could be used during factorization to filter out perfect powers // REF: PARI/GP `Z_ispowerall`, `is_357_power` // FLINT `n_is_perfect_power235`, `fmpz_is_perfect_power` // GMP `mpz_perfect_power_p`