#include <stdio.h>

#include <string.h>

#include <stdlib.h>

#include <assert.h>



typedef unsigned long long bignum;



typedef struct

{

unsigned int *p; /* pointer to array */

int bitsPerByte; /* 8 on normal systems */

int bytesPerInt; /* sizeof(unsigned int) */

int bitsPerInt; /* for bit arithmetic */

int bpiShift; /* 8 bit words = 3, 16=4, etc */

int bpiMask; /* bitsPerInch - 1, for modulus */

bignum bitsInArray; /* how many bits in array */

bignum intsInArray; /* how many uints to give necessary bits */

bignum pageno; /* page 0 is 0-amillion or whatever */

bignum offset; /* ba start at 0, a million, 2 mil? */

} BITARRAY;



void freeBitArray(BITARRAY *ba)

{

free(ba->p);

free(ba);

}



BITARRAY * createBitArray(bignum bits, bignum pageno)

{

BITARRAY *ba = malloc(sizeof(BITARRAY));

assert(ba != NULL);

ba->bitsPerByte = 8;

ba->bytesPerInt = sizeof(unsigned int);

ba->bitsPerInt = ba->bitsPerByte * ba->bytesPerInt;

switch (ba->bitsPerInt)

{

case 8: ba->bpiShift = 3; break;

case 16: ba->bpiShift = 4; break;

case 32: ba->bpiShift = 5; break;

case 64: ba->bpiShift = 6; break;

case 128: ba->bpiShift = 7; break;

case 256: ba->bpiShift = 8; break;

case 512: ba->bpiShift = 9; break;

default: {

perror("ABORTING: Non-standard bits per int

");

exit(1);

break;

}

};

ba->bpiMask = ba->bitsPerInt - 1;

ba->bytesPerInt = sizeof(unsigned int);

ba->bitsInArray = bits;

ba->intsInArray = bits/ba->bitsPerInt + 1;

ba->p = malloc(ba->intsInArray * ba->bytesPerInt);

assert(ba->p != NULL);

ba->pageno = pageno;

ba->offset = ba->pageno * ba->bitsInArray;

return ba;

}



void setBit(BITARRAY *ba, bignum bitSS)

{

unsigned int *pInt = ba->p + (bitSS >> ba->bpiShift);

unsigned int remainder = (bitSS & ba->bpiMask);

*pInt |= (1 << remainder);

}



void clearBit(BITARRAY *ba, bignum bitSS)

{

unsigned int *pInt = ba->p + (bitSS >> ba->bpiShift);

unsigned int remainder = (bitSS & ba->bpiMask);

unsigned int mask = 1 << remainder;

mask = ~mask;

*pInt &= mask;

}



int getBit(BITARRAY *ba, bignum bitSS)

{

unsigned int *pInt = ba->p + (bitSS >> ba->bpiShift);

unsigned int remainder = (bitSS & ba->bpiMask);

unsigned int ret = *pInt;

ret &= (1 << remainder);

return(ret != 0);

}





void clearAll(BITARRAY *ba)

{

bignum intSS;

for(intSS=0; intSS <= ba->intsInArray; intSS++)

{

*(ba->p + intSS) = 0;

}

}



void setAll(BITARRAY *ba)

{

bignum intSS;

for(intSS=0; intSS <= ba->intsInArray; intSS++)

{

*(ba->p + intSS) = ~0;

}

}



void printPrime(bignum bn)

{

printf("%llu

", bn);

}



void fprintPrime(FILE *fp, bignum bn)

{

fprintf(fp, "%llu

", bn);

}



bignum * makePrimeArray(const char *fname)

{

FILE * fin;

char buf[1000], *pbuf;

bignum * factors;

fin = fopen(fname, "r");

assert(fin != NULL);



/* FIND NUMBER OF LINES */

bignum lineCount = 0;

while(fgets(buf, 90, fin) != NULL)

lineCount++;



/* ALLOCATE ARRAY */

factors = malloc(sizeof(bignum) * (lineCount + 3));

assert(factors != NULL);



/* READ THE PRIMES INTO factors */

rewind(fin);

lineCount = 0;

while((pbuf = fgets(buf, 90, fin)) != NULL)

{

*(factors + lineCount) = atoll(pbuf);

lineCount++;

}

fclose(fin);

*(factors + lineCount) = 0;

return factors;

}



void findPrimesInOnePage(bignum topCandidate, bignum pageno, bignum *factors, BITARRAY *ba)

{

bignum ss;

char buf[100];

FILE *fout;



/* SET ba ELEMENTS */

ba->pageno = pageno;

ba->offset = ba->pageno * ba->bitsInArray;



/* SET ALL BUT 0 AND 1 TO PRIME STATUS */

setAll(ba);



/* MARK ALL THE NON-PRIMES */

bignum factorSS = 0;

bignum thisFactor = *(factors + factorSS);

while((thisFactor != 0) && (thisFactor * thisFactor <= topCandidate + ba->offset))

{

/* MARK THE MULTIPLES OF THIS FACTOR */

bignum mark = thisFactor + thisFactor;

if(mark < ba->offset)

{

mark = (ba->offset/thisFactor) * thisFactor;

if(mark < ba->offset)

mark += thisFactor;

}

mark -= ba->offset;

while(mark <= topCandidate)

{

clearBit(ba, mark);

mark += thisFactor;

}



/* SET thisFactor TO NEXT PRIME */

factorSS++;

thisFactor = *(factors + factorSS);

assert(thisFactor <= topCandidate + ba->offset);

}



/* PRINT ALL THE PRIMES */

sprintf(buf, "pri%08llu.pri", ba->pageno);

fout = fopen(buf, "w");

assert(fout != NULL);



for(ss = 0;ss <= topCandidate; ss++)

{

if(getBit(ba, ss)) fprintPrime(fout, ss + ba->offset);

}

fclose(fout);

}



void findPrimes(bignum numbersPerPage, bignum startPageno, bignum pagesToProcess)

{

bignum pageCount;



/* PUT FACTORS IN AN ARRAY */

bignum * factors = makePrimeArray("0.pri");



/* CREATE BITARRAY */

BITARRAY *ba = createBitArray(numbersPerPage, startPageno);

assert(ba != NULL);



/* PROCESS EVERY PAGE */

for(pageCount = 0; pageCount < pagesToProcess; pageCount++)

{

findPrimesInOnePage

(

numbersPerPage,

pageCount+startPageno,

factors,

ba

);

}



/* FREE THE BITARRAY */

freeBitArray(ba);

/* FREE THE FACTOR ARRAY */

}





void test()

{

bignum * factors = makePrimeArray("0.pri");

bignum lineCount = 0;

bignum factor;

FILE * fout = fopen("junky.jnk", "w");

assert(fout != NULL);

for(lineCount = 0; (factor = *(factors + lineCount)) != 0; lineCount++)

fprintPrime(fout, factor);

fprintPrime(fout, 0);

fclose(fout);

}



int main(int argc, char *argv[])

{

bignum numbersPerPage = atoll(argv[1]);

bignum startPageno = atoll(argv[2]);

bignum pagesToProcess = atoll(argv[3]);



printf("numbersPerPage = %llu

", numbersPerPage );

printf("startPageno = %llu

", startPageno );

printf("pagesToProcess = %llu

", pagesToProcess);

/* test(); */



findPrimes(numbersPerPage, startPageno, pagesToProcess) ;

return 0;

}

