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}