/* 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(1200)),
                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: [usize; 8] =  [ 4, 2, 4, 2, 4, 6, 2, 6 ];
// 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_SYNC: [usize; 15] =  [ 4, 0, 2, 4, 0, 2, 4, 0, 6, 0, 0, 2, 6, 0, 0 ];
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 f = j/p;
            let mut wheel_index = WHEEL_SYNC[(j - WHEEL_START)/2%(WHEEL_SYNC.len())/2];
            wheel_index = match p {
                7 => 0,
                11 => 1,
                13 => 2,
                17 => 3,
                19 => 4,
                23 => 5,
                29 => 6,
                31 => 7,
                _ => 0,
            };
            dbg!(wheel_index);
            while p*f <= high {
                dbg!(f, p*f);
                let j = (p*f - low) >> 1;
                unsafe { *sieve_seg.get_unchecked_mut(j) = true };
                f += WHEEL[wheel_index];
                wheel_index = (wheel_index + 1)%WHEEL.len();
            }
            // 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
}