prime.rs
13.7 kB · rust · 440 lines
1use crate::classics::primes;2use crate::factor::factorize_wide;3use crate::series::li;45/// Returns whether the number is prime, by trial division on the six-step wheel.6pub fn is_prime(number: usize) -> bool {7 if number < 2 {8 return false;9 }10 if number < 4 {11 return true;12 }13 if number.is_multiple_of(2) || number.is_multiple_of(3) {14 return false;15 }16 let mut step = 5;17 while step <= number / step {18 if number.is_multiple_of(step) || number.is_multiple_of(step + 2) {19 return false;20 }21 step += 6;22 }23 true24}2526/// Returns the smallest prime at or above the number.27///28/// ```29/// assert_eq!(mrlynum::prime::prime_from(90), 97);30/// ```31pub fn prime_from(number: usize) -> usize {32 let mut n = number.max(2);33 while !is_prime(n) {34 n += 1;35 }36 n37}3839/// Returns every rectangle of the number as a pair of sides, the shorter first, ascending.40///41/// ```42/// assert_eq!(mrlynum::prime::rectangles(6), vec![(1, 6), (2, 3)]);43/// ```44pub fn rectangles(number: usize) -> Vec<(usize, usize)> {45 let mut out = Vec::new();46 if number == 0 {47 return out;48 }49 let mut a = 1;50 while a <= number / a {51 if number.is_multiple_of(a) {52 out.push((a, number / a));53 }54 a += 1;55 }56 out57}5859/// Returns every pair of primes summing to the number, odd numbers included, the smaller first, ascending.60pub fn splits(number: usize) -> Vec<(usize, usize)> {61 let mut out = Vec::new();62 if number < 4 {63 return out;64 }65 for p in primes(number / 2) {66 if is_prime(number - p) {67 out.push((p, number - p));68 }69 }70 out71}7273/// Returns the smallest pair of positive sides whose squares sum to the number, when one exists.74pub fn squares(number: usize) -> Option<(usize, usize)> {75 let mut a = 1;76 while 2 * a * a <= number {77 let rest = number - a * a;78 let b = rest.isqrt();79 if b * b == rest {80 return Some((a, b));81 }82 a += 1;83 }84 None85}8687/// The prime object: one prime with its rank, the step behind it and the shapes it makes.88#[derive(Clone, Copy, Debug, PartialEq, Eq)]89pub struct Prime {90 /// The prime itself.91 pub value: usize,92 /// The one-based rank in the prime sequence, so two has index one.93 pub index: usize,94 /// The distance from the previous prime, zero for two.95 pub gap: usize,96 /// The twin flag: whether a prime sits exactly two away on either side, false for two.97 pub twin: bool,98 /// The two positive sides whose squares sum to the prime, when they exist.99 pub squares: Option<(usize, usize)>,100}101102/// Returns one prime object for every prime up to and including the limit.103pub fn study(limit: usize) -> Vec<Prime> {104 let list = primes(limit);105 list.iter()106 .enumerate()107 .map(|(i, &value)| Prime {108 value,109 index: i + 1,110 gap: if i == 0 { 0 } else { value - list[i - 1] },111 twin: is_prime(value - 2) || is_prime(value + 2),112 squares: squares(value),113 })114 .collect()115}116117/// The sieve of Eratosthenes taken one prime at a time, each number remembering which prime struck it.118#[derive(Clone, Debug, PartialEq, Eq)]119pub struct Sieve {120 types: Vec<u8>,121 at: usize,122 rank: usize,123 struck: usize,124 done: bool,125}126127impl Sieve {128 /// Starts a sieve over zero through the limit with every number untouched; it is done at once when no prime has its square inside.129 pub fn new(limit: usize) -> Sieve {130 let mut sieve = Sieve {131 types: vec![0; limit + 1],132 at: 2,133 rank: 0,134 struck: 0,135 done: false,136 };137 sieve.settle();138 sieve139 }140 fn settle(&mut self) {141 while self.at < self.types.len() && self.types[self.at] != 0 {142 self.at += 1;143 }144 if self.at * self.at >= self.types.len() {145 for number in 2..self.types.len() {146 if self.types[number] == 0 {147 self.types[number] = 1;148 }149 }150 self.done = true;151 }152 }153 /// Uses the next prime: marks it prime, strikes its untouched multiples from its square with its rank plus one, and returns it; zero once done.154 ///155 /// The strike mark saturates at 255, so it is exact through the 254th prime.156 pub fn step(&mut self) -> usize {157 if self.done {158 return 0;159 }160 let prime = self.at;161 self.rank += 1;162 self.types[prime] = 1;163 self.struck = 0;164 let mark = (self.rank + 1).min(255) as u8;165 let mut multiple = prime * prime;166 while multiple < self.types.len() {167 if self.types[multiple] == 0 {168 self.types[multiple] = mark;169 self.struck += 1;170 }171 multiple += prime;172 }173 self.settle();174 prime175 }176 /// Runs the sieve to the end.177 pub fn finish(&mut self) {178 while !self.done {179 self.step();180 }181 }182 /// Returns whether every number is settled.183 pub fn done(&self) -> bool {184 self.done185 }186 /// Returns the type of every number from zero: zero untouched, one prime, and one past the rank of the prime that struck it.187 pub fn types(&self) -> &[u8] {188 &self.types189 }190 /// Returns the count of numbers marked prime so far.191 pub fn count(&self) -> usize {192 self.types.iter().filter(|&&t| t == 1).count()193 }194 /// Returns the count of numbers the last step struck.195 pub fn struck(&self) -> usize {196 self.struck197 }198 /// Returns the count of primes used so far.199 pub fn rank(&self) -> usize {200 self.rank201 }202}203204/// A number as a pile of stones: its prime factors, whether it is prime, and every rectangle the stones make.205#[derive(Clone, Debug, PartialEq, Eq)]206pub struct Pile {207 /// The count of stones.208 pub number: u64,209 /// The prime and exponent pairs, ascending.210 pub factors: Vec<(u64, u32)>,211 /// Whether the stones make a single row and nothing else.212 pub prime: bool,213 /// Every rectangle as a pair of sides, the shorter first, ascending.214 pub rectangles: Vec<(u64, u64)>,215}216217/// Reads a wide number as a pile of stones, its rectangles built from the divisors of its factorization.218///219/// ```220/// let pile = mrlynum::prime::pile(6);221/// assert_eq!(pile.rectangles, vec![(1, 6), (2, 3)]);222/// assert!(!pile.prime);223/// ```224pub fn pile(number: u64) -> Pile {225 let factors = factorize_wide(number);226 let mut sides = vec![1u64];227 for &(prime, power) in &factors {228 let mut next = Vec::with_capacity(sides.len() * (power as usize + 1));229 for &side in &sides {230 let mut value = side;231 next.push(value);232 for _ in 0..power {233 value *= prime;234 next.push(value);235 }236 }237 sides = next;238 }239 sides.sort_unstable();240 let rectangles = sides241 .iter()242 .filter(|&&a| number > 0 && a <= number / a)243 .map(|&a| (a, number / a))244 .collect();245 Pile {246 number,247 prime: factors.len() == 1 && factors[0].1 == 1,248 factors,249 rectangles,250 }251}252253/// One reading of the prime count against its two smooth guesses.254#[derive(Clone, Copy, Debug, PartialEq)]255pub struct Reading {256 /// The point on the number line.257 pub x: usize,258 /// The count of primes up to it.259 pub pi: usize,260 /// The guess x over ln x.261 pub ratio: f64,262 /// The logarithmic integral.263 pub li: f64,264}265266/// Reads the prime count against x over ln x and li at evenly spaced points from two up to the top, at most the given count of them, the top always last.267///268/// ```269/// let readings = mrlynum::prime::chart(100, 10);270/// assert_eq!((readings.len(), readings[9].pi), (10, 25));271/// ```272pub fn chart(top: usize, bins: usize) -> Vec<Reading> {273 let list = primes(top);274 let step = (top / bins.max(1)).max(1);275 let mut out = Vec::new();276 let mut x = step;277 while x <= top {278 if x >= 2 {279 out.push(Reading {280 x,281 pi: list.partition_point(|&p| p <= x),282 ratio: x as f64 / (x as f64).ln(),283 li: li(x as f64),284 });285 }286 if x == top {287 break;288 }289 x = (x + step).min(top);290 }291 out292}293294#[cfg(test)]295mod tests {296 use super::*;297298 #[test]299 fn is_prime_agrees_with_the_sieve() {300 let sieved = primes(100_000);301 for number in 0..=100_000 {302 assert_eq!(303 is_prime(number),304 sieved.binary_search(&number).is_ok(),305 "{number}"306 );307 }308 }309310 #[test]311 fn the_prime_from_a_number_is_the_first_at_or_above_it() {312 assert_eq!(prime_from(90), 97);313 assert_eq!(prime_from(0), 2);314 assert_eq!(prime_from(41), 41);315 for number in 0..=1_000 {316 let next = prime_from(number);317 assert!(is_prime(next) && next >= number, "{number}");318 assert!((number..next).all(|n| !is_prime(n)), "{number}");319 }320 }321322 #[test]323 fn a_single_rectangle_means_prime_above_one() {324 assert!(rectangles(0).is_empty());325 assert_eq!(rectangles(1), vec![(1, 1)]);326 assert_eq!(327 rectangles(36),328 vec![(1, 36), (2, 18), (3, 12), (4, 9), (6, 6)]329 );330 for number in 2..=1_000 {331 assert_eq!(rectangles(number).len() == 1, is_prime(number), "{number}");332 }333 }334335 #[test]336 fn every_even_number_splits_into_two_primes() {337 for number in (4..=2_000).step_by(2) {338 let pairs = splits(number);339 assert!(!pairs.is_empty(), "{number}");340 for pair in pairs.windows(2) {341 assert!(pair[0].0 < pair[1].0, "{number}");342 }343 for (p, q) in pairs {344 assert!(is_prime(p) && is_prime(q), "{number}");345 assert!(p <= q && p + q == number, "{number}");346 }347 }348 }349350 #[test]351 fn an_odd_number_splits_only_through_the_two() {352 assert_eq!(splits(5), vec![(2, 3)]);353 assert_eq!(splits(9), vec![(2, 7)]);354 assert!(splits(27).is_empty());355 }356357 #[test]358 fn a_prime_is_a_sum_of_two_squares_only_when_it_is_two_or_one_past_a_multiple_of_four() {359 for prime in study(10_000) {360 let expected = prime.value == 2 || prime.value % 4 == 1;361 assert_eq!(prime.squares.is_some(), expected, "{}", prime.value);362 }363 }364365 #[test]366 fn study_pins_the_first_primes() {367 let found = study(13);368 let values: Vec<usize> = found.iter().map(|p| p.value).collect();369 let indices: Vec<usize> = found.iter().map(|p| p.index).collect();370 let gaps: Vec<usize> = found.iter().map(|p| p.gap).collect();371 let twins: Vec<bool> = found.iter().map(|p| p.twin).collect();372 assert_eq!(values, vec![2, 3, 5, 7, 11, 13]);373 assert_eq!(indices, vec![1, 2, 3, 4, 5, 6]);374 assert_eq!(gaps, vec![0, 1, 2, 2, 4, 2]);375 assert_eq!(twins, vec![false, true, true, true, true, true]);376 assert_eq!(found[0].squares, Some((1, 1)));377 assert_eq!(found[2].squares, Some((1, 2)));378 assert_eq!(found[3].squares, None);379 assert!(study(17).last().unwrap().twin);380 assert!(!study(23).last().unwrap().twin);381 }382383 #[test]384 fn the_sieve_strikes_one_prime_at_a_time() {385 let mut sieve = Sieve::new(30);386 assert!(!sieve.done());387 assert_eq!((sieve.step(), sieve.struck(), sieve.count()), (2, 14, 1));388 assert_eq!((sieve.step(), sieve.struck(), sieve.count()), (3, 4, 2));389 assert_eq!((sieve.step(), sieve.struck()), (5, 1));390 assert!(sieve.done());391 assert_eq!((sieve.rank(), sieve.count(), sieve.step()), (3, 10, 0));392 assert_eq!(393 &sieve.types()[..13],394 &[0, 0, 1, 1, 2, 1, 2, 1, 2, 3, 2, 1, 2]395 );396 assert_eq!(sieve.types()[25], 4);397 let mut hundred = Sieve::new(100);398 hundred.finish();399 assert_eq!((hundred.rank(), hundred.count()), (4, 25));400 let listed: Vec<usize> = (0..=100).filter(|&n| hundred.types()[n] == 1).collect();401 assert_eq!(listed, primes(100));402 assert_eq!(Sieve::new(3).count(), 2);403 assert_eq!(Sieve::new(0).count(), 0);404 }405406 #[test]407 fn the_pile_agrees_with_the_rectangles_and_the_wheel() {408 let stones = pile(360);409 assert_eq!(stones.factors, vec![(2, 3), (3, 2), (5, 1)]);410 assert_eq!(stones.rectangles.len(), 12);411 assert_eq!(stones.rectangles[11], (18, 20));412 assert!(pile(13).prime && !pile(1).prime && !pile(0).prime);413 assert!(pile(0).rectangles.is_empty());414 assert_eq!(pile(1).rectangles, vec![(1, 1)]);415 for number in 1..=2_000u64 {416 let want: Vec<(u64, u64)> = rectangles(number as usize)417 .iter()418 .map(|&(a, b)| (a as u64, b as u64))419 .collect();420 assert_eq!(pile(number).rectangles, want, "{number}");421 assert_eq!(pile(number).prime, is_prime(number as usize), "{number}");422 }423 assert_eq!(pile(1_000_000_000_000).rectangles.len(), 85);424 assert!(pile(999_999_999_989).prime);425 }426427 #[test]428 fn the_chart_reads_the_prime_count_at_the_top() {429 let readings = chart(10_000, 400);430 assert_eq!(readings.len(), 400);431 let last = readings[399];432 assert_eq!((last.x, last.pi), (10_000, 1229));433 assert!((last.ratio - 1085.7360).abs() < 1e-3);434 assert!((last.li - 1246.1372).abs() < 1e-3);435 assert_eq!(chart(100_000, 100)[99].pi, 9592);436 assert_eq!(chart(100, 400).len(), 99);437 assert_eq!(chart(7, 3)[0].x, 2);438 assert!(chart(1, 5).is_empty());439 }440}