Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /* 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(0)),
- args.next().unwrap_or(Ok(100000)),
- args.next().unwrap_or(Ok(128)),) {
- println!("S15: Finding sexy primes with args {} {} {}", low, high, seg_size_kb, );
- let count_two = if low <= 2 && high >= 2 { 1 } else { 0 };
- 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 mut candidates_per_seg = seg_size_kb*1024*2;
- while (high - low)/candidates_per_seg > 100_000_000 { // increase seg size to reduce mem usage
- candidates_per_seg *= 2;
- }
- dbg!(candidates_per_seg/1024/2);
- 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;
- // print!("{} ", (p << 1) + low);
- 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!("\n{} of {} primes are sexy", sexy_count, nprimes + count_two);
- 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;
- dbg!(high);
- 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
- }
- /*
- 3 5 7
- 1
- 2
- 3
- 5 *
- 7 *
- 9 x
- 11 * 1 *
- 13 2
- 15 x x
- 17 * 3
- 19 4
- 21 x x
- 23 * 5
- 25 x
- 27 x
- 29 * 6
- 31 7
- 33 x
- 35 x x
- 37 0
- 39 x
- 41 1
- 43 2
- 45 x x
- 47 3
- 49 4 x
- 51 x
- 53 5
- 55 x
- 57 x
- 59 6
- 61 7
- 63 x x
- 65 x
- 67 *
- 69 x
- 71 1
- 73 2
- 75 x x
- 77 3 x
- 79 4
- 81 x
- 83 5
- 85 x
- 87 x
- 89 6
- 91 7 x
- 93 x
- 95 x
- 97 *
- 99 x
- 101 1
- 103 2
- 105 x x x
- 107 3
- 109 4
- 111 x
- 113 5
- 115 x
- 117 x
- 119 6 x
- 121 7
- 123 x
- 125 x
- 127 *
- 129 x
- 131 1
- 133 2 x
- 135 x x
- 137 3
- 139 4
- 141 x
- 143 5
- 145 x
- 147 x x
- 149 6
- 151 7
- 153 x
- 155 x
- 157 *
- 159 x
- 161 1 x
- 163 2
- 165 x x
- 167 3
- 169 4
- 171 x
- 173 5
- 175 x x
- 177 x
- 179 6
- 181 7
- 183 x
- 185 x
- 187 *
- 189 x x
- 191
- 193
- 195 x x
- 197
- 199
- 201 x
- 203 x
- 205 x
- 207 x
- 209
- 211
- 213 x
- 215 x
- 217 x
- 219 x
- 221
- 223
- 225 x x
- 227
- 229
- 231 x x
- 233
- 235 x
- 237 x
- 239
- 241
- 243 x
- 245 x x
- 247
- 249 x
- 251
- 253
- 255 x x
- 257
- 259 x
- */
- // const WHEEL: [usize; 2] = [ 1, 2 ];
- // const WHEEL_START: usize = 5;
- // const WHEEL: [usize; 8] = [ 4/2, 2/2, 4/2, 2/2, 4/2, 6/2, 2/2, 6/2 ];
- // const WHEEL: [usize; 15] = [ 4/2, 1, 2/2, 4/2, 1, 2/2, 4/2, 1, 6/2, 2, 1, 2/2, 6/2, 2, 1 ];
- // const WHEEL: [usize; 15] = [ 4/2, 0, 2/2, 4/2, 0, 2/2, 4/2, 0, 6/2, 0, 0, 2/2, 6/2, 0, 0 ];
- // const WHEEL_SYNC: [usize; 15] = [ 0, 1, 1, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 0, 0 ];
- // const WHEEL_SYNC: [usize; 15] = [ 0, 1, 0, 0, 1, 0, 0, 1, 0, 2, 1, 0, 0, 2, 1 ];
- // const WHEEL_SYNC: [(usize, usize); 15] = [ (0,0), (1,1), (0,1), (0,2), (1,3), (0,3), (0,4), (1,5), (0,5), (2,6), (1,6), (0,6), (0,7), (2,0), (1,0) ];
- const WHEEL: [usize; 8] = [ 4, 2, 4, 2, 4, 6, 2, 6 ];
- const WHEEL_SYNC: [usize; 15] = [ 0, 8, 1, 2, 8, 3, 4, 8, 5, 8, 8, 6, 7, 8, 8 ];
- const WHEEL_START: usize = 7;
- // const WHEEL: [usize; 48] = [ 2/2, 4/2, 2/2, 4/2, 6/2, 2/2, 6/2, 4/2, 2/2, 4/2, 6/2, 6/2, 2/2, 6/2, 4/2, 2/2, 6/2, 4/2, 6/2, 8/2, 4/2, 2/2, 4/2, 2/2, 4/2, 8/2, 6/2, 4/2, 6/2, 2/2, 4/2, 6/2, 2/2, 6/2, 6/2, 4/2, 2/2, 4/2, 6/2, 2/2, 6/2, 4/2, 2/2, 4/2, 2/2, 10/2, 2/2, 10/2 ];
- // const WHEEL_START: usize = 11;
- 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];
- // dbg!(low);
- for p in sieve_primes {
- let p = *p as usize;
- let mut j = p*p;
- if j > high {
- break;
- }
- if j < low {
- // j = low - 1 + p - (low - 1)%p;
- // j += (low - 1 - j + p*2)/2/p*p*2;
- j = p*(((low - 1)/p + 1) | 1);
- }
- if p < WHEEL_START {
- j = (j - low) >> 1;
- while j < n {
- unsafe { *sieve_seg.get_unchecked_mut(j) = true };
- j += p;
- }
- } else {
- // dbg!();
- // dbg!(p);
- // let mut wheel_index = (((j << 1) + low - WHEEL_START)%(WHEEL_SYNC.len()*2)) >> 1;
- // wheel_index = WHEEL_SYNC[wheel_index];
- let mut wheel_index = WHEEL_SYNC[(p - WHEEL_START)/2%(WHEEL_SYNC.len())];
- // wheel_index = match p {
- // 7 => 0,
- // 11 => 1,
- // 13 => 2,
- // 17 => 3,
- // 19 => 4,
- // 23 => 5,
- // 29 => 6,
- // 31 => 7,
- // 37 => 0,
- // _ => 0,
- // };
- // dbg!(wheel_index);
- while j <= high {
- unsafe { *sieve_seg.get_unchecked_mut((j - low) >> 1) = true };
- j += p*WHEEL[wheel_index];
- wheel_index = if wheel_index == WHEEL.len() - 1 { 0 } else { wheel_index + 1 };
- }
- // while j < n {
- // dbg!((j << 1) + low, wheel_index, WHEEL[wheel_index]);
- // unsafe { *sieve_seg.get_unchecked_mut(j) = true };
- // j += p*WHEEL[wheel_index];
- // // wheel_index += WHEEL[wheel_index];
- // // wheel_index %= WHEEL.len();
- // // wheel_index = if wheel_index == WHEEL.len() - 1 { 0 } else { wheel_index + 1 };
- // }
- }
- }
- 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
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement