feanor_math/algorithms/
int_factor.rs

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

feanor_math/algorithms/
int_factor.rs
