From bbad5f5dbf3cf264f706f2be1c4b2e3ed8f5e00a Mon Sep 17 00:00:00 2001 From: Jacob Zhong Date: Sun, 1 May 2022 06:54:49 -0400 Subject: [PATCH] Change bound for factorization methods --- examples/profile_factorization.py | 67 ++++++++++++++++++++---------- examples/profile_factorization.rs | 69 ++++++++++++++++++++++++++++--- src/factor.rs | 7 ++++ src/nt_funcs.rs | 50 +++++++++++----------- 4 files changed, 142 insertions(+), 51 deletions(-) diff --git a/examples/profile_factorization.py b/examples/profile_factorization.py index cd32013..60b8b5d 100644 --- a/examples/profile_factorization.py +++ b/examples/profile_factorization.py @@ -1,30 +1,55 @@ # This script plot the result generated by profile_factorization.rs +from re import A 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")) +def plot_n_stats(): + 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 + # 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") + 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") + 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() +def plot_n_min_stats(): + table = pd.read_csv("profile_stats.csv") + table.drop(columns=["n"], inplace=True) + for k in table.columns: # caculate average time + if k.startswith("time_"): + table[k] = table[k] / table[k[5:]] + print(table[k]) + min_table = table.groupby(table['n_bits'] // 4).agg(np.nanmean) + + # MAXITER = 1 << 24 + # table[table >= MAXITER] = np.nan + + ax = min_table.plot("n_bits", ["pollard_rho", "squfof", "one_line"]) + ax.set_yscale("log") + ax.set_ylabel("min iters") + + ax = min_table.plot("n_bits", ["time_pollard_rho", "time_squfof", "time_one_line"]) + ax.set_yscale("log") + ax.set_ylabel("avg time per iter") + +if __name__ == "__main__": + # plot_n_stats() + plot_n_min_stats() + plt.show() diff --git a/examples/profile_factorization.rs b/examples/profile_factorization.rs index 4738734..6a1e6bc 100644 --- a/examples/profile_factorization.rs +++ b/examples/profile_factorization.rs @@ -1,20 +1,24 @@ use std::fs::File; use std::io::{Write, Error}; +use std::time::{Duration, Instant}; use num_prime::factor::{pollard_rho, squfof, one_line, SQUFOF_MULTIPLIERS}; use num_prime::RandPrime; use rand::random; -fn profile_n(n: u128) -> Vec::<(String, usize)> { +/// Collect the the iteration number of each factorization algorithm with different settings +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; + const POLLARD_REPEATS: usize = 2; 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)); + for i in 0..POLLARD_REPEATS { + n_stats.push((format!("pollard_rho{}", i+1), pollard_rho(&n, random(), random(), MAXITER).1)); + } // squfof for &k in &k_squfof { @@ -41,6 +45,58 @@ fn profile_n(n: u128) -> Vec::<(String, usize)> { n_stats } +/// Collect the best case of each factorization algorithm +fn profile_n_min(n: u128) -> Vec::<(String, usize)> { + let k_squfof: Vec = SQUFOF_MULTIPLIERS.iter().cloned().collect(); + let k_oneline: Vec = vec![1, 360, 480]; + const MAXITER: usize = 1 << 24; + const POLLARD_REPEATS: usize = 4; + + let mut n_stats = Vec::new(); + + // pollard rho + let mut pollard_best = (MAXITER, u128::MAX); + for _ in 0..POLLARD_REPEATS { + let tstart = Instant::now(); + let (result, iters) = pollard_rho(&n, random(), random(), pollard_best.0); + if result.is_some() { + pollard_best = pollard_best.min((iters, tstart.elapsed().as_micros())); + } + } + n_stats.push(("pollard_rho".to_string(), pollard_best.0)); + n_stats.push(("time_pollard_rho".to_string(), pollard_best.1 as usize)); + + // squfof + let mut squfof_best = (MAXITER, u128::MAX); + for &k in &k_squfof { + if let Some(kn) = n.checked_mul(k as u128) { + let tstart = Instant::now(); + let (result, iters) = squfof(&n, kn, squfof_best.0); + if result.is_some() { + squfof_best = squfof_best.min((iters, tstart.elapsed().as_micros())); + } + } + } + n_stats.push(("squfof".to_string(), squfof_best.0)); + n_stats.push(("time_squfof".to_string(), squfof_best.1 as usize)); + + // one line + let mut oneline_best = (MAXITER, u128::MAX); + for &k in &k_oneline { + if let Some(kn) = n.checked_mul(k as u128) { + let tstart = Instant::now(); + let (result, iters) = one_line(&n, kn, oneline_best.0); + if result.is_some() { + oneline_best = oneline_best.min((iters, tstart.elapsed().as_micros())); + } + } + } + n_stats.push(("one_line".to_string(), oneline_best.0)); + n_stats.push(("time_one_line".to_string(), squfof_best.1 as usize)); + + 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(); @@ -49,7 +105,7 @@ fn main() -> Result<(), Error> { 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 total_bits in 20..120 { 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); @@ -59,8 +115,9 @@ fn main() -> Result<(), Error> { 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!("Semiprime ({}bits): {} = {} * {}", total_bits, n, p1, p2); + // stats.push(profile_n(n)); + stats.push(profile_n_min(n)); } } diff --git a/src/factor.rs b/src/factor.rs index 625302e..69fb0b6 100644 --- a/src/factor.rs +++ b/src/factor.rs @@ -5,6 +5,9 @@ //! See //! for a detailed comparison between different factorization algorithms + +// XXX: make the factorization method resumable? + use crate::traits::ExactRoots; use num_integer::{Integer, Roots}; use num_modular::{ModularCoreOps, ModularUnaryOps}; @@ -234,6 +237,7 @@ pub const SQUFOF_MULTIPLIERS: [u16; 38] = [ /// where p = next_prime(c^a+d1), p = next_prime(c^b+d2), a and b are close, and c, d1, d2 are small integers. /// /// Reference: Hart, W. B. (2012). A one line factoring algorithm. Journal of the Australian Mathematical Society, 92(1), 61-69. doi:10.1017/S1446788712000146 +// TODO: add multipliers preset for one_line method? pub fn one_line(target: &T, mul_target: T, max_iter: usize) -> (Option, usize) where for<'r> &'r T: RefNum, { @@ -250,6 +254,7 @@ where } } + // prevent overflow ikn = if let Some(n) = (&ikn).checked_add(&mul_target) { n } else { @@ -307,6 +312,8 @@ mod tests { // this case should success at step 276, from https://rosettacode.org/wiki/Talk:Square_form_factorization assert!(matches!(squfof(&4558849u32, 4558849u32, 300).0, Some(_))); + + // TODO(v0.next): add more cases from rosetta code } #[test] diff --git a/src/nt_funcs.rs b/src/nt_funcs.rs index 1f76825..09505bb 100644 --- a/src/nt_funcs.rs +++ b/src/nt_funcs.rs @@ -161,10 +161,6 @@ 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 @@ -269,31 +265,34 @@ 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 << (target.bits() / 4); // empirical lower bound for iterations + let mut max_iter_ratio = 1; // increase max_iter after factorization round let divisor = loop { // try various factorization method iteratively const NMETHODS: usize = 3; match i % NMETHODS { 0 => { - // Pollard's rho + // Pollard's rho (quick check) let start = MontgomeryInt::new(random::(), target); let offset = start.convert(random::()); + let max_iter = max_iter_ratio << (target.bits() / 6); // unoptimized heuristic if let (Some(p), _) = pollard_rho(&Mint::from(target), start.into(), offset.into(), max_iter) { break p.value(); } } 1 => { - // Hart's one-line + // Hart's one-line (quick check) let mul_target = target.checked_mul(480).unwrap_or(target); + let max_iter = max_iter_ratio << (mul_target.bits() / 6); // unoptimized heuristic if let (Some(p), _) = one_line(&target, mul_target, max_iter) { break p; } } 2 => { - // Shanks's squfof + // Shanks's squfof (main power) let mut d = None; for &k in SQUFOF_MULTIPLIERS.iter() { if let Some(mul_target) = target.checked_mul(k as u64) { + let max_iter = max_iter_ratio * 2 * (2 * mul_target.sqrt()).sqrt() as usize; if let (Some(p), _) = squfof(&target, mul_target, max_iter) { d = Some(p); break; @@ -310,7 +309,7 @@ pub(crate) fn factorize64_advanced(cofactors: &[(u64, usize)]) -> Vec<(u64, usiz // increase max iterations after trying all methods if i % NMETHODS == 0 { - max_iter *= 4; + max_iter_ratio *= 2; } }; todo.push((divisor, exp)); @@ -421,23 +420,36 @@ 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 << (target.bits() / 6); // empirical lower bound + let mut max_iter_ratio = 1; + let divisor = loop { - // try various factorization method iteratively + // try various factorization method iteratively, sort by time per iteration const NMETHODS: usize = 3; match i % NMETHODS { 0 => { + // Pollard's rho + let start = MontgomeryInt::new(random::(), target); + let offset = start.convert(random::()); + let max_iter = max_iter_ratio << (target.bits() / 6); // unoptimized heuristic + 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); + let max_iter = max_iter_ratio << (mul_target.bits() / 6); // unoptimized heuristic if let (Some(p), _) = one_line(&target, mul_target, max_iter) { break p; } } - 1 => { + 2 => { // 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(mul_target) = target.checked_mul(k as u128) { + // this bound is from GNU factor + let max_iter = 2*(2 * mul_target.sqrt()).sqrt() as usize; if let (Some(p), _) = squfof(&target, mul_target, max_iter) { d = Some(p); break; @@ -448,23 +460,13 @@ pub(crate) fn factorize128_advanced(cofactors: &[(u128, usize)]) -> Vec<(u128, u break p; } } - 2 => { - // 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(); - } - } - } _ => unreachable!(), } i += 1; // increase max iterations after trying all methods if i % NMETHODS == 0 { - max_iter *= 4; + max_iter_ratio *= 2; } };