diff --git a/.gitignore b/.gitignore index 96ef6c0..2dcd902 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,5 @@ /target Cargo.lock + +# generated profile stats +profile_stats.csv diff --git a/examples/profile_factorization.py b/examples/profile_factorization.py new file mode 100644 index 0000000..cd32013 --- /dev/null +++ b/examples/profile_factorization.py @@ -0,0 +1,30 @@ +# This script plot the result generated by profile_factorization.rs +import pandas as pd +import numpy as np +from matplotlib import pyplot as plt + +table = pd.read_csv("profile_stats.csv") +table.drop(columns=["n"], inplace=True) +pollard_cols = list(k for k in table.columns if k.startswith("pollard")) +squfof_cols = list(k for k in table.columns if k.startswith("squfof")) +oneline_cols = list(k for k in table.columns if k.startswith("one_line")) + +# MAXITER = 1 << 20 +# table[table >= MAXITER] = np.nan + +mean_table = table.groupby(table['n_bits'] // 4).agg(np.nanmean) +fig, ax = plt.subplots() +mean_table.plot("n_bits", pollard_cols, ax=ax) +mean_table.plot("n_bits", squfof_cols, ax=ax) +mean_table.plot("n_bits", oneline_cols, ax=ax) +ax.set_yscale("log") + +min_table = table.groupby(table['n_bits'] // 4).agg(np.nanmin) +fig, ax = plt.subplots() +ax.plot(min_table["n_bits"], np.mean(min_table[pollard_cols], axis=1), label="pollard") +ax.plot(min_table["n_bits"], np.mean(min_table[squfof_cols], axis=1), label="squfof") +ax.plot(min_table["n_bits"], np.mean(min_table[oneline_cols], axis=1), label="one_line") +ax.legend() +ax.set_yscale("log") + +plt.show() diff --git a/examples/profile_factorization.rs b/examples/profile_factorization.rs index f048214..4738734 100644 --- a/examples/profile_factorization.rs +++ b/examples/profile_factorization.rs @@ -1,31 +1,87 @@ +use std::fs::File; +use std::io::{Write, Error}; + use num_prime::factor::{pollard_rho, squfof, one_line, SQUFOF_MULTIPLIERS}; use num_prime::RandPrime; use rand::random; -// TODO: create a plot for iterations needed for factoring random numbers +fn profile_n(n: u128) -> Vec::<(String, usize)> { + let k_squfof: Vec = SQUFOF_MULTIPLIERS.iter().take(10).cloned().collect(); + let k_oneline: Vec = vec![1, 360, 480]; + const MAXITER: usize = 1 << 20; -fn main() { - let mut rng = rand::thread_rng(); + let mut n_stats = Vec::new(); + + // pollard rho + n_stats.push(("pollard_rho1".to_string(), pollard_rho(&n, random(), random(), MAXITER).1)); + n_stats.push(("pollard_rho2".to_string(), pollard_rho(&n, random(), random(), MAXITER).1)); - // let p1: u64 = rng.gen_prime(40, None); - // let p2: u64 = rng.gen_prime(60, None); - - let p1: u64 = 3486784447; - let p2: u64 = 94143178889; - // let p2: u64 = 3486784709; - - println!("Semiprime: {} * {}", p1, p2); - let n = p1 as u128 * p2 as u128; - - // let n: u128 = 133717415095455410877609739380293; - - const MAXITER: usize = 2 << 20; - for &k in SQUFOF_MULTIPLIERS.iter().take(10) { + // squfof + for &k in &k_squfof { + let key = format!("squfof_k{}", k); if let Some(kn) = n.checked_mul(k as u128) { - println!("squfof k={} result: {:?}", k, squfof(&n, kn, MAXITER)); + let n = squfof(&n, kn, MAXITER).1; + n_stats.push((key, n)); + } else { + n_stats.push((key, MAXITER)); + }; + } + + // one line + for &k in &k_oneline { + let key = format!("one_line_k{}", k); + if let Some(kn) = n.checked_mul(k as u128) { + let n = one_line(&n, kn, MAXITER).1; + n_stats.push((key, n)); + } else { + n_stats.push((key, MAXITER)); + }; + } + + n_stats +} + +/// This program try various factorization methods, and log down their iterations number into a csv file +fn main() -> Result<(), Error> { + let mut rng = rand::thread_rng(); + const REPEATS: u32 = 4; + + let mut n_list = Vec::<(u128, f32)>::new(); // n and bits of n + let mut stats: Vec> = Vec::new(); + + for total_bits in 10..80 { + for _ in 0..REPEATS { + let p1: u128 = rng.gen_prime(total_bits / 2, None); + let p2: u128 = rng.gen_prime_exact(total_bits - (128 - p1.leading_zeros()) as usize, None); + if p1 == p2 { + continue; + } + + let n = p1 * p2; + n_list.push((n, (n as f64).log2() as f32)); + println!("Semiprime: {} = {} * {}", n, p1, p2); + stats.push(profile_n(n)); } } - // println!("one_line k=1 result: {:?}", one_line(&n, n, MAXITER)); - // println!("one_line k=480 result: {:?}", one_line(&n, n * 480, MAXITER)); - // println!("pollard rho result: {:?}", pollard_rho(&n, random(), random(), MAXITER)); + + // Log into the CSV file + let mut fout = File::create("profile_stats.csv")?; + fout.write(b"n,n_bits")?; + for k in stats[0].iter().map(|(k, _)| k) { + fout.write(b",")?; + fout.write(k.as_bytes())?; + } + + for ((n, bits), n_stats) in n_list.iter().zip(stats) { + fout.write(b"\n")?; + fout.write(n.to_string().as_bytes())?; + fout.write(b",")?; + fout.write(bits.to_string().as_bytes())?; + for (_, v) in n_stats { + fout.write(b",")?; + fout.write(v.to_string().as_bytes())?; + } + } + + Ok(()) } diff --git a/src/nt_funcs.rs b/src/nt_funcs.rs index 9d03d62..1f76825 100644 --- a/src/nt_funcs.rs +++ b/src/nt_funcs.rs @@ -21,7 +21,7 @@ use crate::tables::{ #[cfg(feature = "big-table")] use crate::tables::{SMALL_PRIMES_INV, ZETA_LOG_TABLE}; use crate::traits::{FactorizationConfig, Primality, PrimalityTestConfig, PrimalityUtils}; -use crate::ExactRoots; +use crate::{ExactRoots, BitTest}; use num_integer::Roots; #[cfg(feature = "num-bigint")] use num_modular::DivExact; @@ -161,6 +161,10 @@ pub fn factorize64(target: u64) -> BTreeMap { // https://github.com/elmomoilanen/prime-factorization // https://github.com/radii/msieve // Pari/GP: ifac_crack + // TODO(v0.next): check the runtime of each factorization and put the fastest first + // TODO(v0.next): add multipliers for one_line method + // TODO(v0.next): quickly increase the limit for squfof, try to match the behavior of gnu factor + // TODO(v0.next): make the factorization method resumable? let mut result = BTreeMap::new(); // quick check on factors of 2 @@ -265,7 +269,7 @@ pub(crate) fn factorize64_advanced(cofactors: &[(u64, usize)]) -> Vec<(u64, usiz // try to find a divisor let mut i = 0usize; - let mut max_iter = 2 << 16; + let mut max_iter = 2 << (target.bits() / 4); // empirical lower bound for iterations let divisor = loop { // try various factorization method iteratively const NMETHODS: usize = 3; @@ -286,11 +290,18 @@ pub(crate) fn factorize64_advanced(cofactors: &[(u64, usize)]) -> Vec<(u64, usiz } } 2 => { - // Shank's squfof - if let Some(mul_target) = target.checked_mul(SQUFOF_MULTIPLIERS[i % SQUFOF_MULTIPLIERS.len()] as u64) { - if let (Some(p), _) = squfof(&target, mul_target, max_iter) { - break p; + // Shanks's squfof + let mut d = None; + for &k in SQUFOF_MULTIPLIERS.iter() { + if let Some(mul_target) = target.checked_mul(k as u64) { + if let (Some(p), _) = squfof(&target, mul_target, max_iter) { + d = Some(p); + break; + } } + }; + if let Some(p) = d { + break p; } } _ => unreachable!(), @@ -410,31 +421,40 @@ pub(crate) fn factorize128_advanced(cofactors: &[(u128, usize)]) -> Vec<(u128, u // try to find a divisor let mut i = 0usize; - let mut max_iter = 2 << 18; // allow more iterations than u64 + let mut max_iter = 2 << (target.bits() / 6); // empirical lower bound let divisor = loop { // try various factorization method iteratively const NMETHODS: usize = 3; match i % NMETHODS { 0 => { - // Pollard's rho - let start = MontgomeryInt::new(random::(), target); - let offset = start.convert(random::()); - if let (Some(p), _) = pollard_rho(&Mint::from(target), start.into(), offset.into(), max_iter) { - break p.value(); - } - } - 1 => { // Hart's one-line let mul_target = target.checked_mul(480).unwrap_or(target); if let (Some(p), _) = one_line(&target, mul_target, max_iter) { break p; } } + 1 => { + // Shanks's squfof, try all mutipliers + let mut d = None; + for &k in SQUFOF_MULTIPLIERS.iter() { + if let Some(mul_target) = target.checked_mul(k as u128) { + if let (Some(p), _) = squfof(&target, mul_target, max_iter) { + d = Some(p); + break; + } + } + }; + if let Some(p) = d { + break p; + } + } 2 => { - // Shanks's squfof - if let Some(mul_target) = target.checked_mul(SQUFOF_MULTIPLIERS[i % NMETHODS] as u128) { - if let (Some(p), _) = squfof(&target, mul_target, max_iter) { - break p; + // Pollard's rho, only twice + if i / NMETHODS < 2 { + let start = MontgomeryInt::new(random::(), target); + let offset = start.convert(random::()); + if let (Some(p), _) = pollard_rho(&Mint::from(target), start.into(), offset.into(), max_iter) { + break p.value(); } } }