Prime Sieve

Cargo.toml

[package]
name = "era"
version = "0.1.0"
authors = ["Lukas Kalbertodt <[email protected]>"]

[dependencies]
time = "0.1"
bit-vec = "*"

main.rs

#![feature(asm)]
extern crate time;
extern crate bit_vec;

use bit_vec::BitVec;

const N: usize = 1_000_000_000;

fn sieve() -> Box<Iterator<Item=usize>> {
    // Configure the wheel:
    // - `circ` is the circonference of the wheel (product of first primes)
    // - `coprimes`: how many numbers we don't skip
    // - the values of `wheel` are added in order to get the new coprime from
    //   the wheel (e.g. starting from 1: 5, 7, 11, 13, ...)
    // - `wheel_val` contains the primes from the innermost wheel (used to
    //   calculate the represented number from vector index)
    // - wheel position lookup: For any (prime % `circ`) this array specifies
    //   the offset within the wheel. With primes greater than `circ`
    //   wheel_pos[prime % `circ`] will never return `inv`.
    let inv = usize::max_value();

    // ------ [2] wheel ------
    // let first_primes = vec![2];
    // let circ = 2;
    // let coprimes = 1;
    // let wheel = [2];
    // let wheel_val = [1];
    // let wheel_pos = [inv, 0];

    // ------ [2, 3] wheel ------
    // let first_primes = vec![2, 3];
    // let circ = 6;
    // let coprimes = 2;
    // let wheel = [4, 2];
    // let wheel_val = [1, 5];
    // let wheel_pos = [inv, 0, inv, inv, inv, 1];

    // ------ [2, 3, 5] wheel ------
    let first_primes = vec![2, 3, 5];
    let circ = 30;
    let coprimes = 8;
    let wheel = [6, 4, 2, 4, 2, 4, 6, 2];
    let wheel_val = [1, 7, 11, 13, 17, 19, 23, 29];
    let wheel_pos = [
        inv, 0,   inv, inv, inv, inv, inv, 1,   inv, inv,
        inv, 2,   inv, 3,   inv, inv, inv, 4,   inv, 5,
        inv, inv, inv, 6,   inv, inv, inv, inv, inv, 7,
    ];

    // ------ [2, 3, 5, 7] wheel ------
    // let first_primes = vec![2, 3, 5, 7];
    // let circ = 210;
    // let coprimes = 48;
    // let wheel = [
    //     10,  2,   4,   2,   4,   6,   2,   6,   4,   2,   4,   6,   6,   2,
    //     6,   4,   2,   6,   4,   6,   8,   4,   2,   4,   2,   4,   8,   6,
    //     4,   6,   2,   4,   6,   2,   6,   6,   4,   2,   4,   6,   2,   6,
    //     4,   2,   4,   2,   10,  2
    // ];
    // let wheel_val = [
    //     1,   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, 121,
    //     127, 131, 137, 139, 143, 149, 151, 157, 163, 167, 169, 173, 179, 181,
    //     187, 191, 193, 197, 199, 209
    // ];
    // let wheel_pos = [
    //     inv, 0,   inv, inv, inv, inv, inv, inv, inv, inv,
    //     inv, 1,   inv, 2,   inv, inv, inv, 3,   inv, 4,
    //     inv, inv, inv, 5,   inv, inv, inv, inv, inv, 6,
    //     inv, 7,   inv, inv, inv, inv, inv, 8,   inv, inv,
    //     inv, 9,   inv, 10,  inv, inv, inv, 11,  inv, inv,
    //     inv, inv, inv, 12,  inv, inv, inv, inv, inv, 13,
    //     inv, 14,  inv, inv, inv, inv, inv, 15,  inv, inv,
    //     inv, 16,  inv, 17,  inv, inv, inv, inv, inv, 18,
    //     inv, inv, inv, 19,  inv, inv, inv, inv, inv, 20,
    //     inv, inv, inv, inv, inv, inv, inv, 21,  inv, inv,
    //     inv, 22,  inv, 23,  inv, inv, inv, 24,  inv, 25,
    //     inv, inv, inv, 26,  inv, inv, inv, inv, inv, inv,
    //     inv, 27,  inv, inv, inv, inv, inv, 28,  inv, inv,
    //     inv, 29,  inv, inv, inv, inv, inv, 30,  inv, 31,
    //     inv, inv, inv, 32,  inv, inv, inv, inv, inv, 33,
    //     inv, 34,  inv, inv, inv, inv, inv, 35,  inv, inv,
    //     inv, inv, inv, 36,  inv, inv, inv, 37,  inv, 38,
    //     inv, inv, inv, 39,  inv, inv, inv, inv, inv, 40,
    //     inv, 41,  inv, inv, inv, inv, inv, 42,  inv, inv,
    //     inv, 43,  inv, 44,  inv, inv, inv, 45,  inv, 46,
    //     inv, inv, inv, inv, inv, inv, inv, inv, inv, 47,
    // ];





    // Prepare the bit vector: By default every number is assumed
    // prime (composite = false). We only save numbers that are coprimes of
    // the first primes that generate the wheel, so we don't need N entries,
    // but only N / ratio many. The first indicies represent: [1, 5, 7, 11, ..]
    // The primes used for wheel generation are added later.
    let len = (coprimes * N) / circ;
    let mut is_composite = BitVec::from_elem(len, false);

    // Iterator that generates the coprimes of the wheel primes
    let coprimes = |index, num, step: usize| (index..).scan(num, move |acc, i| {
        let out = *acc;
        *acc += step * wheel[i % wheel.len()];
        Some(out)  // we add one to match the bitvec index
    });

    // Looping through all coprimes, striking out their multiples
    let sq = (N as f64).sqrt() as usize;
    for (i, n) in coprimes(0, 1, 1).take_while(|&n| n <= sq).enumerate().skip(1) {
        if !is_composite[i] {
            for strike in coprimes(wheel_pos[n % circ], n*n, n).take_while(|&n| n <= N) {
                let idx = wheel_pos[strike % circ] + (strike/circ) * wheel.len();
                is_composite.set(idx, true);
            }
        }
    }

    // Return an iterator that extracts the primes from the bitvector.
    let primes = first_primes.into_iter().chain(is_composite.into_iter()
        .enumerate()
        .skip(1)
        .filter(|&(_, c)| !c)
        .map(move |(i, _)| {
            (i / wheel.len()) * circ + wheel_val[i % wheel.len()]
        })
    );

    Box::new(primes)
}

fn main() {
    let before = time::now();
    let mut res = sieve();

    let after = time::now();
    println!("{}", after - before);

    let first_n = 20;
    println!("{:?}", res.by_ref().take(first_n).collect::<Vec<_>>());
    println!("Prime-Count: {}", res.count() + first_n);

}