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