Seive3

seive3.c

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

#define BLOCK 64

#define SHIFT(X) (1ULL << (X - 1))
#define CHECK(X,Y) (X & SHIFT(Y))

void memset64(void * dest, uint64_t value, uintptr_t size);
uint64_t mask_shift(uint64_t m, uint64_t n);

int main(int argc, char** argv) {
  uint64_t *p, i, j, limit;

  limit = strtoull(argv[1], NULL, 0);

  p = (uint64_t*)calloc(limit, sizeof(uint64_t));

  fprintf(stderr, "\r\e[K%llu", 2ULL);

  memset64(p, 0xaaaaaaaaaaaaaaaa, limit);

  for(i = 3; i < ceil(sqrt(limit * BLOCK)); ++i) {
    if(!CHECK(p[i / BLOCK], i % BLOCK)) {
      fprintf(stderr, "\r\e[K%llu", i);
      if(i < 32) {
        uint64_t *m;

        m = (uint64_t*)calloc(i, sizeof(uint64_t));

        for(j = i; j <= BLOCK; j += i) {
          *m |= SHIFT(j);
        }

        for(j = 1; j < i; ++j) {
          *(m + j) = mask_shift(*(m + j - 1), i);
        }

        for(j = 0; j < limit; ++j) {
          *(p + j) |= *(m + (j % i));
        }

        free(m);
      } else {
        for(j = i * 2; j < limit * BLOCK; j += i) {
          p[j / BLOCK] |= SHIFT(j % BLOCK);
        }
      }
    }
  }

  fprintf(stderr, "\r\e[K");

  free(p);

  return 0;
}

void memset64(void * dest, uint64_t value, uintptr_t size) {
  uintptr_t i;
  for(i = 0; i < (size & (~7)); i += 8) {
    memcpy(((char*)dest) + i, &value, 8);
  }
  for(; i < size; ++i) {
    ((char*)dest)[i] = ((char*)&value)[i&7];
  }
}

uint64_t mask_shift(uint64_t m, uint64_t n) {
  uint64_t i, j, s, t;

  for(i = 0; i < 64; ++i) if(CHECK(m, 64 - i)) break;
  for(j = 0; j < 64; ++j) if(CHECK(m, j + 1)) break;

  if(i != 0) {
    s = i - (n - j - 1);
    t = m;
    m >>= s;
    m |= t << (n - s);
  }

  return m;
}