modify factor impl to eliminate overflow issue

This change does the following:

1. Updates the arithmetic functions in `src/factor/numeric.rs` to
   correctly handle all cases up to 2^64. When numbers are larger
   than 2^63, we fall back to slightly slower routines that check
   for and handle overflow.

2. Since the arithmetic functions will now not overflow, we no longer
   need the safety net trial division implementation. We now always
   use Pollard's rho after eliminating small (<=13 bit) primes.

3. Slight tweak in `src/factor/gen_table.rs` to generate the first
   1027 primes, which means we test every prime of 13 or fewer bits
   before going into Pollard's rho. Includes corresponding update in
   `src/factor/prime_table.rs` and the Makefile to reflect this.

4. Add a new test that generates random numbers with exclusively
   large (14 to 50 bit) prime factors. This exercises the possible
   overflow paths.

5. Add another new test that checks the `is_prime()` function against
   a few dozen 64-bit primes. Again this is to exercise possible
   overflow paths.

src/factor/factor.rs

@@ -19,7 +19,7 @@ extern crate libc;
 extern crate rand;
 
 use numeric::*;
-use prime_table::{P_INVS_U64, NEXT_PRIME};
+use prime_table::P_INVS_U64;
 use std::cmp::{max, min};
 use std::io::{stdin, BufRead, BufReader, Write};
 use std::num::Wrapping;
@@ -36,33 +36,12 @@ mod prime_table;
 static VERSION: &'static str = "1.0.0";
 static NAME: &'static str = "factor";
 
-fn trial_division_slow(mut num: u64, factors: &mut Vec<u64>) {
-    // assumption: this number has already been run through
-    // trial_division, which checks primes from the table.
-    // The first candidate we need to check is NEXT_PRIME.
-    let mut i = NEXT_PRIME;
-    while i * i <= num {
-        while num % i == 0 {
-            num /= i;
-            factors.push(i);
-            if is_prime(num) {
-                factors.push(num);
-                return;
-            }
-            if num < 1 << 63 {
-                // once we're small enough, switch to Pollard's rho
-                return rho_pollard_factor(num, factors);
-            }
-        }
-        i += 2;
-    }
-    if num > 1 {
-        factors.push(num);
-    }
-}
-
 fn rho_pollard_pseudorandom_function(x: u64, a: u64, b: u64, num: u64) -> u64 {
-    (multiply(a, multiply(x, x, num), num) + b) % num
+    if num < 1 << 63 {
+        (sm_mul(a, sm_mul(x, x, num), num) + b) % num
+    } else {
+        big_add(big_mul(a, big_mul(x, x, num), num), b, num)
+    }
 }
 
 fn gcd(mut a: u64, mut b: u64) -> u64 {
@@ -148,13 +127,13 @@ fn table_division(mut num: u64, factors: &mut Vec<u64>) {
     // do we still have more factoring to do?
     // Decide whether to use Pollard Rho or slow divisibility based on
     // number's size:
-    if num >= 1 << 63 {
+    //if num >= 1 << 63 {
         // number is too big to use rho pollard without overflowing
-        trial_division_slow(num, factors);
-    } else if num > 1 {
+        //trial_division_slow(num, factors);
+    //} else if num > 1 {
         // number is still greater than 1, but not so big that we have to worry
         rho_pollard_factor(num, factors);
-    }
+    //}
 }
 
 fn print_factors(num: u64) {

src/factor/gen_table.rs

@@ -8,7 +8,9 @@
 */
 
 //! Generate a table of the multiplicative inverses of p_i mod 2^64
-//! for the first 10000 odd primes.
+//! for the first 1027 odd primes (all 13 bit and smaller primes).
+//! You can supply a commandline argument to override the default
+//! value of 1027 for the number of entries in the table.
 //!
 //! 2 has no multiplicative inverse mode 2^64 because 2 | 2^64,
 //! and in any case divisibility by two is trivial by checking the LSB.
@@ -65,7 +67,7 @@ fn inv_mod_u64(a: u64) -> Option<u64> {
 #[cfg_attr(test, allow(dead_code))]
 fn main() {
     // By default, we print the multiplicative inverses mod 2^64 of the first 1k primes
-    let n = args().skip(1).next().unwrap_or("1000".to_string()).parse::<usize>().ok().unwrap_or(1000);
+    let n = args().skip(1).next().unwrap_or("1027".to_string()).parse::<usize>().ok().unwrap_or(1027);
 
     print!("{}", PREAMBLE);
     let mut cols = 3;
@@ -91,7 +93,7 @@ fn main() {
         x = next;
     }
 
-    print!("\n];\n\npub const NEXT_PRIME: u64 = {};\n", x);
+    print!("\n];\n\n#[allow(dead_code)]\npub const NEXT_PRIME: u64 = {};\n", x);
 }
 
 #[test]

src/factor/numeric.rs

@@ -2,16 +2,36 @@
 * This file is part of the uutils coreutils package.
 *
 * (c) Wiktor Kuropatwa <wiktor.kuropatwa@gmail.com>
+* (c) kwantam <kwantam@gmail.com>
+*     20150507 added big_ routines to prevent overflow when num > 2^63
 *
 * For the full copyright and license information, please view the LICENSE file
 * that was distributed with this source code.
 */
 
+use std::u64::MAX as MAX_U64;
+use std::num::Wrapping;
+
+pub fn big_add(a: u64, b: u64, m: u64) -> u64 {
+    let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
+    let msb_mod_m = msb_mod_m % m;
+
+    let Wrapping(res) = Wrapping(a) + Wrapping(b);
+    let res = if b <= MAX_U64 - a {
+        res
+    } else {
+        (res + msb_mod_m) % m
+    };
+
+    res
+}
+
 // computes (a + b) % m using the russian peasant algorithm
-pub fn multiply(mut a: u64, mut b: u64, m: u64) -> u64 {
+// CAUTION: Will overflow if m >= 2^63
+pub fn sm_mul(mut a: u64, mut b: u64, m: u64) -> u64 {
     let mut result = 0;
     while b > 0 {
-        if b & 1 > 0 {
+        if b & 1 != 0 {
             result = (result + a) % m;
         }
         a = (a << 1) % m;
@@ -20,14 +40,44 @@ pub fn multiply(mut a: u64, mut b: u64, m: u64) -> u64 {
     result
 }
 
+// computes (a + b) % m using the russian peasant algorithm
+// Only necessary when m >= 2^63; otherwise, just wastes time.
+pub fn big_mul(mut a: u64, mut b: u64, m: u64) -> u64 {
+    // precompute 2^64 mod m, since we expect to wrap
+    let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
+    let msb_mod_m = msb_mod_m % m;
+
+    let mut result = 0;
+    while b > 0 {
+        if b & 1 != 0 {
+            let Wrapping(next_res) = Wrapping(result) + Wrapping(a);
+            let next_res = next_res % m;
+            result = if result <= MAX_U64 - a {
+                next_res
+            } else {
+                (next_res + msb_mod_m) % m
+            };
+        }
+        let Wrapping(next_a) = Wrapping(a) << 1;
+        let next_a = next_a % m;
+        a = if a < 1 << 63 {
+            next_a
+        } else {
+            (next_a + msb_mod_m) % m
+        };
+        b >>= 1;
+    }
+    result
+}
+
 // computes a.pow(b) % m
-fn pow(mut a: u64, mut b: u64, m: u64) -> u64 {
+fn pow(mut a: u64, mut b: u64, m: u64, mul: fn(u64, u64, u64) -> u64) -> u64 {
     let mut result = 1;
     while b > 0 {
-        if b & 1 > 0 {
-            result = multiply(result, a, m);
+        if b & 1 != 0 {
+            result = mul(result, a, m);
         }
-        a = multiply(a, a, m);
+        a = mul(a, a, m);
         b >>= 1;
     }
     result
@@ -37,10 +87,17 @@ fn witness(mut a: u64, exponent: u64, m: u64) -> bool {
     if a == 0 {
         return false;
     }
-    if pow(a, m-1, m) != 1 {
+
+    let mul = if m < 1 << 63 {
+        sm_mul as fn(u64, u64, u64) -> u64
+    } else {
+        big_mul as fn(u64, u64, u64) -> u64
+    };
+
+    if pow(a, m-1, m, mul) != 1 {
         return true;
     }
-    a = pow(a, exponent, m);
+    a = pow(a, exponent, m, mul);
     if a == 1 {
         return false;
     }
@@ -51,7 +108,7 @@ fn witness(mut a: u64, exponent: u64, m: u64) -> bool {
         if a == m-1 {
             return false;
         }
-        a = multiply(a, a, m);
+        a = mul(a, a, m);
     }
 }
 
@@ -71,10 +128,5 @@ pub fn is_prime(num: u64) -> bool {
     // These witnesses detect all composites up to at least 2^64.
     // Discovered by Jim Sinclair, according to http://miller-rabin.appspot.com
     let witnesses = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
-    for wit in witnesses.iter() {
-        if witness(*wit % num, exponent, num) {
-            return false;
-        }
-    }
-    true
+    ! witnesses.iter().any(|&wit| witness(wit % num, exponent, num))
 }

src/factor/prime_table.rs

@@ -512,6 +512,21 @@ pub const P_INVS_U64: &'static [(u64, u64, u64)] = &[
     (7879, 15749618908978823927, 2341254483273201), (7883, 13069271299209418467, 2340066481505715),
     (7901, 7641588830939294069, 2334735359284843), (7907, 12420698804657854155, 2332963712369995),
     (7919, 12818844884155027471, 2329428472497733), (7927, 7937661667140567751, 2327077592242910),
+    (7933, 17988656517612138069, 2325317543641693), (7937, 12066964247284867329, 2324145656256715),
+    (7949, 10178314191381317573, 2320637070538376), (7951, 5681810619609444335, 2320053335896057),
+    (7963, 6180574304741565203, 2316557085735219), (7993, 7117447607071219465, 2307862388803897),
+    (8009, 15818733711853814521, 2303251850881452), (8011, 6594866437037093475, 2302676828574404),
+    (8017, 12876135691215997361, 2300953483062186), (8039, 16092426444697734231, 2294656558491050),
+    (8053, 7540876877021214941, 2290667338098789), (8059, 6493784953110063027, 2288961915089905),
+    (8069, 14763796409470353741, 2286125179540160), (8081, 12739918163021037937, 2282730364275405),
+    (8087, 10456272391972868135, 2281036734723575), (8089, 11422888004105716905, 2280472749871374),
+    (8093, 18241602968231501493, 2279345616422779), (8101, 4561020661602299949, 2277094688767997),
+    (8111, 16802434375116035919, 2274287273296702), (8117, 16051417197561976477, 2272606144352538),
+    (8123, 11545395404498259315, 2270927498917832), (8147, 11443457045357563995, 2264237642532165),
+    (8161, 12755174219819017249, 2260353397097114), (8167, 10672323466178233303, 2258692797074758),
+    (8171, 17536936478347301059, 2257587085266130), (8179, 16031233265182478651, 2255378906187743),
+    (8191, 18442239924259250175, 2252074725150720),
 ];
 
-pub const NEXT_PRIME: u64 = 7933;
+#[allow(dead_code)]
+pub const NEXT_PRIME: u64 = 8209;