relp_num/integer/factorization/
size_large.rs

1use std::num::NonZeroU64;
2
3use smallvec::SmallVec;
4
5use crate::{NonZeroFactorizable, NonZeroFactorization, NonZeroSign, Ubig};
6use crate::integer::big::{BITS_PER_WORD, NonZeroUbig};
7use crate::integer::big::ops::div::div_assign_one_word;
8use crate::integer::big::ops::non_zero::{is_one_non_zero, shr};
9use crate::integer::big::ops::normalize::trailing_zeros;
10use crate::integer::factorization::{KEEP_RESIDUAL, NR_SMALL_PRIMES, size_64, start, TRIAL_DIVISION_LIMIT, TRIAL_DIVISION_LIMIT_USIZE};
11use crate::integer::factorization::prime::primes::SMALL_ODD_PRIMES;
12use crate::non_zero::NonZero;
13
14macro_rules! define {
15    ($ty:ident) => {
16        impl<const S: usize> NonZeroFactorizable for $ty<S> {
17            type Factor = usize;
18            type Power = u32;
19
20            fn factorize(&self) -> NonZeroFactorization<Self::Factor, Self::Power> {
21                assert!(self.is_not_zero(), "attempt to factorize zero");
22
23                let factors = factorize::<
24                    NR_SMALL_PRIMES, TRIAL_DIVISION_LIMIT, TRIAL_DIVISION_LIMIT_USIZE, KEEP_RESIDUAL, S
25                >(self.inner());
26                NonZeroFactorization { sign: NonZeroSign::Positive, factors }
27            }
28        }
29    }
30}
31
32define!(Ubig);
33define!(NonZeroUbig);
34
35fn factorize<
36    const NR_SMALL_PRIMES: usize,
37    const TRIAL_DIVISION_LIMIT: u64,
38    const TRIAL_DIVISION_LIMIT_USIZE: usize,
39    const KEEP_RESIDUAL: bool,
40    const S: usize,
41>(values: &[usize]) -> Vec<(usize, u32)> {
42    let (words, bits) = unsafe {
43        // SAFETY: The value is consistent so ends in a nonzero, and is not zero
44        trailing_zeros(values)
45    };
46
47    let mut x = shr::<S>(values, words, bits);
48    let total = words as u32 * BITS_PER_WORD + bits;
49    let mut factors = vec![];
50    if total > 0 {
51        factors.push((2, total));
52    }
53
54    if x.len() == 1 {
55        let as_small = unsafe {
56            // SAFETY: The value is not zero
57            NonZeroU64::new_unchecked(x[0] as u64)
58        };
59        let small = size_64::factorize::<
60            NR_SMALL_PRIMES, TRIAL_DIVISION_LIMIT, 0, KEEP_RESIDUAL,
61        >(as_small);
62        for (factor, power) in small {
63            factors.push((factor as usize, power));
64        }
65    } else {
66        // x.len() > 1
67        'odd: {
68            for divisor in &SMALL_ODD_PRIMES[..NR_SMALL_PRIMES] {
69                let divisor = *divisor as usize;
70
71                let mut counter = 0;
72                loop {
73                    let mut copy = x.clone();
74                    let remainder = unsafe {
75                        // SAFETY: divisor is nonzero and odd, copy is not zero
76                        div_assign_one_word(&mut copy, divisor)
77                    };
78                    if remainder == 0 {
79                        counter += 1;
80                        x = copy;
81                    } else {
82                        break;
83                    }
84                }
85
86                if counter > 0 {
87                    factors.push((divisor, counter));
88                }
89
90                if unsafe { is_one_non_zero(&x) } {
91                    break 'odd;
92                }
93            }
94
95            if TRIAL_DIVISION_LIMIT != 0 {
96                let start = start(NR_SMALL_PRIMES) as usize;
97                trial_division::<TRIAL_DIVISION_LIMIT_USIZE, S>(&mut x, start, &mut factors);
98
99                if unsafe { is_one_non_zero(&x) } {
100                    break 'odd;
101                }
102            }
103
104            if KEEP_RESIDUAL {
105                if let &[value] = x.as_slice() {
106                    factors.push((value, 1));
107                }
108            }
109        }
110    }
111
112    factors
113}
114
115// TODO(PERFORMANCE): Make `start` a const parameter
116fn trial_division<
117    const END: usize,
118    const S: usize,
119>(
120    x: &mut SmallVec<[usize; S]>,
121    start: usize,
122    factors: &mut Vec<(usize, u32)>,
123) {
124    let get_x_bits = |y: &[usize]| y.len() as u32 * BITS_PER_WORD - y[0].leading_zeros();
125
126    let mut divisor = start;
127    let mut x_bits = get_x_bits(x);
128    while {
129        let not_one = unsafe { !is_one_non_zero(x) };
130        let below_limit = divisor <= END;
131        let below_sqrt = {
132            let divisor_bits = BITS_PER_WORD - divisor.leading_zeros();
133            2 * divisor_bits <= x_bits
134        };
135
136        not_one && below_limit && below_sqrt
137    } {
138        let mut counter = 0;
139        loop {
140            let mut copy = x.clone();
141            let remainder = unsafe {
142                // SAFETY: divisor is nonzero and odd, copy is not zero
143                div_assign_one_word(&mut copy, divisor)
144            };
145            if remainder == 0 {
146                counter += 1;
147                *x = copy;
148            } else {
149                break;
150            }
151        }
152
153        if counter > 0 {
154            factors.push((divisor, counter));
155            x_bits = get_x_bits(x);
156        }
157
158        divisor += 2;
159    }
160}
161
162#[cfg(test)]
163mod test {
164    use std::str::FromStr;
165
166    use num_traits::One;
167    use smallvec::smallvec;
168
169    use crate::{NonZeroFactorizable, NonZeroFactorization, NonZeroSign, NonZeroUbig};
170    use crate::integer::factorization::size_large::factorize;
171
172    #[test]
173    fn test_factor() {
174        assert_eq!(NonZeroUbig::<4>::one().factorize(), NonZeroFactorization { sign: NonZeroSign::Positive, factors: vec![]});
175        assert_eq!(
176            NonZeroUbig::<4>::new(4).unwrap().factorize(),
177            NonZeroFactorization { sign: NonZeroSign::Positive, factors: vec![(2, 2)] },
178        );
179        assert_eq!(
180            unsafe { NonZeroUbig::<2>::from_inner_unchecked(smallvec![0, 351684787688]) }.factorize(),
181            NonZeroFactorization {
182                sign: NonZeroSign::Positive,
183                factors: vec![
184                    (2, 64 + 3),
185                    (13, 1),
186                    // (3381584497, 1),
187                ],
188            },
189        );
190        assert_eq!(
191            factorize::<5, 0, 0, true, 2>(&[0, 351684787688]),
192            vec![(2, 64 + 3), (13, 1), (3381584497, 1)],
193        );
194        let composite = NonZeroUbig::<8>::from_str("22429238517634168458101140012627848499653000000").unwrap();
195        assert_eq!(
196            factorize::<256, 100, 100, true, 8>(&composite),
197            vec![(2, 6), (3, 18), (5, 6), (11, 12), (18446744073709551437, 1)],
198        );
199    }
200}
relp_num/integer/factorization/size_large.rs

relp_num/integer/factorization/
size_large.rs
