feanor_math/algorithms/
int_factor.rs

1use crate::algorithms::ec_factor::lenstra_ec_factor;
2use crate::computation::{DontObserve, no_error};
3use crate::divisibility::DivisibilityRingStore;
4use crate::homomorphism::*;
5use crate::integer::*;
6use crate::ordered::{OrderedRing, OrderedRingStore};
7use crate::primitive_int::{StaticRing, StaticRingBase};
8use crate::ring::*;
9use crate::rings::zn::{ZnOperation, ZnRing, ZnRingStore, choose_zn_impl};
10use crate::{DEFAULT_PROBABILISTIC_REPETITIONS, algorithms};
11
12struct ECFactorInt<I>
13where
14    I: RingStore,
15    I::Type: IntegerRing,
16{
17    result_ring: I,
18}
19
20impl<I> ZnOperation for ECFactorInt<I>
21where
22    I: RingStore,
23    I::Type: IntegerRing,
24{
25    type Output<'a>
26        = El<I>
27    where
28        Self: 'a;
29
30    fn call<'a, R>(self, ring: R) -> El<I>
31    where
32        Self: 'a,
33        R: 'a + RingStore + Send + Sync,
34        R::Type: ZnRing,
35        El<R>: Send,
36    {
37        int_cast(
38            lenstra_ec_factor(&ring, DontObserve).unwrap_or_else(no_error),
39            self.result_ring,
40            ring.integer_ring(),
41        )
42    }
43}
44
45/// If the given integer is a power of a prime `p`, returns `Some((p, ln(n)/ln(p)))`.
46pub fn is_prime_power<I>(ZZ: I, n: &El<I>) -> Option<(El<I>, usize)>
47where
48    I: RingStore + Copy,
49    I::Type: IntegerRing,
50{
51    if algorithms::miller_rabin::is_prime(ZZ, n, DEFAULT_PROBABILISTIC_REPETITIONS) {
52        return Some((ZZ.clone_el(n), 1));
53    }
54    let (p, e) = is_power(ZZ, n)?;
55    if algorithms::miller_rabin::is_prime(ZZ, &p, DEFAULT_PROBABILISTIC_REPETITIONS) {
56        return Some((p, e));
57    } else {
58        return None;
59    }
60}
61
62fn is_power<I>(ZZ: I, n: &El<I>) -> Option<(El<I>, usize)>
63where
64    I: RingStore + Copy,
65    I::Type: IntegerRing,
66{
67    assert!(!ZZ.is_zero(n));
68    for i in (2..=ZZ.abs_log2_ceil(n).unwrap()).rev() {
69        let root = algorithms::int_bisect::root_floor(ZZ, ZZ.clone_el(n), i);
70        if ZZ.eq_el(&ZZ.pow(root, i), n) {
71            return Some((algorithms::int_bisect::root_floor(ZZ, ZZ.clone_el(n), i), i));
72        }
73    }
74    return None;
75}
76
77/// Factors the given integer.
78///
79/// Returns a list of all factors with their multipliplicities.
80pub fn factor<I>(ZZ: I, mut n: El<I>) -> Vec<(El<I>, usize)>
81where
82    I: RingStore + Copy,
83    I::Type: IntegerRing + OrderedRing + CanIsoFromTo<BigIntRingBase> + CanIsoFromTo<StaticRingBase<i128>>,
84{
85    const SMALL_PRIME_BOUND: i32 = 10000;
86    let mut result = Vec::new();
87
88    // first make it nonnegative
89    if ZZ.is_neg(&n) {
90        result.push((ZZ.neg_one(), 1));
91        ZZ.negate_inplace(&mut n);
92    }
93
94    // check for special cases
95    if ZZ.is_zero(&n) {
96        result.push((ZZ.zero(), 1));
97        return result;
98    }
99
100    // check if we are done
101    if ZZ.is_one(&n) {
102        return result;
103    } else if algorithms::miller_rabin::is_prime(ZZ, &n, DEFAULT_PROBABILISTIC_REPETITIONS) {
104        result.push((n, 1));
105        return result;
106    }
107
108    // then we remove small factors
109    for p in algorithms::erathostenes::enumerate_primes(StaticRing::<i32>::RING, &SMALL_PRIME_BOUND) {
110        let ZZ_p = ZZ.int_hom().map(p);
111        let mut count = 0;
112        while let Some(quo) = ZZ.checked_div(&n, &ZZ_p) {
113            n = quo;
114            count += 1;
115        }
116        if count >= 1 {
117            result.push((ZZ_p, count));
118        }
119    }
120
121    // check again if we are done
122    if ZZ.is_one(&n) {
123        return result;
124    } else if algorithms::miller_rabin::is_prime(ZZ, &n, DEFAULT_PROBABILISTIC_REPETITIONS) {
125        result.push((n, 1));
126        return result;
127    }
128
129    // then check for powers, as EC factor fails for those
130    if let Some((m, k)) = is_power(ZZ, &n) {
131        let mut power_factors = factor(ZZ, m);
132        for (_, multiplicity) in &mut power_factors {
133            *multiplicity *= k;
134        }
135        result.extend(power_factors);
136        return result;
137    }
138
139    // then we use EC factor to factor the result
140    let m = choose_zn_impl(ZZ, ZZ.clone_el(&n), ECFactorInt { result_ring: ZZ });
141
142    let mut factors1 = factor(ZZ, ZZ.checked_div(&n, &m).unwrap());
143    let mut factors2 = factor(ZZ, m);
144
145    // finally group the prime factors
146    factors1.sort_by(|(a, _), (b, _)| ZZ.cmp(a, b));
147    factors2.sort_by(|(a, _), (b, _)| ZZ.cmp(a, b));
148    let mut iter1 = factors1.into_iter().peekable();
149    let mut iter2 = factors2.into_iter().peekable();
150    loop {
151        match (iter1.peek(), iter2.peek()) {
152            (Some((p1, m1)), Some((p2, m2))) if ZZ.eq_el(p1, p2) => {
153                result.push((ZZ.clone_el(p1), m1 + m2));
154                _ = iter1.next().unwrap();
155                _ = iter2.next().unwrap();
156            }
157            (Some((p1, m1)), Some((p2, _m2))) if ZZ.is_lt(p1, p2) => {
158                result.push((ZZ.clone_el(p1), *m1));
159                _ = iter1.next().unwrap();
160            }
161            (Some((_p1, _m1)), Some((p2, m2))) => {
162                result.push((ZZ.clone_el(p2), *m2));
163                _ = iter2.next().unwrap();
164            }
165            (Some((p1, m1)), None) => {
166                result.push((ZZ.clone_el(p1), *m1));
167                _ = iter1.next().unwrap();
168            }
169            (None, Some((p2, m2))) => {
170                result.push((ZZ.clone_el(p2), *m2));
171                _ = iter2.next().unwrap();
172            }
173            (None, None) => {
174                return result;
175            }
176        }
177    }
178}
179
180#[test]
181fn test_factor() {
182    let ZZbig = BigIntRing::RING;
183    assert_eq!(
184        vec![(3, 2), (5, 1), (29, 1)],
185        factor(&StaticRing::<i64>::RING, 3 * 3 * 5 * 29)
186    );
187    assert_eq!(vec![(2, 8)], factor(&StaticRing::<i64>::RING, 256));
188    assert_eq!(vec![(1009, 2)], factor(&StaticRing::<i64>::RING, 1009 * 1009));
189    assert_eq!(vec![(0, 1)], factor(&StaticRing::<i64>::RING, 0));
190    assert_eq!(Vec::<(i64, usize)>::new(), factor(&StaticRing::<i64>::RING, 1));
191    assert_eq!(vec![(-1, 1)], factor(&StaticRing::<i64>::RING, -1));
192    assert_eq!(
193        vec![(257, 1), (1009, 2)],
194        factor(&StaticRing::<i128>::RING, 257 * 1009 * 1009)
195    );
196
197    let expected = vec![
198        (ZZbig.int_hom().map(-1), 1),
199        (ZZbig.int_hom().map(32771), 1),
200        (ZZbig.int_hom().map(65537), 1),
201    ];
202    let actual = factor(
203        &ZZbig,
204        ZZbig.mul(ZZbig.int_hom().map(-32771), ZZbig.int_hom().map(65537)),
205    );
206    assert_eq!(expected.len(), actual.len());
207    for ((expected_factor, expected_multiplicity), (actual_factor, actual_multiplicity)) in
208        expected.iter().zip(actual.iter())
209    {
210        assert_eq!(expected_multiplicity, actual_multiplicity);
211        assert!(ZZbig.eq_el(expected_factor, actual_factor));
212    }
213
214    let expected = vec![
215        (ZZbig.int_hom().map(257), 2),
216        (ZZbig.int_hom().map(32771), 1),
217        (ZZbig.int_hom().map(65537), 2),
218    ];
219    let actual = factor(
220        &ZZbig,
221        ZZbig.prod(
222            [
223                ZZbig.int_hom().map(257 * 257),
224                ZZbig.int_hom().map(32771),
225                ZZbig.int_hom().map(65537),
226                ZZbig.int_hom().map(65537),
227            ]
228            .into_iter(),
229        ),
230    );
231    assert_eq!(expected.len(), actual.len());
232    for ((expected_factor, expected_multiplicity), (actual_factor, actual_multiplicity)) in
233        expected.iter().zip(actual.iter())
234    {
235        assert_eq!(expected_multiplicity, actual_multiplicity);
236        assert!(ZZbig.eq_el(expected_factor, actual_factor));
237    }
238}
239
240#[test]
241fn test_is_prime_power() {
242    assert_eq!(Some((2, 6)), is_prime_power(&StaticRing::<i64>::RING, &64));
243}
244
245#[test]
246fn test_is_prime_power_large_n() {
247    assert_eq!(
248        Some((5, 25)),
249        is_prime_power(&StaticRing::<i64>::RING, &StaticRing::<i64>::RING.pow(5, 25))
250    );
251}
feanor_math/algorithms/int_factor.rs

feanor_math/algorithms/
int_factor.rs
