Sexy Primes
August 23, 2019
A sexy prime is a prime number p for which p + 6 is also prime. Such primes are called sexy because the Latin word for six is sex.
Your task is to write a program that counts the sexy primes between one billion and two billion. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
In Python. The lazy prime generator is based on a segmented sieve and a little too much code to copy here. Thee deque keeps the tested prime in D[0] and has further 3 next primes to test against.
This takes some 90 secs on my (old) PC and reproduces the answer of ProgrammingPraxis.
from collections import deque start = 10**9 stop = 2 * 10**9 pp = lazy_prime(start) D = deque((next(pp), next(pp), next(pp), next(pp)), maxlen=4) cnt = 0 while D[0] < stop: if D[0] + 6 in D: cnt += 1 D.append(next(pp)) print(cnt)Rust version, takes about 45 secs on my PC:
/* author Bill Wood */ fn sqrt(n: usize) -> usize { let (mut r1, mut r) = (n, n + 1); while r1 < r { r = r1; r1 = (r1 + n/r1) >> 1 } r } fn prime_sieve(begin: usize, end: usize) -> Vec<usize> { let mut sieve = vec![false; end + 1]; let sn = sqrt(end); for i in 2..=sn { if !sieve[i] { let mut j = i*i; while j <= end { sieve[j] = true; j += i; } } } let mut result = vec![]; for i in begin..=end { if !sieve[i] { result.push(i); } } result } fn main() { let primes = prime_sieve(1_000_000_000, 2_000_000_000); let mut count = 0; for (i, x) in primes.iter().enumerate() { for &y in &primes[i + 1..] { if y > x + 6 { break; } if y == x + 6 { count += 1; } } } println!("{} of {} primes are sexy", count, primes.len()); }Output:
Here’s a sieve-based C solution.
#include <stdio.h> #include <stdlib.h> int main(void) { int N = 2000000006; char* array = calloc(N + 1, 1); for (int i = 2; i * i <= N; ++i) { if (array[i]) continue; for (int j = i * i; j <= N; j += i) { array[j] = 1; } } int sum = 0; for (int i = 1000000000; i <= 2000000000; ++i) { sum += !array[i] && !array[i + 6]; } printf("%d\n", sum); return EXIT_SUCCESS; }Output:
I forgot to call free :-(
(my preference is to explicitly clean up, even if the memory is reclaimed by the OS)
Here’s the C++ version of my code above. It uses vector instead of a char array, to reduce the program’s memory usage from 2GB to 250MB, since it can use a single bit for each boolean instead of a byte. Additionally, the execution time on my laptop decreases from 22 seconds for the C version to 17 seconds for the C++ version (both compiled with -O3).
#include <cstdlib> #include <iostream> #include <vector> using namespace std; int main(void) { int N = 2000000006; vector<bool> array(N + 1, false); for (int i = 2; i * i <= N; ++i) { if (array[i]) continue; for (int j = i * i; j <= N; j += i) { array[j] = 1; } } int sum = 0; for (int i = 1000000000; i <= 2000000000; ++i) { sum += !array[i] && !array[i + 6]; } printf("%d\n", sum); return EXIT_SUCCESS; }Output:
Here’s a solution in C that may not quite be in the spirit of the
question since it uses the excellent primesieve library to do all the
heavy lifting but, hey, it works and it’s fast, and it isn’t my
homework. Runs in about half a second on my very modest
computer. Bam!
#include <primesieve.h> #include <inttypes.h> #include <stdio.h> static int sexy_count = 0; void record_prime(uint64_t p) { static uint64_t buf[] = {0, 0, 0}; static short cur = 0; uint64_t q = p - 6; if (q == buf[0] || q == buf[1] || q == buf[2]) sexy_count++; buf[cur] = p; cur = (cur + 1) % 3; } int main() { primesieve_callback_primes(1000000000, 2000000000, &record_prime); printf("%u\n", sexy_count); return 0; }<
pre class="src src-sh">time ./sexy-primes
Output:
Daniel inspired me to re-write the Rust version with a bit vector (home-rolled) instead of a bool vector. His C++ version and this Rust version both run in 25 secs on my home PC. The Rust bool version above runs 3 secs slower (28 secs) on my home PC (timing of 45 secs above was on my work PC).
/* author Bill Wood */ fn sqrt(n: usize) -> usize { let (mut r1, mut r) = (n, n + 1); while r1 < r { r = r1; r1 = (r1 + n/r1) >> 1 } r } struct Bvec { v: Vec<u8> } impl Bvec { fn new(s: usize, val: u8) -> Self { Bvec { v: vec![val; s/8 + if s%8 == 0 { 0 } else { 1 } ] } } #[inline(always)] fn get(&self, index: usize) -> u8 { self.v[index/8] & 1 << index%8 } #[inline(always)] fn set(&mut self, index: usize, val: u8) { self.v[index/8] |= val << index%8; } } fn prime_sieve(begin: usize, end: usize) -> Vec<usize> { let mut sieve = Bvec::new(end + 1, 0); let sn = sqrt(end); for i in 2..=sn { if sieve.get(i) == 0 { let mut j = i*i; while j <= end { sieve.set(j, 1); j += i; } } } let mut result = vec![]; for i in begin..=end { if sieve.get(i) == 0 { result.push(i); } } result } fn main() { let primes = prime_sieve(1_000_000_000, 2_000_000_000); let mut count = 0; for (i, x) in primes.iter().enumerate() { for &y in &primes[i + 1..] { if y > x + 6 { break; } else if y == x + 6 { count += 1; } } } println!("{} of {} primes are sexy", count, primes.len()); }Rust output:
Daniel’s C++ version:
Here’s a basic segmented sieve solution, using a vector (ie. a bitmap representation) and working in segments of 128k entries. Takes a little under 3 seconds on my laptop, versus 28 seconds for Daniel’s C++ solution above. Not competitive with Chaw’s impressively fast primesieve solution yet though, I must have a look at what that is doing:
#include <vector> #include <iostream> #include <algorithm> #include <math.h> std::vector<int> getoddprimes(int max) { std::vector<int> primes; max /= 2; std::vector<bool> a(max); for (int i = 1; i < max; i++) { if (a[i] == 0) { int p = 2*i+1; primes.push_back(p); for (int j = p*p/2; j < max; j += p) { a[j] = 1; } } } return primes; } int main() { int start = 1000000000, end = 2*start+6, block = 128*1024; int primemax = sqrt(end); // Maybe should use isqrt. std::vector<int> primes = getoddprimes(primemax+1); int nprimes = primes.size(); start /= 2; end /= 2; // Just look at odd numbers std::vector<bool> a(end-start); std::vector<int> offsets(nprimes); // Keeps track of sieve position for (int i = 0; i < nprimes; i++) { int p = primes[i]; int offset = p*p/2; if (offset < start) { offset += (start-offset+p-1)/p*p; // Move into region of interest } offsets[i] = offset; } int blockstart = start, sexystart = blockstart+3, count = 0; while (blockstart != end) { int blockend = std::min(end, blockstart + block); for (int i = 0; i < nprimes; i++) { int p = primes[i]; int offset = offsets[i]; // Faster to cache in register while (offset < blockend) { a[offset-start] = 1; // Cross off multiple of p offset += p; } offsets[i] = offset; } for (int i = sexystart; i < blockend; i++) { count += a[i-start-3] == 0 && a[i-start] == 0; // Look backwards } sexystart = blockstart = blockend; } std::cout << count << "\n"; }This Rust one runs in under 0.6 seconds on my Google Pixelbook with Crostini Linux. It’s faster than the previous Rust one above because:
1) it only looks at odd numbers
2) it calculates the primes in segments which is faster because the segments fit in the CPU cache
3) removed bit vectors as ultimately they were slower
4) it uses rayon to parallelize computing the segment sieves and find the sexy primes
/* author Bill Wood */ use rayon::prelude::*; use std::io::Write; fn main() { let mut args = std::env::args() .skip(1) .map(|a| a .split('_') .fold("".to_string(), |a, s| a + s) .parse::<usize>()); if args.len() <= 3 { if let (Ok(mut low), Ok(mut high), Ok(seg_size_kb), ) = ( args.next().unwrap_or(Ok(1_000_000_000)), args.next().unwrap_or(Ok(2_000_000_000)), args.next().unwrap_or(Ok(128)),) { println!("S15: Finding sexy primes with args {} {} {}", low, high, seg_size_kb, ); low = if low < 3 { 3 } else { low | 1 }; // lowest odd >= low high = (high - 1) | 1; // highest odd <= high // find sieving primes let sieve_primes = simple_sieve(high); // prepare for parallel segment sieve let candidates_per_seg = seg_size_kb*1024*2; // generate a Vec of (low, high) segment tuples let segs: Vec<(usize, usize)> = { low..=high } .step_by(candidates_per_seg) .map(|seg_low| { let seg_high = seg_low + candidates_per_seg - 2; (seg_low, if seg_high > high { high } else { seg_high } + 6) }) .collect(); // do parallel segment sieve let (sexy_count, nprimes) = segs.into_par_iter() .map(|(low, high)| { // sieve this segment let sieve_seg = seg_sieve(low, high, &sieve_primes); // find sexy primes in this segment let (mut sexy_count, mut nprimes) = (0_usize, 0_usize); for p in 0..sieve_seg.len() - 6/2 { if unsafe { !sieve_seg.get_unchecked(p) } { nprimes += 1; if unsafe { !sieve_seg.get_unchecked(p + 6/2) } { sexy_count += 1; } } } (sexy_count, nprimes) }) .reduce(|| (0, 0), |a, b| (a.0 + b.0, a.1 + b.1)); println!("{} of {} primes are sexy", sexy_count, nprimes); return; } } writeln!(std::io::stderr(), "Usage: sexy_primes [low [high [seg_size_kb]]]").unwrap(); std::process::exit(1); } fn simple_sieve(mut high: usize) -> Vec<u32> { // x where x: odd and (x+2)^2 > high and x^2 < high high = (sqrt(high) - 1) | 1; let mut pvec = vec![false; (high + 1)/2]; let mut primes = vec![]; let n = high >> 1; let mut i = 3; while i*i <= high { if unsafe { !pvec.get_unchecked(i >> 1) } { primes.push(i as u32); let mut j = (i*i) >> 1; while j <= n { unsafe { *pvec.get_unchecked_mut(j) = true }; j += i; } } i += 2; } i >>= 1; while i <= n { if unsafe { !pvec.get_unchecked(i) } { primes.push(((i << 1) | 1) as u32); } i += 1; } primes } fn seg_sieve(low: usize, high: usize, sieve_primes: &Vec<u32>) -> Vec<bool> { let n = (high + 2 - low) >> 1; let mut sieve_seg = vec![false; n]; for p in sieve_primes { let p = *p as usize; let mut j = p*p; if j > high { break; } if j < low { // lowest odd multiple of p which is >= low j = p*(((low - 1)/p + 1) | 1); } j = (j - low) >> 1; while j < n { unsafe { *sieve_seg.get_unchecked_mut(j) = true }; j += p; } } sieve_seg } fn sqrt(n: usize) -> usize { let (mut r1, mut r) = (n, n + 1); while r1 < r { r = r1; r1 = (r1 + n/r1) >> 1 } r }Playground link: https://play.rust-lang.org/?version=stable&mode=release&edition=2018&gist=63d60dbd91efa0106ee96bba0c6983ed