num_prime/
factor.rs

1//! Implementations for various factorization algorithms.
2//!
3//! Note general prime number field sieve is not planned to be implemented, since it's too complex
4//!
5//! See <https://web.archive.org/web/20110331180514/https://diamond.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf>
6//! for a detailed comparison between different factorization algorithms
7
8// XXX: make the factorization method resumable? Maybe let all of them returns a Future
9
10use crate::traits::ExactRoots;
11use num_integer::{Integer, Roots};
12use num_modular::{ModularCoreOps, ModularUnaryOps};
13use num_traits::{CheckedAdd, FromPrimitive, NumRef, RefNum};
14use std::collections::BTreeMap;
15
16/// Find factors by trial division, returns a tuple of the found factors and the residual.
17///
18/// The target is guaranteed fully factored only if bound * bound > target, where bound = max(primes).
19/// The parameter limit additionally sets the maximum of primes to be tried.
20/// The residual will be Ok(1) or Ok(p) if fully factored.
21///
22/// TODO: implement fast check for small primes with BigInts in the precomputed table, and skip them in this function
23pub fn trial_division<
24    I: Iterator<Item = u64>,
25    T: Integer + Clone + Roots + NumRef + FromPrimitive,
26>(
27    primes: I,
28    target: T,
29    limit: Option<u64>,
30) -> (BTreeMap<u64, usize>, Result<T, T>)
31where
32    for<'r> &'r T: RefNum<T>,
33{
34    let tsqrt: T = Roots::sqrt(&target) + T::one();
35    let limit = if let Some(l) = limit {
36        tsqrt.clone().min(T::from_u64(l).unwrap())
37    } else {
38        tsqrt.clone()
39    };
40
41    let mut residual = target;
42    let mut result = BTreeMap::new();
43    let mut factored = false;
44    for (p, pt) in primes.map(|p| (p, T::from_u64(p).unwrap())) {
45        if &pt > &tsqrt {
46            factored = true;
47        }
48        if &pt > &limit {
49            break;
50        }
51
52        while residual.is_multiple_of(&pt) {
53            residual = residual / &pt;
54            *result.entry(p).or_insert(0) += 1;
55        }
56        if residual.is_one() {
57            factored = true;
58            break;
59        }
60    }
61
62    if factored {
63        (result, Ok(residual))
64    } else {
65        (result, Err(residual))
66    }
67}
68
69/// Find factors using Pollard's rho algorithm with Brent's loop detection algorithm
70///
71/// The returned values are the factor and the count of passed iterations.
72pub fn pollard_rho<
73    T: Integer
74        + FromPrimitive
75        + NumRef
76        + Clone
77        + for<'r> ModularCoreOps<&'r T, &'r T, Output = T>
78        + for<'r> ModularUnaryOps<&'r T, Output = T>,
79>(
80    target: &T,
81    start: T,
82    offset: T,
83    max_iter: usize,
84) -> (Option<T>, usize)
85where
86    for<'r> &'r T: RefNum<T>,
87{
88    let mut a = start.clone();
89    let mut b = start.clone();
90    let mut z = T::one() % target; // accumulator for gcd
91
92    // using Brent's loop detection, i = tortoise, j = hare
93    let (mut i, mut j) = (0usize, 1usize);
94
95    // backtracing states
96    let mut s = start;
97    let mut backtrace = false;
98
99    while i < max_iter {
100        i += 1;
101        a = a.sqm(&target).addm(&offset, &target);
102        if a == b {
103            return (None, i);
104        }
105
106        // FIXME: optimize abs_diff for montgomery form if we are going to use the abs_diff in the std lib
107        let diff = if b > a { &b - &a } else { &a - &b }; // abs_diff
108        z = z.mulm(&diff, &target);
109        if z.is_zero() {
110            // the factor is missed by a combined GCD, do backtracing
111            if backtrace {
112                // ultimately failed
113                return (None, i);
114            } else {
115                backtrace = true;
116                a = std::mem::replace(&mut s, T::one()); // s is discarded
117                z = T::one() % target; // clear gcd
118                continue;
119            }
120        }
121
122        // here we check gcd every 2^k steps or 128 steps
123        // larger batch size leads to large overhead when backtracing.
124        // reference: https://www.cnblogs.com/812-xiao-wen/p/10544546.html
125        if i == j || i & 127 == 0 || backtrace {
126            let d = z.gcd(target);
127            if !d.is_one() && &d != target {
128                return (Some(d), i);
129            }
130
131            // save state
132            s = a.clone();
133        }
134
135        // when tortoise catches up with hare, let hare jump to the next stop
136        if i == j {
137            b = a.clone();
138            j <<= 1;
139        }
140    }
141
142    (None, i)
143}
144
145/// This function implements Shanks's square forms factorization (SQUFOF).
146///
147/// The input is usually multiplied by a multiplier, and the multiplied integer should be put in
148/// the `mul_target` argument. The multiplier can be choosen from SQUFOF_MULTIPLIERS, or other square-free odd numbers.
149/// The returned values are the factor and the count of passed iterations.
150///
151/// The max iteration can be choosed as 2*n^(1/4), based on Theorem 4.22 from [1].
152///
153/// Reference: Gower, J., & Wagstaff Jr, S. (2008). Square form factorization.
154/// In [1] [Mathematics of Computation](https://homes.cerias.purdue.edu/~ssw/gowerthesis804/wthe.pdf)
155/// or [2] [his thesis](https://homes.cerias.purdue.edu/~ssw/gowerthesis804/wthe.pdf)
156/// The code is from [3] [Rosetta code](https://rosettacode.org/wiki/Square_form_factorization)
157pub fn squfof<T: Integer + NumRef + Clone + ExactRoots + std::fmt::Debug>(
158    target: &T,
159    mul_target: T,
160    max_iter: usize,
161) -> (Option<T>, usize)
162where
163    for<'r> &'r T: RefNum<T>,
164{
165    assert!(
166        &mul_target.is_multiple_of(&target),
167        "mul_target should be multiples of target"
168    );
169    let rd = Roots::sqrt(&mul_target); // root of k*N
170
171    /// Reduction operator for binary quadratic forms. It's equivalent to
172    /// the one used in the `num-irrational` crate, in a little different form.
173    ///
174    /// This function reduces (a, b, c) = (qm1, p, q), updates qm1 and q, returns new p.
175    #[inline]
176    fn rho<T: Integer + Clone + NumRef>(rd: &T, p: &T, q: &mut T, qm1: &mut T) -> T
177    where
178        for<'r> &'r T: RefNum<T>,
179    {
180        let b = (rd + p).div_floor(&*q);
181        let new_p = &b * &*q - p;
182        let new_q = if p > &new_p {
183            &*qm1 + b * (p - &new_p)
184        } else {
185            &*qm1 - b * (&new_p - p)
186        };
187
188        *qm1 = std::mem::replace(q, new_q);
189        new_p
190    }
191
192    // forward loop, search principal cycle
193    let (mut p, mut q, mut qm1) = (rd.clone(), &mul_target - &rd * &rd, T::one());
194    if q.is_zero() {
195        // shortcut for perfect square
196        return (Some(rd), 0);
197    }
198
199    for i in 1..max_iter {
200        p = rho(&rd, &p, &mut q, &mut qm1);
201        if i.is_odd() {
202            if let Some(rq) = q.sqrt_exact() {
203                let b = (&rd - &p) / &rq;
204                let mut u = b * &rq + &p;
205                let (mut v, mut vm1) = ((&mul_target - &u * &u) / &rq, rq);
206
207                // backward loop, search ambiguous cycle
208                loop {
209                    let new_u = rho(&rd, &u, &mut v, &mut vm1);
210                    if new_u == u {
211                        break;
212                    } else {
213                        u = new_u
214                    }
215                }
216
217                let d = target.gcd(&u);
218                if d > T::one() && &d < target {
219                    return (Some(d), i);
220                }
221            }
222        }
223    }
224    (None, max_iter)
225}
226
227/// Good squfof multipliers sorted by efficiency descendingly, from Dana Jacobsen.
228///
229/// Note: square-free odd numbers are suitable as SQUFOF multipliers
230pub const SQUFOF_MULTIPLIERS: [u16; 38] = [
231    3 * 5 * 7 * 11,
232    3 * 5 * 7,
233    3 * 5 * 7 * 11 * 13,
234    3 * 5 * 7 * 13,
235    3 * 5 * 7 * 11 * 17,
236    3 * 5 * 11,
237    3 * 5 * 7 * 17,
238    3 * 5,
239    3 * 5 * 7 * 11 * 19,
240    3 * 5 * 11 * 13,
241    3 * 5 * 7 * 19,
242    3 * 5 * 7 * 13 * 17,
243    3 * 5 * 13,
244    3 * 7 * 11,
245    3 * 7,
246    5 * 7 * 11,
247    3 * 7 * 13,
248    5 * 7,
249    3 * 5 * 17,
250    5 * 7 * 13,
251    3 * 5 * 19,
252    3 * 11,
253    3 * 7 * 17,
254    3,
255    3 * 11 * 13,
256    5 * 11,
257    3 * 7 * 19,
258    3 * 13,
259    5,
260    5 * 11 * 13,
261    5 * 7 * 19,
262    5 * 13,
263    7 * 11,
264    7,
265    3 * 17,
266    7 * 13,
267    11,
268    1,
269];
270
271/// William Hart's one line factorization algorithm for 64 bit integers.
272///
273/// The number to be factored could be multiplied by a smooth number (coprime to the target)
274/// to speed up, put the multiplied number in the `mul_target` argument. A good multiplier given by Hart is 480.
275/// `iters` determine the range for iterating the inner multiplier itself. The returned values are the factor
276/// and the count of passed iterations.
277///
278///
279/// The one line factorization algorithm is especially good at factoring semiprimes with form pq,
280/// 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.
281///
282/// Reference: Hart, W. B. (2012). A one line factoring algorithm. Journal of the Australian Mathematical Society, 92(1), 61-69. doi:10.1017/S1446788712000146
283// TODO: add multipliers preset for one_line method?
284pub fn one_line<T: Integer + NumRef + FromPrimitive + ExactRoots + CheckedAdd>(
285    target: &T,
286    mul_target: T,
287    max_iter: usize,
288) -> (Option<T>, usize)
289where
290    for<'r> &'r T: RefNum<T>,
291{
292    assert!(
293        &mul_target.is_multiple_of(&target),
294        "mul_target should be multiples of target"
295    );
296
297    let mut ikn = mul_target.clone();
298    for i in 1..max_iter {
299        let s = ikn.sqrt() + T::one(); // assuming target is not perfect square
300        let m = &s * &s - &ikn;
301        if let Some(t) = m.sqrt_exact() {
302            let g = target.gcd(&(s - t));
303            if !g.is_one() && &g != target {
304                return (Some(g), i);
305            }
306        }
307
308        // prevent overflow
309        ikn = if let Some(n) = (&ikn).checked_add(&mul_target) {
310            n
311        } else {
312            return (None, i);
313        }
314    }
315    return (None, max_iter);
316}
317
318// TODO: ECM, (self initialize) Quadratic sieve, Lehman's Fermat(https://en.wikipedia.org/wiki/Fermat%27s_factorization_method, n_factor_lehman)
319// REF: https://pypi.org/project/primefac/
320//      http://flintlib.org/doc/ulong_extras.html#factorisation
321//      https://github.com/zademn/facto-rs/
322//      https://github.com/elmomoilanen/prime-factorization
323//      https://cseweb.ucsd.edu/~ethome/teaching/2022-cse-291-14/
324fn pollard_pp1() {}
325fn williams_pp1() {}
326
327#[cfg(test)]
328mod tests {
329    use super::*;
330    use crate::mint::SmallMint;
331    use num_modular::MontgomeryInt;
332    use rand::random;
333
334    #[test]
335    fn pollard_rho_test() {
336        assert_eq!(pollard_rho(&8051u16, 2, 1, 100).0, Some(97));
337        assert!(matches!(pollard_rho(&8051u16, random(), 1, 100).0, Some(i) if i == 97 || i == 83));
338        assert_eq!(pollard_rho(&455459u32, 2, 1, 100).0, Some(743));
339
340        // Mint test
341        for _ in 0..10 {
342            let target = random::<u16>() | 1;
343            let start = random::<u16>() % target;
344            let offset = random::<u16>() % target;
345
346            let expect = pollard_rho(&target, start, offset, 65536);
347            let mint_result = pollard_rho(
348                &SmallMint::from(target),
349                MontgomeryInt::new(start, &target).into(),
350                MontgomeryInt::new(offset, &target).into(),
351                65536,
352            );
353            assert_eq!(expect.0, mint_result.0.map(|v| v.value()));
354        }
355    }
356
357    #[test]
358    fn squfof_test() {
359        // case from wikipedia
360        assert_eq!(squfof(&11111u32, 11111u32, 100).0, Some(41));
361
362        // cases from https://rosettacode.org/wiki/Square_form_factorization
363        let cases: Vec<u64> = vec![
364            2501,
365            12851,
366            13289,
367            75301,
368            120787,
369            967009,
370            997417,
371            7091569,
372            5214317,
373            20834839,
374            23515517,
375            33409583,
376            44524219,
377            13290059,
378            223553581,
379            2027651281,
380            11111111111,
381            100895598169,
382            60012462237239,
383            287129523414791,
384            9007199254740931,
385            11111111111111111,
386            314159265358979323,
387            384307168202281507,
388            419244183493398773,
389            658812288346769681,
390            922337203685477563,
391            1000000000000000127,
392            1152921505680588799,
393            1537228672809128917,
394            // this case should success at step 276, from https://rosettacode.org/wiki/Talk:Square_form_factorization
395            4558849,
396        ];
397        for n in cases {
398            let d = squfof(&n, n, 40000)
399                .0
400                .or(squfof(&n, 3 * n, 40000).0)
401                .or(squfof(&n, 5 * n, 40000).0)
402                .or(squfof(&n, 7 * n, 40000).0)
403                .or(squfof(&n, 11 * n, 40000).0);
404            assert!(matches!(d, Some(_)), "{}", n);
405        }
406    }
407
408    #[test]
409    fn one_line_test() {
410        assert_eq!(one_line(&11111u32, 11111u32, 100).0, Some(271));
411    }
412}
num_prime/factor.rs

num_prime/
factor.rs
