primeutils/
lib.rs

1use std::sync::{Arc, Mutex};
2use std::thread;
3
4mod bits;
5mod cpu;
6
7// Check if a number is prime
8pub fn is_prime(num: u64) -> bool {
9
10  // If num is 2, then it is prime
11  if num == 2 { return true }
12  // If num is less than 2 or is even, then it isn't prime
13  if num < 2 || num % 2 == 0 { return false }
14
15  // For each odd number between 3 and the square root of num,
16  // if it's a divisor of num then num isn't prime
17  let mut i: u64 = 3;
18  while i < (num as f64).sqrt() as u64 {
19    
20    if num % i == 0 {
21      return false;
22    }
23
24    i += 2;
25  }
26
27  // Else num is prime
28  true
29}
30
31// Split a number into its prime factors
32pub fn split_into_factors(num: u64) -> Vec<u64> {
33
34  // Duplicate num as mutable
35  let mut num: u64 = num;
36  // Create a vector to store the factors
37  let mut factors: Vec<u64> = Vec::new();
38
39  if num <= 1 { return factors }
40
41  // First, while the number is divisible by two,
42  // divide it by two and push 2 to the vector
43  while num % 2 == 0 {
44    factors.push(2);
45    num /= 2;
46  }
47
48  // Then, for each odd number between 3 and the square root of the current num,
49  // if it's a divisor of num divide num by that number and add the number to the vector
50  let mut i: u64 = 3;
51  let mut num_sqrt: u64 = (num as f64).sqrt() as u64;
52  while i <= num_sqrt {
53
54    while num % i == 0 {
55      factors.push(i);
56      num /= i;
57      num_sqrt = (num as f64).sqrt() as u64;
58    }
59
60    i += 2;
61  }
62
63  // If num is a prime number add itself to the factors
64  if num > 1 {
65    factors.push(num);
66  }
67
68  // Return the vector
69  factors
70}
71
72// Greatest Common Divisor
73pub fn gcd(x: u64, y: u64) -> u64 {
74
75  // If any of x or y is 0, return the other one
76  if x == 0 || y == 0 { return x + y }
77  // If the numbers are equal, return any of them
78  if x == y { return x }
79
80  // Declare x and y as mutable
81  let (mut x, mut y) = (x, y);
82
83  // While y is not 0; x = y, y = x mod y
84  while y != 0 {
85    (x, y) = (y, x % y);
86  }
87  
88  // Return x
89  x
90}
91
92// Least Common Multiple
93pub fn lcm(x: u64, y: u64) -> u128 {
94  if x == 0 && y == 0 { return 0 }
95  x as u128 * ((y as u128) / (gcd(x, y) as u128))
96}
97
98// Initial sieve
99fn simple_sieve(size: u32) -> Vec<u32> {
100
101  // Create a vector to store all prime numbers found
102  let mut primes: Vec<u32> = Vec::new();
103
104  // Handle specific scenarios 
105  if size >= 2 { primes.push(2) }
106  if size <= 2 { return primes }
107
108  // Create the sieve
109  let mut sieve: Vec<u8> = vec![0xff; ((size + 13) / 16) as usize];
110  // Set the last bits that doesn't have to be sieved to not prime
111  if sieve.len() != 0 {
112    bits::unset_last_bits(sieve.last_mut().unwrap(), (7 - ((size - 3) / 2) % 8) as u8);
113  }
114
115  for i in 0..sieve.len() {
116    for bit in 0..8 {
117
118      // If the bit corresponding to this number is unset,
119      // it is composite, so continue.
120      if bits::is_bit_unset(&sieve[i as usize], bit) { continue }
121
122      // If it is prime, add it to the vector
123      let i: u32 = (i*16) as u32 + (bit*2) as u32 + 3;
124      primes.push(i);
125
126      // If i to the power of 2 is greater than the last element,
127      // i can't be divisor of any of the remaining numbers, so continue.
128      if (i as u64).pow(2) > size as u64 { continue }
129
130      // Unset the bit corresponding to all multiples of i
131      let mut multiple: u32 = i.pow(2);
132      while multiple <= size {
133        bits::unset_bit(&mut sieve[((multiple - 3) / 16) as usize], (((multiple - 3) / 2) % 8) as u8);
134        // Add 2 times i because even numbers are not on the sieve
135        multiple += 2*i;
136      }
137
138    }
139  }
140
141  primes
142}
143
144// Sieve a segment
145fn segment_sieve(sieve: &mut Vec<u8>, primes: &Vec<u32>, low: usize, high: usize) -> u32 {
146
147  // Handle specific scenarios
148  if low == 2 && high == 2 { return 1 }
149  else if low > high { return 0 }
150  else if low == high && low % 2 == 0 { return 0 }
151
152  // Update low and high to be odd numbers
153  let size: usize = std::cmp::max((high - low).div_ceil(16), 1);
154  let low: usize = if low < 3 { 3 } else if low == high { low } else if low % 2 == 0 { low + 1 } else { low };
155  let high: usize = if low == high { high } else if high % 2 == 0 { high - 1 } else { high };
156
157  // Fill the sieve and reset the count
158  assert!(size <= sieve.len());
159  sieve.fill(0xff);
160  let mut count: u32 = 0;
161
162  // Set the last bits that doesn't have to be sieved to not prime
163  bits::unset_last_bits(&mut sieve[size - 1], (7 - ((high - low) / 2) % 8) as u8);
164
165  // For each prime (skipping the 2, which is not needed on odd numbers)
166  for prime in primes.iter().skip(1) {
167    // Get the first multiple
168    let mut multiple: usize = low.div_ceil(*prime as usize) * (*prime as usize);
169    if multiple % 2 == 0 { multiple += *prime as usize }
170    
171    // Mark all multiples as not primes
172    while multiple <= high {
173      bits::unset_bit(&mut sieve[(multiple - low) / 16], (((multiple - low) / 2) % 8) as u8);
174      multiple += 2 * (*prime as usize);
175    }
176  }
177
178  // Count how many primes are there in this segment
179  for byte in sieve.iter().take(size) {
180    count += bits::count_set_bits(&byte) as u32;
181  }
182
183  count
184}
185
186// Count the number of prime numbers below or equal to limit
187pub fn count_primes(limit: usize, start: Option<usize>, threads: Option<usize>, cache: Option<usize>) -> usize {
188
189  if limit < 2 { return 0 }
190
191  let threads: usize = threads.unwrap_or(cpu::get_cores());
192  let cache: usize = cache.unwrap_or(cpu::get_cache_size());
193  let start: usize = std::cmp::max(start.unwrap_or(2), 2);
194  
195  let sqrt: u32 = (limit as f64).sqrt() as u32;
196  let segment_size: usize = std::cmp::min(std::cmp::max(sqrt as usize, cache * 16), limit - std::cmp::max(sqrt as usize, start - 1)).div_ceil(16) * 16;
197  let segment_sieve_size: usize = segment_size.div_ceil(16);
198  
199  let small_primes: Arc<Vec<u32>> = Arc::new(simple_sieve(sqrt));
200  let count: Arc<Mutex<usize>> = Arc::new(Mutex::new(small_primes.len()));
201  let iter: Arc<Mutex<Option<u32>>> = Arc::new(Mutex::new(Some(0)));
202
203  if start > 2 {
204    let mut num = count.lock().unwrap();
205    *num -= if start > sqrt as usize {
206      small_primes.len()
207    }
208    else {
209      small_primes.iter().filter(|&&x| (x as usize) < start).count()
210    };
211  }
212
213  let mut handles = vec![];
214
215  for _ in 0..threads {
216
217    let mut sieve: Vec<u8> = vec![0xff; segment_sieve_size];
218    let count: Arc<Mutex<usize>> = Arc::clone(&count);
219    let iter: Arc<Mutex<Option<u32>>> = Arc::clone(&iter);
220    let small_primes: Arc<Vec<u32>> = Arc::clone(&small_primes);
221
222    let handle = thread::spawn(move || {
223
224      let mut low: usize;
225      let mut high: usize;
226
227      loop {
228        {
229          let mut iter = iter.lock().unwrap();
230          if let Option::None = *iter {
231            break;
232          }
233
234          low = start + (iter.unwrap() as usize * segment_size);
235          high = std::cmp::min(low + segment_size - 1, limit);
236          if high >= limit {
237            *iter = None;
238          }
239          else {
240            *iter = Some(iter.unwrap() + 1);
241          }
242        }
243
244        let current_count = segment_sieve(&mut sieve, &small_primes, low, high) as usize;
245
246        {
247          let mut num = count.lock().unwrap();
248          *num += current_count;
249        }
250      }
251    });
252
253    handles.push(handle);
254  }
255
256  for handle in handles {
257    handle.join().unwrap();
258  }
259
260  let num = count.lock().unwrap();
261  num.clone()
262}
263
264#[cfg(test)]
265mod tests {
266  use crate::*;
267
268  #[test]
269  fn test_is_prime() {
270    assert_eq!(is_prime(0), false);
271    assert_eq!(is_prime(1), false);
272    assert_eq!(is_prime(2), true);
273    assert_eq!(is_prime(3), true);
274    assert_eq!(is_prime(4), false);
275    assert_eq!(is_prime(5), true);
276    assert_eq!(is_prime(6), false);
277    assert_eq!(is_prime(7), true);
278    assert_eq!(is_prime(99), false);
279    assert_eq!(is_prime(100), false);
280    assert_eq!(is_prime(101), true);
281    assert_eq!(is_prime(102), false);
282    assert_eq!(is_prime(103), true);
283    assert_eq!(is_prime(4_294_967_290), false);
284    assert_eq!(is_prime(4_294_967_291), true);
285    assert_eq!(is_prime(4_294_967_292), false);
286    assert_eq!(is_prime(4_294_967_295), false);
287    assert_eq!(is_prime(4_294_967_296), false);
288    assert_eq!(is_prime(18_446_744_073_709_551_614), false);
289    assert_eq!(is_prime(18_446_744_073_709_551_615), false);
290  }
291
292  #[test]
293  fn test_factors() {
294    assert_eq!(split_into_factors(0), vec![]);
295    assert_eq!(split_into_factors(1), vec![]);
296    assert_eq!(split_into_factors(2), vec![2]);
297    assert_eq!(split_into_factors(3), vec![3]);
298    assert_eq!(split_into_factors(4), vec![2, 2]);
299    assert_eq!(split_into_factors(5), vec![5]);
300    assert_eq!(split_into_factors(6), vec![2, 3]);
301    assert_eq!(split_into_factors(7), vec![7]);
302    assert_eq!(split_into_factors(99), vec![3, 3, 11]);
303    assert_eq!(split_into_factors(100), vec![2, 2, 5, 5]);
304    assert_eq!(split_into_factors(101), vec![101]);
305    assert_eq!(split_into_factors(102), vec![2, 3, 17]);
306    assert_eq!(split_into_factors(103), vec![103]);
307    assert_eq!(split_into_factors(4_294_967_290), vec![2, 5, 19, 22_605_091]);
308    assert_eq!(split_into_factors(4_294_967_291), vec![4_294_967_291]);
309    assert_eq!(split_into_factors(4_294_967_292), vec![2, 2, 3, 3, 7, 11, 31, 151, 331]);
310    assert_eq!(split_into_factors(4_294_967_295), vec![3, 5, 17, 257, 65_537]);
311    assert_eq!(split_into_factors(4_294_967_296), vec![2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]);
312    assert_eq!(split_into_factors(18_446_744_073_709_551_614), vec![2, 7, 7, 73, 127, 337, 92_737, 649_657]);
313    assert_eq!(split_into_factors(18_446_744_073_709_551_615), vec![3, 5, 17, 257, 641, 65_537, 6_700_417]);
314  }
315
316  #[test]
317  fn test_gcd() {
318    assert_eq!(gcd(0, 0), 0);
319    assert_eq!(gcd(0, 1), 1);
320    assert_eq!(gcd(1, 0), 1);
321    assert_eq!(gcd(1, 1), 1);
322    assert_eq!(gcd(2, 4), 2);
323    assert_eq!(gcd(4, 2), 2);
324    assert_eq!(gcd(3, 9), 3);
325    assert_eq!(gcd(9, 3), 3);
326    assert_eq!(gcd(6, 8), 2);
327    assert_eq!(gcd(7, 13), 1);
328    assert_eq!(gcd(99, 121), 11);
329    assert_eq!(gcd(100, 80), 20);
330    assert_eq!(gcd(101, 103), 1);
331    assert_eq!(gcd(102, 170), 34);
332    assert_eq!(gcd(103, 206), 103);
333    assert_eq!(gcd(1_234_567_890, 987_654_321), 9);
334    assert_eq!(gcd(4_294_967_295, 65_536), 1);
335    assert_eq!(gcd(4_294_967_296, 65_536), 65_536);
336    assert_eq!(gcd(4_294_967_290, 4_294_967_295), 5);
337    assert_eq!(gcd(4_294_967_291, 4_294_967_292), 1);
338    assert_eq!(gcd(4_294_967_292, 4_294_967_296), 4);
339    assert_eq!(gcd(18_446_744_073_709_551_614, 18_446_744_073_709_551_615), 1);
340    assert_eq!(gcd(18_446_744_073_709_551_615, 18_446_744_073_709_551_615), 18_446_744_073_709_551_615);
341  }
342
343  #[test]
344  fn test_lcm() {
345    assert_eq!(lcm(0, 0), 0);
346    assert_eq!(lcm(0, 1), 0);
347    assert_eq!(lcm(1, 0), 0);
348    assert_eq!(lcm(1, 1), 1);
349    assert_eq!(lcm(2, 4), 4);
350    assert_eq!(lcm(4, 2), 4);
351    assert_eq!(lcm(3, 9), 9);
352    assert_eq!(lcm(6, 8), 24);
353    assert_eq!(lcm(7, 13), 91);
354    assert_eq!(lcm(99, 121), 1089);
355    assert_eq!(lcm(100, 80), 400);
356    assert_eq!(lcm(101, 103), 10403);
357    assert_eq!(lcm(102, 170), 510);
358    assert_eq!(lcm(103, 206), 206);
359    assert_eq!(lcm(4_294_967_295, 65_536), 281_474_976_645_120);
360    assert_eq!(lcm(4_294_967_296, 65_536), 4_294_967_296);
361    assert_eq!(lcm(1_234_567_890, 987_654_321), 135_480_701_236_261_410);
362    assert_eq!(lcm(4_294_967_290, 4_294_967_295), 3_689_348_808_728_956_110);
363    assert_eq!(lcm(4_294_967_291, 4_294_967_292), 18_446_744_035_054_845_972);
364    assert_eq!(lcm(4_294_967_292, 4_294_967_296), 4_611_686_014_132_420_608);
365    assert_eq!(lcm(18_446_744_073_709_551_614, 18_446_744_073_709_551_615), 340_282_366_920_938_463_408_034_375_210_639_556_610);
366  }
367
368  #[test]
369  fn test_simple_sieve() {
370    assert_eq!(simple_sieve(0), vec![]);
371    assert_eq!(simple_sieve(1), vec![]);
372    assert_eq!(simple_sieve(2), vec![2]);
373    assert_eq!(simple_sieve(3), vec![2, 3]);
374    assert_eq!(simple_sieve(4), vec![2, 3]);
375    assert_eq!(simple_sieve(10), vec![2, 3, 5, 7]);
376    assert_eq!(simple_sieve(20), vec![2, 3, 5, 7, 11, 13, 17, 19]);
377    assert_eq!(simple_sieve(30), vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29]);
378    assert_eq!(simple_sieve(50), vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]);
379    assert_eq!(simple_sieve(100), vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]);
380    assert_eq!(simple_sieve(101), vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]);
381    assert_eq!(simple_sieve(102), vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]);
382    assert_eq!(simple_sieve(199), vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199]);
383    assert_eq!(simple_sieve(1_000), vec![
384      2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
385      101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
386      211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
387      331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
388      449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577,
389      587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
390      709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
391      853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997
392    ]);
393    assert_eq!(simple_sieve(10_000).len(), 1229);
394    assert_eq!(simple_sieve(100_000).len(), 9592);
395    assert_eq!(simple_sieve(1_000_000).len(), 78498);
396    assert_eq!(simple_sieve(10_000_000).len(), 664579);
397  }
398  
399  #[test]
400  fn test_count_primes() {
401    assert_eq!(count_primes(10, None, None, None), 4);
402    assert_eq!(count_primes(10, Some(2), None, None), 4);
403    assert_eq!(count_primes(10, Some(3), None, None), 3);
404    assert_eq!(count_primes(10, Some(5), None, None), 2);
405    assert_eq!(count_primes(10, Some(11), None, None), 0);
406    assert_eq!(count_primes(100, None, None, None), 25);
407    assert_eq!(count_primes(1000, None, None, None), 168);
408    assert_eq!(count_primes(10_000, None, None, None), 1229);
409    assert_eq!(count_primes(100_000, None, None, None), 9592);
410    assert_eq!(count_primes(1_000_000, None, None, None), 78498);
411    assert_eq!(count_primes(10_000_000, None, None, None), 664579);
412    assert_eq!(count_primes(100, Some(50), None, None), 10);
413    assert_eq!(count_primes(100, Some(97), None, None), 1);
414    assert_eq!(count_primes(100, Some(98), None, None), 0);
415    assert_eq!(count_primes(2, None, None, None), 1);
416    assert_eq!(count_primes(1, None, None, None), 0);
417    assert_eq!(count_primes(0, None, None, None), 0);
418    // Test with explicit threads and cache
419    assert_eq!(count_primes(100, None, Some(1), Some(1)), 25);
420    assert_eq!(count_primes(100, None, Some(4), Some(2)), 25);
421  }
422}
primeutils/lib.rs

primeutils/
lib.rs
