SHARE
TWEET

Untitled

a guest Sep 23rd, 2019 85 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /* author Bill Wood */
  2.  
  3. use rayon::prelude::*;
  4. use std::io::Write;
  5.  
  6. fn main() {
  7.     let mut args = std::env::args()
  8.         .skip(1)
  9.         .map(|a| a
  10.             .split('_')
  11.             .fold("".to_string(), |a, s| a + s)
  12.             .parse::<usize>());
  13.  
  14.     if args.len() <= 3 {
  15.         if let (Ok(mut low), Ok(mut high), Ok(seg_size_kb), ) = (
  16.                 args.next().unwrap_or(Ok(0)),
  17.                 args.next().unwrap_or(Ok(100000)),
  18.                 args.next().unwrap_or(Ok(128)),) {
  19.             println!("S15: Finding sexy primes with args {} {} {}", low, high, seg_size_kb, );
  20.  
  21.             let count_two = if low <= 2 && high >= 2 { 1 } else { 0 };
  22.             low = if low < 3 { 3 } else { low | 1 };   // lowest odd >= low
  23.             high = (high - 1) | 1;  // highest odd <= high
  24.  
  25.             // find sieving primes
  26.             let sieve_primes = simple_sieve(high);
  27.  
  28.             // prepare for parallel segment sieve
  29.             let mut candidates_per_seg = seg_size_kb*1024*2;
  30.             while (high - low)/candidates_per_seg > 100_000_000 {   // increase seg size to reduce mem usage
  31.                 candidates_per_seg *= 2;
  32.             }
  33.             dbg!(candidates_per_seg/1024/2);
  34.             let segs: Vec<(usize, usize)> = { low..=high }
  35.                     .step_by(candidates_per_seg)
  36.                     .map(|seg_low| {
  37.                             let seg_high = seg_low + candidates_per_seg - 2;
  38.                             (seg_low, if seg_high > high { high } else { seg_high } + 6)
  39.                         })
  40.                     .collect();
  41.  
  42.             // do parallel segment sieve
  43.             let (sexy_count, nprimes) =
  44.                 segs.into_par_iter()
  45.                     .map(|(low, high)| {
  46.                         // sieve this segment
  47.                         let sieve_seg = seg_sieve(low, high, &sieve_primes);
  48.                         // find sexy primes in this segment
  49.                         let (mut sexy_count, mut nprimes) = (0_usize, 0_usize);
  50.                         for p in 0..sieve_seg.len() - 6/2 {
  51.                             if unsafe { !sieve_seg.get_unchecked(p) } {
  52.                                 nprimes += 1;
  53.                                 // print!("{} ", (p << 1) + low);
  54.                                 if unsafe { !sieve_seg.get_unchecked(p + 6/2) } {
  55.                                     sexy_count += 1;
  56.                                 }
  57.                             }
  58.                         }
  59.                         (sexy_count, nprimes)
  60.                     })
  61.                     .reduce(|| (0, 0),
  62.                             |a, b| (a.0 + b.0, a.1 + b.1));
  63.  
  64.             println!("\n{} of {} primes are sexy", sexy_count, nprimes + count_two);
  65.             return;
  66.         }
  67.     }
  68.  
  69.     writeln!(std::io::stderr(), "Usage: sexy_primes [low [high [seg_size_kb]]]").unwrap();
  70.     std::process::exit(1);
  71. }
  72.  
  73. fn simple_sieve(mut high: usize) -> Vec<u32> {
  74.     // x where x: odd and (x+2)^2 > high and x^2 < high
  75.     high = (sqrt(high) - 1) | 1;
  76.     dbg!(high);
  77.     let mut pvec = vec![false; (high + 1)/2];
  78.     let mut primes = vec![];
  79.     let n = high >> 1;
  80.     let mut i = 3;
  81.     while i*i <= high {
  82.         if unsafe { !pvec.get_unchecked(i >> 1) } {
  83.             primes.push(i as u32);
  84.             let mut j = (i*i) >> 1;
  85.             while j <= n {
  86.                 unsafe { *pvec.get_unchecked_mut(j) = true };
  87.                 j += i;
  88.             }
  89.         }
  90.         i += 2;
  91.     }
  92.     i >>= 1;
  93.     while i <= n {
  94.         if unsafe { !pvec.get_unchecked(i) } {
  95.             primes.push(((i << 1) | 1) as u32);
  96.         }
  97.         i += 1;
  98.     }
  99.     primes
  100. }
  101.  
  102. /*
  103.     3   5   7
  104. 1
  105. 2
  106. 3
  107. 5   *
  108. 7       *
  109. 9   x
  110. 11  *   1   *
  111. 13      2
  112. 15  x   x
  113. 17  *   3
  114. 19      4
  115. 21  x       x
  116. 23  *   5
  117. 25      x
  118. 27  x
  119. 29  *   6
  120. 31      7
  121. 33  x
  122. 35      x   x
  123. 37      0
  124. 39  x
  125. 41      1
  126. 43      2
  127. 45  x   x
  128. 47      3
  129. 49      4    x
  130. 51  x
  131. 53      5
  132. 55      x
  133. 57  x
  134. 59      6
  135. 61      7
  136. 63  x       x
  137. 65      x
  138. 67      *
  139. 69  x
  140. 71      1
  141. 73      2
  142. 75  x   x
  143. 77      3   x
  144. 79      4
  145. 81  x
  146. 83      5
  147. 85      x
  148. 87  x
  149. 89      6
  150. 91      7   x
  151. 93  x
  152. 95      x
  153. 97      *
  154. 99  x
  155. 101     1
  156. 103     2
  157. 105 x   x   x
  158. 107     3
  159. 109     4
  160. 111 x
  161. 113     5
  162. 115     x
  163. 117 x
  164. 119     6   x
  165. 121     7
  166. 123 x
  167. 125     x
  168. 127     *
  169. 129 x
  170. 131     1
  171. 133     2   x
  172. 135 x   x
  173. 137     3
  174. 139     4
  175. 141 x
  176. 143     5
  177. 145     x
  178. 147 x       x
  179. 149     6
  180. 151     7
  181. 153 x
  182. 155     x
  183. 157     *
  184. 159 x
  185. 161     1   x
  186. 163     2
  187. 165 x   x
  188. 167     3
  189. 169     4
  190. 171 x
  191. 173     5
  192. 175     x   x
  193. 177 x
  194. 179     6
  195. 181     7
  196. 183 x
  197. 185     x
  198. 187     *
  199. 189 x       x
  200. 191
  201. 193
  202. 195 x   x
  203. 197
  204. 199
  205. 201 x
  206. 203         x
  207. 205     x
  208. 207 x
  209. 209
  210. 211
  211. 213 x
  212. 215     x
  213. 217         x
  214. 219 x
  215. 221
  216. 223
  217. 225 x   x
  218. 227
  219. 229
  220. 231 x       x
  221. 233
  222. 235     x
  223. 237 x
  224. 239
  225. 241
  226. 243 x
  227. 245     x   x
  228. 247
  229. 249 x
  230. 251
  231. 253
  232. 255 x   x
  233. 257
  234. 259         x
  235. */
  236. // const WHEEL: [usize; 2] = [ 1, 2 ];
  237. // const WHEEL_START: usize = 5;
  238. // const WHEEL: [usize; 8] =  [ 4/2, 2/2, 4/2, 2/2, 4/2, 6/2, 2/2, 6/2 ];
  239. // 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 ];
  240. // 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 ];
  241. // const WHEEL_SYNC: [usize; 15] =  [ 0, 1, 1, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 0, 0 ];
  242. // const WHEEL_SYNC: [usize; 15] = [ 0, 1, 0, 0, 1, 0, 0, 1, 0, 2, 1, 0, 0, 2, 1 ];
  243. // 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) ];
  244. const WHEEL: [usize; 8] =  [ 4, 2, 4, 2, 4, 6, 2, 6 ];
  245. const WHEEL_SYNC: [usize; 15] =  [ 0, 8, 1, 2, 8, 3, 4, 8, 5, 8, 8, 6, 7, 8, 8 ];
  246. const WHEEL_START: usize = 7;
  247. // 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 ];
  248. // const WHEEL_START: usize = 11;
  249. fn seg_sieve(low: usize, high: usize, sieve_primes: &Vec<u32>) -> Vec<bool> {
  250.     let n = (high + 2 - low) >> 1;
  251.     let mut sieve_seg = vec![false; n];
  252.     // dbg!(low);
  253.     for p in sieve_primes {
  254.         let p = *p as usize;
  255.         let mut j = p*p;
  256.         if j > high {
  257.             break;
  258.         }
  259.         if j < low {
  260.             // j = low - 1 + p - (low - 1)%p;
  261.             // j += (low - 1 - j + p*2)/2/p*p*2;
  262.             j = p*(((low - 1)/p + 1) | 1);
  263.         }
  264.         if p < WHEEL_START {
  265.             j = (j - low) >> 1;
  266.             while j < n {
  267.                 unsafe { *sieve_seg.get_unchecked_mut(j) = true };
  268.                 j += p;
  269.             }
  270.         } else {
  271.             // dbg!();
  272.             // dbg!(p);
  273.             // let mut wheel_index = (((j << 1) + low - WHEEL_START)%(WHEEL_SYNC.len()*2)) >> 1;
  274.             // wheel_index = WHEEL_SYNC[wheel_index];
  275.             let mut wheel_index = WHEEL_SYNC[(p - WHEEL_START)/2%(WHEEL_SYNC.len())];
  276.             // wheel_index = match p {
  277.             //     7 => 0,
  278.             //     11 => 1,
  279.             //     13 => 2,
  280.             //     17 => 3,
  281.             //     19 => 4,
  282.             //     23 => 5,
  283.             //     29 => 6,
  284.             //     31 => 7,
  285.             //     37 => 0,
  286.             //     _ => 0,
  287.             // };
  288.             // dbg!(wheel_index);
  289.             while j <= high {
  290.                 unsafe { *sieve_seg.get_unchecked_mut((j - low) >> 1) = true };
  291.                 j += p*WHEEL[wheel_index];
  292.                 wheel_index = if wheel_index == WHEEL.len() - 1 { 0 } else { wheel_index + 1 };
  293.             }
  294.             // while j < n {
  295.             //     dbg!((j << 1) + low, wheel_index, WHEEL[wheel_index]);
  296.             //     unsafe { *sieve_seg.get_unchecked_mut(j) = true };
  297.             //     j += p*WHEEL[wheel_index];
  298.             //     // wheel_index += WHEEL[wheel_index];
  299.             //     // wheel_index %= WHEEL.len();
  300.             //     // wheel_index = if wheel_index == WHEEL.len() - 1 { 0 } else { wheel_index + 1 };
  301.             // }
  302.         }
  303.     }
  304.     sieve_seg
  305. }
  306.  
  307. fn sqrt(n: usize) -> usize {
  308.     let (mut r1, mut r) = (n, n + 1);
  309.     while r1 < r {
  310.         r = r1;
  311.         r1 = (r1 + n/r1) >> 1
  312.     }
  313.     r
  314. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top