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};
11use crate::integer::factorization::prime::primes::SMALL_ODD_PRIMES;
12use crate::non_zero::NonZero;
13use crate::Prime;
14use crate::traits::factorization::FactorizationResidual;
15
16/// A big integer that was not decomposed entirely leaves a big integer behind.
17impl<const S: usize> FactorizationResidual for Ubig<S> {
18    #[inline]
19    fn one() -> Self {
20        <Self as num_traits::One>::one()
21    }
22
23    #[inline]
24    fn is_one(&self) -> bool {
25        <Self as num_traits::One>::is_one(self)
26    }
27}
28
29macro_rules! define {
30    ($ty:ident) => {
31        impl<const S: usize> NonZeroFactorizable for $ty<S> {
32            type Factor = usize;
33            type Power = u32;
34            /// A cofactor that could not be decomposed can be larger than a machine word, so it
35            /// doesn't fit in a `Factor`.
36            ///
37            /// Both `Ubig` and `NonZeroUbig` use this same type, such that a `Ratio` of the two has
38            /// a residual with two halves of the same type.
39            type Residual = Ubig<S>;
40
41            fn factorize(&self) -> NonZeroFactorization<Self::Factor, Self::Power, Self::Residual> {
42                assert!(self.is_not_zero(), "attempt to factorize zero");
43
44                let (factors, residual) = factorize::<
45                    NR_SMALL_PRIMES, TRIAL_DIVISION_LIMIT, KEEP_RESIDUAL, S
46                >(self.inner());
47                let residual = unsafe {
48                    // SAFETY: The residual is a divisor of a nonzero value, so it is nonzero, and
49                    // it was normalized by the routines that produced it
50                    Ubig::from_inner_unchecked(residual)
51                };
52
53                NonZeroFactorization { sign: NonZeroSign::Positive, factors, residual }
54            }
55        }
56    }
57}
58
59define!(Ubig);
60define!(NonZeroUbig);
61
62/// Approximately factorize a big integer.
63///
64/// # Return value
65///
66/// The `(factor, power)` tuples and the residual: the part of the value that was not decomposed.
67/// The residual is one when the factorization is complete, and
68/// `value == residual * product(factor ^ power)` always holds.
69fn factorize<
70    const NR_SMALL_PRIMES: usize,
71    const TRIAL_DIVISION_LIMIT: u64,
72    const KEEP_RESIDUAL: bool,
73    const S: usize,
74>(values: &[usize]) -> (Vec<(usize, u32)>, SmallVec<[usize; S]>) {
75    let (words, bits) = unsafe {
76        // SAFETY: The value is consistent so ends in a nonzero, and is not zero
77        trailing_zeros(values)
78    };
79
80    let mut x = shr::<S>(values, words, bits);
81    let total = words as u32 * BITS_PER_WORD + bits;
82    let mut factors = vec![];
83    if total > 0 {
84        factors.push((2, total));
85    }
86
87    if x.len() == 1 {
88        let as_small = unsafe {
89            // SAFETY: The value is not zero
90            NonZeroU64::new_unchecked(x[0] as u64)
91        };
92        let (small, residual) = size_64::factorize::<
93            NR_SMALL_PRIMES, TRIAL_DIVISION_LIMIT, 0, KEEP_RESIDUAL,
94        >(as_small);
95        for (factor, power) in small {
96            factors.push((factor as usize, power));
97        }
98        // Propagate the residual of the single word routine; it is nonzero, so `x` stays well
99        // formed. Note that `x.len() == 1`.
100        x[0] = residual as usize;
101    } else {
102        // x.len() > 1
103        'odd: {
104            for divisor in &SMALL_ODD_PRIMES[..NR_SMALL_PRIMES] {
105                let divisor = *divisor as usize;
106
107                let mut counter = 0;
108                loop {
109                    let mut copy = x.clone();
110                    let remainder = unsafe {
111                        // SAFETY: divisor is nonzero and odd, copy is not zero
112                        div_assign_one_word(&mut copy, divisor)
113                    };
114                    if remainder == 0 {
115                        counter += 1;
116                        x = copy;
117                    } else {
118                        break;
119                    }
120                }
121
122                if counter > 0 {
123                    factors.push((divisor, counter));
124                }
125
126                if unsafe { is_one_non_zero(&x) } {
127                    break 'odd;
128                }
129            }
130
131            if TRIAL_DIVISION_LIMIT != 0 {
132                let start = start(NR_SMALL_PRIMES) as usize;
133                trial_division::<TRIAL_DIVISION_LIMIT, S>(&mut x, start, &mut factors);
134
135                if unsafe { is_one_non_zero(&x) } {
136                    break 'odd;
137                }
138            }
139
140            // A remainder that fits in a single word and is proven prime is a genuine factor, so it
141            // is moved out of the residual. Note that `x` is set to one, such that it is not
142            // counted twice. Anything else stays in the residual: a larger remainder doesn't even
143            // fit in a `Factor`, and one that isn't proven prime would be a lie as a factor.
144            if KEEP_RESIDUAL && let &[value] = x.as_slice() && (value as u64).is_prime() {
145                factors.push((value, 1));
146                x[0] = 1;
147            }
148        }
149    }
150
151    (factors, x)
152}
153
154// TODO(PERFORMANCE): Make `start` a const parameter
155fn trial_division<
156    const END: u64,
157    const S: usize,
158>(
159    x: &mut SmallVec<[usize; S]>,
160    start: usize,
161    factors: &mut Vec<(usize, u32)>,
162) {
163    // The words are stored little endian, so the most significant word is the last one.
164    let get_x_bits = |y: &[usize]| {
165        let most_significant = *y.last().expect("the value is nonzero, so it has a word");
166        (y.len() as u32 - 1) * BITS_PER_WORD + most_significant.bit_width()
167    };
168
169    let mut divisor = start;
170    let mut x_bits = get_x_bits(x);
171    while {
172        let not_one = unsafe { !is_one_non_zero(x) };
173        let below_limit = divisor as u64 <= END;
174        let below_sqrt = {
175            let divisor_bits = divisor.bit_width();
176            2 * divisor_bits <= x_bits
177        };
178
179        not_one && below_limit && below_sqrt
180    } {
181        let mut counter = 0;
182        loop {
183            let mut copy = x.clone();
184            let remainder = unsafe {
185                // SAFETY: divisor is nonzero and odd, copy is not zero
186                div_assign_one_word(&mut copy, divisor)
187            };
188            if remainder == 0 {
189                counter += 1;
190                *x = copy;
191            } else {
192                break;
193            }
194        }
195
196        if counter > 0 {
197            factors.push((divisor, counter));
198            x_bits = get_x_bits(x);
199        }
200
201        divisor += 2;
202    }
203}
204
205#[cfg(test)]
206mod test {
207    use std::str::FromStr;
208
209    use num_traits::One;
210    use smallvec::smallvec;
211
212    use crate::{NonZeroFactorizable, NonZeroFactorization, NonZeroSign, NonZeroUbig, Ubig};
213    use crate::integer::factorization::size_large::factorize;
214
215    #[test]
216    fn test_factor() {
217        assert_eq!(
218            NonZeroUbig::<4>::one().factorize(),
219            NonZeroFactorization { sign: NonZeroSign::Positive, factors: vec![], residual: Ubig::one() },
220        );
221        assert_eq!(
222            NonZeroUbig::<4>::new(4).unwrap().factorize(),
223            NonZeroFactorization { sign: NonZeroSign::Positive, factors: vec![(2, 2)], residual: Ubig::one() },
224        );
225        // The crate wide constants don't attempt to split the cofactor, so it ends up in the
226        // residual rather than being dropped.
227        assert_eq!(
228            unsafe { NonZeroUbig::<2>::from_inner_unchecked(smallvec![0, 351684787688]) }.factorize(),
229            NonZeroFactorization {
230                sign: NonZeroSign::Positive,
231                factors: vec![
232                    (2, 64 + 3),
233                    (13, 1),
234                ],
235                residual: Ubig::new(3381584497),
236            },
237        );
238        assert_eq!(
239            factorize::<5, 0, true, 2>(&[0, 351684787688]),
240            (vec![(2, 64 + 3), (13, 1), (3381584497, 1)], smallvec![1]),
241        );
242        let composite = NonZeroUbig::<8>::from_str("22429238517634168458101140012627848499653000000").unwrap();
243        assert_eq!(
244            factorize::<256, 100, true, 8>(&composite),
245            (vec![(2, 6), (3, 18), (5, 6), (11, 12), (18446744073709551437, 1)], smallvec![1]),
246        );
247    }
248
249    /// A cofactor that spans more than one word must survive in the residual.
250    #[test]
251    fn test_multi_word_residual() {
252        // Two primes just below `2 ** 64`, so their product needs two words.
253        let left = 2_u128.pow(64) - 59;
254        let right = 2_u128.pow(64) - 83;
255        let product = left * right;
256
257        let value = NonZeroUbig::<2>::new_u128(product).unwrap();
258        let NonZeroFactorization { sign, factors, residual } = value.factorize();
259        assert_eq!(sign, NonZeroSign::Positive);
260        assert_eq!(factors, vec![]);
261        assert_eq!(residual, Ubig::<2>::new_u128(product));
262
263        // Even with the residual explicitly kept as a factor, a multi word cofactor can't be one,
264        // so it stays in the residual.
265        let (factors, residual) = factorize::<256, 0, true, 2>(value.inner());
266        assert_eq!(factors, vec![]);
267        assert_eq!(residual.as_slice(), value.inner().as_slice());
268    }
269}
relp_num/integer/factorization/size_large.rs

relp_num/integer/factorization/
size_large.rs
