#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <getopt.h>
#include <sys/mman.h>

/*---------------------------------------------------------------------------*/

#ifdef MAP_ANON
# define MAP (MAP_PRIVATE | MAP_ANON)
#else
# define MAP (MAP_PRIVATE)
#endif

static int fd = -1;

/*---------------------------------------------------------------------------*/

typedef unsigned long   dword;
typedef unsigned short  word;
typedef unsigned char   byte;

/*---------------------------------------------------------------------------*/

static dword  start;                                    /* ab start ausgeben */
static dword  ziel;                              /* von 0 bis ziel berechnen */
static dword  mem;                          /* Benoetigter Speicher in Bytes */
static byte *bitmap;                       /* Bitmap, fuer jede Zahl ein Bit */

static dword est,est_i;
static dword loope;
static dword *loopt;
static dword *loopc;

/*---------------------------------------------------------------------------*/

dword
estimate_max_prim()
{
    unsigned char *data = NULL;
    int  high = 1023;
    int  size;
    int  low  = 0;
    
#ifndef MAP_ANON
    if (-1 == (fd = open("/dev/zero",O_RDONLY))) {
	perror("open /dev/zero");
	exit(1);
    }
#endif

    /* check how much we can map */
    for (;;) {
	if (high - low < 2)
	    break;
	size = (high+low)/2;
	data = mmap(NULL, size*1024*1024,
		    PROT_READ | PROT_WRITE, MAP, fd, 0);
	if ((void*)-1 == data) {
	    high = size;
	} else {
	    munmap(data,size*1024*1024);
	    low = size;
	}
	fflush(stdout);
    }
    low = low*3/4; /* use max. 75% of mem */

    return (low*1024*1024-1)*16;
}

/*---------------------------------------------------------------------------*/

static void inline
clearzahl(dword zahl)
{
    if (zahl & 1) {
	zahl >>= 1;
	bitmap[zahl>>3] |= 1 << (zahl&7);
    }
}

static int inline
isprim(dword zahl)
{
    if (zahl & 1) {
	zahl >>= 1;
	return !(bitmap[zahl>>3] & (1 << (zahl&7)));
    } else
	return zahl==2;
}

/*---------------------------------------------------------------------------*/

void calc_init(int max)
{
    dword  count,testing,i;

    for (i = 0, testing = 3; testing*testing <= ziel; testing++) {
	if (!isprim(testing))
	    continue;
	for (count = 3*testing; count <= max; count += testing, count += testing)
	    clearzahl(count);
	loopt[i] = testing;
	loopc[i] = count;
	i++;
    }
    loope = i;
    return;
}

void calc_step(int max)
{
    register dword  count,testing;
    dword           i;

    for (i = 0; i < loope; i++) {
	testing = loopt[i]<<1;
	for (count = loopc[i]; count <= max; count += testing)
	    clearzahl(count);
	loopc[i] = count;
    }
    return;
}

dword
calculate1(dword max)
{
    ziel = max;
    mem = ziel/16+1;
    bitmap = mmap(NULL, mem, PROT_READ | PROT_WRITE, MAP, fd, 0);
    if (-1 == (int)bitmap) {
	fprintf(stderr,"out of memory\n");
	exit(1);
    }
    clearzahl(0);
    clearzahl(1);

    for (est = 1000; est*est <= ziel; est += 1000);
    est *= 4; /* tuneable */
    if (est > ziel)
	est = ziel;
    loopt = malloc((est+1)*sizeof(dword));
    loopc = malloc((est+1)*sizeof(dword));
    if (NULL == loopt || NULL == loopc) {
	fprintf(stderr,"out of memory\n");
	exit(1);
    }

    est_i = est;
    calc_init(est_i);
    return est_i;
}

dword
calculate2()
{
    est_i += est;
    if (est_i > ziel)
	est_i = ziel;
    calc_step(est_i);
    return est_i;
}

/*---------------------------------------------------------------------------*/

void usage(char *name)
{
    fprintf(stderr,
	    "\n"
	    "I'll compute prime numbers up to max\n"
	    "usage:\n"
	    "  %s [ -d | -q ] [ -s min ] max\n"
	    "\n"
	    "options:\n"
	    "  -d       print only twins (two primes, like 17 and 19)\n"
	    "  -q       print only groups of four primes (191,193,197,199 for example)\n"
	    "  -s min   don't print primes smaller than min\n"
	    "\n"
	    "-- \n"
	    "(c) 1998 Gerd Knorr <kraxel@goldbach.in-berlin.de>\n",
	    name);
    exit(1);
}

int
main(int argc, char *argv[])
{
    dword count,max;
    int   dbl = 0, qua = 0, c;

    max = estimate_max_prim();
    for (;;) {
	c = getopt(argc, argv, "dqs:");
	if (c == -1)
	    break;
	switch (c) {
	case 'd':
	    dbl = 1;
	    break;
	case 'q':
	    qua = 1;
	    break;
	case 's':
	    start = atoi(optarg);
	    break;
	default:
	    usage(argv[0]);
	}
    }

    if (optind+1 != argc)
	usage(argv[0]);
    ziel = atol(argv[optind]);
#if 0
    if (ziel > max) {
	fprintf(stderr,"%ld: too big\n",ziel);
	exit(1);
    }
#endif

    fprintf(stderr,"processing... ");
    calculate1(ziel);
    while (calculate2() < ziel);
    fprintf(stderr,"ok\n");

    if (dbl) {
	if (start < 6) {
	    printf("3 5\n");
	} /* endif */
	for (count = start+6-start%6; count < ziel; count+= 6) {
	    if (isprim(count-1) && isprim(count+1)) {
		printf("%lu %lu\n",count-1,count+1);
	    } /* endif */
	} /* endfor */
    } else if (qua) {
	if (start < 10) {
	    printf("5 7 11 13\n");
	} /* endif */
	for (count = start-start%10+5; count < ziel-3; count += 10) {
	    if (isprim(count-4) && isprim(count-2) && isprim(count+2) && isprim(count+4)) {
		printf("%lu %lu %lu %lu\n",count-4,count-2,count+2,count+4);
	    } /* endif */
	} /* endif */
    } else {
	for (count = start; count <= ziel; count++) {
	    if (isprim(count)) {
		printf("%lu ",count);
	    } /* endif */
	} /* endfor */
    } /* endif */
    printf("\n");
    munmap(bitmap,mem);
    return(0);
}