uutils/num-prime
0
0
Fork 0
mirror of https://github.com/uutils/num-prime.git synced 2026-06-10 16:12:35 -07:00
Code Issues Packages Projects Releases Wiki Activity

Some refactor

Browse Source
This commit is contained in:
Jacob Zhong
2022-04-27 11:30:01 -04:00
parent 0f0fa88c44
commit 2c5d8aa6d5
4 changed files with 115 additions and 60 deletions
Download Patch File Download Diff File Expand all files Collapse all files
README.md
+1 -1
View File
@@ -17,7 +17,7 @@ This crate provides utilities for prime number related functionalities:
- Trial division
- Pollard's rho algorithm
- Shanks's square forms factorization (SQUFOF)
- Fast factorization of `u64` integers
- Fast factorization of `u64` and `u128` integers
- Number theoretic functions
- Prime Pi function (number of primes under limit), its estimation and its bounds
- Nth prime, its estimation and its bounds
examples/divisor_sigma.rs
+1 -1
View File
@@ -1,5 +1,5 @@
// https://en.wikipedia.org/wiki/Divisor_function
// TODO(v0.3.2): implement divisor sigma as example for factorization
// TODO(v0.3.3): implement divisor sigma as example for factorization
fn main() {
println!("nothing here")
src/factor.rs
+33 -3
View File
@@ -1,4 +1,7 @@
//! Implementations for various factorization algorithms
//!
//! See <https://web.archive.org/web/20110331180514/https://diamond.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf>
//! for a detailed comparison between different factorization algorithms
use crate::traits::ExactRoots;
use num_integer::{Integer, Roots};
@@ -106,8 +109,11 @@ 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]`
/// The multiplier can be choosen from SQUFOF_MULTIPLIERS, or other square-free odd numbers.
///
/// 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<T: Integer + NumRef + Clone + ExactRoots>(target: &T, multiplier: T) -> Option<T>
where
for<'r> &'r T: RefNum<T>,
@@ -181,7 +187,31 @@ where
}
}
// TODO: ECM, One line, Quadratic sieve / Prime field sieve, Fermat(https://en.wikipedia.org/wiki/Fermat%27s_factorization_method)
// Square-free even numbers are suitable as SQUFOF multipliers
pub const SQUFOF_MULTIPLIERS: [u16; 16] = [
1,
3,
5,
7,
11,
3 * 5,
3 * 7,
3 * 11,
5 * 7,
5 * 11,
7 * 11,
3 * 5 * 7,
3 * 5 * 11,
3 * 7 * 11,
5 * 7 * 11,
3 * 5 * 7 * 11,
];
// TODO(v0.3.3): implement one line factorization and its optimization
// REF: doi:10.1017/S1446788712000146
// https://math.mit.edu/research/highschool/primes/materials/2019/Gopalakrishna.pdf
// TODO: ECM, 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/
src/nt_funcs.rs
+80 -55
View File
@@ -13,7 +13,7 @@
use crate::buffer::{NaiveBuffer, PrimeBufferExt};
use crate::mint::Mint;
use crate::factor::{pollard_rho, squfof};
use crate::factor::{pollard_rho, squfof, SQUFOF_MULTIPLIERS};
use crate::primality::{PrimalityBase, PrimalityRefBase};
use crate::tables::{MOEBIUS_ODD, SMALL_PRIMES, SMALL_PRIMES_NEXT, WHEEL_NEXT, WHEEL_PREV, WHEEL_SIZE};
#[cfg(feature = "big-table")]
@@ -29,7 +29,7 @@ use num_modular::DivExact;
use num_traits::{CheckedAdd, FromPrimitive, Num, RefNum, ToPrimitive};
use rand::{random, Rng};
use std::collections::BTreeMap;
use std::convert::{TryFrom, TryInto};
use std::convert::TryFrom;
#[cfg(feature = "big-table")]
use crate::tables::{MILLER_RABIN_BASE32, MILLER_RABIN_BASE64};
@@ -146,8 +146,8 @@ fn is_prime64_miller(target: u64) -> bool {
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.
/// if target is larger than 2^64 or more controlled primality tests are desired, please use [factors()].
/// Fast integer factorization on a u64 target. It's based on a selection of factorization methods.
/// if target is larger than 2^128 or more controlled primality tests are desired, please use [factors()].
///
/// The factorization can be quite faster under 2^64 because: 1) faster and deterministic primality check,
/// 2) efficient montgomery multiplication implementation of u64
@@ -234,42 +234,24 @@ pub fn factorize64(target: u64) -> BTreeMap<u64, usize> {
}
// then try advanced methods to find a divisor util fully factored
for (p, exp) in factorize64_advance(vec![(residual, 1usize)]).into_iter() {
for (p, exp) in factorize64_advanced(&[(residual, 1usize)]).into_iter() {
*result.entry(p).or_insert(0) += exp;
}
result
}
// This function factorize all cofactors after some trivial division steps
pub(crate) fn factorize64_advance(cofactors: Vec<(u64, usize)>) -> Vec<(u64, usize)> {
let mut todo = cofactors;
pub(crate) fn factorize64_advanced(cofactors: &[(u64, usize)]) -> Vec<(u64, usize)> {
let mut todo: Vec<_> = cofactors.iter().cloned().collect();
let mut factored: Vec<(u64, usize)> = Vec::new(); // prime factor, exponent
const SQUFOF_MULTIPLIERS: [u16; 16] = [
1,
3,
5,
7,
11,
3 * 5,
3 * 7,
3 * 11,
5 * 7,
5 * 11,
7 * 11,
3 * 5 * 7,
3 * 5 * 11,
3 * 7 * 11,
5 * 7 * 11,
3 * 5 * 7 * 11,
];
while let Some((target, exp)) = todo.pop() {
if is_prime64_miller(target) {
factored.push((target, exp));
continue;
}
// check perfect powers before other methods
// check perfect powers before other methods, this is required for SQUFOF
// it suffices to check square and cubic if big-table is enabled, since fifth power of
// the smallest prime that haven't been checked is 8167^5 > 2^64
if let Some(d) = target.sqrt_exact() {
@@ -282,19 +264,29 @@ pub(crate) fn factorize64_advance(cofactors: Vec<(u64, usize)>) -> Vec<(u64, usi
}
// try to find a divisor
let mut i = 1usize;
let mut i = 0usize;
let divisor = loop {
// try SQUFOF after 4 failed pollard rho trials
if i % 5 == 0 && (i / 5) < SQUFOF_MULTIPLIERS.len() {
if let Some(p) = squfof(&target, SQUFOF_MULTIPLIERS[i / 5] as u64) {
break p;
}
} else {
let start = MontgomeryInt::new(random::<u64>(), target);
let offset = start.convert(random::<u64>());
if let Some(p) = pollard_rho(&Mint::from(target), start.into(), offset.into()) {
break p.value();
}
// try various factorization method iteratively
const NMETHODS: usize = 2;
match i % NMETHODS {
0 => {
// pollard rho
let start = MontgomeryInt::new(random::<u64>(), target);
let offset = start.convert(random::<u64>());
if let Some(p) = pollard_rho(&Mint::from(target), start.into(), offset.into()) {
break p.value();
}
},
1 => {
// squfof
if i / NMETHODS >= SQUFOF_MULTIPLIERS.len() {
continue;
}
if let Some(p) = squfof(&target, SQUFOF_MULTIPLIERS[i / NMETHODS] as u64) {
break p;
}
},
_ => unreachable!()
}
i += 1;
};
@@ -304,6 +296,9 @@ pub(crate) fn factorize64_advance(cofactors: Vec<(u64, usize)>) -> Vec<(u64, usi
factored
}
/// Fast integer factorization on a u128 target. It's based on a selection of factorization methods.
/// if target is larger than 2^128 or more controlled primality tests are desired, please use [factors()].
// TODO(v0.next): return BTreeMap<u64, u8> instead of BTreeMap<u64, usize>
pub fn factorize128(target: u128) -> BTreeMap<u128, usize> {
// shortcut for u64
if target < (1u128 << 64) {
@@ -323,7 +318,7 @@ pub fn factorize128(target: u128) -> BTreeMap<u128, usize> {
let mut residual = target >> f2;
// trial division using primes in the table
// note that p^2 is never larger than target, so we don't need to shortcut trial division
// note that p^2 is never larger than target (at least 64 bits), so we don't need to shortcut trial division
#[cfg(not(feature = "big-table"))]
for p in SMALL_PRIMES.iter().skip(1).map(|&v| v as u128) {
while residual % p == 0 {
@@ -357,15 +352,26 @@ pub fn factorize128(target: u128) -> BTreeMap<u128, usize> {
}
// then try advanced methods to find a divisor util fully factored
let (mut todo128, mut todo64) = if let Ok(r64) = u64::try_from(residual) {
(Vec::new(), vec![(r64, 1usize)])
} else {
(vec![(residual, 1usize)], Vec::new())
}; // cofactors to be processed
for (p, exp) in factorize128_advanced(&[(residual, 1usize)]).into_iter() {
*result.entry(p).or_insert(0) += exp;
}
result
}
pub(crate) fn factorize128_advanced(cofactors: &[(u128, usize)]) -> Vec<(u128, usize)> {
let (mut todo128, mut todo64) = (Vec::new(), Vec::new()); // cofactors to be processed
let mut factored: Vec<(u128, usize)> = Vec::new(); // prime factor, exponent
for &(co, e) in cofactors.iter() {
if let Ok(co64) = u64::try_from(co) {
todo64.push((co64, e));
} else {
todo128.push((co, e));
};
}
while let Some((target, exp)) = todo128.pop() {
if is_prime(&Mint::from(target), None).probably() {
*result.entry(target).or_insert(0) += 1;
if is_prime(&Mint::from(target), Some(PrimalityTestConfig::bpsw())).probably() {
factored.push((target, exp));
continue;
}
@@ -388,14 +394,34 @@ pub fn factorize128(target: u128) -> BTreeMap<u128, usize> {
}
continue;
}
// TODO: check 5-th, 7-th power
// try to find a divisor
let mut i = 0usize;
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();
// try various factorization method iteratively
const NMETHODS: usize = 2;
match i % NMETHODS {
0 => {
// pollard rho
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();
}
},
1 => {
// squfof
if i / NMETHODS >= SQUFOF_MULTIPLIERS.len() {
continue;
}
if let Some(p) = squfof(&target, SQUFOF_MULTIPLIERS[i / NMETHODS] as u128) {
break p;
}
},
_ => unreachable!()
}
i += 1;
};
if let Ok(d64) = u64::try_from(divisor) {
@@ -411,10 +437,9 @@ pub fn factorize128(target: u128) -> BTreeMap<u128, usize> {
}
}
for (p, exp) in factorize64_advance(todo64).into_iter() {
*result.entry(p as u128).or_insert(0) += exp;
}
result
// forward 64 bit cofactors
factored.extend(factorize64_advanced(&todo64).into_iter().map(|(p, exp)| (p as u128, exp)));
factored
}
/// This function re-exports [PrimeBufferExt::is_prime()][crate::buffer::PrimeBufferExt::is_prime()] with a default buffer distance