//
// Terje Mathisen wrote this beautiful sieve algorithm in Feb, 1998
//
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

unsigned char  *buf;

unsigned char   mask_tab[30] = {
    0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4, 0, 8, 0,
    0, 0, 16, 0, 32, 0, 0, 0, 64, 0, 0, 0, 0, 0, 128
};

unsigned        max_prime_number;

int      test_dm(unsigned d, unsigned m)
{
    return mask_tab[m] && ((buf[d] & mask_tab[m]) == 0);
}

void     set_dm(unsigned d, unsigned m)
{
    buf[d] |= mask_tab[m];
}

unsigned divmod30(unsigned index, unsigned *pm)
{
//  d = index / 30;
    unsigned        d;
    __int64         t = (__int64) index * 2290649225;
    d = *(((unsigned *) &t) + 1);
    d >>= 4;
    *pm = index - d * 30;
    return d;
}

unsigned next_prime(unsigned p)
{
    unsigned        d,
                    m;

    if (p & 1)
        p++;
    p++;
    d = divmod30(p, &m);
    while (m < 30) {
        if (p > max_prime_number)
            return 0;
        if (test_dm(d, m))
            return p;
        m += 2;
        p += 2;
    }
    m = 1;
    d++;
    while (buf[d] == 0xff) {
        p += 30;
        d++;
        if (p > max_prime_number)
            return 0;
    }
    do {
        if (test_dm(d, m))
            return p;
        m += 2;
        p += 2;
    } while (p <= max_prime_number);
    return 0;
}

unsigned        sieve(unsigned alloc, unsigned max)
{
    unsigned        maxtest = (unsigned) sqrt((double) max);
    unsigned        curr,
                    step,
                    dstep[30],
                    next_m[30],
                    dcurr,
                    mcurr;
    unsigned        prime,
                    i;
    unsigned        count;
    unsigned        zero_bits[256];

    max_prime_number = max;

    prime = 5;
    for (;;) {
        prime = next_prime(prime);
        if (prime > maxtest)
            break;

        step = prime * 2;
        curr = prime * prime;
        dcurr = divmod30(curr, &mcurr);

        for (i = 1; i < 30; i += 2) {
            unsigned        d,
                            m,
                            s = i;
            do {
                s += step;
                d = divmod30(s, &m);
            } while (mask_tab[m] == 0);
            dstep[i] = d;
            next_m[i] = m;
        }

        do {
            if ((buf[dcurr] & mask_tab[mcurr]) == 0)
                set_dm(dcurr, mcurr);

            dcurr += dstep[mcurr];
            mcurr = next_m[mcurr];
        } while (dcurr < alloc);
    }
    // Count the primes!
    printf("...Counting!\n");


    for (curr = 0; curr < 256; curr++) {
        unsigned        c = curr;
        for (count = 8; c; count--) {
            c &= c - 1;
        }
        zero_bits[curr] = count;
    }

    for (count = 2, curr = 0; curr < (max / 30); curr++) {
        count += zero_bits[buf[curr]];
    }
    for (prime = curr * 30; prime < max;) {
        prime = next_prime(prime);
        if (!prime)
            break;
        count++;
    }

    return count;
}

int             main(int argc, char *argv[])
{
    unsigned        count,
                    alloc,
                    s,
                    max,
                    show;
    time_t          begin;
    if (argc > 1)
        max = (unsigned) atof(argv[1]);
    else
        max = 1000L;

    if (argc > 2)
        show = (unsigned) atol(argv[2]);
    else
        show = 100;

    alloc = (max + 29) / 30 + 4;
    buf = (unsigned char *) calloc(alloc, 1);
    if (buf == NULL) {
        printf("Cannot handle this many numbers!\n");
        return EXIT_FAILURE;
    }
    printf("Searching prime numbers to : %ld", max);
    begin = time(NULL);

    count = sieve(alloc, max);

    printf(" %ld prime numbers found in %ld secs\n",
           count, time(NULL) - begin);

    s = 1;
    for (;;) {
        if (s >= 5) {
            s = next_prime(s);
            if (!s || (s >= show))
                break;
        } else if (s >= 3)
            s = 5;
        else if (s >= 2)
            s = 3;
        else
            s = 2;
        printf("%8d", s);
    }
    return 0;
}