Layout prime (powers) into fixed-size blocks.

primes.rs

// use std::fmt;
// use std::iter::Iterator;
use std::collections::BinaryHeap;
use std::cmp::Ordering;

const BITS : u32 = 50;
const MAX_PRIME : u32 = 1 << 25;
const MAX_NUM : u64 = 1 << 50;

const ALL_PP : u32 = 256;
const COUNTERED : u64 = 512;

const QSIZE : usize = 4095;

// 56 = 31 + 25 (with 25 + 4 + 1 = 30)
//const MAX_NORMAL : u64 = 31_786_648;
const MAX_NORMAL : u64 = 1 << 36;
//const NUM_BUCKETS : usize = 4;
const NUM_BUCKETS : usize = 16;
//const EVEN_SMALL_LIMIT : u64 = 1 << 18;
const EVEN_SMALL_LIMIT : u64 = !0;

type Primes = Vec<u32>;
type Powers = (usize, [u64; BITS as usize]);

fn get_primes() -> Primes {
    let mut sieve = [true; MAX_PRIME as usize];
    sieve[0] = false;
    sieve[1] = false;

    for p in 0..MAX_PRIME {
        if p*p > MAX_PRIME {
            break
        }
        if sieve[p as usize] {
            for x in (p*p .. MAX_PRIME).step_by(p as usize) {
                sieve[x as usize] = false
            }
        }
    }

    let mut primes = Vec::new();

    for (n, &p) in sieve.iter().enumerate() {
        if p {
            primes.push(n as u32)
        }
    }
    primes
}


fn get_powersof(p : u32) -> Powers {
    let mut n = 0;
    let mut pp = [0u64; BITS as usize];

    let mut i = 1u64;
    while i <= MAX_NUM / p as u64 {
        i *= p as u64;
        pp[n as usize] = i;
        n += 1;
    }
    (n, pp)
}


fn get_usedpowers(primes : &Primes) -> Vec<u64> {
    let mut usedpowers = Vec::new();
    for &p in primes {
        let (n, pp) = get_powersof(p);
        // If p is small then use all powers...
        assert!(n > 1);
        let n = if p < ALL_PP { n } else { n - 1 };
        for i in 0..n {
            usedpowers.push(pp[i])
        }
    }
    usedpowers.sort();
    usedpowers
}


fn get_split(used : &Vec<u64>) -> (Vec<u64>, Vec<u64>, Vec<u64>) {
    let mut small  = Vec::new();
    let mut normal = Vec::new();
    let mut big    = Vec::new();

    for &p in used {
        if p < COUNTERED || (p % 2 == 0 && p < EVEN_SMALL_LIMIT) {
            small.push(p)
        } else if p > MAX_NORMAL || p % 2 == 0 {
            big.push(p)
        } else {
            normal.push(p)
        }
    }
    (small, normal, big)
}


fn recip_sum(pp : &Vec<u64>) -> f64 {
    pp.iter().rev().map(|&p| 1.0 / p as f64).sum()
}


fn bucket_split(pp : &Vec<u64>) -> Vec<Vec<u64>> {
    let mut res : Vec<_> = (0..NUM_BUCKETS).map(|_| Vec::new()).collect();
    // let mut res : [Vec<u64>; NUM_BUCKETS] = Default::default();
    for &p in pp {
        assert!(p & 1 == 1);
        res[p as usize / 2 % NUM_BUCKETS].push(p);
    }
    res
}


struct Queue {
    items : Vec<u64>,
    weight : f64,
}


impl Queue {
    fn new() -> Queue {
        Queue{items: Vec::new(), weight: 0.0}
    }
    fn weight(&self) -> f64 {
        self.weight + (QSIZE - self.items.len()) as f64 / MAX_NORMAL as f64
    }
}


impl Ord for Queue {
    fn cmp(&self, that: &Self) -> Ordering {
        let a_space = self.items.len() < QSIZE;
        let b_space = that.items.len() < QSIZE;
        if a_space != b_space {
            Ord::cmp(&a_space, &b_space)
        }
        else {
            PartialOrd::partial_cmp(&that.weight(), &self.weight()).unwrap()
        }
    }
}


impl PartialOrd for Queue {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}


impl PartialEq for Queue {
    fn eq(&self, other: &Self) -> bool {
        self.weight() == other.weight() && self.items == other.items
    }
}


impl Eq for Queue { }


fn queue_split(pp : &Vec<u64>) -> Vec<Queue> {
    let num_queues = (pp.len() + QSIZE - 1) / QSIZE;
    let mut heap = BinaryHeap::new();
    for _ in 0..num_queues {
        heap.push(Queue::new())
    }
    for &p in pp {
        let mut h = heap.pop().unwrap();
        h.items.push(p);
        h.weight += 1.0 / p as f64;
        heap.push(h);
    }
    heap.into_sorted_vec()
}


fn main() {
    let primes : Primes = get_primes();
    let used = get_usedpowers(&primes);

    println!("Num primes = {}", primes.len());
    println!("Num ppused = {}", used.len());

    let (small, normal, big) = get_split(&used);
    println!("Num small, normal, big = {}, {}, {}",
             small.len(), normal.len(), big.len());
    println!("Total rate {:.8}", recip_sum(&used));
    println!("Small rate {:.8}", recip_sum(&small));
    println!("Norml rate {:.8}", recip_sum(&normal));
    println!("Big 1/rate {:.8}", 1.0 / recip_sum(&big));

    let buckets = bucket_split(&normal);
    for (b, pp) in buckets.iter().enumerate() {
        let queues = queue_split(pp);
        print!("Bucket {:2} has {} ({}):", b, pp.len(), queues.len());
        println!(" {} {} {:.6} {:.6}",
                 queues.iter().map(|q| q.items.len()).max().unwrap(),
                 queues.iter().map(|q| q.items.len()).min().unwrap(),
                 1.0 / queues.iter().map(|q| q.weight()).min_by(
                     |a,b| a.partial_cmp(b).unwrap()).unwrap(),
                 1.0 / queues.iter().map(|q| q.weight()).max_by(
                     |a,b| a.partial_cmp(b).unwrap()).unwrap());
    }
}