Home | History | Annotate | Line # | Download | only in primes
primes.c revision 1.14
      1  1.14    matt /*	$NetBSD: primes.c,v 1.14 2008/02/02 17:45:05 matt Exp $	*/
      2   1.4     cgd 
      3   1.1     cgd /*
      4   1.4     cgd  * Copyright (c) 1989, 1993
      5   1.4     cgd  *	The Regents of the University of California.  All rights reserved.
      6   1.1     cgd  *
      7   1.1     cgd  * This code is derived from software contributed to Berkeley by
      8   1.1     cgd  * Landon Curt Noll.
      9   1.1     cgd  *
     10   1.1     cgd  * Redistribution and use in source and binary forms, with or without
     11   1.1     cgd  * modification, are permitted provided that the following conditions
     12   1.1     cgd  * are met:
     13   1.1     cgd  * 1. Redistributions of source code must retain the above copyright
     14   1.1     cgd  *    notice, this list of conditions and the following disclaimer.
     15   1.1     cgd  * 2. Redistributions in binary form must reproduce the above copyright
     16   1.1     cgd  *    notice, this list of conditions and the following disclaimer in the
     17   1.1     cgd  *    documentation and/or other materials provided with the distribution.
     18  1.11     agc  * 3. Neither the name of the University nor the names of its contributors
     19   1.1     cgd  *    may be used to endorse or promote products derived from this software
     20   1.1     cgd  *    without specific prior written permission.
     21   1.1     cgd  *
     22   1.1     cgd  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
     23   1.1     cgd  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     24   1.1     cgd  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     25   1.1     cgd  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
     26   1.1     cgd  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     27   1.1     cgd  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     28   1.1     cgd  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     29   1.1     cgd  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     30   1.1     cgd  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     31   1.1     cgd  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     32   1.1     cgd  * SUCH DAMAGE.
     33   1.1     cgd  */
     34   1.1     cgd 
     35   1.7   lukem #include <sys/cdefs.h>
     36   1.1     cgd #ifndef lint
     37   1.7   lukem __COPYRIGHT("@(#) Copyright (c) 1989, 1993\n\
     38   1.7   lukem 	The Regents of the University of California.  All rights reserved.\n");
     39   1.1     cgd #endif /* not lint */
     40   1.1     cgd 
     41   1.1     cgd #ifndef lint
     42   1.4     cgd #if 0
     43   1.6     tls static char sccsid[] = "@(#)primes.c	8.5 (Berkeley) 5/10/95";
     44   1.4     cgd #else
     45  1.14    matt __RCSID("$NetBSD: primes.c,v 1.14 2008/02/02 17:45:05 matt Exp $");
     46   1.4     cgd #endif
     47   1.1     cgd #endif /* not lint */
     48   1.1     cgd 
     49   1.1     cgd /*
     50   1.1     cgd  * primes - generate a table of primes between two values
     51   1.1     cgd  *
     52   1.4     cgd  * By: Landon Curt Noll chongo (at) toad.com, ...!{sun,tolsoft}!hoptoad!chongo
     53   1.1     cgd  *
     54   1.4     cgd  * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
     55   1.1     cgd  *
     56   1.1     cgd  * usage:
     57   1.1     cgd  *	primes [start [stop]]
     58   1.1     cgd  *
     59   1.1     cgd  *	Print primes >= start and < stop.  If stop is omitted,
     60   1.1     cgd  *	the value 4294967295 (2^32-1) is assumed.  If start is
     61   1.1     cgd  *	omitted, start is read from standard input.
     62   1.1     cgd  *
     63   1.1     cgd  * validation check: there are 664579 primes between 0 and 10^7
     64   1.1     cgd  */
     65   1.1     cgd 
     66   1.4     cgd #include <ctype.h>
     67   1.4     cgd #include <err.h>
     68   1.4     cgd #include <errno.h>
     69   1.4     cgd #include <limits.h>
     70   1.1     cgd #include <math.h>
     71   1.1     cgd #include <memory.h>
     72   1.4     cgd #include <stdio.h>
     73   1.4     cgd #include <stdlib.h>
     74   1.6     tls #include <unistd.h>
     75   1.4     cgd 
     76   1.1     cgd #include "primes.h"
     77   1.1     cgd 
     78   1.1     cgd /*
     79   1.1     cgd  * Eratosthenes sieve table
     80   1.1     cgd  *
     81   1.1     cgd  * We only sieve the odd numbers.  The base of our sieve windows are always
     82   1.1     cgd  * odd.  If the base of table is 1, table[i] represents 2*i-1.  After the
     83   1.1     cgd  * sieve, table[i] == 1 if and only iff 2*i-1 is prime.
     84   1.1     cgd  *
     85   1.1     cgd  * We make TABSIZE large to reduce the overhead of inner loop setup.
     86   1.1     cgd  */
     87   1.1     cgd char table[TABSIZE];	 /* Eratosthenes sieve of odd numbers */
     88   1.1     cgd 
     89   1.1     cgd /*
     90   1.1     cgd  * prime[i] is the (i-1)th prime.
     91   1.1     cgd  *
     92   1.1     cgd  * We are able to sieve 2^32-1 because this byte table yields all primes
     93   1.1     cgd  * up to 65537 and 65537^2 > 2^32-1.
     94   1.1     cgd  */
     95   1.9     jsm extern const ubig prime[];
     96   1.9     jsm extern const ubig *pr_limit;		/* largest prime in the prime array */
     97   1.1     cgd 
     98   1.1     cgd /*
     99   1.1     cgd  * To avoid excessive sieves for small factors, we use the table below to
    100   1.1     cgd  * setup our sieve blocks.  Each element represents a odd number starting
    101   1.1     cgd  * with 1.  All non-zero elements are factors of 3, 5, 7, 11 and 13.
    102   1.1     cgd  */
    103   1.9     jsm extern const char pattern[];
    104   1.9     jsm extern const int pattern_size;	/* length of pattern array */
    105   1.1     cgd 
    106  1.12     jsm int	main(int, char *[]);
    107  1.12     jsm void	primes(ubig, ubig);
    108  1.12     jsm ubig	read_num_buf(void);
    109  1.13   perry void	usage(void) __dead;
    110   1.1     cgd 
    111   1.4     cgd int
    112   1.1     cgd main(argc, argv)
    113   1.4     cgd 	int argc;
    114   1.4     cgd 	char *argv[];
    115   1.1     cgd {
    116   1.4     cgd 	ubig start;		/* where to start generating */
    117   1.4     cgd 	ubig stop;		/* don't generate at or above this value */
    118   1.4     cgd 	int ch;
    119   1.4     cgd 	char *p;
    120   1.4     cgd 
    121   1.7   lukem 	while ((ch = getopt(argc, argv, "")) != -1)
    122   1.4     cgd 		switch (ch) {
    123   1.4     cgd 		case '?':
    124   1.4     cgd 		default:
    125   1.4     cgd 			usage();
    126   1.4     cgd 		}
    127   1.4     cgd 	argc -= optind;
    128   1.4     cgd 	argv += optind;
    129   1.1     cgd 
    130   1.1     cgd 	start = 0;
    131   1.1     cgd 	stop = BIG;
    132   1.1     cgd 
    133   1.4     cgd 	/*
    134   1.4     cgd 	 * Convert low and high args.  Strtoul(3) sets errno to
    135   1.4     cgd 	 * ERANGE if the number is too large, but, if there's
    136   1.4     cgd 	 * a leading minus sign it returns the negation of the
    137   1.4     cgd 	 * result of the conversion, which we'd rather disallow.
    138   1.4     cgd 	 */
    139   1.4     cgd 	switch (argc) {
    140   1.4     cgd 	case 2:
    141   1.4     cgd 		/* Start and stop supplied on the command line. */
    142   1.4     cgd 		if (argv[0][0] == '-' || argv[1][0] == '-')
    143   1.4     cgd 			errx(1, "negative numbers aren't permitted.");
    144   1.4     cgd 
    145   1.4     cgd 		errno = 0;
    146   1.4     cgd 		start = strtoul(argv[0], &p, 10);
    147   1.4     cgd 		if (errno)
    148   1.4     cgd 			err(1, "%s", argv[0]);
    149   1.4     cgd 		if (*p != '\0')
    150   1.4     cgd 			errx(1, "%s: illegal numeric format.", argv[0]);
    151   1.4     cgd 
    152   1.4     cgd 		errno = 0;
    153   1.4     cgd 		stop = strtoul(argv[1], &p, 10);
    154   1.4     cgd 		if (errno)
    155   1.4     cgd 			err(1, "%s", argv[1]);
    156   1.4     cgd 		if (*p != '\0')
    157   1.4     cgd 			errx(1, "%s: illegal numeric format.", argv[1]);
    158   1.4     cgd 		break;
    159   1.4     cgd 	case 1:
    160   1.4     cgd 		/* Start on the command line. */
    161   1.4     cgd 		if (argv[0][0] == '-')
    162   1.4     cgd 			errx(1, "negative numbers aren't permitted.");
    163   1.4     cgd 
    164   1.4     cgd 		errno = 0;
    165   1.4     cgd 		start = strtoul(argv[0], &p, 10);
    166   1.4     cgd 		if (errno)
    167   1.4     cgd 			err(1, "%s", argv[0]);
    168   1.4     cgd 		if (*p != '\0')
    169   1.4     cgd 			errx(1, "%s: illegal numeric format.", argv[0]);
    170   1.4     cgd 		break;
    171   1.4     cgd 	case 0:
    172   1.4     cgd 		start = read_num_buf();
    173   1.4     cgd 		break;
    174   1.4     cgd 	default:
    175   1.4     cgd 		usage();
    176   1.4     cgd 	}
    177   1.1     cgd 
    178   1.4     cgd 	if (start > stop)
    179   1.4     cgd 		errx(1, "start value must be less than stop value.");
    180   1.1     cgd 	primes(start, stop);
    181   1.1     cgd 	exit(0);
    182   1.1     cgd }
    183   1.1     cgd 
    184   1.1     cgd /*
    185   1.4     cgd  * read_num_buf --
    186   1.4     cgd  *	This routine returns a number n, where 0 <= n && n <= BIG.
    187   1.1     cgd  */
    188   1.4     cgd ubig
    189   1.4     cgd read_num_buf()
    190   1.1     cgd {
    191   1.4     cgd 	ubig val;
    192   1.4     cgd 	char *p, buf[100];		/* > max number of digits. */
    193   1.1     cgd 
    194   1.4     cgd 	for (;;) {
    195   1.4     cgd 		if (fgets(buf, sizeof(buf), stdin) == NULL) {
    196   1.4     cgd 			if (ferror(stdin))
    197   1.4     cgd 				err(1, "stdin");
    198   1.4     cgd 			exit(0);
    199   1.1     cgd 		}
    200   1.4     cgd 		for (p = buf; isblank(*p); ++p);
    201   1.4     cgd 		if (*p == '\n' || *p == '\0')
    202   1.1     cgd 			continue;
    203   1.4     cgd 		if (*p == '-')
    204   1.4     cgd 			errx(1, "negative numbers aren't permitted.");
    205   1.4     cgd 		errno = 0;
    206   1.4     cgd 		val = strtoul(buf, &p, 10);
    207   1.4     cgd 		if (errno)
    208   1.4     cgd 			err(1, "%s", buf);
    209   1.4     cgd 		if (*p != '\n')
    210   1.4     cgd 			errx(1, "%s: illegal numeric format.", buf);
    211   1.4     cgd 		return (val);
    212   1.4     cgd 	}
    213   1.1     cgd }
    214   1.1     cgd 
    215   1.1     cgd /*
    216   1.1     cgd  * primes - sieve and print primes from start up to and but not including stop
    217   1.1     cgd  */
    218   1.1     cgd void
    219   1.1     cgd primes(start, stop)
    220   1.1     cgd 	ubig start;	/* where to start generating */
    221   1.1     cgd 	ubig stop;	/* don't generate at or above this value */
    222   1.1     cgd {
    223   1.7   lukem 	char *q;		/* sieve spot */
    224   1.7   lukem 	ubig factor;		/* index and factor */
    225   1.7   lukem 	char *tab_lim;		/* the limit to sieve on the table */
    226   1.9     jsm 	const ubig *p;		/* prime table pointer */
    227   1.7   lukem 	ubig fact_lim;		/* highest prime for current block */
    228  1.10  itojun 	ubig mod;		/* temp storage for mod */
    229   1.1     cgd 
    230   1.1     cgd 	/*
    231   1.4     cgd 	 * A number of systems can not convert double values into unsigned
    232   1.4     cgd 	 * longs when the values are larger than the largest signed value.
    233   1.4     cgd 	 * We don't have this problem, so we can go all the way to BIG.
    234   1.1     cgd 	 */
    235   1.1     cgd 	if (start < 3) {
    236   1.1     cgd 		start = (ubig)2;
    237   1.1     cgd 	}
    238   1.1     cgd 	if (stop < 3) {
    239   1.1     cgd 		stop = (ubig)2;
    240   1.1     cgd 	}
    241   1.1     cgd 	if (stop <= start) {
    242   1.1     cgd 		return;
    243   1.1     cgd 	}
    244   1.1     cgd 
    245   1.1     cgd 	/*
    246   1.1     cgd 	 * be sure that the values are odd, or 2
    247   1.1     cgd 	 */
    248   1.1     cgd 	if (start != 2 && (start&0x1) == 0) {
    249   1.1     cgd 		++start;
    250   1.1     cgd 	}
    251   1.1     cgd 	if (stop != 2 && (stop&0x1) == 0) {
    252   1.1     cgd 		++stop;
    253   1.1     cgd 	}
    254   1.1     cgd 
    255   1.1     cgd 	/*
    256   1.1     cgd 	 * quick list of primes <= pr_limit
    257   1.1     cgd 	 */
    258   1.1     cgd 	if (start <= *pr_limit) {
    259   1.1     cgd 		/* skip primes up to the start value */
    260   1.1     cgd 		for (p = &prime[0], factor = prime[0];
    261   1.4     cgd 		    factor < stop && p <= pr_limit; factor = *(++p)) {
    262   1.1     cgd 			if (factor >= start) {
    263   1.7   lukem 				printf("%lu\n", (unsigned long) factor);
    264   1.1     cgd 			}
    265   1.1     cgd 		}
    266   1.1     cgd 		/* return early if we are done */
    267   1.1     cgd 		if (p <= pr_limit) {
    268   1.1     cgd 			return;
    269   1.1     cgd 		}
    270   1.1     cgd 		start = *pr_limit+2;
    271   1.1     cgd 	}
    272   1.1     cgd 
    273   1.1     cgd 	/*
    274   1.1     cgd 	 * we shall sieve a bytemap window, note primes and move the window
    275   1.1     cgd 	 * upward until we pass the stop point
    276   1.1     cgd 	 */
    277   1.1     cgd 	while (start < stop) {
    278   1.1     cgd 		/*
    279   1.1     cgd 		 * factor out 3, 5, 7, 11 and 13
    280   1.1     cgd 		 */
    281   1.1     cgd 		/* initial pattern copy */
    282   1.1     cgd 		factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
    283   1.1     cgd 		memcpy(table, &pattern[factor], pattern_size-factor);
    284   1.1     cgd 		/* main block pattern copies */
    285   1.1     cgd 		for (fact_lim=pattern_size-factor;
    286   1.4     cgd 		    fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) {
    287   1.1     cgd 			memcpy(&table[fact_lim], pattern, pattern_size);
    288   1.1     cgd 		}
    289   1.1     cgd 		/* final block pattern copy */
    290   1.1     cgd 		memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
    291   1.1     cgd 
    292   1.1     cgd 		/*
    293   1.1     cgd 		 * sieve for primes 17 and higher
    294   1.1     cgd 		 */
    295   1.1     cgd 		/* note highest useful factor and sieve spot */
    296   1.1     cgd 		if (stop-start > TABSIZE+TABSIZE) {
    297   1.1     cgd 			tab_lim = &table[TABSIZE]; /* sieve it all */
    298   1.1     cgd 			fact_lim = (int)sqrt(
    299   1.1     cgd 					(double)(start)+TABSIZE+TABSIZE+1.0);
    300   1.1     cgd 		} else {
    301   1.1     cgd 			tab_lim = &table[(stop-start)/2]; /* partial sieve */
    302   1.1     cgd 			fact_lim = (int)sqrt((double)(stop)+1.0);
    303   1.1     cgd 		}
    304   1.1     cgd 		/* sieve for factors >= 17 */
    305   1.1     cgd 		factor = 17;	/* 17 is first prime to use */
    306   1.1     cgd 		p = &prime[7];	/* 19 is next prime, pi(19)=7 */
    307   1.1     cgd 		do {
    308   1.1     cgd 			/* determine the factor's initial sieve point */
    309  1.10  itojun 			mod = start%factor;
    310  1.10  itojun 			if (mod & 0x1) {
    311  1.10  itojun 				q = &table[(factor-mod)/2];
    312   1.1     cgd 			} else {
    313  1.10  itojun 				q = &table[mod ? factor-(mod/2) : 0];
    314   1.1     cgd 			}
    315  1.14    matt 			/* sieve for our current factor */
    316   1.1     cgd 			for ( ; q < tab_lim; q += factor) {
    317   1.1     cgd 				*q = '\0'; /* sieve out a spot */
    318   1.1     cgd 			}
    319   1.1     cgd 		} while ((factor=(ubig)(*(p++))) <= fact_lim);
    320   1.1     cgd 
    321   1.1     cgd 		/*
    322   1.1     cgd 		 * print generated primes
    323   1.1     cgd 		 */
    324   1.1     cgd 		for (q = table; q < tab_lim; ++q, start+=2) {
    325   1.1     cgd 			if (*q) {
    326   1.7   lukem 				printf("%lu\n", (unsigned long) start);
    327   1.1     cgd 			}
    328   1.1     cgd 		}
    329   1.1     cgd 	}
    330   1.4     cgd }
    331   1.4     cgd 
    332   1.4     cgd void
    333   1.4     cgd usage()
    334   1.4     cgd {
    335   1.4     cgd 	(void)fprintf(stderr, "usage: primes [start [stop]]\n");
    336   1.4     cgd 	exit(1);
    337   1.1     cgd }
    338