## Wednesday, June 17, 2015

### Compressed Sieve of Eratosthenes

In my previous post, we used "2-3-5-7" wheel factorization to produce a faster Sieve of Eratosthenes with a compression factor of 2 (by storing only the odd numbers). Samuel Tardieu reminded me that the sieve variation used by Factor uses less memory with a higher compression factor of 3.75.

I thought would be a fun implementation to show below.

## Version 6

The "2-3-5" wheel with the first 30 numbers that are not divisible by 2, 3, or 5.

```CONSTANT: wheel-2-3-5 \$[
30 [1,b] [
{ 2 3 5 } [ divisor? ] with any? not
] B{ } filter-as
]```
Note: We use 30 here because the "`{ 2 3 5 } product`" cycle is 30.

If you look at the wheel, you see that there are 8 candidates in every 30 numbers:

```IN: scratchpad wheel-2-3-5 .
{ 1 7 11 13 17 19 23 29 }```

The number `8` is interesting because it suggests that we can use a single byte to store the 8 candidates in each block of 30 numbers, using a bitmask for possible primes or `f`:

```CONSTANT: masks \$[
30 [0,b) [
wheel-2-3-5 index [ 7 swap - 2^ ] [ f ] if*
] map
]```

We can index into this array to find the byte/mask for every number:

```: byte-mask ( n -- byte mask/f )

Using the byte/mask we can check if a number is marked in the sieve byte-array.

```:: marked? ( n sieve -- ? )
mask [ byte sieve nth mask bitand zero? not ] [ f ] if ;```

And we can unmark a number in the sieve in a similar fashion:

```:: unmark ( n sieve -- )
mask [ byte sieve [ mask bitnot bitand ] change-nth ] when ;```

Using that, we can unmark multiples of prime numbers:

```:: unmark-multiples ( i upper sieve -- )
i sq upper i 2 * <range> [ sieve unmark ] each ;```

That's all we need to build our sieve, starting from a byte-array initialized with all numbers marked and then progressively unmarking multiples of primes:

```:: sieve6 ( n -- sieve )
n 30 /i 1 + [ 0xff ] B{ } replicate-as :> sieve
sieve length 30 * 1 - :> upper
3 upper sqrt 2 <range> [| i |
i sieve marked? [
i upper sieve unmark-multiples
] when
] each sieve ;```

As you can see, storing prime numbers up to 10 million takes only 333 KB!

```IN: scratchpad 10,000,000 sieve6 length .
333334```

This approach is used by the math.primes vocabulary to quickly check primality of any number less than 9 million using only 300 KB of memory.

Neat!