Black Key sieve algorithm implementation

blackkeysieve.c

/*
 * Black Key sieve implementation
 *
 * This is a C99 implementation of the Black Key Sieve algorithm,
 * as described by Larry Deering in [1]. It relies on the fact that all
 * primes above 5 will be equal to {1, 7, 11, 13, 17, 19, 23, 29} modulo
 * 30, allowing us to care only about the numbers satisfying this property,
 * and to pack the information on 30 integers in a single byte.
 *
 * This code is released in the public domain.
 *
 * [1] http://www.qsl.net/w2gl/blackkey.html
 */

#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>

#define MAX_ENTRY 10000000
#define SIEVE_SIZE (MAX_ENTRY / 3)
#define PRIMES_BOUND (MAX_ENTRY / 10)

#define PACK(n, off, bn) do {off = n / 30; bn = invlookup(n % 30);} while (0)
#define INCR(off, bn) if (bn == 7) {off++; bn = 0;} else {bn++;}
#define UNPACK(off, bn) (30 * off + lookup[bn])
#define MARK(res, off, bn) (res[off] |= (1 << bn))
#define MARKED(res, off, bn) (res[off] & (1 << bn))

static const uint8_t lookup[] = {1, 7, 11, 13, 17, 19, 23, 29};
static inline int invlookup(int n)
{
  switch (n) {
    case 1:  return  0;
    case 7:  return  1;
    case 11: return  2;
    case 13: return  3;
    case 17: return  4;
    case 19: return  5;
    case 23: return  6;
    case 29: return  7;
    default: return -1;
  }
}

static int primes[PRIMES_BOUND];

size_t sieve(size_t max)
{
  assert(max <= MAX_ENTRY);

  size_t nprimes = 0;

  primes[nprimes++] = 2;
  primes[nprimes++] = 3;
  primes[nprimes++] = 5;

  static uint8_t res[SIEVE_SIZE]; // Initially null.

  for (size_t off = 0, bn = 1; UNPACK(off, bn) <= max;) {
    if (!MARKED(res, off, bn)) {

      int prime = UNPACK(off,bn);

      primes[nprimes++] = prime;

      size_t i, o, b;
      for (i = 0, o = off, b = bn; i < 8; i++) {
        int p = UNPACK(o, b);
        p *= prime;
        size_t o1, b1;
        PACK(p, o1, b1);
        for (; UNPACK(o1,b1) <= max; o1 += prime) {
          MARK(res, o1, b1);
        }

        INCR(o, b);
      }
    }

    INCR(off, bn);
  }

  return nprimes;
}