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
43pub fn is_prime_power<I>(ZZ: I, n: &El<I>) -> Option<(El<I>, usize)>
44    where I: RingStore + Copy,
45        I::Type: IntegerRing
46{
47    if algorithms::miller_rabin::is_prime(ZZ, n, DEFAULT_PROBABILISTIC_REPETITIONS) {
48        return Some((ZZ.clone_el(n), 1));
49    }
50    let (p, e) = is_power(ZZ, n)?;
51    if algorithms::miller_rabin::is_prime(ZZ, &p, DEFAULT_PROBABILISTIC_REPETITIONS) {
52        return Some((p, e));
53    } else {
54        return None;
55    }
56}
57
58fn is_power<I>(ZZ: I, n: &El<I>) -> Option<(El<I>, usize)>
59    where I: RingStore + Copy,
60        I::Type: IntegerRing
61{
62    assert!(!ZZ.is_zero(n));
63    for i in (2..=ZZ.abs_log2_ceil(n).unwrap()).rev() {
64        let root = algorithms::int_bisect::root_floor(ZZ, ZZ.clone_el(n), i);
65        if ZZ.eq_el(&ZZ.pow(root, i), n) {
66            return Some((algorithms::int_bisect::root_floor(ZZ, ZZ.clone_el(n), i), i));
67        }
68    }
69    return None;
70}
71
72pub fn factor<I>(ZZ: I, mut n: El<I>) -> Vec<(El<I>, usize)> 
73    where I: RingStore + Copy, 
74        I::Type: IntegerRing + OrderedRing + CanIsoFromTo<BigIntRingBase> + CanIsoFromTo<StaticRingBase<i128>>
75{
76    const SMALL_PRIME_BOUND: i32 = 10000;
77    let mut result = Vec::new();
78
79    // first make it nonnegative
80    if ZZ.is_neg(&n) {
81        result.push((ZZ.neg_one(), 1));
82        ZZ.negate_inplace(&mut n);
83    }
84
85    // check for special cases
86    if ZZ.is_zero(&n) {
87        result.push((ZZ.zero(), 1));
88        return result;
89    }
90
91    // check if we are done
92    if ZZ.is_one(&n) {
93        return result;
94    } else if algorithms::miller_rabin::is_prime(ZZ, &n, DEFAULT_PROBABILISTIC_REPETITIONS) {
95        result.push((n, 1));
96        return result;
97    }
98
99    // then we remove small factors
100    for p in algorithms::erathostenes::enumerate_primes(StaticRing::<i32>::RING, &SMALL_PRIME_BOUND) {
101        let ZZ_p = ZZ.int_hom().map(p);
102        let mut count = 0;
103        while let Some(quo) = ZZ.checked_div(&n, &ZZ_p) {
104            n = quo;
105            count += 1;
106        }
107        if count >= 1 {
108            result.push((ZZ_p, count));
109        }
110    }
111
112    // check again if we are done
113    if ZZ.is_one(&n) {
114        return result;
115    } else if algorithms::miller_rabin::is_prime(ZZ, &n, DEFAULT_PROBABILISTIC_REPETITIONS) {
116        result.push((n, 1));
117        return result;
118    }
119
120    // then check for powers, as EC factor fails for those
121    if let Some((m, k)) = is_power(ZZ, &n) {
122        let mut power_factors = factor(ZZ, m);
123        for (_, multiplicity) in &mut power_factors {
124            *multiplicity *= k;
125        }
126        result.extend(power_factors.into_iter());
127        return result;
128    }
129
130    // then we use EC factor to factor the result
131    let m = choose_zn_impl(ZZ, ZZ.clone_el(&n), ECFactorInt { result_ring: ZZ });
132
133    let mut factors1 = factor(ZZ, ZZ.checked_div(&n, &m).unwrap());
134    let mut factors2 = factor(ZZ, m);
135
136    // finally group the prime factors
137    factors1.sort_by(|(a, _), (b, _)| ZZ.cmp(a, b));
138    factors2.sort_by(|(a, _), (b, _)| ZZ.cmp(a, b));
139    let mut iter1 = factors1.into_iter().peekable();
140    let mut iter2 = factors2.into_iter().peekable();
141    loop {
142        match (iter1.peek(), iter2.peek()) {
143            (Some((p1, m1)), Some((p2, m2))) if ZZ.eq_el(p1, p2) => {
144                result.push((ZZ.clone_el(p1), m1 + m2));
145                _ = iter1.next().unwrap();
146                _ = iter2.next().unwrap();
147            },
148            (Some((p1, m1)), Some((p2, _m2))) if ZZ.is_lt(p1, p2) => {
149                result.push((ZZ.clone_el(p1), *m1));
150                _ = iter1.next().unwrap();
151            },
152            (Some((_p1, _m1)), Some((p2, m2))) => {
153                result.push((ZZ.clone_el(p2), *m2));
154                _ = iter2.next().unwrap();
155            },
156            (Some((p1, m1)), None) => {
157                result.push((ZZ.clone_el(p1), *m1));
158                _ = iter1.next().unwrap();
159            },
160            (None, Some((p2, m2))) => {
161                result.push((ZZ.clone_el(p2), *m2));
162                _ = iter2.next().unwrap();
163            },
164            (None, None) => {
165                return result;
166            }
167        }
168    }
169}
170
171#[test]
172fn test_factor() {
173    let ZZbig = BigIntRing::RING;
174    assert_eq!(vec![(3, 2), (5, 1), (29, 1)], factor(&StaticRing::<i64>::RING, 3 * 3 * 5 * 29));
175    assert_eq!(vec![(2, 8)], factor(&StaticRing::<i64>::RING, 256));
176    assert_eq!(vec![(1009, 2)], factor(&StaticRing::<i64>::RING, 1009 * 1009));
177    assert_eq!(vec![(0, 1)], factor(&StaticRing::<i64>::RING, 0));
178    assert_eq!(Vec::<(i64, usize)>::new(), factor(&StaticRing::<i64>::RING, 1));
179    assert_eq!(vec![(-1, 1)], factor(&StaticRing::<i64>::RING, -1));
180    assert_eq!(vec![(257, 1), (1009, 2)], factor(&StaticRing::<i128>::RING, 257 * 1009 * 1009));
181
182    let expected = vec![(ZZbig.int_hom().map(-1), 1), (ZZbig.int_hom().map(32771), 1), (ZZbig.int_hom().map(65537), 1)];
183    let actual = factor(&ZZbig, ZZbig.mul(ZZbig.int_hom().map(-32771), ZZbig.int_hom().map(65537)));
184    assert_eq!(expected.len(), actual.len());
185    for ((expected_factor, expected_multiplicity), (actual_factor, actual_multiplicity)) in expected.iter().zip(actual.iter()) {
186        assert_eq!(expected_multiplicity, actual_multiplicity);
187        assert!(ZZbig.eq_el(expected_factor, actual_factor));
188    }
189
190    let expected = vec![(ZZbig.int_hom().map(257), 2), (ZZbig.int_hom().map(32771), 1), (ZZbig.int_hom().map(65537), 2)];
191    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()));
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
199#[test]
200fn test_is_prime_power() {
201    assert_eq!(Some((2, 6)), is_prime_power(&StaticRing::<i64>::RING, &64));
202}
203
204#[test]
205fn test_is_prime_power_large_n() {
206    assert_eq!(Some((5, 25)), is_prime_power(&StaticRing::<i64>::RING, &StaticRing::<i64>::RING.pow(5, 25)));
207}
feanor_math/algorithms/int_factor.rs

feanor_math/algorithms/
int_factor.rs
