1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

// * This file is part of the uutils coreutils package.
// *
// * (c) 2015 kwantam <kwantam@gmail.com>
// * (c) 2020 nicoo   <nicoo@debian.org>
// *
// * For the full copyright and license information, please view the LICENSE file
// * that was distributed with this source code.

// spell-checker: ignore (ToDO) INVS

use crate::Factors;

include!(concat!(env!("OUT_DIR"), "/prime_table.rs"));

pub fn factor(num: &mut u64, factors: &mut Factors) {
    for &(prime, inv, ceil) in PRIME_INVERSIONS_U64 {
        if *num == 1 {
            break;
        }

        // inv = prime^-1 mod 2^64
        // ceil = floor((2^64-1) / prime)
        // if (num * inv) mod 2^64 <= ceil, then prime divides num
        // See https://math.stackexchange.com/questions/1251327/
        // for a nice explanation.
        let mut k = 0;
        loop {
            let x = num.wrapping_mul(inv);

            // While prime divides num
            if x <= ceil {
                *num = x;
                k += 1;
                #[cfg(feature = "coz")]
                coz::progress!("factor found");
            } else {
                if k > 0 {
                    factors.add(prime, k);
                }
                break;
            }
        }
    }
}

pub const CHUNK_SIZE: usize = 8;
pub fn factor_chunk(n_s: &mut [u64; CHUNK_SIZE], f_s: &mut [Factors; CHUNK_SIZE]) {
    for &(prime, inv, ceil) in PRIME_INVERSIONS_U64 {
        if n_s[0] == 1 && n_s[1] == 1 && n_s[2] == 1 && n_s[3] == 1 {
            break;
        }

        for (num, factors) in n_s.iter_mut().zip(f_s.iter_mut()) {
            if *num == 1 {
                continue;
            }
            let mut k = 0;
            loop {
                let x = num.wrapping_mul(inv);

                // While prime divides num
                if x <= ceil {
                    *num = x;
                    k += 1;
                } else {
                    if k > 0 {
                        factors.add(prime, k);
                    }
                    break;
                }
            }
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::Factors;
    use quickcheck::quickcheck;
    use rand::{rngs::SmallRng, Rng, SeedableRng};

    quickcheck! {
        fn chunk_vs_iter(seed: u64) -> () {
            let mut rng = SmallRng::seed_from_u64(seed);
            let mut n_c: [u64; CHUNK_SIZE] = rng.gen();
            let mut f_c: [Factors; CHUNK_SIZE] = rng.gen();

            let mut n_i = n_c;
            let mut f_i = f_c.clone();
            for (n, f) in n_i.iter_mut().zip(f_i.iter_mut()) {
                factor(n, f);
            }

            factor_chunk(&mut n_c, &mut f_c);

            assert_eq!(n_i, n_c);
            assert_eq!(f_i, f_c);
        }
    }
}