diff --git a/src/factor.rs b/src/factor.rs index cff1d9c..17aeb8f 100644 --- a/src/factor.rs +++ b/src/factor.rs @@ -85,6 +85,7 @@ where return None; } + // FIXME: optimize abs_diff for montgomery form let diff = if b > a { &b - &a } else { &a - &b }; // abs_diff let d = diff.gcd(target); if d > T::one() && &d < target { @@ -183,13 +184,27 @@ fn williams_pp1() {} #[cfg(test)] mod tests { use super::*; + use crate::mint::Mint; + use num_modular::MontgomeryInt; use rand::random; #[test] fn pollard_rho_test() { - assert!(matches!(pollard_rho(&8051u16, 2, 1), Some(97))); + assert_eq!(pollard_rho(&8051u16, 2, 1), Some(97)); assert!(matches!(pollard_rho(&8051u16, random(), 1), Some(i) if i == 97 || i == 83)); - assert!(matches!(pollard_rho(&455459u32, 2, 1), Some(743))) + assert_eq!(pollard_rho(&455459u32, 2, 1), Some(743)); + + // Mint test + for _ in 0..10 { + let target = random::() | 1; + let start = random::() % target; + let offset = random::() % target; + assert_eq!(pollard_rho(&target, start, offset), + pollard_rho(&Mint::from(target), + MontgomeryInt::new(start, target).into(), + MontgomeryInt::new(offset, target).into() + ).map(|v| v.value())); + } } #[test] diff --git a/src/mint.rs b/src/mint.rs index d2ff61e..1973e1a 100644 --- a/src/mint.rs +++ b/src/mint.rs @@ -1,17 +1,19 @@ //! Wrapper of integer to makes it efficient in modular arithmetics but still have the same //! API of normal integers. -use core::{ops::*, panic}; +use core::ops::*; use either::*; use num_traits::{Num, Zero, One, FromPrimitive, ToPrimitive, Pow}; use num_integer::{Integer, Roots}; -use num_modular::{ModularInteger, Montgomery, MontgomeryInt}; +use num_modular::{ModularInteger, Montgomery, MontgomeryInt, ModularCoreOps, ModularUnaryOps, ModularSymbols, ModularPow}; use crate::{ExactRoots, BitTest}; -// TODO (v0.3.x): Implement PrimalityBase for Mint - -/// Integer with fast modular arithmetics support, based on MontgomeryInt +/// Integer with fast modular arithmetics support, based on [MontgomeryInt] under the hood +/// +/// This struct only designed to be working with this crate. Most binary operators assume that +/// the modulus of two operands (when in montgomery form) are the same, and most implicit conversions +/// between conventional form and montgomery form will be forbidden pub struct Mint(Either>); // it seems that auto derivation struggles to provide an implementation for Copy, Clone and Debug with proper trait bounds @@ -35,12 +37,19 @@ impl From for Mint { Self(Left(v)) } } +impl From> for Mint { + #[inline(always)] + fn from(v: MontgomeryInt) -> Self { + Self(Right(v)) + } +} + #[inline(always)] fn left_only(lhs: Mint, rhs: Mint) -> (T, T) { match(lhs.0, rhs.0) { (Left(v1), Left(v2)) => (v1, v2), - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } @@ -48,7 +57,7 @@ fn left_only(lhs: Mint, rhs: Mint) -> (T, T) { fn left_ref_only<'a, T: Integer + Montgomery>(lhs: &'a Mint, rhs: &'a Mint) -> (&'a T, &'a T) { match(&lhs.0, &rhs.0) { (Left(v1), Left(v2)) => (v1, v2), - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } @@ -88,12 +97,8 @@ where fn eq(&self, other: &Self) -> bool { match (&self.0, &other.0) { (Left(v1), Left(v2)) => v1 == v2, - (Left(v1), Right(v2)) => v1 == &v2.residue(), - (Right(v1), Left(v2)) => &v1.residue() == v2, - (Right(v1), Right(v2)) => { - debug_assert!(v1.modulus() == v2.modulus()); - v1.residue() == v2.residue() - }, + (Right(v1), Right(v2)) => v1 == v2, + (_, _) => unreachable!() // force optimization of equality test } } } @@ -211,14 +216,14 @@ macro_rules! forward_binops_right { }) } } - impl $imp<&'r T, Output = T>> $imp for &Mint + impl<'a, 'b, T: Integer + Montgomery + Clone + for<'r> $imp<&'r T, Output = T>> $imp<&'b Mint> for &'a Mint where T::Double: From, T::Inv: Clone, { type Output = Mint; #[inline] - fn $method(self, rhs: Self) -> Self::Output { + fn $method(self, rhs: &Mint) -> Self::Output { Mint(match (&self.0, &rhs.0) { (Left(v1), Left(v2)) => Left(v1.clone().$method(v2)), (Left(v1), Right(v2)) => Right(v2.convert(v1.clone()).$method(v2)), @@ -261,7 +266,7 @@ where fn div(self, rhs: &Self) -> Self::Output { match(self.0, &rhs.0) { (Left(v1), Left(v2)) => Self(Left(v1.div(v2))), - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } } @@ -276,22 +281,21 @@ where fn div(self, rhs: Mint) -> Self::Output { match(&self.0, rhs.0) { (Left(v1), Left(v2)) => Mint(Left(v1.clone().div(v2))), - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } } -impl Div<&'r T, Output=T>> Div for &Mint +impl<'a, 'b, T: Integer + Montgomery + Clone + for<'r> Div<&'r T, Output=T>> Div<&'b Mint> for &'a Mint where T::Double: From, T::Inv: Clone, { type Output = Mint; - #[inline] - fn div(self, rhs: Self) -> Self::Output { + fn div(self, rhs: &Mint) -> Self::Output { match(&self.0, &rhs.0) { (Left(v1), Left(v2)) => Mint(Left(v1.clone().div(v2))), - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } } @@ -311,7 +315,7 @@ where debug_assert!(v1.modulus() == &v2); Self(Right(v1)) } - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } } @@ -330,7 +334,7 @@ where debug_assert!(v1.modulus() == v2); Self(Right(v1)) } - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } } @@ -349,11 +353,11 @@ where debug_assert!(v1.modulus() == &v2); Mint(Right(v1.clone())) } - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } } -impl Rem<&'r T, Output=T>> Rem for &Mint +impl<'a, 'b, T: Integer + Montgomery + Clone + for<'r> Rem<&'r T, Output=T>> Rem<&'b Mint> for &'a Mint where T::Double: From, T::Inv: Clone, @@ -361,14 +365,14 @@ where type Output = Mint; #[inline] - fn rem(self, rhs: Self) -> Self::Output { + fn rem(self, rhs: &Mint) -> Self::Output { match(&self.0, &rhs.0) { (Left(v1), Left(v2)) => Mint(Right(MontgomeryInt::new(v1.clone(), v2.clone()))), (Right(v1), Left(v2)) => { debug_assert!(v1.modulus() == v2); Mint(Right(v1.clone())) } - (_, _) => panic!("not supported"), + (_, _) => unreachable!(), } } } @@ -423,7 +427,6 @@ where { forward_binops_left_ref_only!(div_floor); forward_binops_left_ref_only!(mod_floor); - forward_binops_left_ref_only!(gcd); forward_binops_left_ref_only!(lcm); forward_binops_left_ref_only!(divides => bool); forward_binops_left_ref_only!(is_multiple_of => bool); @@ -436,16 +439,26 @@ where let (q, r) = v1.div_rem(v2); (Self(Left(q)), Self(Left(r))) } + #[inline(always)] + fn gcd(&self, other: &Self) -> Self { + Self(Left(match (&self.0, &other.0) { + (Left(v1), Left(v2)) => v1.gcd(v2), + (Right(v1), Left(v2)) => v1.residue().gcd(v2), + (Left(v1), Right(v2)) => v1.gcd(&v2.residue()), + (Right(v1), Right(v2)) => v1.residue().gcd(&v2.residue()) + })) + } } impl Roots for Mint where T::Double: From, T::Inv: Clone, { + #[inline] fn nth_root(&self, n: u32) -> Self { match &self.0 { Left(v) => Self(Left(v.nth_root(n))), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } } @@ -454,12 +467,15 @@ impl FromPrimitive for Mint where T::Double: From, T::Inv: Clone, { + #[inline] fn from_f64(n: f64) -> Option { T::from_f64(n).map(|v| Self(Left(v))) } + #[inline] fn from_i64(n: i64) -> Option { T::from_i64(n).map(|v| Self(Left(v))) } + #[inline] fn from_u64(n: u64) -> Option { T::from_u64(n).map(|v| Self(Left(v))) } @@ -469,18 +485,21 @@ impl ToPrimitive for Mint where T::Double: From, T::Inv: Clone, { + #[inline] fn to_f64(&self) -> Option { match &self.0 { Left(v) => v.to_f64(), Right(m) => m.residue().to_f64() } } + #[inline] fn to_i64(&self) -> Option { match &self.0 { Left(v) => v.to_i64(), Right(m) => m.residue().to_i64() } } + #[inline] fn to_u64(&self) -> Option { match &self.0 { Left(v) => v.to_u64(), @@ -494,10 +513,11 @@ where T::Double: From, T::Inv: Clone, { type Output = Self; + #[inline] fn pow(self, rhs: u32) -> Self::Output { match self.0 { Left(v) => Self(Left(v.pow(rhs))), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } } @@ -506,10 +526,11 @@ impl ExactRoots for Mint where T::Double: From, T::Inv: Clone, { + #[inline] fn nth_root_exact(&self, n: u32) -> Option { match &self.0 { Left(v) => v.nth_root_exact(n).map(|v| Self(Left(v))), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } } @@ -518,53 +539,229 @@ impl BitTest for Mint where T::Double: From, T::Inv: Clone, { + #[inline] fn bit(&self, position: usize) -> bool { match &self.0 { Left(v) => v.bit(position), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } + #[inline] fn bits(&self) -> usize { match &self.0 { Left(v) => v.bits(), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } + #[inline] fn trailing_zeros(&self) -> usize { match &self.0 { Left(v) => v.trailing_zeros(), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } } -impl> Shr for Mint +impl> Shr for Mint where T::Double: From, T::Inv: Clone, { type Output = Self; - fn shr(self, rhs: u32) -> Self::Output { + #[inline] + fn shr(self, rhs: usize) -> Self::Output { match self.0 { Left(v) => Self(Left(v >> rhs)), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } } -impl> Shr for &Mint +impl> Shr for &Mint where T::Double: From, T::Inv: Clone, { type Output = Mint; - fn shr(self, rhs: u32) -> Self::Output { + #[inline] + fn shr(self, rhs: usize) -> Self::Output { match &self.0 { Left(v) => Mint(Left(v.clone() >> rhs)), - Right(_) => panic!("not supported") + Right(_) => unreachable!() } } } +impl ModularCoreOps<&Self, &Self> for Mint where + T::Double: From, + T::Inv: Clone, { + type Output = Self; + #[inline] + fn addm(self, rhs: &Self, m: &Self) -> Self::Output { + match(self.0, &rhs.0, &m.0) { + (Right(v1), Right(v2), Left(m)) => { + debug_assert!(v1.modulus() == m && v2.modulus() == m); + Self(Right(v1 + v2)) + }, + (_, _, _) => unreachable!() + } + } + #[inline] + fn subm(self, rhs: &Self, m: &Self) -> Self::Output { + match(self.0, &rhs.0, &m.0) { + (Right(v1), Right(v2), Left(m)) => { + debug_assert!(v1.modulus() == m && v2.modulus() == m); + Self(Right(v1 - v2)) + }, + (_, _, _) => unreachable!() + } + } + #[inline] + fn mulm(self, rhs: &Self, m: &Self) -> Self::Output { + match(self.0, &rhs.0, &m.0) { + (Right(v1), Right(v2), Left(m)) => { + debug_assert!(v1.modulus() == m && v2.modulus() == m); + Self(Right(v1 * v2)) + }, + (_, _, _) => unreachable!() + } + } +} +impl<'a, 'b, T: Integer + Montgomery + Clone> ModularCoreOps<&'b Mint, &'b Mint> for &'a Mint where + T::Double: From, + T::Inv: Clone, { + type Output = Mint; + #[inline] + fn addm(self, rhs: &Mint, m: &Mint) -> Self::Output { + match(&self.0, &rhs.0, &m.0) { + (Right(v1), Right(v2), Left(m)) => { + debug_assert!(v1.modulus() == m && v2.modulus() == m); + Mint(Right(v1 + v2)) + }, + (_, _, _) => unreachable!() + } + } + #[inline] + fn subm(self, rhs: &Mint, m: &Mint) -> Self::Output { + match(&self.0, &rhs.0, &m.0) { + (Right(v1), Right(v2), Left(m)) => { + debug_assert!(v1.modulus() == m && v2.modulus() == m); + Mint(Right(v1 - v2)) + }, + (_, _, _) => unreachable!() + } + } + #[inline] + fn mulm(self, rhs: &Mint, m: &Mint) -> Self::Output { + match(&self.0, &rhs.0, &m.0) { + (Right(v1), Right(v2), Left(m)) => { + debug_assert!(v1.modulus() == m && v2.modulus() == m); + Mint(Right(v1 * v2)) + }, + (_, _, _) => unreachable!() + } + } +} +impl ModularUnaryOps<&Self> for Mint +where + T::Double: From, + T::Inv: Clone, { + type Output = Self; + #[inline] + fn negm(self, m: &Self) -> Self::Output { + Self(Right(match (self.0, &m.0) { + (Left(v), Left(m)) => MontgomeryInt::new(v, m.clone()).neg(), + (Right(v), Left(m)) => { debug_assert!(v.modulus() == m); v.neg() } + (_, Right(_)) => unreachable!() + })) + } + fn invm(self, _: &Self) -> Option { + unreachable!() // not used in this crate + } + #[inline] + fn dblm(self, m: &Self) -> Self::Output { + Self(Right(match (self.0, &m.0) { + (Left(v), Left(m)) => MontgomeryInt::new(v, m.clone()).double(), + (Right(v), Left(m)) => { debug_assert!(v.modulus() == m); v.double() } + (_, Right(_)) => unreachable!() + })) + } + #[inline] + fn sqm(self, m: &Self) -> Self::Output { + Self(Right(match (self.0, &m.0) { + (Left(v), Left(m)) => MontgomeryInt::new(v, m.clone()).square(), + (Right(v), Left(m)) => { debug_assert!(v.modulus() == m); v.square() } + (_, Right(_)) => unreachable!() + })) + } +} +impl<'a, 'b, T: Integer + Montgomery + Clone> ModularUnaryOps<&'b Mint> for &'a Mint +where + T::Double: From, + T::Inv: Clone, { + type Output = Mint; + #[inline] + fn negm(self, m: &Mint) -> Self::Output { + Mint(Right(match (&self.0, &m.0) { + (Left(v), Left(m)) => MontgomeryInt::new(v.clone(), m.clone()).neg(), + (Right(v), Left(m)) => { debug_assert!(v.modulus() == m); v.clone().neg() } + (_, Right(_)) => unreachable!() + })) + } + fn invm(self, _: &Mint) -> Option { + unreachable!() // not used in this crate + } + #[inline] + fn dblm(self, m: &Mint) -> Self::Output { + Mint(Right(match (&self.0, &m.0) { + (Left(v), Left(m)) => MontgomeryInt::new(v.clone(), m.clone()).double(), + (Right(v), Left(m)) => { debug_assert!(v.modulus() == m); v.clone().double() } + (_, Right(_)) => unreachable!() + })) + } + #[inline] + fn sqm(self, m: &Mint) -> Self::Output { + Mint(Right(match (&self.0, &m.0) { + (Left(v), Left(m)) => MontgomeryInt::new(v.clone(), m.clone()).square(), + (Right(v), Left(m)) => { debug_assert!(v.modulus() == m); v.clone().square() } + (_, Right(_)) => unreachable!() + })) + } +} +impl ModularSymbols<&'r T>> ModularSymbols<&Self> for Mint +where + T::Double: From, + T::Inv: Clone, { + #[inline] + fn checked_jacobi(&self, n: &Self) -> Option { + let (a, n) = left_ref_only(self, n); + a.checked_jacobi(n) + } + #[inline] + fn checked_legendre(&self, n: &Self) -> Option { + let (a, n) = left_ref_only(self, n); + a.checked_legendre(n) + } + #[inline] + fn kronecker(&self, n: &Self) -> i8 { + let (a, n) = left_ref_only(self, n); + a.kronecker(n) + } +} + +impl ModularPow<&Self, &Self> for Mint +where + T::Double: From, + T::Inv: Clone, { + type Output = Self; + #[inline] + fn powm(self, exp: &Self, m: &Self) -> Self::Output { + Self(Right(match (self.0, &exp.0, &m.0) { + (Left(v), Left(e), Left(m)) => MontgomeryInt::new(v, m.clone()).pow(e.clone()), + (Right(v), Left(e), Left(m)) => { debug_assert!(v.modulus() == m); v.pow(e.clone()) } + (_, _, _) => unreachable!() + })) + } +} // TODO: implement ModularRefOps #[cfg(test)] @@ -573,10 +770,8 @@ mod tests { #[test] fn test_basics() { - let a: Mint = 12.into(); + let a: Mint = 19.into(); let b: Mint = 8.into(); - assert_eq!(a + b, 20.into()); - - dbg!(crate::nt_funcs::is_prime(&a, None)); + assert_eq!(a + b, 27.into()); } } diff --git a/src/nt_funcs.rs b/src/nt_funcs.rs index 76a392d..b96534d 100644 --- a/src/nt_funcs.rs +++ b/src/nt_funcs.rs @@ -12,9 +12,10 @@ //! use crate::buffer::{NaiveBuffer, PrimeBufferExt}; +use crate::mint::Mint; use crate::factor::{pollard_rho, squfof}; use crate::primality::{PrimalityBase, PrimalityRefBase}; -use crate::tables::{MOEBIUS_ODD, SMALL_PRIMES, WHEEL_NEXT, WHEEL_PREV, WHEEL_SIZE}; +use crate::tables::{MOEBIUS_ODD, SMALL_PRIMES, SMALL_PRIMES_NEXT, WHEEL_NEXT, WHEEL_PREV, WHEEL_SIZE}; #[cfg(feature = "big-table")] use crate::tables::{SMALL_PRIMES_INV, SMALL_PRIMES_INVLIM, ZETA_LOG_TABLE}; use crate::traits::{FactorizationConfig, Primality, PrimalityTestConfig, PrimalityUtils}; @@ -22,7 +23,7 @@ use crate::RandPrime; #[cfg(feature = "num-bigint")] use num_bigint::{BigUint, RandBigInt}; use num_integer::Roots; -use num_modular::ModularCoreOps; +use num_modular::{ModularCoreOps, MontgomeryInt, ModularInteger}; use num_traits::{CheckedAdd, FromPrimitive, Num, RefNum, ToPrimitive}; use rand::{random, Rng}; use std::collections::BTreeMap; @@ -37,32 +38,40 @@ use crate::tables::{MILLER_RABIN_BASE32, MILLER_RABIN_BASE64}; #[cfg(not(feature = "big-table"))] pub fn is_prime64(target: u64) -> bool { // shortcuts - if target < 1 { + if target < 2 { return false; } if target & 1 == 0 { return target == 2; - } - - // first find in the prime list + } if let Ok(u) = u8::try_from(target) { + // find in the prime list if the target is small enough return SMALL_PRIMES.binary_search(&u).is_ok(); + } else { + // check remainder against the wheel table + let pos = (target % WHEEL_SIZE as u64) as usize; + if pos == 0 || WHEEL_NEXT[pos] < WHEEL_NEXT[pos-1] { + return false; + } } // Then do a deterministic Miller-rabin test // The collection of witnesses are from http://miller-rabin.appspot.com/ if let Ok(u) = u16::try_from(target) { // 2, 3 for u16 range - return u.is_sprp(2) && u.is_sprp(3); + let u = Mint::from(u); + return u.is_sprp(Mint::from(2)) && u.is_sprp(Mint::from(3)); } if let Ok(u) = u32::try_from(target) { // 2, 7, 61 for u32 range - return u.is_sprp(2) && u.is_sprp(7) && u.is_sprp(61); + let u = Mint::from(u); + return u.is_sprp(Mint::from(2)) && u.is_sprp(Mint::from(7)) && u.is_sprp(Mint::from(61)); } // 2, 325, 9375, 28178, 450775, 9780504, 1795265022 for u64 range const WITNESS64: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]; - WITNESS64.iter().all(|&x| target.is_sprp(x)) + let u = Mint::from(target); + WITNESS64.iter().all(|&x| u.is_sprp(Mint::from(x))) } /// Very fast primality test on a u64 integer is a prime number. It's based on @@ -71,32 +80,40 @@ pub fn is_prime64(target: u64) -> bool { #[cfg(feature = "big-table")] pub fn is_prime64(target: u64) -> bool { // shortcuts - if target < 1 { + if target < 2 { return false; } if target & 1 == 0 { return target == 2; } - - // first find in the prime list - if target < 8167 { + 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(); + } else { + // check remainder against the wheel table + let pos = (target % WHEEL_SIZE as u64) as usize; + if pos == 0 || WHEEL_NEXT[pos] < WHEEL_NEXT[pos-1] { + return false; + } } // 32bit test const MAGIC: u32 = 0xAD625B89; if let Ok(u) = u32::try_from(target) { let base = u.wrapping_mul(MAGIC) >> 24; - return u.is_sprp(MILLER_RABIN_BASE32[base as usize] as u32); + let u = Mint::from(u); + return u.is_sprp(Mint::from(MILLER_RABIN_BASE32[base as usize] as u32)); } // 49bit test - if !target.is_sprp(2) { + let mt = Mint::from(target); + if !mt.is_sprp(2.into()) { return false; } let u = target as u32; // truncate + let base = u.wrapping_mul(MAGIC) >> 18; - if !target.is_sprp(MILLER_RABIN_BASE64[base as usize] as u64) { + if !mt.is_sprp(Mint::from(MILLER_RABIN_BASE64[base as usize] as u64)) { return false; } if target < (1u64 << 49) { @@ -106,7 +123,7 @@ pub fn is_prime64(target: u64) -> bool { // 64bit test const SECOND_BASES: [u64; 8] = [15, 135, 13, 60, 15, 117, 65, 29]; let base = base >> 13; - target.is_sprp(SECOND_BASES[base as usize]) + mt.is_sprp(Mint::from(SECOND_BASES[base as usize])) } /// Fast integer factorization on a u64 target. It's based on pollard's rho method and SQUFOF. @@ -222,6 +239,8 @@ 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) { *result.entry(target).or_insert(0) += 1; } else { @@ -229,15 +248,15 @@ pub fn factorize64(target: u64) -> BTreeMap { let divisor = loop { // try SQUFOF after 4 failed pollard rho trials if i % 5 == 0 && (i / 5) < SQUFOF_MULTIPLIERS.len() { - // TODO: check if the residual is a sqaure number before SQUFOF + // TODO: check if the residual is a sqaure number before SQUFOF (and also pollard_rho?) if let Some(p) = squfof(&target, SQUFOF_MULTIPLIERS[i / 5] as u64) { break p; } } else { - let start = random::() % target; - let offset = random::() % target; - if let Some(p) = pollard_rho(&target, start, offset) { - break p; + 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(); } } i += 1; @@ -249,7 +268,7 @@ pub fn factorize64(target: u64) -> BTreeMap { result } -// XXX: support factorize128, as we have efficient modular arithmetic for u128 +// TODO: support factorize128, as we have efficient modular arithmetic for u128 /// 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 diff --git a/src/primality.rs b/src/primality.rs index 695eec6..222fa71 100644 --- a/src/primality.rs +++ b/src/primality.rs @@ -148,6 +148,7 @@ where #[inline] fn is_sprp(&self, base: Self) -> bool { + // TODO(v0.3.1): test performance without either dispatch self.test_sprp(base).either(|v| v, |_| false) } @@ -161,23 +162,27 @@ where let shift = tm1.trailing_zeros(); let u = &tm1 >> shift; + // prevent reduction if the input is in montgomery form + let m1 = T::one() % self; + let mm1 = (&m1).negm(self); + let mut x = base.powm(&u, self); - if x == T::one() || x == tm1 { + if x == m1 || x == mm1 { return Either::Left(true); } for _ in 0..shift { let y = (&x).sqm(self); - if y.is_one() { + if y == m1 { return Either::Right(self.gcd(&(x - T::one()))); } - if y == tm1 { + if y == mm1 { return Either::Left(true); } x = y; } - Either::Left(x == T::one()) + Either::Left(x == m1) } fn is_lprp(&self, p: Option, q: Option) -> bool { @@ -386,6 +391,7 @@ impl PrimalityRefBase for T where #[cfg(test)] mod tests { use super::*; + use crate::mint::Mint; use num_modular::{ModularAbs, ModularSymbols}; use rand::random; @@ -411,10 +417,12 @@ mod tests { let spsp: [u16; 5] = [2047, 3277, 4033, 4681, 8321]; for psp in spsp { assert!(psp.is_sprp(2)); + assert!(Mint::from(psp).is_sprp(2.into())); // test Mint execution } // test cofactor return - assert!(matches!(341u16.test_sprp(2), Either::Right(31))); + assert_eq!(341u16.test_sprp(2), Either::Right(31)); + assert_eq!(Mint::from(341u16).test_sprp(2.into()), Either::Right(31.into())); } #[test]