This page is sponsored by
http://www.A-Wee-Bit-of-Ireland.com/
Your source in the USA for the finest Irish wool sweaters, celtic ruanas,
capes, scarves and hats.
Here are photos of Ireland from our trips there to find the best Irish products.
/* 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
