diff --git a/.gitignore b/.gitignore index 4fffb2f..2ebc5ea 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,2 @@ /target -/Cargo.lock +/Cargo.lock \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml index 0544736..d8f4fbf 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,17 +4,19 @@ authors = ["Johanna Sörngård (jsorngard@gmail.com)"] version = "0.4.8" edition = "2021" license = "MIT OR Apache-2.0" +keywords = ["primes", "const"] +categories = ["mathematics", "no-std::no-alloc"] description = "Generate and work with prime numbers in const contexts" repository = "https://github.com/JSorngard/const-primes" -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - -[dependencies] - [dev-dependencies] -criterion = {version = "0.5", features = ["html_reports"]} +criterion = { version = "0.5", features = ["html_reports"] } rand = "0.8" +[features] +std = [ "alloc"] +alloc = [] + [[bench]] name = "prime_benches" -harness = false \ No newline at end of file +harness = false diff --git a/README.md b/README.md index ce4b1cd..d7b73dd 100644 --- a/README.md +++ b/README.md @@ -9,13 +9,13 @@ A crate for generating and working with prime numbers in const contexts. `#![no_std]` compatible. ## Examples -Generate arrays of prime numbers with the function `primes` which uses a [segmented sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve). +Generate arrays of prime numbers with the function `primes` which uses a [segmented sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve): ```rust const PRIMES: [u32; 10] = primes(); assert_eq!(PRIMES[5], 13); assert_eq!(PRIMES, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]); ``` -or with the type `Primes` which ensures that a non-zero number of primes are generated: +or with the wrapping type [`Primes`]: ```rust const PRIMES: Primes<10> = Primes::new(); assert_eq!(PRIMES[5], 13); @@ -40,30 +40,54 @@ assert_eq!(PRIMES_LEQ_100, Some(25)); assert!(CACHE.is_prime(1000).is_none()); assert!(CACHE.count_primes_leq(1000).is_none()); ``` -Creating a `Primes<0>` is a compile fail in const contexts and a panic otherwise. - -### Other functionality -Use `is_prime` to test whether a given number is prime -```rust -const CHECK: bool = is_prime(18_446_744_073_709_551_557); -assert!(CHECK); -``` -or `are_prime` to compute the prime status of the `N` first integers, +Sieve a range of numbers for their prime status with `sieve`: ```rust const N: usize = 10; -const PRIME_STATUS: [bool; N] = are_prime(); +const PRIME_STATUS: [bool; N] = sieve(); // 0 1 2 3 4 5 6 7 8 9 assert_eq!(PRIME_STATUS, [false, false, true, true, false, true, false, true, false, false]); +``` + +## Arbitrary ranges +The crate provides prime generation and sieving functions with suffixes, e.g. `primes_geq` and `sieve_lt` +that can be used to work with ranges that don't start at zero. They take two generics: the number of elements +to store in the binary and the size of the sieve used during evaluation. The sieve size must be the cieling +of the square root of the largest encountered value: +```rust +// ceil(sqrt(5_000_000_063)) = 70_711 +const PRIMES_GEQ: const_primes::Result<3> = primes_geq::<3, 70_711>(5_000_000_031); +assert_eq!(PRIMES_GEQ?, [5_000_000_039, 5_000_000_059, 5_000_000_063]); ``` -or `are_prime_below` to compute the prime status of the `N` largest integers below a given value, ```rust const N: usize = 70711; -const BIG_PRIME_STATUS: [bool; N] = are_prime_below(5_000_000_031); -// 5_000_000_028 5_000_000_029 5_000_000_030 -assert_eq!(BIG_PRIME_STATUS[N - 3..], [false, true, false]); +const PRIME_STATUS_LT: [bool; N] = sieve_lt(5_000_000_031); +// 5_000_000_028 5_000_000_029 5_000_000_030 +assert_eq!(PRIME_STATUS_LT[N - 3..], [false, true, false]); +``` +The sieving functions have yet to be modified for two generics, and must save the entire sieve in the binary. +## Other functionality +Use `is_prime` to test whether a given number is prime: +```rust +const CHECK: bool = is_prime(18_446_744_073_709_551_557); +assert!(CHECK); +``` +Find the next or previous prime numbers with `next_prime` and `previous_prime` if they exist: +```rust +const NEXT: Option = next_prime(25); +const PREV: Option = previous_prime(25); +const NOSUCH: Option = previous_prime(2); + +assert_eq!(NEXT, Some(29)); +assert_eq!(PREV, Some(23)); +assert_eq!(NOSUCH, None); ``` and more! +## Features + +`std`: derives the `Error` trait for the error types. +`alloc`: enables conversion of the type returned by `primes_geq` and `primes_lt` into `Vec`s and `Box`ed slices. + ## License Licensed under either of diff --git a/benches/prime_benches.rs b/benches/prime_benches.rs index 29eaa95..3794bb1 100644 --- a/benches/prime_benches.rs +++ b/benches/prime_benches.rs @@ -1,37 +1,54 @@ -use const_primes::{are_prime, is_prime, moebius, primes}; -use criterion::{criterion_group, criterion_main, BatchSize, Criterion}; +use const_primes::{is_prime, primes, primes_geq, primes_lt, sieve, sieve_geq, sieve_lt}; +use criterion::{criterion_group, criterion_main, BatchSize, Criterion, Throughput}; use rand::prelude::*; use std::hint::black_box; fn benchmarks(c: &mut Criterion) { - c.bench_function("generate 10000 primes", |b| { - b.iter(|| black_box(primes::<10_000>())) - }); + { + const N: usize = 10_000; + let mut prime_generation = c.benchmark_group("prime generation"); + prime_generation.bench_function(format!("first {N} primes"), |b| { + b.iter(|| black_box(primes::())) + }); + prime_generation.bench_function(format!("{N} primes < 100000000"), |b| { + b.iter(|| black_box(primes_lt::(100000000))) + }); + prime_generation.bench_function(format!("{N} primes >= 99990000"), |b| { + b.iter(|| black_box(primes_geq::(99990000))) + }); + } - c.bench_function("is_prime on random numbers", |b| { - b.iter_batched( - || (0..10_000).map(|_| random()).collect::>(), - |data| { - for number in data.iter() { - black_box(is_prime(*number)); - } - }, - BatchSize::SmallInput, - ) - }); + { + const N: u64 = 10_000; + let mut rng = StdRng::seed_from_u64(1234567890); + let mut primality_testing = c.benchmark_group("primality testing"); + primality_testing.throughput(Throughput::Elements(N)); + primality_testing.bench_function(format!("is_prime on {N} random numbers"), |b| { + b.iter_batched( + || (0..N).map(|_| rng.gen()).collect::>(), + |data| { + for number in data.iter() { + black_box(is_prime(*number)); + } + }, + BatchSize::SmallInput, + ) + }); + } - c.bench_function("sieve 10000 integers", |b| { - b.iter(|| black_box(are_prime::<10_000>())) - }); - - let ints: Vec<_> = (1..1_000_000).map(|n| n).collect(); - c.bench_function("möbius of first 1e6 integers", |b| { - b.iter(|| { - for &i in &ints { - black_box(moebius(i)); - } - }) - }); + { + const N: usize = 10_000; + let mut sieving = c.benchmark_group("prime sieving"); + sieving.bench_function(format!("first {N} integers"), |b| { + b.iter(|| black_box(sieve::())) + }); + sieving.bench_function(format!("{N} integers < 100000000"), |b| { + b.iter(|| black_box(sieve_lt::(100000000))) + }); + sieving.bench_function(format!("{N} integers >= 99990000"), |b| { + b.iter(|| black_box(sieve_geq::(99990000))) + }); + } } criterion_group!(benches, benchmarks); diff --git a/rust-toolchain.toml b/rust-toolchain.toml new file mode 100644 index 0000000..271800c --- /dev/null +++ b/rust-toolchain.toml @@ -0,0 +1,2 @@ +[toolchain] +channel = "nightly" \ No newline at end of file diff --git a/src/generation.rs b/src/generation.rs new file mode 100644 index 0000000..f773702 --- /dev/null +++ b/src/generation.rs @@ -0,0 +1,435 @@ +use core::fmt; + +use crate::{sieve, sieving::sieve_segment, Underlying}; + +/// Returns the `N` first prime numbers. +/// Fails to compile if `N` is 0. +/// +/// [`Primes`](crate::Primes) might be relevant for you if you intend to later use these prime numbers for related computations. +/// +/// Uses a [segmented sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve). +/// +/// # Example +/// +/// ``` +/// # use const_primes::primes; +/// const PRIMES: [u32; 10] = primes(); +/// assert_eq!(PRIMES, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]); +/// ``` +/// Fails to compile if `N = 0`: +/// ```compile_fail +/// # use const_primes::primes; +/// let primes: [u32; 0] = primes(); +/// ``` +/// +#[must_use = "the function only returns a new value"] +pub const fn primes() -> [Underlying; N] { + const { assert!(N > 0, "`N` must be at least 1") } + + if N == 1 { + return [2; N]; + } else if N == 2 { + let mut primes = [0; N]; + primes[0] = 2; + primes[1] = 3; + return primes; + } + + // This is a segmented sieve that runs until it has found enough primes. + + // This array is the output in the end + let mut primes = [0; N]; + // This keeps track of how many primes we've found so far. + let mut prime_count = 0; + + // Sieve the first primes below N + let mut sieve: [bool; N] = sieve(); + + // Count how many primes we found + // and store them in the final array + let mut number = 0; + while number < N { + if sieve[number] { + primes[prime_count] = number as Underlying; + prime_count += 1; + } + + number += 1; + } + + // For every segment of N numbers + let mut low = N - 1; + let mut high = 2 * N - 1; + 'generate: while prime_count < N { + // reset the sieve for the segment + sieve = [true; N]; + let mut i = 0; + + // and repeat for each prime found so far: + while i < prime_count { + let prime = primes[i] as usize; + + // Find the smallest composite in the current segment, + let mut composite = (low / prime) * prime; + if composite < low { + composite += prime; + } + + // and sieve all numbers in the segment that are multiples of the prime. + while composite < high { + sieve[composite - low] = false; + composite += prime; + } + + i += 1; + } + + // Move the found primes into the final array + i = low; + while i < high { + if sieve[i - low] { + primes[prime_count] = i as Underlying; + prime_count += 1; + // and stop the generation of primes if we're done. + if prime_count >= N { + break 'generate; + } + } + i += 1; + } + + // Update low and high for the next segment + low += N; + high += N; + } + + primes +} + +/// Returns the `N` largest primes less than `upper_limit`. +/// +/// This function uses a segmented sieve of size `MEM` for computation, +/// but only saves the `N` requested primes in the binary. +/// +/// Set `MEM` such that `MEM*MEM >= upper_limit`. +/// +/// Fails to compile if `N` or `MEM` is 0, if `MEM < N` or if `MEM`^2 does not fit in a u64. +/// +/// If you want to compute primes that are larger than some limit, take a look at [`primes_geq`]. +/// +/// # Example +/// +/// Basic usage: +/// ``` +/// # use const_primes::{primes_lt, Error}; +/// // Sieving up to 100 means the sieve needs to be of size ceil(sqrt(100)) = 10. +/// // However, we only save the 4 largest primes in the constant. +/// const PRIMES: [u64;4] = match primes_lt::<4, 10>(100) {Ok(ps) => ps, Err(_) => panic!()}; +/// assert_eq!(PRIMES, [79, 83, 89, 97]); +/// ``` +/// Compute larger primes: +/// ``` +/// # use const_primes::{primes_lt, Error}; +/// use const_primes::isqrt; +/// const LIMIT: u64 = 5_000_000_030; +/// # #[allow(long_running_const_eval)] +/// const BIG_PRIMES: Result<[u64; 3], Error> = primes_lt::<3, {isqrt(LIMIT) as usize + 1}>(LIMIT); +/// +/// assert_eq!(BIG_PRIMES, Ok([4_999_999_903, 4_999_999_937, 5_000_000_029])); +/// ``` +/// # Errors +/// +/// If the number of primes requested, `N`, is larger than +/// the number of primes that exists below the `upper_limit` we +/// will get an error: +/// ``` +/// # use const_primes::{primes_lt, Error}; +/// const N: usize = 9; +/// const PRIMES: Result<[u64; N], Error> = primes_lt::(10); +/// assert_eq!(PRIMES, Err(Error::OutOfPrimes)); +/// ``` +/// +/// We will also get an error if `upper_limit` is larger than `MEM`^2 or if `upper_limit` is smaller than or equal to 2. +/// ``` +/// # use const_primes::{primes_lt, Error}; +/// const TOO_LARGE_LIMIT: Result<[u64; 3], Error> = primes_lt::<3, 5>(26); +/// const TOO_SMALL_LIMIT: Result<[u64 ;1], Error> = primes_lt::<1, 1>(1); +/// assert!(matches!(TOO_LARGE_LIMIT, Err(Error::TooLargeLimit(_, _)))); +/// assert!(matches!(TOO_SMALL_LIMIT, Err(Error::TooSmallLimit(_)))); +/// ``` +#[must_use = "the function only returns a new value and does not modify its input"] +pub const fn primes_lt( + mut upper_limit: u64, +) -> Result<[u64; N], Error> { + const { + assert!(N > 0, "`N` must be at least 1"); + assert!(MEM >= N, "`MEM` must be at least as large as `N`"); + } + let mem_sqr = const { + let mem64 = MEM as u64; + match mem64.checked_mul(mem64) { + Some(prod) => prod, + None => panic!("`MEM`^2 must fit in a u64"), + } + }; + + if upper_limit <= 2 { + return Err(Error::TooSmallLimit(upper_limit)); + } + + if upper_limit > mem_sqr { + return Err(Error::TooLargeLimit(upper_limit, mem_sqr)); + } + + let mut primes: [u64; N] = [0; N]; + + // This will be used to sieve all upper ranges. + let base_sieve: [bool; MEM] = sieve(); + + let mut total_primes_found: usize = 0; + 'generate: while total_primes_found < N { + // This is the smallest prime we have found so far. + let mut smallest_found_prime = primes[N - 1 - total_primes_found]; + // Sieve for primes in the segment. + let upper_sieve: [bool; MEM] = sieve_segment(&base_sieve, upper_limit); + + let mut i: usize = 0; + while i < MEM { + // Iterate backwards through the upper sieve. + if upper_sieve[MEM - 1 - i] { + smallest_found_prime = upper_limit - 1 - i as u64; + // Write every found prime to the primes array. + primes[N - 1 - total_primes_found] = smallest_found_prime; + total_primes_found += 1; + if total_primes_found >= N { + // If we have found enough primes we stop sieving. + break 'generate; + } + } + i += 1; + } + upper_limit = smallest_found_prime; + if upper_limit <= 2 && total_primes_found < N { + return Err(Error::OutOfPrimes); + } + } + + Ok(primes) +} + +/// Call [`primes_geq`] or [`primes_lt`], and automatically compute the memory requirement. +/// +/// # Example +/// +/// ``` +/// # use const_primes::{primes_segment, Error}; +/// const LIMIT: u64 = 5_000_000_031; +/// const PRIMES_GEQ: Result<[u64; 3], Error> = primes_segment!(3; >= LIMIT); +/// const PRIMES_LT: Result<[u64; 3], Error> = primes_segment!(3; < LIMIT); +/// // Can also be used at runtime: +/// let primes_geq = primes_segment!(3; >= LIMIT); +/// +/// assert_eq!(PRIMES_GEQ, primes_geq); +/// assert_eq!(PRIMES_GEQ, Ok([5000000039, 5000000059, 5000000063])); +/// assert_eq!(PRIMES_LT, Ok([4999999903, 4999999937, 5000000029])); +/// ``` +#[macro_export] +macro_rules! primes_segment { + ($n:expr; < $lim:expr) => { + $crate::primes_lt::< + $n, + { + let mem = { $lim }; + $crate::isqrt(mem) as ::core::primitive::usize + 1 + }, + >($lim) + }; + ($n:expr; >= $lim:expr) => { + $crate::primes_geq::< + $n, + { + let mem = { $lim }; + $crate::isqrt(mem) as ::core::primitive::usize + 1 + { $n } + }, + >($lim) + }; +} + +/// Returns the `N` smallest primes greater than or equal to `lower_limit`. +/// Fails to compile if `N` is 0. If `lower_limit` is less than 2 this functions assumes that it is 2, +/// since there are no primes smaller than 2. +/// +/// If you want to compute primes smaller than some limit, take a look at [`primes_lt`]. +/// +/// # Examples +/// +/// Basic usage: +/// ``` +/// use const_primes::{primes_geq, Error}; +/// // Compute 5 primes larger than 40. The largest will be 59, so `MEM` needs to be at least 8. +/// const PRIMES: [u64; 5] = match primes_geq::<5, 8>(40) {Ok(ps) => ps, Err(_) => panic!()}; +/// assert_eq!(PRIMES, [41, 43, 47, 53, 59]); +/// ``` +/// Estimate the size of `MEM` using the square root of the limit (and some extra, proportional to `N`): +/// ``` +/// # use const_primes::{primes_geq, Error}; +/// use const_primes::isqrt; +/// const N: usize = 3; +/// const LIMIT: u64 = 5_000_000_030; +/// # #[allow(long_running_const_eval)] +/// const PRIMES_GEQ: Result<[u64; N], Error> = primes_geq::(LIMIT); +/// assert_eq!(PRIMES_GEQ, Ok([5_000_000_039, 5_000_000_059, 5_000_000_063])); +/// # Ok::<(), Error>(()) +/// ``` +/// # Errors +/// +/// Only primes smaller than `MEM^2` can be generated, so if the sieve +/// encounters a number larger than that it results in an error: +/// ``` +/// # use const_primes::{primes_geq, Error}; +/// const PRIMES: Result<[u64; 3], Error> = primes_geq::<3, 3>(5); +/// assert!(matches!(PRIMES, Err(Error::SieveOverrun(_)))); +/// ``` +/// +/// Also returns an error if `lower_limit` is larger than or equal to `MEM^2`: +/// ``` +/// # use const_primes::{primes_geq, Error}; +/// const PRIMES: Result<[u64; 5], Error> = primes_geq::<5, 5>(26); +/// assert!(matches!(PRIMES, Err(Error::TooLargeLimit(_, _)))); +/// ``` +#[must_use = "the function only returns a new value and does not modify its input"] +pub const fn primes_geq( + lower_limit: u64, +) -> Result<[u64; N], Error> { + const { + assert!(N > 0, "`N` must be at least 1"); + assert!(MEM >= N, "`MEM` must be at least as large as `N`"); + } + let (mem64, mem_sqr) = const { + let mem64 = MEM as u64; + let Some(mem_sqr) = mem64.checked_mul(mem64) else { + panic!("`MEM`^2 must fit in a `u64`") + }; + (mem64, mem_sqr) + }; + + // There are no primes smaller than 2, so we will always start looking at 2. + let new_lower_limit = if lower_limit >= 2 { lower_limit } else { 2 }; + + if new_lower_limit >= mem_sqr { + return Err(Error::TooLargeLimit(lower_limit, mem_sqr)); + } + + let lower_limit = new_lower_limit; + + let mut primes = [0; N]; + let mut total_found_primes = 0; + let mut largest_found_prime = 0; + let base_sieve: [bool; MEM] = sieve(); + let mut sieve_limit = lower_limit; + 'generate: while total_found_primes < N { + let upper_sieve = sieve_segment(&base_sieve, sieve_limit + mem64); + + let mut i = 0; + while i < MEM { + if upper_sieve[i] { + largest_found_prime = sieve_limit + i as u64; + + // We can not know whether this is actually a prime since + // the base sieve contains no information + // about numbers larger than or equal to `MEM`^2. + if largest_found_prime >= mem_sqr { + return Err(Error::SieveOverrun(largest_found_prime)); + } + + if largest_found_prime >= lower_limit { + primes[total_found_primes] = largest_found_prime; + total_found_primes += 1; + if total_found_primes >= N { + // We've found enough primes. + break 'generate; + } + } + } + i += 1; + } + sieve_limit = largest_found_prime + 1; + } + + Ok(primes) +} + +/// The error returned by [`primes_lt`] and [`primes_geq`] if the input +/// is invalid or does not work to produce the requested primes. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub enum Error { + /// The limit was larger than `MEM^2`. + TooLargeLimit(u64, u64), + /// The limit was smaller than or equal to 2. + TooSmallLimit(u64), + /// Encountered a number larger than `MEM`^2. + SieveOverrun(u64), + /// Ran out of primes. + OutOfPrimes, +} + +impl fmt::Display for Error { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + match self { + Self::TooLargeLimit(limit, mem_sqr) => write!( + f, + "the limit ({limit}) was larger than `MEM`^2 ({mem_sqr})" + ), + Self::TooSmallLimit(limit) => write!( + f, + "the limit was {limit}, which is smaller than or equal to 2" + ), + Self::SieveOverrun(number) => write!( + f, + "encountered the number {number} which would have needed `MEM` to be at least {} to sieve", crate::imath::isqrt(*number) + 1 + ), + Self::OutOfPrimes => write!(f, "ran out of primes before the array was filled"), + } + } +} + +#[cfg(feature = "std")] +impl std::error::Error for Error {} + +#[cfg(test)] +mod test { + use crate::is_prime; + + use super::*; + + #[test] + fn sanity_check_primes_geq() { + { + const P: Result<[u64; 5], Error> = primes_geq::<5, 5>(10); + assert_eq!(P, Ok([11, 13, 17, 19, 23])); + } + { + const P: Result<[u64; 5], Error> = primes_geq::<5, 5>(0); + assert_eq!(P, Ok([2, 3, 5, 7, 11])); + } + { + const P: Result<[u64; 1], Error> = primes_geq::<1, 1>(0); + assert_eq!(P, Err(Error::TooLargeLimit(0, 1))); + } + for &prime in primes_geq::<2_000, 2_008>(3_998_000).unwrap().as_slice() { + assert!(is_prime(prime)); + } + } + + #[test] + fn check_primes_segment() { + const P_GEQ: Result<[u64; 10], Error> = primes_segment!(10; >= 1000); + const P_LT: Result<[u64; 10], Error> = primes_segment!(10; < 1000); + + assert_eq!( + P_GEQ, + Ok([1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061]) + ); + assert_eq!(P_LT, Ok([937, 941, 947, 953, 967, 971, 977, 983, 991, 997])); + } +} diff --git a/src/imath.rs b/src/imath.rs index b600d2a..d8c11ec 100644 --- a/src/imath.rs +++ b/src/imath.rs @@ -19,7 +19,7 @@ pub const fn isqrt(n: u64) -> u64 { } } -/// Calculates `(base ^ exp) % modulo` without overflow. +/// Calculates `(base ^ exp) mod modulo` without overflow. #[must_use] pub const fn mod_pow(mut base: u64, mut exp: u64, modulo: u64) -> u64 { let mut res = 1; @@ -37,7 +37,7 @@ pub const fn mod_pow(mut base: u64, mut exp: u64, modulo: u64) -> u64 { res } -/// Calculates `(a * b) % modulo` without overflow. +/// Calculates `(a * b) mod modulo` without overflow. #[must_use] pub const fn mod_mul(a: u64, b: u64, modulo: u64) -> u64 { ((a as u128 * b as u128) % modulo as u128) as u64 diff --git a/src/lib.rs b/src/lib.rs index 9ec4b9c..6c0df6c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -4,21 +4,21 @@ //! //! # Examples //! Generate arrays of prime numbers with the function [`primes`] which uses a -//! [segmented sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve). +//! [segmented sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve): //! ``` //! use const_primes::primes; //! const PRIMES: [u32; 10] = primes(); //! assert_eq!(PRIMES[5], 13); //! assert_eq!(PRIMES, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]); //! ``` -//! or with the type [`Primes`], which ensures that a non-zero number of primes are generated +//! or with the wrapping type [`Primes`]: //! ``` //! use const_primes::Primes; //! const PRIMES: Primes<10> = Primes::new(); //! assert_eq!(PRIMES[5], 13); //! assert_eq!(PRIMES, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]); //! ``` -//! It also lets you reuse it as a cache of primes for related computations +//! which lets you reuse it as a cache of primes for related computations: //! ``` //! # use const_primes::Primes; //! const CACHE: Primes<100> = Primes::new(); @@ -37,442 +37,104 @@ //! assert!(CACHE.is_prime(1000).is_none()); //! assert!(CACHE.count_primes_leq(1000).is_none()); //! ``` -//! Creating a `Primes<0>` is a compile fail in const contexts and a panic otherwise. -//! ```compile_fail -//! # use const_primes::Primes; -//! const PRIMES: Primes<0> = Primes::new(); +//! Sieve a range of numbers for their prime status with [`sieve`]: +//! ``` +//! # use const_primes::sieve; +//! const N: usize = 10; +//! const PRIME_STATUS: [bool; N] = sieve(); +//! // 0 1 2 3 4 5 6 7 8 9 +//! assert_eq!(PRIME_STATUS, [false, false, true, true, false, true, false, true, false, false]); +//! ``` +//! ## Arbitrary ranges +//! +//! The crate provides prime generation and sieving functions with suffixes, e.g. [`primes_lt`] and [`sieve_geq`] +//! that can be used to work with ranges that don't start at zero. They take two generics: the number of primes +//! to store in the binary, and the size of the sieve used during evaluation +//! (which must be at least the ceiling of the square root of the largest encountered number). +//! This means that one can sieve to large numbers, but doesn't need to store the entire sieve in the binary. +//! ``` +//! use const_primes::{primes_lt, Error, isqrt}; +//! const LIMIT: u64 = 5_000_000_031; +//! const N: usize = 3; +//! const PRIMES_LT: Result<[u64; N], Error> = primes_lt::(LIMIT); +//! +//! assert_eq!(PRIMES_LT, Ok([4_999_999_903, 4_999_999_937, 5_000_000_029])); +//! ``` +//! If you do not wish to compute the required sieve size yourself, +//! you can use the provided macro [`primes_segment!`]: +//! ``` +//! # use const_primes::{primes_segment, Error}; +//! const PRIMES_OVER_100: Result<[u64; 3], Error> = primes_segment!(3; >= 100); +//! const PRIMES_UNDER_100: Result<[u64; 3], Error> = primes_segment!(3; < 100); +//! +//! assert_eq!(PRIMES_OVER_100, Ok([101, 103, 107])); +//! assert_eq!(PRIMES_UNDER_100, Ok([83, 89, 97])); //! ``` +//! it may, however, overestimate the required sieve size. +//! +//! ``` +//! # use const_primes::sieve_lt; +//! const N: usize = 70711; +//! const PRIME_STATUS_LT: [bool; N] = sieve_lt(5_000_000_031); +//! // 5_000_000_028 5_000_000_029 5_000_000_030 +//! assert_eq!(PRIME_STATUS_LT[N - 3..], [false, true, false]); +//! ``` +//! Unfortunately the output array must be large enough to contain the prime sieve, which scales with +//! the square root of largest relavant number, which is why the examples use a size of over 70000 even though +//! they're only interested in three numbers. +//! //! ## Other functionality //! -//! Use [`is_prime`] to test whether a given number is prime, +//! Use [`is_prime`] to test whether a given number is prime: //! ``` //! # use const_primes::is_prime; //! const CHECK: bool = is_prime(18_446_744_073_709_551_557); //! assert!(CHECK); //! ``` -//! or [`are_prime`] to compute the prime status of the `N` first integers, +//! Find the next or previous prime numbers with [`next_prime`] and [`previous_prime`] if they exist: //! ``` -//! # use const_primes::are_prime; -//! const N: usize = 10; -//! const PRIME_STATUS: [bool; N] = are_prime(); -//! // 0 1 2 3 4 5 6 7 8 9 -//! assert_eq!(PRIME_STATUS, [false, false, true, true, false, true, false, true, false, false]); -//! ``` -//! or [`are_prime_below`] to compute the prime status of the `N` largest integers below a given value, -//! ``` -//! # use const_primes::are_prime_below; -//! const N: usize = 70711; -//! const BIG_PRIME_STATUS: [bool; N] = are_prime_below(5_000_000_031); -//! // 5_000_000_028 5_000_000_029 5_000_000_030 -//! assert_eq!(BIG_PRIME_STATUS[N - 3..], [false, true, false]); +//! # use const_primes::{next_prime, previous_prime}; +//! const NEXT: Option = next_prime(25); +//! const PREV: Option = previous_prime(25); +//! const NOSUCH: Option = previous_prime(2); +//! +//! assert_eq!(NEXT, Some(29)); +//! assert_eq!(PREV, Some(23)); +//! assert_eq!(NOSUCH, None); //! ``` //! and more! +//! +//! # Features +//! +//! `std`: derives the [`Error`](std::error::Error) trait for the error types. +//! `alloc`: enables conversion of the type returned by [`primes_geq`] and [`primes_lt`] into [`Vec`]s and [`Box`]ed slices. #![forbid(unsafe_code)] -#![cfg_attr(not(test), no_std)] +#![cfg_attr(all(not(test), not(feature = "std")), no_std)] /// The type that `Primes` stores, and `primes::()`` returns. Currently `u32`. // Just change this to whatever unsigned primitive integer type you want and it should work as long as it has enough bits for your purposes. // This is used since there is currently no way to be generic over types that can do arithmetic at compile time. type Underlying = u32; +mod generation; mod imath; mod miller_rabin; +mod other_prime; +mod sieving; mod wrapper; -use imath::isqrt; + +pub use generation::{primes, primes_geq, primes_lt, Error}; +pub use imath::isqrt; pub use miller_rabin::is_prime; +pub use other_prime::{next_prime, previous_prime}; +pub use sieving::{sieve, sieve_geq, sieve_lt}; pub use wrapper::Primes; -/// Returns the `N` first prime numbers. -/// -/// [`Primes`] might be relevant for you if you intend to later use these prime numbers for related computations. -/// -/// Uses a [segmented sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve). -/// -/// # Example -/// ``` -/// # use const_primes::primes; -/// const PRIMES: [u32; 10] = primes(); -/// assert_eq!(PRIMES, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]); -/// ``` -/// # Panics -/// Panics if a computed prime overflows a `u32`. This will result in a compile error in a const context. -#[must_use = "the function only returns a new value"] -pub const fn primes() -> [Underlying; N] { - if N == 0 { - return [0; N]; - } else if N == 1 { - return [2; N]; - } else if N == 2 { - let mut primes = [0; N]; - primes[0] = 2; - primes[1] = 3; - return primes; - } - - // This is a segmented sieve that runs until it has found enough primes. - - // This array is the output in the end - let mut primes = [0; N]; - // This keeps track of how many primes we've found so far. - let mut prime_count = 0; - - // Sieve the first primes below N - let mut sieve: [bool; N] = are_prime(); - - // Count how many primes we found - // and store them in the final array - let mut number = 0; - while number < N { - if sieve[number] { - primes[prime_count] = number as Underlying; - prime_count += 1; - } - - number += 1; - } - - // For every segment of N numbers - let mut low = N - 1; - let mut high = 2 * N - 1; - 'generate: while prime_count < N { - // reset the sieve for the segment - sieve = [true; N]; - let mut i = 0; - - // and repeat for each prime found so far: - while i < prime_count { - let prime = primes[i] as usize; - - // Find the smallest composite in the current segment, - let mut composite = (low / prime) * prime; - if composite < low { - composite += prime; - } - - // and sieve all numbers in the segment that are multiples of the prime. - while composite < high { - sieve[composite - low] = false; - composite += prime; - } - - i += 1; - } - - // Move the found primes into the final array - i = low; - while i < high { - if sieve[i - low] { - primes[prime_count] = i as Underlying; - prime_count += 1; - // and stop the generation of primes if we're done. - if prime_count == N { - break 'generate; - } - } - i += 1; - } - - // Update low and high for the next segment - low += N; - high += N; - } - - primes -} - -/// Returns an array of size `N` that indicates which of the integers in `[upper_limit - N, upper_limit)` are prime, -/// or in other words: the value at a given index represents whether `index + upper_limit - N` is prime. -/// -/// If you just want the prime status of the first `N` integers, see [`are_prime`]. -/// -/// Uses a sieve of Eratosthenes to sieve the first `N` integers -/// and then uses the result to sieve the output range if needed. -/// -/// # Examples -/// Basic usage -/// ``` -/// # use const_primes::are_prime_below; -/// const PRIME_STATUSES: [bool; 10] = are_prime_below(30); -/// -/// assert_eq!( -/// PRIME_STATUSES, -/// // 20 21 22 23 24 25 26 27 28 29 -/// [false, false, false, true, false, false, false, false, false, true], -/// ); -/// ``` -/// Sieve limited ranges of very large values -/// ``` -/// # use const_primes::are_prime_below; -/// const BIG_NUMBER: u64 = 5_000_000_031; -/// const CEIL_SQRT_BIG_NUMBER: usize = 70711; -/// const BIG_PRIME_STATUSES: [bool; CEIL_SQRT_BIG_NUMBER] = are_prime_below(BIG_NUMBER); -/// assert_eq!( -/// BIG_PRIME_STATUSES[CEIL_SQRT_BIG_NUMBER - 3..], -/// // 5_000_000_028 5_000_000_029 5_000_000_030 -/// [false, true, false], -/// ); -/// ``` -/// -/// # Panics -/// Panics if `upper_limit` is not in the range `[N, N^2]`, or if `N^2` overflows a `u64`. -/// These are compile errors in const contexts. -/// ```compile_fail -/// # use const_primes::are_prime_below; -/// const PRIME_STATUSES: [bool; 5] = are_prime_below(26); -/// ``` -/// ```compile_fail -/// # use const_primes::are_prime_below; -/// const PRIME_STATUSES: [bool; 5] = are_prime_below(4); -/// ``` -#[must_use = "the function returns a new value and does not modify its input"] -pub const fn are_prime_below(upper_limit: u64) -> [bool; N] { - let n64 = N as u64; - - // Since panics are compile time errors in const contexts - // we check all the preconditions here and panic early. - match n64.checked_mul(n64) { - Some(prod) => assert!( - upper_limit <= prod, - "`upper_limit` must be smaller than or equal to `N`^2" - ), - None => panic!("`N`^2 must fit in a `u64`"), - } - assert!(upper_limit >= n64, "`upper_limit` must be at least `N`"); - - // Use a normal sieve of Eratosthenes for the first N numbers. - let base_sieve: [bool; N] = are_prime(); - - if upper_limit == n64 { - // If we are not interested in sieving a larger range we can just return early. - return base_sieve; - } - - // Otherwise we use the base sieve to sieve the upper range - let mut segment_sieve = [true; N]; - - let lower_limit = upper_limit - n64; - - // In case 0 and/or 1 are included in the upper sieve we need to treat them as a special case - // since they are not multiples of any prime in `base_sieve` even though they are not primes. - if lower_limit == 0 && N > 1 { - segment_sieve[0] = false; - segment_sieve[1] = false; - } else if lower_limit == 1 && N > 0 { - segment_sieve[0] = false; - } - - // For all the found primes - let mut i = 0; - while i < N { - if base_sieve[i] { - let prime = i as u64; - - // Find the smallest multiple of the prime larger than or equal to `lower_limit`. - let mut composite = (lower_limit / prime) * prime; - if composite < lower_limit { - composite += prime; - } - if composite == prime { - composite += prime; - } - - // Sieve all numbers in the segment that are multiples of the prime. - while composite < upper_limit { - segment_sieve[(composite - lower_limit) as usize] = false; - composite += prime; - } - } - i += 1; - } - - segment_sieve -} - -/// Returns an array of size `N` where the value at a given index indicates whether the index is prime. -/// -/// Uses a sieve of Eratosthenes. -/// -/// # Example -/// ``` -/// # use const_primes::are_prime; -/// const PRIMALITY: [bool; 10] = are_prime(); -/// // 0 1 2 3 4 5 6 7 8 9 -/// assert_eq!(PRIMALITY, [false, false, true, true, false, true, false, true, false, false]); -/// ``` -#[must_use = "the function only returns a new value"] -pub const fn are_prime() -> [bool; N] { - let mut sieve = [true; N]; - if N > 0 { - sieve[0] = false; - } - if N > 1 { - sieve[1] = false; - } - - let mut number: usize = 2; - let bound = isqrt(N as u64); - // For all numbers up to and including sqrt(n): - while (number as u64) <= bound { - if sieve[number] { - // If a number is prime we enumerate all multiples of it - // starting from its square, - let Some(mut composite) = number.checked_mul(number) else { - break; - }; - - // and mark them as not prime. - while composite < N { - sieve[composite] = false; - composite = match composite.checked_add(number) { - Some(sum) => sum, - None => break, - }; - } - } - number += 1; - } - - sieve -} - -/// Returns the value of the [Möbius function](https://en.wikipedia.org/wiki/M%C3%B6bius_function). -/// -/// This function is -/// - 1 if `x` is a square-free integer with an even number of prime factors, -/// - -1 if `x` is a square-free integer with an odd number of prime factors, -/// - 0 if `x` has a squared prime factor. -/// -/// Uses a small wheel to check prime factors up to `√x` and exits early if -/// there is a squared factor. -/// -/// # Example -/// ``` -/// # use const_primes::moebius; -/// const N: u64 = 1001; -/// const MÖBIUS1001: i8 = moebius(N); -/// assert_eq!(MÖBIUS1001, -1); -/// ``` -// I have added this code to Rosettacode: https://www.rosettacode.org/wiki/M%C3%B6bius_function#Rust -// as of the writing of this comment. -pub const fn moebius(mut x: u64) -> i8 { - let mut prime_count: u64 = 0; - - /// If x is divisible by the any of the given factors this macro counts the factor and divides it out. - /// It then returns zero if x is still divisible by the factor. - macro_rules! handle_factors { - ($($factor:expr),+) => { - $( - if x % $factor == 0 { - x /= $factor; - prime_count += 1; - // If x is still divisible by the factor that means x has a - // square or more prime factor, and we return 0. - if x % $factor == 0 { return 0; } - } - )+ - }; - } - - // Handle 2 and 3 separately, since the wheel will not find these factors. - handle_factors!(2, 3); - - // Handle the remaining factors <= √x with the wheel. - let mut i = 5; - let bound = isqrt(x); - while i <= bound { - handle_factors!(i, i + 2); - i += 6; - } - - // There can exist one prime factor larger than √x, - // in that case we can check if x is still larger than one, and then count it. - if x > 1 { - prime_count += 1; - } - - if prime_count % 2 == 0 { - 1 - } else { - -1 - } -} - -/// Returns the largest prime smaller than or equal to `n` if there is one. -/// -/// Scans for primes downwards from the input with [`is_prime`]. -/// -/// # Examples -/// ``` -/// # use const_primes::largest_prime_leq; -/// const LPLEQ: Option = largest_prime_leq(400); -/// assert_eq!(LPLEQ, Some(397)); -/// ``` -/// There's no prime smaller than or equal to one -/// ``` -/// # use const_primes::largest_prime_leq; -/// const NOSUCH: Option = largest_prime_leq(1); -/// assert!(NOSUCH.is_none()); -/// ``` -#[must_use = "the function only returns a new value and does not modify its input"] -pub const fn largest_prime_leq(mut n: u64) -> Option { - if n == 0 || n == 1 { - None - } else if n == 2 { - Some(2) - } else { - if n % 2 == 0 { - n -= 1; - } - - while !is_prime(n) { - n -= 2; - } - - Some(n) - } -} - -/// Returns the smallest prime greater than or equal to `n` if there is one that -/// can be represented by a `u64`. -/// -/// Scans for primes upwards from the input with [`is_prime`]. -/// -/// # Example -/// ``` -/// # use const_primes::smallest_prime_geq; -/// const SPGEQ: Option = smallest_prime_geq(400); -/// assert_eq!(SPGEQ, Some(401)); -/// ``` -/// Primes larger than 18446744073709551557 can not be represented by a `u64` -/// ``` -/// # use const_primes::smallest_prime_geq; -/// const NOSUCH: Option = smallest_prime_geq(18_446_744_073_709_551_558); -/// assert!(NOSUCH.is_none()); -/// ``` -#[must_use = "the function only returns a new value and does not modify its input"] -pub const fn smallest_prime_geq(mut n: u64) -> Option { - // The largest prime smaller than 2^64 - if n > 18_446_744_073_709_551_557 { - None - } else if n <= 2 { - Some(2) - } else { - if n % 2 == 0 { - n += 1; - } - - while !is_prime(n) { - n += 2; - } - - Some(n) - } -} - /// Returns an array of size `N` where the value at a given index is how many primes are less than or equal to the index. +/// Fails to compile if `N` is 0. /// -/// Sieves primes with [`are_prime`] and then counts them. +/// Sieves primes with [`sieve`] and then counts them. /// /// # Example /// Basic usage @@ -483,12 +145,13 @@ pub const fn smallest_prime_geq(mut n: u64) -> Option { /// ``` #[must_use = "the function only returns a new value"] pub const fn prime_counts() -> [usize; N] { + const { assert!(N > 0, "`N` must be at least 1") } let mut counts = [0; N]; if N <= 2 { return counts; } counts[2] = 1; - let prime_status: [bool; N] = are_prime(); + let prime_status: [bool; N] = sieve(); let mut count = 1; let mut i = 3; while i < N { @@ -521,21 +184,24 @@ mod test { } #[test] - fn check_are_prime() { + fn check_sieve() { macro_rules! test_to { ($($n:expr),+) => { $( - println!("{}", $n); - assert_eq!(&PRIMALITIES[..$n], are_prime::<$n>()); + { + const P: [bool; $n] = sieve(); + assert_eq!(&PRIMALITIES[..$n], P); + assert_eq!(&PRIMALITIES[..$n], sieve::<$n>()); + } )+ }; } test_to!( - 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, - 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, - 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, - 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, - 90, 91, 92, 93, 94, 95, 96, 97, 98, 99 + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, + 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, + 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, + 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, + 91, 92, 93, 94, 95, 96, 97, 98, 99 ); } @@ -555,38 +221,69 @@ mod test { i += 1; } } - assert_eq!(PRIME_COUNTS, counts); + assert_eq!(counts, prime_counts()); } )+ }; } - test_prime_counts_up_to!(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100); + test_prime_counts_up_to!(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100); } #[test] - fn check_möbius() { - #[rustfmt::skip] - const TEST_CASES: [i8; 51] = [0, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, 0, 1, 1, -1, 0, 0, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1, 0, -1, -1, -1, 0, 0, 1, -1, 0, 0, 0]; - for (n, ans) in TEST_CASES.into_iter().enumerate() { - assert_eq!(moebius(n as u64), ans); + fn verify_primes() { + macro_rules! test_to_n { + ($($n:expr),+) => { + $( + { + const PRIMES: [Underlying; $n] = primes(); + assert_eq!(PRIMES, PRECOMPUTED_PRIMES[..$n]); + assert_eq!(PRECOMPUTED_PRIMES[..$n], primes::<$n>()); + } + )+ + }; } - #[rustfmt::skip] - const BIG_TEST_CASES: [i8; 51] = [0, -1, -1, 1, 0, -1, 1, 1, 0, -1, -1, 1, 0, -1, 0, -1, 0, 0, 1, -1, 0, -1, -1, -1, 0, 0, 0, 1, 0, 0, -1, -1, 0, -1, -1, 0, 0, 1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1, 0, -1, 0]; - for (n, ans) in BIG_TEST_CASES.into_iter().enumerate() { - assert_eq!(moebius(n as u64 + 1000), ans); + + test_to_n!(1, 2, 3, 4, 5, 10, 100, 1000, 10000); + } + + #[test] + fn check_primes_lt() { + macro_rules! test_n_below_100 { + ($($n:expr),+) => { + $( + { + const P: Result<[u64; $n], Error> = primes_lt::<$n, $n>(100); + for (i, prime) in P.unwrap().as_slice().into_iter().enumerate() { + assert_eq!(PRECOMPUTED_PRIMES[25-$n..25][i], *prime as u32); + } + assert_eq!( + PRECOMPUTED_PRIMES[25-$n..25], + primes_lt::<$n, $n>(100).unwrap().as_slice().into_iter().map(|i| *i as u32).collect::>() + ); + } + )+ + }; } - assert_eq!(moebius(u32::MAX.into()), -1); + + test_n_below_100!(10, 15, 20); + + assert_eq!([2, 3, 5, 7], primes_lt::<4, 9>(10).unwrap().as_slice()); + + assert_eq!(primes_lt::<1, 2>(3).unwrap().as_slice(), [2]); } #[test] - fn check_are_prime_below() { + fn check_sieve_lt() { macro_rules! test_n_below_100 { ($($n:expr),+) => { $( - println!("{}", $n); - assert_eq!(&PRIMALITIES[100-$n..], are_prime_below::<$n>(100)); + { + const P: [bool; $n] = sieve_lt(100); + assert_eq!(&PRIMALITIES[100-$n..], P); + assert_eq!(&PRIMALITIES[100-$n..], sieve_lt::<$n>(100)); + } )+ }; } @@ -600,40 +297,42 @@ mod test { } #[test] - fn verify_primes() { - macro_rules! test_to_n { + fn check_sieve_geq() { + macro_rules! test_n_geq_10 { ($($n:expr),+) => { $( { - const PRIMES: [Underlying; $n] = $crate::primes(); - assert_eq!(PRIMES, PRECOMPUTED_PRIMES[..$n]); + const P: [bool; $n] = sieve_geq(10); + assert_eq!(&PRIMALITIES[10..10+$n], P); + assert_eq!(&PRIMALITIES[10..10+$n], sieve_geq::<$n>(10)); } )+ }; } - - test_to_n!(0, 1, 2, 3, 4, 5, 10, 100, 1000, 10000); + test_n_geq_10!(4, 5, 10, 20, 30, 40, 50); } #[test] fn check_next_prime() { for i in 1..PRECOMPUTED_PRIMES.len() - 1 { assert_eq!( - smallest_prime_geq(PRECOMPUTED_PRIMES[i] as u64 + 1), + next_prime(PRECOMPUTED_PRIMES[i] as u64), Some(PRECOMPUTED_PRIMES[i + 1] as u64) ); assert_eq!( - largest_prime_leq(PRECOMPUTED_PRIMES[i] as u64 - 1), + previous_prime(PRECOMPUTED_PRIMES[i] as u64), Some(PRECOMPUTED_PRIMES[i - 1] as u64) ); } - assert_eq!(smallest_prime_geq(0), Some(2)); - assert_eq!(smallest_prime_geq(1), Some(2)); - assert_eq!(smallest_prime_geq(2), Some(2)); - assert_eq!(largest_prime_leq(0), None); - assert_eq!(largest_prime_leq(1), None); - assert_eq!(largest_prime_leq(2), Some(2)); + assert_eq!(next_prime(18_446_744_073_709_551_558), None); + assert_eq!(next_prime(0), Some(2)); + assert_eq!(next_prime(1), Some(2)); + assert_eq!(next_prime(2), Some(3)); + assert_eq!(previous_prime(0), None); + assert_eq!(previous_prime(1), None); + assert_eq!(previous_prime(2), None); + assert_eq!(previous_prime(3), Some(2)); } #[rustfmt::skip] diff --git a/src/miller_rabin.rs b/src/miller_rabin.rs index c916fce..7f28ada 100644 --- a/src/miller_rabin.rs +++ b/src/miller_rabin.rs @@ -23,7 +23,8 @@ pub const fn is_prime(n: u64) -> bool { false } else { let mut candidate_factor = 5; - while candidate_factor <= n.ilog2() as u64 { + let trial_limit = n.ilog2() as u64; + while candidate_factor <= trial_limit { if n % candidate_factor == 0 || n % (candidate_factor + 2) == 0 { return false; } diff --git a/src/other_prime.rs b/src/other_prime.rs new file mode 100644 index 0000000..b37690a --- /dev/null +++ b/src/other_prime.rs @@ -0,0 +1,79 @@ +use crate::is_prime; + +/// Returns the largest prime smaller than `n` if there is one. +/// +/// Scans for primes downwards from the input with [`is_prime`]. +/// +/// # Examples +/// Basic usage: +/// ``` +/// # use const_primes::previous_prime; +/// const LPLEQ: Option = previous_prime(400); +/// assert_eq!(LPLEQ, Some(397)); +/// ``` +/// There's no prime smaller than two: +/// ``` +/// # use const_primes::previous_prime; +/// const NOSUCH: Option = previous_prime(2); +/// assert!(NOSUCH.is_none()); +/// ``` +#[must_use = "the function only returns a new value and does not modify its input"] +pub const fn previous_prime(mut n: u64) -> Option { + if n <= 2 { + None + } else if n == 3 { + Some(2) + } else { + n -= 1; + + if n % 2 == 0 { + n -= 1; + } + + while !is_prime(n) { + n -= 2; + } + + Some(n) + } +} + +/// Returns the smallest prime greater than `n` if there is one that +/// can be represented by a `u64`. +/// +/// Scans for primes upwards from the input with [`is_prime`]. +/// +/// # Example +/// Basic usage: +/// ``` +/// # use const_primes::next_prime; +/// const SPGEQ: Option = next_prime(400); +/// assert_eq!(SPGEQ, Some(401)); +/// ``` +/// Primes larger than 18446744073709551557 can not be represented by a `u64`: +/// ``` +/// # use const_primes::next_prime; +/// const NOSUCH: Option = next_prime(18_446_744_073_709_551_557); +/// assert!(NOSUCH.is_none()); +/// ``` +#[must_use = "the function only returns a new value and does not modify its input"] +pub const fn next_prime(mut n: u64) -> Option { + // The largest prime smaller than u64::MAX + if n >= 18_446_744_073_709_551_557 { + None + } else if n <= 1 { + Some(2) + } else { + n += 1; + + if n % 2 == 0 { + n += 1; + } + + while !is_prime(n) { + n += 2; + } + + Some(n) + } +} diff --git a/src/sieving.rs b/src/sieving.rs new file mode 100644 index 0000000..aa0e861 --- /dev/null +++ b/src/sieving.rs @@ -0,0 +1,241 @@ +use crate::isqrt; + +/// Uses the primalities of the first `N` integers in `base_sieve` to sieve the numbers in the range `[upper_limit - N, upper_limit)`. +/// Assumes that the base sieve contains the prime status of the `N` fist integers. The output is only meaningful +/// for the numbers below `N^2`. Fails to compile if `N` is 0. +#[must_use = "the function only returns a new value and does not modify its inputs"] +pub(crate) const fn sieve_segment( + base_sieve: &[bool; N], + upper_limit: u64, +) -> [bool; N] { + const { assert!(N > 0, "`N` must be at least 1") } + + let mut segment_sieve = [true; N]; + + let lower_limit = upper_limit - N as u64; + + if lower_limit == 0 && N > 1 { + // If the lower limit is 0 we can just return the base sieve. + return *base_sieve; + } else if lower_limit == 1 && N > 0 { + // In case 1 is included in the upper sieve we need to treat it as a special case + // since it's not a multiple of any prime in `base_sieve` even though it's not prime. + segment_sieve[0] = false; + } + + let mut i = 0; + while i < N { + if base_sieve[i] { + let prime = i as u64; + + // Find the smallest multiple of the prime larger than or equal to `lower_limit`. + let mut composite = (lower_limit / prime) * prime; + if composite < lower_limit { + composite += prime; + } + if composite == prime { + composite += prime; + } + + // Sieve all numbers in the segment that are multiples of the prime. + while composite < upper_limit { + segment_sieve[(composite - lower_limit) as usize] = false; + composite += prime; + } + } + i += 1; + } + + segment_sieve +} + +/// Returns an array of size `N` that indicates which of the integers in `[upper_limit - N, upper_limit)` are prime, +/// or in other words: the value at a given index represents whether `index + upper_limit - N` is prime. +/// +/// If you just want the prime status of the first `N` integers, see [`sieve`]. +/// +/// Uses a sieve of Eratosthenes to sieve the first `N` integers +/// and then uses the result to sieve the output range if needed. +/// +/// Fails to compile if `N` is 0. +/// +/// # Examples +/// +/// Basic usage +/// ``` +/// # use const_primes::sieve_lt; +/// const PRIME_STATUSES: [bool; 10] = sieve_lt(30); +/// +/// assert_eq!( +/// PRIME_STATUSES, +/// // 20 21 22 23 24 25 26 27 28 29 +/// [false, false, false, true, false, false, false, false, false, true], +/// ); +/// ``` +/// Sieve limited ranges of very large values +/// ``` +/// # use const_primes::sieve_lt; +/// const BIG_NUMBER: u64 = 5_000_000_031; +/// const CEIL_SQRT_BIG_NUMBER: usize = 70711; +/// const BIG_PRIME_STATUSES: [bool; CEIL_SQRT_BIG_NUMBER] = sieve_lt(BIG_NUMBER); +/// assert_eq!( +/// BIG_PRIME_STATUSES[CEIL_SQRT_BIG_NUMBER - 3..], +/// // 5_000_000_028 5_000_000_029 5_000_000_030 +/// [false, true, false], +/// ); +/// ``` +/// +/// # Panics +/// +/// Panics if `upper_limit` is not in the range `[N, N^2]`. In const contexts these are compile errors: +/// ```compile_fail +/// # use const_primes::sieve_lt; +/// const PRIME_STATUSES: [bool; 5] = sieve_lt(26); +/// ``` +/// ```compile_fail +/// # use const_primes::sieve_lt; +/// const PRIME_STATUSES: [bool; 5] = sieve_lt(4); +/// ``` +#[must_use = "the function only returns a new value and does not modify its input"] +pub const fn sieve_lt(upper_limit: u64) -> [bool; N] { + const { assert!(N > 0, "`N` must be at least 1") } + + let n64 = N as u64; + + // Since panics are compile time errors in const contexts + // we check all the preconditions here and panic early. + match n64.checked_mul(n64) { + Some(prod) => assert!( + upper_limit <= prod, + "`upper_limit` must be smaller than or equal to `N^2`" + ), + None => panic!("`N^2` must fit in a `u64`"), + } + assert!(upper_limit >= n64, "`upper_limit` must be at least `N`"); + + // Use a normal sieve of Eratosthenes for the first N numbers. + let base_sieve: [bool; N] = sieve(); + + if upper_limit == n64 { + // If we are not interested in sieving a larger range we can just return early. + return base_sieve; + } + + sieve_segment(&base_sieve, upper_limit) +} + +/// Returns an array of size `N` where the value at a given index indicates whether the index is prime. +/// Fails to compile if `N` is 0. +/// +/// Uses a sieve of Eratosthenes. +/// +/// # Example +/// +/// ``` +/// # use const_primes::sieve; +/// const PRIMALITY: [bool; 10] = sieve(); +/// // 0 1 2 3 4 5 6 7 8 9 +/// assert_eq!(PRIMALITY, [false, false, true, true, false, true, false, true, false, false]); +/// ``` +#[must_use = "the function only returns a new value"] +pub const fn sieve() -> [bool; N] { + const { assert!(N > 0, "`N` must be at least 1") } + + let mut sieve = [true; N]; + if N > 0 { + sieve[0] = false; + } + if N > 1 { + sieve[1] = false; + } + + let mut number: usize = 2; + let bound = isqrt(N as u64); + // For all numbers up to and including sqrt(n): + while (number as u64) <= bound { + if sieve[number] { + // If a number is prime we enumerate all multiples of it + // starting from its square, + let Some(mut composite) = number.checked_mul(number) else { + break; + }; + + // and mark them as not prime. + while composite < N { + sieve[composite] = false; + composite = match composite.checked_add(number) { + Some(sum) => sum, + None => break, + }; + } + } + number += 1; + } + + sieve +} + +/// Returns the prime status of the `N` smallest integers greater than or equal to `lower_limit`. +/// Fails to compile if `N` is 0. +/// +/// # Example +/// +/// Basic usage: +/// ``` +/// # use const_primes::sieve_geq; +/// const PRIME_STATUS: [bool; 5] = sieve_geq(10); +/// // 10 11 12 13 14 +/// assert_eq!(PRIME_STATUS, [false, true, false, true, false]); +/// ``` +/// # Panics +/// +/// Panics if `N + lower_limit` is larger than to `N^2`. In const contexts this is a compile error: +/// ```compile_fail +/// # use const_primes::sieve_geq; +/// const P: [bool; 5] = sieve_geq(21); +/// ``` +#[must_use = "the function only returns a new value and does not modify its input"] +pub const fn sieve_geq(lower_limit: u64) -> [bool; N] { + const { assert!(N > 0, "`N` must be at least 1") } + + let n64 = N as u64; + + // Since panics are compile time errors in const contexts + // we check all the preconditions here and panic early. + let upper_limit = if let Some(sum) = n64.checked_add(lower_limit) { + sum + } else { + panic!("`N + lower_limit` must fit in a `u64`") + }; + if let Some(n_sqr) = n64.checked_mul(n64) { + assert!( + upper_limit <= n_sqr, + "`lower_limit + N` must be less than or equal to `N^2`" + ); + } else { + panic!("`N^2` must fit in a `u64`") + } + + let base_sieve: [bool; N] = sieve(); + + // If `lower_limit` is zero the upper range is the same as what we already sieved, + // so we return early. + if lower_limit == 0 { + return base_sieve; + } + + sieve_segment(&base_sieve, upper_limit) +} + +#[cfg(test)] +mod test { + use super::{sieve, sieve_segment}; + + #[test] + fn test_consistency_of_sieve_segment() { + const P: [bool; 10] = sieve_segment(&sieve(), 10); + const PP: [bool; 10] = sieve_segment(&sieve(), 11); + assert_eq!(P, sieve()); + assert_eq!(PP, sieve::<11>()[1..]); + } +} diff --git a/src/wrapper.rs b/src/wrapper.rs index f588261..8781115 100644 --- a/src/wrapper.rs +++ b/src/wrapper.rs @@ -3,9 +3,10 @@ use crate::{primes, Underlying}; // region: Primes /// A wrapper around an array that consists of the first `N` primes. -/// Can be created in const contexts, and if so it ensures that `N` is non-zero at compile time. +/// Can be created in const contexts, and ensures that `N` is non-zero at compile time. /// /// # Examples +/// /// Basic usage /// ``` /// # use const_primes::Primes; @@ -36,6 +37,7 @@ impl Primes { /// Uses a [segmented sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve). /// /// # Examples + /// /// Basic usage /// ``` /// # use const_primes::Primes; @@ -54,19 +56,18 @@ impl Primes { /// assert_eq!(primes, [2, 3, 5, 7, 11]); /// ``` /// - /// # Panics - /// - /// Panics if `N` is zero. In const contexts this will fail to compile + /// Fails to compile if `N` is zero. /// ```compile_fail /// # use const_primes::Primes; /// const NO_PRIMES: Primes<0> = Primes::new(); /// ``` - /// In other contexts it may panic at runtime instead. + /// + /// # Panics + /// /// If any of the primes overflow a `u32` it will panic in const contexts or debug mode. #[must_use = "the associated method only returns a new value"] pub const fn new() -> Self { - assert!(N >= 1, "`N` must be at least 1"); - + const { assert!(N > 0, "`N` must be at least 1") } Self { primes: primes() } } @@ -127,45 +128,57 @@ impl Primes { // region: Next prime - /// Returns the largest prime less than or equal to `n`. - /// If `n` is 0, 1, or larger than the largest prime in `self` this returns `None`. + /// Returns the largest prime less than `n`. + /// If `n` is 0, 1, 2, or larger than the largest prime in `self` this returns `None`. /// /// Uses a binary search. /// # Example /// ``` /// # use const_primes::Primes; /// const CACHE: Primes<100> = Primes::new(); - /// const LPLEQ400: Option = CACHE.largest_prime_leq(400); - /// assert_eq!(LPLEQ400, Some(397)); + /// const PREV400: Option = CACHE.previous_prime(400); + /// assert_eq!(PREV400, Some(397)); /// ``` #[must_use = "the method only returns a new value and does not modify `self`"] - pub const fn largest_prime_leq(&self, n: Underlying) -> Option { - if n <= 1 { + pub const fn previous_prime(&self, n: Underlying) -> Option { + if n <= 2 { None } else { match self.binary_search(n) { - Ok(i) => Some(self.primes[i]), - Err(Some(i)) => Some(self.primes[i - 1]), + Ok(i) | Err(Some(i)) => { + if i > 0 { + Some(self.primes[i - 1]) + } else { + None + } + } Err(None) => None, } } } - /// Returns the smallest prime greater than or equal to `n`. - /// If `n` is larger than the largest prime in `self` this returns `None`. + /// Returns the smallest prime greater than `n`. + /// If `n` is larger than or equal to the largest prime in `self` this returns `None`. /// /// Uses a binary search. /// # Example /// ``` /// # use const_primes::Primes; /// const CACHE: Primes<100> = Primes::new(); - /// const SPGEQ: Option = CACHE.smallest_prime_geq(400); - /// assert_eq!(SPGEQ, Some(401)); + /// const NEXT: Option = CACHE.next_prime(400); + /// assert_eq!(NEXT, Some(401)); /// ``` #[must_use = "the method only returns a new value and does not modify `self`"] - pub const fn smallest_prime_geq(&self, n: Underlying) -> Option { + pub const fn next_prime(&self, n: Underlying) -> Option { match self.binary_search(n) { - Ok(i) | Err(Some(i)) => Some(self.primes[i]), + Ok(i) => { + if i + 1 < self.len() { + Some(self.primes[i + 1]) + } else { + None + } + } + Err(Some(i)) => Some(self.primes[i]), Err(None) => None, } } @@ -175,10 +188,10 @@ impl Primes { /// Searches the underlying array of primes for the target integer. /// If the target is found it returns a [`Result::Ok`] that contains the index of the matching element. /// If the target is not found in the array a [`Result::Err`] is returned that contains an [`Option`]. - /// If the target could be inserted into the array while maintaining the sorted order, the [`Some`](Option::Some) - /// variant contains the index of that location. + /// If the target could be inserted into the array while maintaining the sorted order, the [`Option::Some`] + /// variant is returned and contains the index of that location. /// If the target is larger than the largest prime in the array no information about where it might fit is available, - /// and a [`None`](Option::None) is returned. + /// and an [`Option::None`] is returned. #[must_use = "the method only returns a new value and does not modify `self`"] pub const fn binary_search(&self, target: Underlying) -> Result> { if target > *self.last() { @@ -208,6 +221,7 @@ impl Primes { /// Converts `self` into an array of size `N`. /// /// # Example + /// /// Basic usage /// ``` /// # use const_primes::Primes; @@ -296,6 +310,14 @@ impl Primes { } } +/// This statically asserts that N > 0. +impl Default for Primes { + fn default() -> Self { + const { assert!(N > 0, "`N` must be at least 1") } + Self::new() + } +} + impl core::ops::Index for Primes where I: core::slice::SliceIndex<[Underlying]>, @@ -303,7 +325,7 @@ where type Output = I::Output; #[inline] fn index(&self, index: I) -> &Self::Output { - &self.primes[index] + self.primes.index(index) } } @@ -487,40 +509,47 @@ mod test { } #[test] - fn check_largest_prime_leq() { + fn check_previous_prime() { const CACHE: Primes<100> = Primes::new(); - const LPLEQ0: Option = CACHE.largest_prime_leq(0); - const LPLEQ400: Option = CACHE.largest_prime_leq(400); - const LPLEQ541: Option = CACHE.largest_prime_leq(541); - const LPLEQ542: Option = CACHE.largest_prime_leq(542); - assert_eq!(LPLEQ0, None); - assert_eq!(LPLEQ400, Some(397)); - assert_eq!(LPLEQ541, Some(541)); - assert_eq!(LPLEQ542, None); + const PREV0: Option = CACHE.previous_prime(0); + const PREV400: Option = CACHE.previous_prime(400); + const PREV541: Option = CACHE.previous_prime(541); + const PREV542: Option = CACHE.previous_prime(542); + const PREVS: [Underlying; 18] = [ + 2, 3, 3, 5, 5, 7, 7, 7, 7, 11, 11, 13, 13, 13, 13, 17, 17, 19, + ]; + for (i, prev) in PREVS.into_iter().enumerate() { + println!("n = {}", i + 3); + assert_eq!(Some(prev), CACHE.previous_prime(i as u32 + 3)); + } + assert_eq!(PREV0, None); + assert_eq!(PREV400, Some(397)); + assert_eq!(PREV541, Some(523)); + assert_eq!(PREV542, None); } #[test] - fn check_smallest_prime_geq() { + fn check_next_prime() { const CACHE: Primes<100> = Primes::new(); - const SPGEQ0: Option = CACHE.smallest_prime_geq(0); - const SPGEQ400: Option = CACHE.smallest_prime_geq(400); - const SPGEQ541: Option = CACHE.smallest_prime_geq(541); - const SPGEQ542: Option = CACHE.smallest_prime_geq(542); + const SPGEQ0: Option = CACHE.next_prime(0); + const SPGEQ400: Option = CACHE.next_prime(400); + const SPGEQ541: Option = CACHE.next_prime(540); + const SPGEQ542: Option = CACHE.next_prime(541); assert_eq!(SPGEQ0, Some(2)); assert_eq!(SPGEQ400, Some(401)); assert_eq!(SPGEQ541, Some(541)); assert_eq!(SPGEQ542, None); - const N: usize = 32; + const N: usize = 31; const NEXT_PRIME: [u32; N] = [ - 2, 2, 2, 3, 5, 5, 7, 7, 11, 11, 11, 11, 13, 13, 17, 17, 17, 17, 19, 19, 23, 23, 23, 23, + 2, 2, 3, 5, 5, 7, 7, 11, 11, 11, 11, 13, 13, 17, 17, 17, 17, 19, 19, 23, 23, 23, 23, 29, 29, 29, 29, 29, 29, 31, 31, ]; const P: Primes = Primes::new(); for n in 0..N { println!("{n}"); - assert_eq!(P.smallest_prime_geq(n as u32), Some(NEXT_PRIME[n])); + assert_eq!(P.next_prime(n as u32), Some(NEXT_PRIME[n])); } } @@ -533,7 +562,7 @@ mod test { } #[test] - fn verity_as_slice() { + fn verify_as_slice() { const N: usize = 10; const P: Primes = Primes::new(); const A: [Underlying; N] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29];