Add factorize128

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<T: PrimalityBase>(
         &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!(

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<T: Integer + NumRef + Clone + ExactRoots>(target: &T, multiplier: T) -> Option<T>
 where
     for<'r> &'r T: RefNum<T>,
@@ -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() {}
 

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;

src/mint.rs

@@ -733,18 +733,27 @@ where
     T::Inv: Clone, {
     #[inline]
     fn checked_jacobi(&self, n: &Self) -> Option<i8> {
-        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<i8> {
-        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 {

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<u64, usize> {
     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<u64, usize> {
         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<u64, usize> {
     result
 }
 
-// TODO: support factorize128, as we have efficient modular arithmetic for u128
+pub fn factorize128(target: u128) -> BTreeMap<u128, usize> {
+    // 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::<u128>(), target);
+                let offset = start.convert(random::<u128>());
+                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<T: PrimalityBase>(target: &T, config: Option<PrimalityTestConfig>) -> 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);
         }
     }
 

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) {

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<Primality> 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<u32, Output = Self> + 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`