diff --git a/examples/profile_factorization.rs b/examples/profile_factorization.rs
index 3947024..f048214 100644
--- a/examples/profile_factorization.rs
+++ b/examples/profile_factorization.rs
@@ -2,6 +2,8 @@ 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 main() {
     let mut rng = rand::thread_rng();
 
@@ -18,7 +20,7 @@ fn main() {
     // let n: u128 = 133717415095455410877609739380293;
 
     const MAXITER: usize = 2 << 20;
-    for k in SQUFOF_MULTIPLIERS {
+    for &k in SQUFOF_MULTIPLIERS.iter().take(10) {
         if let Some(kn) = n.checked_mul(k as u128) {
             println!("squfof k={} result: {:?}", k, squfof(&n, kn, MAXITER));
         }
diff --git a/src/factor.rs b/src/factor.rs
index efeacca..625302e 100644
--- a/src/factor.rs
+++ b/src/factor.rs
@@ -209,25 +209,17 @@ where
     (None, max_iter)
 }
 
-// Square-free even numbers are suitable as SQUFOF multipliers
-// TODO(v0.next): which multiplier is more efficient?
-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,
+/// Good squfof multipliers sorted by efficiency descendingly, from Dana Jacobsen.
+/// 
+/// Note: square-free odd numbers are suitable as SQUFOF multipliers
+pub const SQUFOF_MULTIPLIERS: [u16; 38] = [
+    3 * 5 * 7 * 11, 3 * 5 * 7, 3 * 5 * 7 * 11 * 13, 3 * 5 * 7 * 13, 3 * 5 * 7 * 11 * 17, 3 * 5 * 11,
+    3 * 5 * 7 * 17, 3 * 5, 3 * 5 * 7 * 11 * 19, 3 * 5 * 11 * 13, 3 * 5 * 7 * 19, 3 * 5 * 7 * 13 * 17,
+    3 * 5 * 13, 3 * 7 * 11, 3 * 7, 5 * 7 * 11, 3 * 7 * 13, 5 * 7,
+    3 * 5 * 17, 5 * 7 * 13, 3 * 5 * 19, 3 * 11, 3 * 7 * 17, 3,
+    3 * 11 * 13, 5 * 11, 3 * 7 * 19, 3 * 13, 5, 5 * 11 * 13,
+    5 * 7 * 19, 5 * 13, 7 * 11, 7, 3 * 17, 7 * 13,
+    11, 1
 ];
 
 /// William Hart's one line factorization algorithm for 64 bit integers.
diff --git a/src/nt_funcs.rs b/src/nt_funcs.rs
index a53dbbc..dc3d990 100644
--- a/src/nt_funcs.rs
+++ b/src/nt_funcs.rs
@@ -267,11 +267,10 @@ 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 divisor = loop {
             // try various factorization method iteratively
             const NMETHODS: usize = 3;
-            let mut max_iter = 2 << 16;
-
             match i % NMETHODS {
                 0 => {
                     // Pollard's rho
@@ -301,9 +300,8 @@ pub(crate) fn factorize64_advanced(cofactors: &[(u64, usize)]) -> Vec<(u64, usiz
             i += 1;
 
             // increase max iterations after trying all methods
-            #[allow(unused_assignments)]
             if i % NMETHODS == 0 {
-                max_iter *= 2;
+                max_iter *= 4;
             }
         };
         todo.push((divisor, exp));
@@ -414,11 +412,10 @@ 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 divisor = loop {
             // try various factorization method iteratively
             const NMETHODS: usize = 3;
-            let mut max_iter = 2 << 18; // allow more iterations than u64
-
             match i % NMETHODS {
                 0 => {
                     // Pollard's rho
@@ -448,7 +445,6 @@ pub(crate) fn factorize128_advanced(cofactors: &[(u128, usize)]) -> Vec<(u128, u
             i += 1;
 
             // increase max iterations after trying all methods
-            #[allow(unused_assignments)]
             if i % NMETHODS == 0 {
                 max_iter *= 4;
             }
@@ -775,8 +771,6 @@ where
 {
     let buf = NaiveBuffer::new();
     let config = Some(PrimalityTestConfig::strict());
-    // XXX: use miller-rabin for large numbers (more than 256 bits?), as BPSW could be too slow (need check)
-    //      the NIST recommends 5 rounds for 512 and 1024 bits. For 1536 bits, the recommendation is 4 rounds.
 
     // test (n-1)/2 first since its smaller
     let sophie_p = buf.is_prime(&(target >> 1), config);
diff --git a/src/traits.rs b/src/traits.rs
index 9dfdab4..170fdba 100644
--- a/src/traits.rs
+++ b/src/traits.rs
@@ -109,9 +109,12 @@ impl Default for PrimalityTestConfig {
 }
 
 impl PrimalityTestConfig {
-    /// Create a configuration with the **stongest deterministic** primality test available
+    /// Create a configuration with a very strong primality check.
+    /// Currently the configuration is BPSW test + SPRP test with 1 random base
     pub fn strict() -> Self {
-        Self::bpsw() // TODO: change to 2-base SPRP + VPRP
+        let mut config = Self::bpsw();
+        config.sprp_random_trials = 1;
+        config
     }
 
     /// Create a configuration for Baillie-PSW test (base 2 SPRP test + SLPRP test)