use crate::factor::{pollard_rho, trial_division};
use crate::nt_funcs::{
factorize128, is_prime64, next_prime, nth_prime_bounds, nth_prime_est, prev_prime,
};
use crate::primality::{PrimalityBase, PrimalityRefBase};
use crate::tables::{SMALL_PRIMES, SMALL_PRIMES_NEXT};
use crate::traits::{
FactorizationConfig, Primality, PrimalityTestConfig, PrimalityUtils, PrimeBuffer,
};
use bitvec::{bitvec, prelude::Msb0};
use lru::LruCache;
use num_integer::Roots;
use rand::random;
use std::collections::BTreeMap;
use std::num::NonZeroUsize;
pub trait PrimeBufferExt: for<'a> PrimeBuffer<'a> {
fn is_prime<T: PrimalityBase>(
&self,
target: &T,
config: Option<PrimalityTestConfig>,
) -> Primality
where
for<'r> &'r T: PrimalityRefBase<T>,
{
if target.is_even() {
return if target == &T::from_u8(2u8).unwrap() {
Primality::Yes
} else {
Primality::No
};
}
if let Some(x) = target.to_u64() {
return match is_prime64(x) {
true => Primality::Yes,
false => Primality::No,
};
}
let config = config.unwrap_or(PrimalityTestConfig::default());
let mut probability = 1.;
let mut witness_list: Vec<u64> = Vec::new();
if config.sprp_trials > 0 {
witness_list.extend(self.iter().take(config.sprp_trials));
probability *= 1. - 0.25f32.powi(config.sprp_trials as i32);
}
if config.sprp_random_trials > 0 {
for _ in 0..config.sprp_random_trials {
let mut w: u64 = rand::random();
while witness_list.iter().find(|&x| x == &w).is_some() {
w = rand::random();
}
witness_list.push(w);
}
probability *= 1. - 0.25f32.powi(config.sprp_random_trials as i32);
}
if !witness_list
.into_iter()
.all(|x| target.is_sprp(T::from_u64(x).unwrap()))
{
return Primality::No;
}
if config.slprp_test {
probability *= 1. - 4f32 / 15f32;
if !target.is_slprp(None, None) {
return Primality::No;
}
}
if config.eslprp_test {
probability *= 1. - 4f32 / 15f32;
if !target.is_eslprp(None) {
return Primality::No;
}
}
Primality::Probable(probability)
}
fn factors<T: PrimalityBase>(
&self,
target: T,
config: Option<FactorizationConfig>,
) -> (BTreeMap<T, usize>, Option<Vec<T>>)
where
for<'r> &'r T: PrimalityRefBase<T>,
{
if let Some(x) = target.to_u128() {
let factors = factorize128(x)
.into_iter()
.map(|(k, v)| (T::from_u128(k).unwrap(), v))
.collect();
return (factors, None);
}
let config = config.unwrap_or(FactorizationConfig::default());
let (result, factored) = trial_division(self.iter().cloned(), target, config.td_limit);
let mut result: BTreeMap<T, usize> = result
.into_iter()
.map(|(k, v)| (T::from_u64(k).unwrap(), v))
.collect();
let mut failed = Vec::new();
let mut config = config;
config.td_limit = Some(0); match factored {
Ok(res) => {
if !res.is_one() {
result.insert(res, 1);
}
}
Err(res) => {
let mut todo = vec![res];
while let Some(target) = todo.pop() {
if self
.is_prime(&target, Some(config.primality_config))
.probably()
{
*result.entry(target).or_insert(0) += 1;
} else {
if let Some(divisor) = self.divisor(&target, &mut config) {
todo.push(divisor.clone());
todo.push(target / divisor);
} else {
failed.push(target);
}
}
}
}
};
if failed.is_empty() {
(result, None)
} else {
(result, Some(failed))
}
}
fn factorize<T: PrimalityBase>(&self, target: T) -> BTreeMap<T, usize>
where
for<'r> &'r T: PrimalityRefBase<T>,
{
loop {
let (result, remainder) = self.factors(target.clone(), None);
if remainder.is_none() {
break result;
}
}
}
fn divisor<T: PrimalityBase>(&self, target: &T, config: &mut FactorizationConfig) -> Option<T>
where
for<'r> &'r T: PrimalityRefBase<T>,
{
if matches!(config.td_limit, Some(0)) {
let tsqrt: T = Roots::sqrt(target) + T::one();
let limit = if let Some(l) = config.td_limit {
tsqrt.clone().min(T::from_u64(l).unwrap())
} else {
tsqrt.clone()
};
for p in self.iter().map(|p| T::from_u64(*p).unwrap()) {
if &p > &tsqrt {
return None; }
if &p > &limit {
break;
}
if target.is_multiple_of(&p) {
return Some(p);
}
}
}
let below64 = target.to_u64().is_some();
while config.rho_trials > 0 {
let (start, offset) = if below64 {
(
T::from_u8(random::<u8>()).unwrap() % target,
T::from_u8(random::<u8>()).unwrap() % target,
)
} else {
(
T::from_u64(random::<u64>()).unwrap() % target,
T::from_u64(random::<u64>()).unwrap() % target,
)
};
config.rho_trials -= 1;
if let (Some(p), _) = pollard_rho(target, start, offset, 1048576) {
return Some(p);
}
}
None
}
}
impl<T> PrimeBufferExt for T where for<'a> T: PrimeBuffer<'a> {}
pub struct NaiveBuffer {
list: Vec<u64>, next: u64, }
impl NaiveBuffer {
#[inline]
pub fn new() -> Self {
let list = SMALL_PRIMES.iter().map(|&p| p as u64).collect();
NaiveBuffer {
list,
next: SMALL_PRIMES_NEXT,
}
}
}
impl<'a> PrimeBuffer<'a> for NaiveBuffer {
type PrimeIter = std::slice::Iter<'a, u64>;
fn contains(&self, num: u64) -> bool {
self.list.binary_search(&num).is_ok()
}
fn clear(&mut self) {
self.list.truncate(16);
self.list.shrink_to_fit();
self.next = 55; }
fn iter(&'a self) -> Self::PrimeIter {
self.list.iter()
}
fn bound(&self) -> u64 {
*self.list.last().unwrap()
}
fn reserve(&mut self, limit: u64) {
let sieve_limit = (limit | 1) + 2; let current = self.next; debug_assert!(current % 2 == 1);
if sieve_limit < current {
return;
}
let mut sieve = bitvec![usize, Msb0; 0; ((sieve_limit - current) / 2) as usize];
for p in self.list.iter().skip(1) {
let start = if p * p < current {
p * ((current / p) | 1) } else {
p * p
};
for multi in (start..sieve_limit).step_by(2 * (*p as usize)) {
if multi >= current {
sieve.set(((multi - current) / 2) as usize, true);
}
}
}
for p in (current..Roots::sqrt(&sieve_limit) + 1).step_by(2) {
for multi in (p * p..sieve_limit).step_by(2 * (p as usize)) {
if multi >= current {
sieve.set(((multi - current) / 2) as usize, true);
}
}
}
self.list
.extend(sieve.iter_zeros().map(|x| (x as u64) * 2 + current));
self.next = sieve_limit;
}
}
impl NaiveBuffer {
pub fn primorial<T: PrimalityBase + std::iter::Product>(&mut self, n: usize) -> T {
self.nprimes(n).map(|&p| T::from_u64(p).unwrap()).product()
}
pub fn primes(&mut self, limit: u64) -> std::iter::Take<<Self as PrimeBuffer>::PrimeIter> {
self.reserve(limit);
let position = match self.list.binary_search(&limit) {
Ok(p) => p + 1,
Err(p) => p,
}; return self.list.iter().take(position);
}
pub fn into_primes(mut self, limit: u64) -> std::vec::IntoIter<u64> {
self.reserve(limit);
let position = match self.list.binary_search(&limit) {
Ok(p) => p + 1,
Err(p) => p,
}; self.list.truncate(position);
return self.list.into_iter();
}
pub fn nprimes(&mut self, count: usize) -> std::iter::Take<<Self as PrimeBuffer>::PrimeIter> {
let (_, bound) = nth_prime_bounds(&(count as u64))
.expect("Estimated size of the largest prime will be larger than u64 limit");
self.reserve(bound);
debug_assert!(self.list.len() >= count);
self.list.iter().take(count)
}
pub fn into_nprimes(mut self, count: usize) -> std::vec::IntoIter<u64> {
let (_, bound) = nth_prime_bounds(&(count as u64))
.expect("Estimated size of the largest prime will be larger than u64 limit");
self.reserve(bound);
debug_assert!(self.list.len() >= count);
self.list.truncate(count);
return self.list.into_iter();
}
pub fn nth_prime(&mut self, n: u64) -> u64 {
if n < self.list.len() as u64 {
return self.list[n as usize - 1];
}
const THRESHOLD_NTH_PRIME_SIEVE: u64 = 4096;
if n <= THRESHOLD_NTH_PRIME_SIEVE {
return *self.nprimes(n as usize).last().unwrap();
}
let mut x = prev_prime(&nth_prime_est(&n).unwrap(), None).unwrap();
let mut pi = self.prime_pi(x);
while pi > n {
x = prev_prime(&x, None).unwrap();
pi -= 1;
}
while pi < n {
x = next_prime(&x, None).unwrap();
pi += 1;
}
x
}
pub fn prime_phi(&mut self, x: u64, a: usize, cache: &mut LruCache<(u64, usize), u64>) -> u64 {
if a == 1 {
return (x + 1) / 2;
}
if let Some(v) = cache.get(&(x, a)) {
return *v;
}
let t1 = self.prime_phi(x, a - 1, cache);
let pa = self.nth_prime(a as u64);
let t2 = self.prime_phi(x / pa, a - 1, cache);
let t = t1 - t2;
cache.put((x, a), t);
t
}
pub fn prime_pi(&mut self, limit: u64) -> u64 {
const THRESHOLD_PRIME_PI_SIEVE: u64 = 38873; if &limit <= self.list.last().unwrap() || limit <= THRESHOLD_PRIME_PI_SIEVE {
self.reserve(limit);
return match self.list.binary_search(&limit) {
Ok(pos) => (pos + 1) as u64,
Err(pos) => pos as u64,
};
}
let b = limit.sqrt();
let a = b.sqrt();
let c = limit.cbrt();
self.reserve(b);
let a = self.prime_pi(a);
let b = self.prime_pi(b);
let c = self.prime_pi(c);
let cache_cap = NonZeroUsize::new(a as usize).unwrap(); let mut phi_cache = LruCache::new(cache_cap);
let mut sum =
self.prime_phi(limit, a as usize, &mut phi_cache) + (b + a - 2) * (b - a + 1) / 2;
for i in a + 1..b + 1 {
let w = limit / self.nth_prime(i);
sum -= self.prime_pi(w);
if i <= c {
let l = self.prime_pi(w.sqrt());
sum += (l * (l - 1) - i * (i - 3)) / 2 - 1;
for j in i..(l + 1) {
let pj = self.nth_prime(j);
sum -= self.prime_pi(w / pj);
}
}
}
return sum;
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::mint::SmallMint;
#[cfg(feature = "num-bigint")]
use core::str::FromStr;
#[cfg(feature = "num-bigint")]
use num_bigint::BigUint;
use rand::random;
#[test]
fn prime_generation_test() {
const PRIME50: [u64; 15] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47];
const PRIME300: [u64; 62] = [
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271,
277, 281, 283, 293,
];
let mut pb = NaiveBuffer::new();
assert_eq!(pb.primes(50).cloned().collect::<Vec<_>>(), PRIME50);
assert_eq!(pb.primes(300).cloned().collect::<Vec<_>>(), PRIME300);
pb.clear();
assert_eq!(pb.primes(293).cloned().collect::<Vec<_>>(), PRIME300);
pb = NaiveBuffer::new();
assert_eq!(*pb.primes(257).last().unwrap(), 257); pb = NaiveBuffer::new();
assert_eq!(*pb.primes(8167).last().unwrap(), 8167); }
#[test]
fn nth_prime_test() {
let mut pb = NaiveBuffer::new();
assert_eq!(pb.nth_prime(10000), 104729);
assert_eq!(pb.nth_prime(20000), 224737);
assert_eq!(pb.nth_prime(10000), 104729); assert_eq!(pb.nth_prime(10u64.pow(4)), 104729);
assert_eq!(pb.nth_prime(10u64.pow(5)), 1299709);
assert_eq!(pb.nth_prime(10u64.pow(6)), 15485863);
assert_eq!(pb.nth_prime(10u64.pow(7)), 179424673);
}
#[test]
fn prime_pi_test() {
let mut pb = NaiveBuffer::new();
assert_eq!(pb.prime_pi(8161), 1024); assert_eq!(pb.prime_pi(10000), 1229);
assert_eq!(pb.prime_pi(20000), 2262);
assert_eq!(pb.prime_pi(38873), 4096);
pb.clear(); assert_eq!(pb.prime_pi(8161), 1024);
assert_eq!(pb.prime_pi(10000), 1229);
assert_eq!(pb.prime_pi(20000), 2262);
assert_eq!(pb.prime_pi(10000), 1229); assert_eq!(pb.prime_pi(10u64.pow(5)), 9592);
assert_eq!(pb.prime_pi(10u64.pow(6)), 78498);
assert_eq!(pb.prime_pi(10u64.pow(7)), 664579);
assert_eq!(pb.prime_pi(10u64.pow(8)), 5761455);
}
#[test]
fn is_prime_test() {
let pb = NaiveBuffer::new();
assert_eq!(pb.is_prime(&(2u32.pow(19) - 1), None), Primality::Yes);
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(&SmallMint::from(2u128.pow(89) - 1), None),
Primality::Probable(_)
));
for _ in 0..100 {
let target = random::<u64>();
assert_eq!(
!is_prime64(target),
matches!(
pb.is_prime(&target, Some(PrimalityTestConfig::bpsw())),
Primality::No
)
);
}
const P: u128 = 18699199384836356663; 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(&SmallMint::from(P), None),
Primality::Probable(_)
));
assert!(matches!(
pb.is_prime(&SmallMint::from(P), Some(PrimalityTestConfig::bpsw())),
Primality::Probable(_)
));
const P2: u128 = 2019922777445599503530083;
assert!(matches!(pb.is_prime(&P2, None), Primality::Probable(_)));
assert!(matches!(
pb.is_prime(&P2, Some(PrimalityTestConfig::bpsw())),
Primality::Probable(_)
));
assert!(matches!(
pb.is_prime(&SmallMint::from(P2), None),
Primality::Probable(_)
));
assert!(matches!(
pb.is_prime(&SmallMint::from(P2), Some(PrimalityTestConfig::bpsw())),
Primality::Probable(_)
));
#[cfg(feature = "num-bigint")]
{
let large_primes = [
"98920366548084643601728869055592650835572950932266967461790948584315647051443",
"94560208308847015747498523884063394671606671904944666360068158221458669711639",
"449417999055441493994709297093108513015373787049558499205492347871729927573118262811508386655998299074566974373711472560655026288668094291699357843464363003144674940345912431129144354948751003607115263071543163",
"230975859993204150666423538988557839555560243929065415434980904258310530753006723857139742334640122533598517597674807096648905501653461687601339782814316124971547968912893214002992086353183070342498989426570593",
"5521712099665906221540423207019333379125265462121169655563495403888449493493629943498064604536961775110765377745550377067893607246020694972959780839151452457728855382113555867743022746090187341871655890805971735385789993",
"203956878356401977405765866929034577280193993314348263094772646453283062722701277632936616063144088173312372882677123879538709400158306567338328279154499698366071906766440037074217117805690872792848149112022286332144876183376326512083574821647933992961249917319836219304274280243803104015000563790123",
"3618502788666131106986593281521497120414687020801267626233049500247285301239", "57896044618658097711785492504343953926634992332820282019728792003956564819949", "9850501549098619803069760025035903451269934817616361666987073351061430442874302652853566563721228910201656997576599", "42307582002575910332922579714097346549017899709713998034217522897561970639123926132812109468141778230245837569601494931472367", "6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151", "169511182982703321453314585423962898651587669459838234386506572286328885534468792292646838949809616446341407457141008401355628947670484184607678853094537849610289912805960069455687743151708433319901176932959509872662610091644590437761688516626993416011399330087939042347256922771590903190536793274742859624657"
];
for pstr in large_primes {
assert!(
pb.is_prime(&BigUint::from_str(pstr).unwrap(), None)
.probably(),
"false negative prime {}",
pstr
)
}
let large_composites = [
"21284175091214687912771199898307297748211672914763848041968395774954376176754",
"6084766654921918907427900243509372380954290099172559290432744450051395395951",
"84594350493221918389213352992032324280367711247940675652888030554255915464401",
"82793403787388584738507275144194252681",
"1195068768795265792518361315725116351898245581", "8038374574536394912570796143419421081388376882875581458374889175222974273765333652186502336163960045457915042023603208766569966760987284043965408232928738791850869166857328267761771029389697739470167082304286871099974399765441448453411558724506334092790222752962294149842306881685404326457534018329786111298960644845216191652872597534901",
];
for cstr in large_composites {
assert!(
!pb.is_prime(
&BigUint::from_str(cstr).unwrap(),
Some(PrimalityTestConfig::bpsw())
)
.probably(),
"false positive prime {}",
cstr
)
}
}
}
#[test]
fn pb_factors_test() {
let pb = NaiveBuffer::new();
#[cfg(feature = "num-bigint")]
{
let m131 = BigUint::from(2u8).pow(131) - 1u8; let (fac, r) = pb.factors(m131, None);
assert!(fac.len() == 2 && r.is_none());
}
}
}