C Source Code for a Sieve Program

/* Copyright 1999-2000 (C) John Moyer */
/* assume 32 bit ints */
/* mailto:jrm@rsok.com */
#include <stdlib.h>
#include <stddef.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#define STEP_CONST 128

int small_primes[343] = {
2,3,5,7,11,13,17,19,23,29,
31,37,41,43,47,53,59,61,67,71,
73,79,83,89,97,101,103,107,109,113,
127,131,137,139,149,151,157,163,167,173,
179,181,191,193,197,199,211,223,227,229,
233,239,241,251,257,263,269,271,277,281,
283,293,307,311,313,317,331,337,347,349,
353,359,367,373,379,383,389,397,401,409,
419,421,431,433,439,443,449,457,461,463,
467,479,487,491,499,503,509,521,523,541,
547,557,563,569,571,577,587,593,599,601,
607,613,617,619,631,641,643,647,653,659,
661,673,677,683,691,701,709,719,727,733,
739,743,751,757,761,769,773,787,797,809,
811,821,823,827,829,839,853,857,859,863,
877,881,883,887,907,911,919,929,937,941,
947,953,967,971,977,983,991,997,1009,1013,
1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,
1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,
1153,1163,1171,1181,1187,1193,1201,1213,1217,1223,
1229,1231,1237,1249,1259,1277,1279,1283,1289,1291,
1297,1301,1303,1307,1319,1321,1327,1361,1367,1373,
1381,1399,1409,1423,1427,1429,1433,1439,1447,1451,
1453,1459,1471,1481,1483,1487,1489,1493,1499,1511,
1523,1531,1543,1549,1553,1559,1567,1571,1579,1583,
1597,1601,1607,1609,1613,1619,1621,1627,1637,1657,
1663,1667,1669,1693,1697,1699,1709,1721,1723,1733,
1741,1747,1753,1759,1777,1783,1787,1789,1801,1811,
1823,1831,1847,1861,1867,1871,1873,1877,1879,1889,
1901,1907,1913,1931,1933,1949,1951,1973,1979,1987,
1993,1997,1999,2003,2011,2017,2027,2029,2039,2053,
2063,2069,2081,2083,2087,2089,2099,2111,2113,2129,
2131,2137,2141,2143,2153,2161,2179,2203,2207,2213,
2221,2237,2239,2243,2251,2267,2269,2273,2281,2287,
2293,2297,2309
};

/* expand scheme from 30 to 2310
   30 == 2*3*5 and 8 == (2-1)*(3-1)*(5-1)
   2*3*5*7*11 == 2310
   1*2*4*6*10 ==  480 bits == 60 bytes ==> prime or not prime for 
   38.5 integers per byte

In each 2310 integers, for N >= 1, the numbers that might be prime are:
N*2310+1,
N*2310+13,
N*2310+17,
.
.
.
N*2310+13*13,
N*2310+13*17,
.
.
.
N*2310+2309
*/


unsigned long maskbits[32] = 
{
0x00000001, 0x00000002, 0x00000004, 0x00000008, 
0x00000010, 0x00000020, 0x00000040, 0x00000080,
0x00000100, 0x00000200, 0x00000400, 0x00000800, 
0x00001000, 0x00002000, 0x00004000, 0x00008000,
0x00010000, 0x00020000, 0x00040000, 0x00080000, 
0x00100000, 0x00200000, 0x00400000, 0x00800000,
0x01000000, 0x02000000, 0x04000000, 0x08000000, 
0x10000000, 0x20000000, 0x40000000, 0x80000000
};

/*
   value of mask_bit_index is bit corresponding to integer modulo 2310
   mask_bit_index[i] is zero if this is integer that could not possibly
   be prime or bit number from 1 to 480 if it could be prime.
*/
unsigned int mask_bit_index[2310];

/* bit is one to 480 */
void clear_one_mask_bit(unsigned long *sieve_array, int bit)
{
int k;

k = (bit - 1) >> 5 ;				/* divide by 32 */
sieve_array[k] &= ~(maskbits[(bit - 1) & 31]);	/* modulo 32 */
#ifdef DEBUG
printf("sieve_array[%d]=0x%08x, ~(maskbits[(%d - 1) & 31]=0x%08x\n",
	k, sieve_array[k], bit, ~(maskbits[(bit - 1) & 31]));
#endif
}

/* bit is one to 480 */
int test_one_mask_bit(unsigned long *sieve_array, int bit)
{
int k;

k = (bit - 1) >> 5 ;				/* divide by 32 */
return (sieve_array[k] & (maskbits[(bit - 1) & 31]));	/* modulo 32 */
}


extern char *optarg;	/* used by getopt */
extern int optind, opterr, optopt;

int main(int argc, char *argv[])
{
unsigned b, j, ii, k;
unsigned long i;
unsigned int bit_index = 0;
unsigned long s, indx;
unsigned long start, stop;
int c;
unsigned long *list;
FILE *fp;
unsigned long stop_here, t_stop_here;
char ifnam[2048];
int errflag = 0;
int write_flag = 0;

strcpy(ifnam, "primes.dat");

while ((c = getopt(argc, argv, "s:e:f:n")) != EOF)
  {
  switch(c)
    {
    case 's':
	start = strtoul(optarg,NULL,0);
	break;
    case 'e':
	stop = strtoul(optarg,NULL,0);
	s = (strtoul(optarg,NULL, 0) +2309)/2310 * 60 +60; /* bytes required */
	break;
    case 'f':
	strncpy(ifnam, optarg, sizeof(ifnam) -1);
	ifnam[sizeof(ifnam) -1] = '\0';
	write_flag = 1;
	break;
    case 'n':
	write_flag = 0;
	break;
    default:
	errflag++;
	break;
    }
  }

if ( (errflag != 0) || ( argc > optind) || (argc < 5) || (stop < start))
  {
  fprintf(stderr,"Usage: %s -s start -e end [-f file]\n", argv[0]);
  fprintf(stderr,"Start and end must be positive integers\n");
  fprintf(stderr,
"If a file name is given, primes will be saved to file as 480bits/2310integers.\n");
  if (stop < start)
    fprintf(stderr,"end must be greater than start\n");
  return -1;
  }

/* print the first few small primes if they were requested */
for ( i = 0 ; i < 343 ; i ++)
  {
  if ( small_primes[i] > stop )
    return 0;	/* return to OS if nothing more to do */
  if ( small_primes[i] >= start )
    fprintf(stdout,"%10ld\n",small_primes[i]);
  }


fprintf(stderr,"attempting to malloc %d bytes\n", s);
list = malloc(s);
if ( list == NULL )
  {
  fprintf(stderr, "Could not allocate %lu bytes of memory\n", s);
  exit(1);
  }

memset(list,0xff,s);


j = 5;

/* create array mapping 2310 integers to 480 bits */
mask_bit_index[0] = bit_index++;
mask_bit_index[1] = bit_index++;

for ( ii = 2; ii < 2310 ; ii++ )
  {
  while (small_primes[j] < ii )
    j++;
  /* if modulo any of the prime factors of 2310, then could not be prime */
  if ( ii == small_primes[j] || 
	((ii%2)!=0 && (ii%3)!=0 && (ii%5)!=0 && (ii%7)!=0 && (ii%11)!=0))
    {
    mask_bit_index[ii] = bit_index++;
    }
  else
    mask_bit_index[ii] = 0;
#ifdef DEBUG
  printf("mask_bit_index[%d]=%u, j=%d\n", ii, mask_bit_index[ii], j);
#endif
  }


indx = 0L;
stop_here = (unsigned long) (sqrt((stop+2309)/2310 *2310) + 0.5);

for ( k = 0 ; k*2310 < stop ; k+=STEP_CONST )
  {
  /* no need to sieve numbers larger than this range */
  t_stop_here = (k+STEP_CONST)*2310UL;
  if ( stop_here < t_stop_here )
    t_stop_here = stop_here;
  for( i = 13 ; i <= t_stop_here ; i+=2 )
    {
/* start with 13 since multiples of 2,3,5,7,11 are handled by the storage 
   method */
/* increment by 2 for odd numbers only */
    b = i % 2310;

  /* i could not possibly be prime if remainder is 2,3,4,7,11 */
    if ( mask_bit_index[b] == 0
	|| (i < k*2310
	&& (test_one_mask_bit(&list[i/2310 *15], mask_bit_index[b]))==0 )
	)
      continue;	/* or this one already marked so it is not a prime */
/* */
    if ( k == 0 )
      indx = i*i;
    else
      indx = (k*2310) - (k*2310)%i +i;
    if ( (indx & 1) == 0 )
      indx += i;
/* start with i*i since any integer < i has already been sieved */
/* add 2 * i to avoid even numbers and mark all multiples of this prime */
    for ( ; indx < (k+STEP_CONST)*2310 && indx <= (stop+2309)/2310*2310
		; indx +=(i+i))
      {
      b = indx % 2310;		/* modulo 2310 */
      if ( mask_bit_index[b] != 0 )
        {
        clear_one_mask_bit(&list[indx/2310 *15],mask_bit_index[b]);
#ifdef DEBUG
        printf("indx = %lu, mask_bit_index[%d]=%d\n", indx, b, mask_bit_index[b]);
#endif
        }
      } /* for indx */
    } /* for i */

  if ( start < (k+STEP_CONST)*2310 ) /* are there some to print now? */
    {
    if ( k*2310 < start )
      i = start;
    else
      i = k*2310;
    if ( (i & 1) == 0 )
      i++;	/* force it to be odd */
    if ( 2311 >= i )
      i = 2311;
    for ( ; i <= stop && i < (k+STEP_CONST)*2310; i+=2 )
      {
      b = i % 2310;
      if ( mask_bit_index[b] != 0 && i >= start &&
  	(test_one_mask_bit(&list[i/2310 *15], mask_bit_index[b]))!=0 )
        fprintf(stdout,"%10lu\n",i);
      } /* for i */
    }
  } /* for k */

if ( write_flag )
  {
  if ( (fp = fopen(ifnam,"wb")) == NULL )
    {
    perror(ifnam);
    return 1;
    }
  i = fwrite( list, 1L, s, fp );
  if ( i != s )
    fprintf(stderr,"write error: i=%lu, s=%lu\n",i,s);
  }
return 0;
}


send email to John Moyer, jrm@rsok.com

Copyright 2000, John Moyer. All rights reserved.