ognlib/algorithm/
prime.rs

1//! Primality tests. These tests are divided into 2 major groups: first group of tests gives exact
2//! results; second group is for probabilistic tests, which can only suppose whether number is prime
3//! or not.
4
5#[cfg(feature = "std")] use rayon::prelude::*;
6
7use {
8    crate::num::power::modpow, fastrand::Rng, integer_sqrt::IntegerSquareRoot, num_bigint::BigUint,
9    snafu::Snafu,
10};
11
12/// If number is less than 2, we can't say that number is either prime or composite.
13#[non_exhaustive]
14#[derive(Debug, Snafu, PartialEq, Eq, Clone, Copy)]
15#[snafu(display("This number is neither prime nor composite"))]
16pub struct PrimeStatusError;
17
18#[non_exhaustive]
19#[derive(Debug, PartialEq, Eq, Clone, Copy)]
20pub enum PrimeStatus {
21    Prime,
22    Composite,
23    ProbablyPrime,
24}
25
26impl PrimeStatus {
27    /// Returns `true` if [`PrimeStatus`] is [`Prime`].
28    /// # Examples
29    ///
30    /// ```
31    /// use ognlib::algorithm::prime::{PrimeStatus, wilson_th};
32    ///
33    /// assert!(wilson_th(13usize).unwrap().is_prime());
34    /// assert!(!wilson_th(455usize).unwrap().is_prime());
35    /// ```
36    ///
37    /// [`Prime`]: PrimeStatus::Prime
38    #[inline]
39    #[must_use]
40    pub fn is_prime(self) -> bool { self == Self::Prime }
41
42    /// Returns `true` if [`PrimeStatus`] is not [`Composite`].
43    /// # Examples
44    ///
45    /// ```
46    /// use ognlib::algorithm::prime::{PrimeStatus, miller_rabin};
47    ///
48    /// assert!(miller_rabin(13usize).unwrap().is_probably_prime());
49    /// assert!(miller_rabin(7usize).unwrap().is_probably_prime());
50    /// ```
51    ///
52    /// [`Composite`]: PrimeStatus::Composite
53    #[inline]
54    #[must_use]
55    pub fn is_probably_prime(self) -> bool { self != Self::Composite }
56
57    /// Returns `true` if [`PrimeStatus`] is [`Composite`].
58    /// # Examples
59    ///
60    /// ```
61    /// use ognlib::algorithm::prime::{PrimeStatus, wilson_th};
62    ///
63    /// assert!(!wilson_th(13usize).unwrap().is_composite());
64    /// assert!(wilson_th(455usize).unwrap().is_composite());
65    /// ```
66    ///
67    /// [`Composite`]: PrimeStatus::Composite
68    #[inline]
69    #[must_use]
70    pub fn is_composite(self) -> bool { self == Self::Composite }
71}
72
73/// Methods to check prime status.
74pub trait Prime {
75    /// Returns `true` if number is prime.
76    fn is_prime(&self) -> bool;
77
78    /// Returns `true` if number is either prime or probably prime.
79    fn is_probably_prime(&self) -> bool;
80
81    /// Returns `true` if number is composite.
82    fn is_composite(&self) -> bool;
83}
84
85impl Prime for usize {
86    /// Returns `true` if number is prime.
87    /// # Examples
88    ///
89    /// ```
90    /// use ognlib::algorithm::prime::Prime;
91    ///
92    /// assert!(13usize.is_prime());
93    /// assert!(!455usize.is_prime());
94    /// ```
95    #[inline]
96    fn is_prime(&self) -> bool { wilson_th(*self) == Ok(PrimeStatus::Prime) }
97
98    /// Returns `true` if number is either prime or probably prime.
99    /// # Examples
100    ///
101    /// ```
102    /// use ognlib::algorithm::prime::Prime;
103    ///
104    /// assert!(13usize.is_probably_prime());
105    /// assert!(7usize.is_probably_prime());
106    /// ```
107    #[inline]
108    fn is_probably_prime(&self) -> bool { miller_rabin(*self) != Ok(PrimeStatus::Composite) }
109
110    /// Returns `true` if number is composite.
111    /// # Examples
112    ///
113    /// ```
114    /// use ognlib::algorithm::prime::Prime;
115    ///
116    /// assert!(!13usize.is_composite());
117    /// assert!(455usize.is_composite());
118    /// ```
119    #[inline]
120    fn is_composite(&self) -> bool { wilson_th(*self) == Ok(PrimeStatus::Composite) }
121}
122
123impl Prime for Result<PrimeStatus, PrimeStatusError> {
124    /// Applied to result of prime test.
125    ///
126    /// Returns true if test is successful and value under result is [`PrimeStatus::Prime`].
127    /// # Examples
128    ///
129    /// ```
130    /// use ognlib::algorithm::prime::{Prime, sqrtest};
131    ///
132    /// assert!(sqrtest(13).is_prime());
133    /// assert!(!sqrtest(455).is_prime());
134    /// ```
135    #[inline]
136    fn is_prime(&self) -> bool { self.is_ok_and(|st| st == PrimeStatus::Prime) }
137
138    /// Applied to result of prime test.
139    ///
140    /// Returns true if test is successful and value under result is [`PrimeStatus::ProbablyPrime`].
141    /// # Examples
142    ///
143    /// ```
144    /// use ognlib::algorithm::prime::{Prime, miller_rabin};
145    ///
146    /// assert!(miller_rabin(13).is_probably_prime());
147    /// assert!(miller_rabin(7).is_probably_prime());
148    /// ```
149    #[inline]
150    fn is_probably_prime(&self) -> bool { self.is_ok_and(|st| st == PrimeStatus::ProbablyPrime) }
151
152    /// Applied to result of prime test.
153    ///
154    /// Returns true if test is successful and value under result is [`PrimeStatus::Composite`].
155    /// # Examples
156    ///
157    /// ```
158    /// use ognlib::algorithm::prime::{Prime, sqrtest};
159    ///
160    /// assert!(!sqrtest(13).is_composite());
161    /// assert!(sqrtest(455).is_composite());
162    /// ```
163    #[inline]
164    fn is_composite(&self) -> bool { self.is_ok_and(|st| st == PrimeStatus::Composite) }
165}
166
167/// Simple prime test.
168///
169/// Prime test that takes ceil of sqrt(n) as upper bound and checks if there is any divisor from 3
170/// to ceil with step 2.
171///
172/// # Errors
173/// Returns a [`PrimeStatusError`] if number <= 1
174///
175/// # Examples
176/// ```
177/// use ognlib::algorithm::prime::{PrimeStatus, PrimeStatusError, sqrtest};
178///
179/// assert_eq!(
180///     sqrtest(1usize).unwrap_err().to_string(),
181///     "This number is neither prime nor composite",
182/// );
183/// assert_eq!(sqrtest(13usize).ok(), Some(PrimeStatus::Prime));
184/// assert_eq!(sqrtest(455usize).ok(), Some(PrimeStatus::Composite));
185/// ```
186pub fn sqrtest(num: usize) -> Result<PrimeStatus, PrimeStatusError> {
187    if num < 2 {
188        Err(PrimeStatusError)
189    } else {
190        // FIXME: https://github.com/rust-lang/rust/issues/116226
191        // let sqrt_res = num.isqrt() + 1;
192        let sqrt_res = num.integer_sqrt() + 1;
193        cfg_if::cfg_if! {
194            if #[cfg(feature = "std")] {
195                let cond = (3..=sqrt_res).into_par_iter().find_any(|&i| num % i == 0).is_some();
196            } else {
197                let cond = (3..=sqrt_res).any(|i| num & i == 0);
198            }
199        };
200        if cond { Ok(PrimeStatus::Composite) } else { Ok(PrimeStatus::Prime) }
201    }
202}
203
204/// Wilson's theory.
205///
206/// From [Wikipedia](https://en.wikipedia.org/wiki/Wilson%27s_theorem): "Wilson's theorem states
207/// that a natural number n > 1 is a prime number if and only if the product of all the positive
208/// integers less than n is one less than a multiple of n. That is the factorial (n - 1)! satisfies
209/// (n - 1)! % n == -1".
210///
211/// # Errors
212/// Returns a [`PrimeStatusError`] if number <= 1
213///
214/// # Examples
215/// ```
216/// use ognlib::algorithm::prime::{PrimeStatus, PrimeStatusError, wilson_th};
217///
218/// assert_eq!(
219///     wilson_th(1usize).unwrap_err().to_string(),
220///     "This number is neither prime nor composite",
221/// );
222/// assert_eq!(wilson_th(13usize).ok(), Some(PrimeStatus::Prime));
223/// assert_eq!(wilson_th(455usize).ok(), Some(PrimeStatus::Composite));
224/// ```
225pub fn wilson_th(num: usize) -> Result<PrimeStatus, PrimeStatusError> {
226    if num < 2 {
227        return Err(PrimeStatusError);
228    }
229
230    let mut fact = BigUint::from(1u8);
231    for i in 2..num {
232        fact = (fact * i) % num;
233    }
234
235    if fact + 1u8 == BigUint::from(num) {
236        Ok(PrimeStatus::Prime)
237    } else {
238        Ok(PrimeStatus::Composite)
239    }
240}
241
242/// Miller-Rabin's prime test.
243///
244/// From
245/// [Wikipedia](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test): the Miller–Rabin
246/// primality test or Rabin–Miller primality test is a probabilistic primality test: an algorithm
247/// which determines whether a given number is likely to be prime.
248///
249/// NOTE:
250/// This function is *not* cryptographically safe
251///
252/// # Errors
253/// Returns a [`PrimeStatusError`] if number <= 1
254///
255/// # Examples
256/// ```
257/// use ognlib::algorithm::prime::{PrimeStatus, PrimeStatusError, miller_rabin};
258///
259/// assert_eq!(
260///     miller_rabin(1usize).unwrap_err().to_string(),
261///     "This number is neither prime nor composite",
262/// );
263/// assert_eq!(miller_rabin(13usize).ok(), Some(PrimeStatus::ProbablyPrime));
264/// assert_eq!(miller_rabin(455usize).ok(), Some(PrimeStatus::Composite));
265/// ```
266pub fn miller_rabin(num: usize) -> Result<PrimeStatus, PrimeStatusError> {
267    match num {
268        0 | 1 => Err(PrimeStatusError),
269        5 => Ok(PrimeStatus::Prime),
270        _ if num % 2 == 0 || num % 3 == 0 => Ok(PrimeStatus::Composite),
271        _ => {
272            let log_sqr = (num.ilog2() * num.ilog2()).into();
273            let (mut temp, mut su) = (num - 1, 0);
274            while temp % 2 == 0 {
275                temp /= 2;
276                su += 1;
277            }
278            'outer: for i in 0..log_sqr {
279                let rand_num = Rng::with_seed(i).usize(2..num - 1);
280                let mut x_num = modpow(rand_num, temp, num);
281                if x_num == 1 || x_num == num - 1 {
282                    continue;
283                }
284                for _j in 0..su - 1 {
285                    x_num = modpow(x_num, 2, num);
286                    if x_num == 1 {
287                        return Ok(PrimeStatus::Composite);
288                    }
289                    if x_num == num - 1 {
290                        continue 'outer;
291                    }
292                }
293                return Ok(PrimeStatus::Composite);
294            }
295            Ok(PrimeStatus::ProbablyPrime)
296        },
297    }
298}