/*
* erato: sieve of eratosthenes for up to 2**32.
* Two steps: first find all primes less than 2**16, then use them to sift
* the remaining numbers in blocks.
*
* Jim Gillogly, 24 Apr 1999
*
* Hacked by Peter Gunn 26 Apr 1999
*
* NOTE: Switch off Gloabal Optimsation & Full Optimisation when using
*       MS VC++ 5, coz its broken!!!
*/

#include <stdio.h>
#include <memory.h>
#include <time.h>

#define MAXSMALL        6542    /* 6542 for primes up to 1<<16. */
#define TOP             (1 << 16)       /* Top of each block to be sifted. */
#define PRFILE          "primes"/* Save 'em in a file. */
#define BLOCKS          (1 << 16)       /* How many blocks to do total? */

unsigned long   sieve[TOP / 32];/* Turned off if composite. */
#define bit(x) (sieve[(x)>>5]&(1<<(x&31)))
#define cbit(x) sieve[(x)>>5]&=~(1<<(x&31))

/* bitmasks for all primes < 32 in all rotations */
unsigned long   sieve2[32][32][TOP / 32];
#define csbit(y,r,x) sieve2[y][r][(x)>>5]&=~(1<<(x&31))

//char *  stop = &sieve[TOP];

unsigned long   smalls[MAXSMALL];       /* All primes under 2**16. */
unsigned long   smallnum = 0;

main()
{
    unsigned long   b,
                    b1,
                    i,
                    j,
                    k,
                    nprimes,
                    start,
                    step,
                   *s;
    char            xxx[40];
    unsigned long   start_time,
                    stop_time;

    printf("hit enter...\n");
    gets(xxx);
    start_time = time(NULL);
    memset(sieve, 0xff, sizeof(sieve));
    cbit(0);
    cbit(1);
    b = 0;
    do {
        b++;
        while (b < TOP && !bit(b))
            b++;                /* Scan to next prime. */
        b1 = b;
        while ((b1 += b) < TOP)
            cbit(b1);
    } while (b < TOP);
    for (i = 0, s = smalls; i < TOP; i++) {
        if (bit(i)) {
            *s = i;

            if (i < 32) {
                unsigned long   rot;
                for (rot = 0; rot < 32; rot++) {
                    memset(&sieve2[s - smalls][rot], 0xff, sizeof(sieve2[0][0]));
                    b1 = rot % i;
                    csbit(s - smalls, rot, b1);
                    while ((b1 += i) < TOP)
                        csbit(s - smalls, rot, b1);
                }
            }
            s++;
        }
    }
    smallnum = nprimes = (s - smalls);
    printf("%d small primes.\n", smallnum);
    /* Now sift each 2**16-number block with these small primes. */
    start = 0;
    for (k = 1; k < BLOCKS; k++) {
        start = k * TOP;
        memset(sieve, 0xff, sizeof(sieve));
        cbit(0);
        for (i = 0, s = smalls; i < smallnum; i++) {    /* small primes */
            if ((step = *s++) < 32) {
                unsigned long   rot = step - (start % step);
                unsigned long  *s1 = &sieve[TOP / 32 - 1],
                               *s2 = &sieve2[s - smalls - 1][rot][TOP / 32 - 1];
                while (s1 >= sieve)
                    *s1-- &= *s2--;
            } else {
                for (b = step - (start % step); b < TOP; b += step)
                    cbit(b);
            }
        }
        for (i = 0; i < TOP; i++) {
            if (bit(i)) {
                j = i + start;
                nprimes++;
                if (!(nprimes & 0xfffff))
                    printf("%ld primes in %ld secs (k=%ld)\n", nprimes,
                           time(NULL) - start_time, k);
            }
        }
    }
    stop_time = time(NULL);
    printf("%d primes in %ld secs \n", nprimes, stop_time - start_time);
    printf("hit enter\n");
    gets(xxx);
    return 0;
}