primal_check/
perfect_power.rs

1use num_integer::Integer;
2
3fn wrapping_pow(mut base: u64, mut exp: u32) -> u64 {
4    let mut acc: u64 = 1;
5    while exp > 0 {
6        if exp % 2 == 1 {
7            acc = acc.wrapping_mul(base)
8        }
9        base = base.wrapping_mul(base);
10        exp /= 2;
11    }
12    acc
13}
14
15/// Returns integers `(y, k)` such that `x = y^k` with `k` maximised
16/// (other than for `x = 0, 1`, in which case `y = x`, `k = 1`).
17///
18/// # Examples
19///
20/// ```rust
21/// # use primal_check as primal;
22/// assert_eq!(primal::as_perfect_power(2), (2, 1));
23/// assert_eq!(primal::as_perfect_power(4), (2, 2));
24/// assert_eq!(primal::as_perfect_power(8), (2, 3));
25/// assert_eq!(primal::as_perfect_power(1024), (2, 10));
26///
27/// assert_eq!(primal::as_perfect_power(1000), (10, 3));
28///
29/// assert_eq!(primal::as_perfect_power(15), (15, 1));
30/// ```
31pub fn as_perfect_power(x: u64) -> (u64, u8) {
32    if x == 0 || x == 1 {
33        return (x, 1)
34    }
35
36    let floor_log_2 = 64 - x.leading_zeros() as u32 - 1;
37
38    let x_ = x as f64;
39    let mut last = (x, 1);
40    // TODO: we could be smarter about this: we know all the possible
41    // primes that can divide the exponent (since we have a list up to
42    // 251 >= 64), so we really only need to check them.
43    let mut expn: u32 = 2;
44    let mut step = 1;
45    while expn <= floor_log_2 {
46        let factor = x_.powf(1.0/expn as f64).round() as u64;
47        // the only case this will wrap is if x is close to 2^64 and
48        // the round() rounds up, pushing this calculation over the
49        // edge, however, the overflow will be well away from x, so we
50        // still correctly don't take this branch. (x can't be a
51        // perfect power if the result rounds away.)
52        if wrapping_pow(factor, expn) == x {
53            last = (factor, expn as u8);
54            // if x is a 2nd and 5th power, it's going to be a 10th
55            // power too, meaning we can search faster.
56            // TODO: check if this is actually saving work
57            step = step.lcm(&expn);
58        }
59
60        expn += step;
61    }
62    last
63}
64
65/// Return `Some((p, k))` if `x = p^k` for some prime `p` and `k >= 1`
66/// (that is, including when `x` is itself a prime).
67///
68/// Returns `None` if `x` not a perfect power.
69///
70/// # Examples
71///
72/// ```rust
73/// # use primal_check as primal;
74/// assert_eq!(primal::as_prime_power(2), Some((2, 1)));
75/// assert_eq!(primal::as_prime_power(4), Some((2, 2)));
76/// assert_eq!(primal::as_prime_power(8), Some((2, 3)));
77/// assert_eq!(primal::as_prime_power(1024), Some((2, 10)));
78///
79/// assert_eq!(primal::as_prime_power(1000), None);
80///
81/// assert_eq!(primal::as_prime_power(15), None);
82/// ```
83pub fn as_prime_power(x: u64) -> Option<(u64, u8)> {
84    let (y, k) = as_perfect_power(x);
85    if crate::miller_rabin(y) {
86        Some((y, k))
87    } else {
88        None
89    }
90}
91
92#[cfg(test)]
93mod tests {
94    use primal::Sieve;
95
96    use super::{as_perfect_power, as_prime_power};
97
98    #[test]
99    fn perfect_and_prime_power() {
100        let tests = [
101            (0, (0, 1), false),
102            (1, (1, 1), false),
103            (2, (2, 1), true),
104            (3, (3, 1), true),
105            (4, (2, 2), true),
106            (5, (5, 1), true),
107            (6, (6, 1), false),
108            (8, (2, 3), true),
109            (9, (3, 2), true),
110            (16, (2, 4), true),
111            (25, (5, 2), true),
112            (32, (2, 5), true),
113            (36, (6, 2), false),
114            (100, (10, 2), false),
115            (1000, (10, 3), false),
116                ];
117
118        for &(x, expected, is_prime) in tests.iter() {
119            assert_eq!(as_perfect_power(x), expected);
120            assert_eq!(as_prime_power(x),
121                       if is_prime { Some(expected)} else { None })
122        }
123
124        let sieve = Sieve::new(200);
125        let mut primes = sieve.primes_from(0);
126        const MAX: f64 = 0xFFFF_FFFF_FFFF_FFFFu64 as f64;
127        // test a whole pile of (semi)primes
128        loop {
129            let p = match primes.next() {
130                Some(p) => p as u64,
131                None => break
132            };
133
134            let subprimes = primes.clone().map(|x| (x, false));
135            // include 1 to test p itself.
136            for (q, is_prime) in Some((1, true)).into_iter().chain(subprimes) {
137                let pq = p * q as u64;
138                for n in 1..(MAX.log(pq as f64) as u32) {
139                    let x = pq.pow(n);
140
141                    let expected = (pq, n as u8);
142                    assert_eq!(as_perfect_power(x), expected);
143                    assert_eq!(as_prime_power(x),
144                               if is_prime { Some(expected) } else { None });
145                }
146            }
147        }
148    }
149}
primal_check/perfect_power.rs

primal_check/
perfect_power.rs
