Category : Recently Uploaded Files
Archive : PSIEVE.ZIP
Filename : PSIEVE.C
Output of file : PSIEVE.C contained in archive : PSIEVE.ZIP
* psieve.c
*
* A sieve program using bitmaps. Based on a program written back
* in April of '88.
*
* Copyright 1991 by John M. Gamble.
* email to: [email protected]
*
* 19 Nov 1991
* Psieve version 1.0 is on the air.
*/
#include
#include
#include
#include
#include
/*
* Make sure that MAPSTYLE has a value of 1, 2, or 3.
* Default value is 3.
*/
#if MAPSTYLE < 1 || MAPSTYLE > 3
#define MAPSTYLE 3
#endif
/*
* I've yet to find a machine/compiler combination for which
* the following do not hold true, but there's always a first time...
*/
typedef short int int16;
typedef long int int32;
typedef unsigned short int uint16;
typedef unsigned long int uint32;
char *strrev(char *str);
int map_alloc(uint32 highest);
int primetest(uint32 c);
int ultostr(char buf[], unsigned long int num, int base);
uint16 pfile(uint16 fcount, uint32 highest);
uint32 isqrt(uint32 n);
void dex_out(uint32 x, FILE *fp);
void map_print(uint16 mapcount, uint32 highest);
void sift_to(uint32 highest);
void sieve(uint32 pn, uint32 upto);
void space_out(int a, FILE *fp);
void txt_print(uint32 highest);
void unprime(uint32 c);
/*
* The global variable bpar holds our bitmaps, which
* is what we sift through when we toss our composite
* numbers out.
*/
char **bpar;
#define TRUE 1
#define FALSE 0
#define BIT(x) (1 << (x))
#if MAPSTYLE == 3
#define REDUCE(x) (x /= 3)
#else
#if MAPSTYLE == 2
#define REDUCE(x) (x >>= 1)
#else
#define REDUCE(x)
#endif
#endif
#define MAP_SEG(x) ((x) / mapbits)
#define MAP_BIT(x) ((x) % mapbits)
/*
* mapbits represents the number of integers (each represented by a bit)
* that can be represented in our bitmap.
* Our bitmap is SHRT_MAX in size.
*/
#ifdef DEBUG
#undef SHRT_MAX
#define SHRT_MAX 3
#endif
const uint16 mapsize = SHRT_MAX;
const uint32 mapbits = (((uint32) SHRT_MAX) << 3);
int flag_txtfile = TRUE;
int flag_verbose = FALSE;
int flag_index = FALSE;
int flag_time = FALSE;
int flag_diff = FALSE;
int flag_mem = FALSE;
int siftto = 0;
int base_n = 10;
int fieldwidth = 8;
int paperwidth = 80;
uint32 chapsize = 0x10000L; /* 2**16, the default sift limit.*/
/*
* psieve [-many options] [upto]
*/
main(int argc, char **argv)
{
char *subarg;
int mapcount;
uint32 highest = 0xffffL;/* 2**16 - 1: Default value to sift up to.*/
while (--argc && (*(subarg = *++argv) == '-'))
{
switch (*++subarg)
{
case 'b': flag_txtfile = FALSE; break;
case 'c': chapsize = atol(++subarg); break;
case 'd': flag_diff = TRUE; break;
case 'f': fieldwidth = atoi(++subarg); break;
case 'i': flag_index = TRUE; break;
case 'm': flag_mem = flag_verbose = TRUE; break;
case 's': siftto = atoi(++subarg); break;
case 't': flag_time = TRUE; break;
case 'v': flag_verbose = TRUE;
fputs("Psieve version 1.0.\nCompiled with MAPSTYLE = ", stdout);
fputc(MAPSTYLE + '0', stdout);
fputc('\n', stdout); break;
case 'w': paperwidth = atoi(++subarg); break;
case 'x': base_n = 16; break;
default: exit(0);
}
}
if (argc)
highest = atol(*argv);
if (flag_verbose)
{
fputs("Highest = ", stdout);
dex_out(highest, stdout);
}
mapcount = map_alloc(highest);
sift_to(highest);
if (flag_txtfile)
txt_print(highest);
else
map_print(mapcount, highest);
if (flag_verbose)
fputs("psieve: finished writing file(s).\n", stdout);
}
/*
* int map_alloc(uint32 highest)
*
*/
int map_alloc(uint32 highest)
{
char buf[16];
int idx;
int mapcount;
REDUCE(highest);
mapcount = MAP_SEG(highest) + (MAP_BIT(highest) != 0);
if (flag_verbose)
{
fputs("Mapbits = ", stdout); dex_out(mapbits, stdout);
fputs("Mapsize = ", stdout); dex_out(mapsize, stdout);
fputs("Mapcount = ", stdout); dex_out(mapcount, stdout);
}
if ((bpar = (char **)calloc(mapcount, sizeof(char *))) == NULL)
{
fputs("Out of memory allocating array of pointers\n", stderr);
exit(0);
}
for (idx = 0; idx < mapcount; idx++)
{
if ((*(bpar + idx) = (char *) calloc(mapsize, sizeof(char))) == NULL)
{
fputs("Null pointer at ", stderr);
ultostr(buf, idx + 1, 10);
fputs(buf, stderr);
fputs(" calloc call.\n", stderr);
exit(0);
}
}
if (flag_mem) exit(0);
return (mapcount);
}
/*
* void sift_to(uint32 highest)
*
* Search for prime numbers. Since multiples of 2 are not stored in
* our bitmaps, the sifting starts with the number 3.
*
* Our increment is set to toggle between 2 and 4, so as to
* automatically skip multiples of 2 and 3. Therefore, one special
* sieving is carried out to handle the prime number 3, before we
* enter the loop.
*
* We could set up an array of increments to skip the next prime
* numbers, 5, 7, 11, and so forth, but i decided the speed gain was
* not worth the code complexity and space lost.
*/
void sift_to(uint32 highest)
{
char buf[16];
time_t secs;
register int incr;
uint32 lj;
uint32 hiroot;
hiroot = isqrt(highest);
if (siftto > 1 && siftto < hiroot)
hiroot = siftto;
fputs("Beginning sift.\n", stdout);
if (flag_time) secs = time(NULL);
unprime(1L); /* One is not a prime number.*/
/*
* If our maps need it, sift out 2 and/or 3 specially.
* The sieve loop will be skipping multiples of 2 and 3.
*/
#if MAPSTYLE == 1
sieve(2L, highest);
#endif
#if MAPSTYLE != 3
sieve(3L, highest);
#endif
incr = 4;
for (lj = 5; lj <= hiroot; lj += (incr ^= 6))
sieve(lj, highest);
if (flag_time)
secs = time(NULL) - secs;
fputs("End of sift.\n", stdout);
if (flag_time)
{
fputs("Elapsed time in seconds: ", stdout);
ultostr(buf, secs, 10);
strcat(buf, "\n");
fputs(buf, stdout);
}
}
/*
* void sieve(uint32 pn, uint32 upto)
*
* Sift out the potentially prime pn from the array of bitmaps.
*/
void sieve(uint32 pn, uint32 upto)
{
char buf[16];
uint32 v;
#if MAPSTYLE == 3
int toggle3;
#endif
/*
* If the number in question is indeed prime, sift it out.
* Depending upon the map style, even-numbered multiples of
* pn and/or multiples of three will not be sifted, since
* we may not storing these numbers in the bitmap. Our
* starting point is the number squared.
*/
if (primetest(pn) == 0)
{
if (flag_verbose)
{
ultostr(buf, pn, 10);
strcat(buf, "\n");
fputs(buf, stdout);
}
#if MAPSTYLE == 3
toggle3 = (pn % 3 != 1); /* Toggle to skip multiples of 3.*/
#endif
v = pn * pn;
#if MAPSTYLE > 1
pn <<= 1; /* Skip even-numbered multiples.*/
#endif
upto -= pn;
for (;;)
{
unprime(v);
if (v > upto) break;
v += pn;
#if MAPSTYLE == 3
if (toggle3 ^= 1)
{
if (v > upto) break;
v += pn;
}
#endif
}
}
}
/*
* map_print(uint16 mapcount, uint32 highest)
*/
void map_print(uint16 mapcount, uint32 highest)
{
const char *fname = "primes.map";
uint16 mcount;
uint32 map_off;
int nthbyte;
FILE *fp;
/*
* Store the bit maps.
*/
if ((fp = fopen(fname, "wb")) == NULL)
{
fputs("Cannot open file ", stdout);
fputs(fname, stdout);
}
else
{
fputs("Writing to file ", stdout);
fputs(fname, stdout);
}
fputc('\n', stdout);
for (mcount = 0; mcount < mapcount - 1; mcount++)
fwrite(*(bpar + mcount), sizeof(char), mapsize, fp);
/*
* Write out the last map, which may be only partially filled.
*/
REDUCE(highest);
map_off = MAP_BIT(highest);
nthbyte = map_off >> 3;
if (map_off & 7) nthbyte++;
fwrite(*(bpar + mcount), sizeof(char), nthbyte, fp);
fclose(fp);
}
/*
* void txt_print(uint32 highest)
*/
void txt_print(uint32 highest)
{
const char *pr_index = "primes.idx";
char buf[16];
FILE *fp = NULL;
int len;
uint16 fcount;
uint16 prime_count = 0;
uint32 prime_countsum = 0L;
uint32 chapters;
chapters = (highest/chapsize) + ((highest % chapsize) != 0);
/*
* Print out the prime numbers, chapter by chapter.
*/
if (flag_index)
if ((fp = fopen(pr_index, "w")) == NULL)
fputs("Index file cannot be opened.\n", stderr);
for (fcount = 0; fcount < chapters; fcount++)
{
prime_countsum += prime_count = pfile(fcount, highest);
if (fp != NULL)
{
len = ultostr(buf, fcount, 10);
strcat(buf, ":");
space_out(8 - len, fp);
fputs(buf, fp);
len = ultostr(buf, prime_count, 10);
strcat(buf, ":");
space_out(8 - len, fp);
len = ultostr(buf, prime_countsum, 10);
space_out(11 - len, fp);
fputs(buf, fp);
fputc('\n', fp);
}
}
if (fp != NULL)
fclose(fp);
}
/*
* uint16 pfile(uint16 fcount, uint32 highest)
*
* Print out the primes stored in the chapter.
*/
uint16 pfile(uint16 fcount, uint32 highest)
{
static uint16 indent = 0;
static uint32 last_prime = 3;
FILE *fp;
uint32 ival;
uint32 upto;
uint16 prime_count;
int len, ind;
int incr;
char buf[16];
char fname[16];
strcpy(fname, "00000000");
if (fcount > 0)
{
len = ultostr(buf, fcount, 16);
fname[8 - len] = '\0';
strcat(fname, buf);
}
strcat(fname, flag_diff? ".pdf": ".ptx");
if ((fp = fopen(fname, "w")) == NULL)
{
fputs("Cannot open file ", stderr);
fputs(fname, stderr);
fputc('\n', stderr);
return(0);
}
else
{
fputs("Writing file ", stdout);
fputs(fname, stdout);
fputc('\n', stdout);
}
/*
* Special case: in the first chapter, take care of 2, which is not
* stored, and print out 3, so that we can set up our increment to
* skip multiples of 2 and 3.
*
* For the other cases, ival is the value at the start of the chapter,
* adjusted for the increment.
*/
if (fcount == 0)
{
ind = fieldwidth;
space_out(fieldwidth - 1, fp);
if (flag_diff)
fputc('1', fp);
else
{
fputc('2', fp);
ind += fieldwidth;
space_out(fieldwidth - 1, fp);
fputc('3', fp);
}
prime_count = 2;
ival = 5L;
incr = 4;
upto = chapsize;
}
else
{
ind = indent * fieldwidth;
space_out(ind, fp);
prime_count = 0;
ival = 1 + chapsize * fcount;
incr = (ival % 3 == 1)? 2: 4;
if (ival % 3 == 0) ival += 2;
upto = chapsize * (fcount + 1);
}
/*
* Adjust upto so that the loop will quit when
* ival reaches its value.
*/
if (upto > highest) upto = highest;
if ((upto & 1) == 0) upto--;
if ((upto % 3) == 0) upto -= 2;
/*
* Go through the bitmap, looking for unset bits.
*/
for (;;)
{
if (primetest(ival) == 0)
{
ind += fieldwidth;
if (ind > paperwidth)
{
ind = fieldwidth;
fputc('\n', fp);
}
len = ultostr(buf, flag_diff? ival - last_prime: ival,
base_n);
space_out(fieldwidth - len, fp);
fputs(buf, fp);
prime_count++;
last_prime = ival;
}
if (ival >= upto) break;
ival += (incr ^= 6);
}
if (ind < paperwidth)
fputc('\n', fp);
fclose(fp);
indent = ind/fieldwidth;
return(prime_count);
}
/*
* int primetest(uint32 c)
*
* Reduce the number c to its bitmap location, and return the bit
* value of zero or one.
*/
int primetest(uint32 c)
{
register char *bptr;
#ifdef DEBUG
if (flag_verbose)
fprintf(stdout, "Testing %lu\n", c);
#endif
REDUCE(c);
#ifdef DEBUG
if (flag_verbose)
fprintf(stdout, "MAP_SEG(c) = %lu, MAP_BIT(c) = %lu\n",
MAP_SEG(c), MAP_BIT(c));
#endif
bptr = *(bpar + MAP_SEG(c));
c = MAP_BIT(c);
bptr += (uint16)(c >> 3);
return ((*bptr >> ((uint16)c & 7)) & 1);
}
/*
* void unprime(uint32 c)
*
* Reduce the number c to its bitmap location, and set the bit value
* to one.
*/
void unprime(uint32 c)
{
register char *bptr;
#ifdef DEBUG
if (flag_verbose)
fprintf(stdout, "\t%lu sifted\n", c);
#endif
REDUCE(c);
bptr = *(bpar + MAP_SEG(c));
c = MAP_BIT(c);
bptr += (uint16)(c >> 3);
*bptr |= BIT((uint16)c & 7);
}
/*
* uint32 isqrt(uint32 n)
*
* Take the integer square root. Returns an integer value less than
* or equal to the square root of n.
* From Richard Campbell, March 1986 DDJ.
*/
uint32 isqrt(uint32 n)
{
register uint32 guess1, guess2;
if (n < 2)
return (n);
guess1 = 1;
guess2 = n;
/*
* Find reasonable starting guesses.
*/
while ((guess1 << 1) < guess2)
{
guess1 <<= 1;
guess2 >>= 1;
}
/*
* Take our starting guesses, and find their mean, which
* we store into guess1. This will be larger than the
* root (if it isn't already the root). Divide it into
* n, and recalculate until guess1 dips below or equals
* guess2. This means we have closed in on the root,
* which we return.
*/
do
{
guess1 += guess2;
guess1 >>= 1;
guess2 = n/guess1;
} while (guess1 > guess2);
return (guess1);
}
/*
* int ultostr(char buf[], unsigned long int num, int base)
*
* Convert an unsigned long integer into a string using base
* as the radix. Returns the length of the string in the buffer.
*/
int ultostr(char buf[], unsigned long int num, int base)
{
unsigned int len, m;
len = 0;
if (base > 1)
{
do
{
m = num % base;
if (m > 9)
buf[len] = (char) (m - 10) + 'a';
else
buf[len] = (char) m + '0';
len++;
} while ((num /= base) > 0);
}
buf[len] = '\0';
strrev(buf);
return (len);
} /* end of ultostr.*/
/*
* void dex_out(uint32 x, FILE *fp)
*
* Put out a number in "%ld (0x%lx)" format.
*/
void dex_out(uint32 x, FILE *fp)
{
char buf[32];
int len;
len = ultostr(buf, x, 10);
strcat(buf, " (0x");
ultostr(&buf[len + 4], x, 16);
strcat(buf, ")\n");
fputs(buf, fp);
}
/*
* void space_out(int a, FILE *fp)
*
*/
void space_out(int a, FILE *fp)
{
register int j;
for (j = 0; j < a; j++)
fputc(' ', fp);
}
#if __ZTC__ == 0 && __TURBOC__ == 0
/*
* char *strrev(char *str)
*
* Reverse the contents of the string in str.
*/
char *strrev(char *str)
{
int s;
char *ptr1, *ptr2;
char c;
s = strlen(ptr1 = str);
ptr2 = ptr1 + (s - 1);
s >>= 1;
while (s--)
{
c = *ptr1;
*ptr1++ = *ptr2;
*ptr2-- = c;
}
return (str);
}
#endif
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/