BSD 4_4 release

[unix-history] / usr / src / games / primes / primes.c

diff --git a/usr/src/games/primes/primes.c b/usr/src/games/primes/primes.c
index fe737d1..998e8c2 100644 (file)
--- a/usr/src/games/primes/primes.c
+++ b/usr/src/games/primes/primes.c
@@ -1,151 +1,410 @@
 /*
 /*

- * Copyright (c) 1980 Regents of the University of California.
- * All rights reserved.  The Berkeley software License Agreement
- * specifies the terms and conditions for redistribution.
+ * Copyright (c) 1989, 1993
+ *     The Regents of the University of California.  All rights reserved.
+ *
+ * This code is derived from software contributed to Berkeley by
+ * Landon Curt Noll.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *     This product includes software developed by the University of
+ *     California, Berkeley and its contributors.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.

  */
 
 #ifndef lint
  */
 
 #ifndef lint

-char copyright[] =
-"@(#) Copyright (c) 1980 Regents of the University of California.\n\
- All rights reserved.\n";
-#endif not lint
+static char copyright[] =
+"@(#) Copyright (c) 1989, 1993\n\
+       The Regents of the University of California.  All rights reserved.\n";
+#endif /* not lint */

 
 #ifndef lint
 
 #ifndef lint

-static char sccsid[] = "@(#)primes.c   5.1 (Wollongong) %G%";
-#endif not lint
+static char sccsid[] = "@(#)primes.c   8.1 (Berkeley) 5/31/93";
+#endif /* not lint */

 
 /*
 
 /*

- *     primes [ number ]
+ * primes - generate a table of primes between two values
+ *
+ * By: Landon Curt Noll   chongo@toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
+ *
+ *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
+ *
+ * usage:
+ *     primes [start [stop]]

  *
  *

- *     Print all primes greater than argument (or number read from stdin).
+ *     Print primes >= start and < stop.  If stop is omitted,
+ *     the value 4294967295 (2^32-1) is assumed.  If start is
+ *     omitted, start is read from standard input.

  *
  *

- *     A free translation of 'primes.s'
+ *     Prints "ouch" if start or stop is bogus.

  *
  *

+ * validation check: there are 664579 primes between 0 and 10^7

  */
 
 #include <stdio.h>
 #include <math.h>
  */
 
 #include <stdio.h>
 #include <math.h>

+#include <memory.h>
+#include <ctype.h>
+#include "primes.h"

 
 

-#define        TABSIZE 1000            /* size of sieve table */
-#define        BIG     4294967296.     /* largest unsigned int */
+/*
+ * Eratosthenes sieve table
+ *
+ * We only sieve the odd numbers.  The base of our sieve windows are always
+ * odd.  If the base of table is 1, table[i] represents 2*i-1.  After the
+ * sieve, table[i] == 1 if and only iff 2*i-1 is prime.
+ *
+ * We make TABSIZE large to reduce the overhead of inner loop setup.
+ */
+char table[TABSIZE];    /* Eratosthenes sieve of odd numbers */

 
 

-char   table[TABSIZE];         /* table for sieve of Eratosthenes */
-int    tabbits = 8*TABSIZE;    /* number of bits in table */
+/*
+ * prime[i] is the (i-1)th prime.
+ *
+ * We are able to sieve 2^32-1 because this byte table yields all primes 
+ * up to 65537 and 65537^2 > 2^32-1.
+ */
+extern ubig prime[];
+extern ubig *pr_limit;         /* largest prime in the prime array */

 
 

-float  fstart;
-unsigned       start;                  /* lowest number to test for prime */
-char   bittab[] = {            /* bit positions (to save shifting) */
-       01, 02, 04, 010, 020, 040, 0100, 0200
-};
+/*
+ * To avoid excessive sieves for small factors, we use the table below to 
+ * setup our sieve blocks.  Each element represents a odd number starting 
+ * with 1.  All non-zero elements are factors of 3, 5, 7, 11 and 13.
+ */
+extern char pattern[];
+extern int pattern_size;       /* length of pattern array */

 
 

-unsigned pt[] =        {               /* primes < 100 */
-       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
-};
+#define MAX_LINE 255    /* max line allowed on stdin */

 
 

-unsigned factab[] = {          /* difference between succesive trial factors */
-       10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4,
-       2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4,
-       6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2
-};
+char *read_num_buf();   /* read a number buffer */
+void primes();          /* print the primes in range */
+char *program;          /* our name */

 
 main(argc, argv)
 
 main(argc, argv)

-int    argc;
-char   **argv;
+       int argc;       /* arg count */
+       char *argv[];   /* args */

 {
 {

-       register unsigned       *fp;
-       register char   *p;
-       register int    i;
-       unsigned        quot;
-       unsigned        factor, v;
-
-       if (argc >= 2) {                /* get starting no. from arg */
-               if (sscanf(argv[1], "%f", &fstart) != 1
-                   || fstart < 0.0 || fstart >= BIG) {
-                       ouch();
+       char buf[MAX_LINE+1];   /* input buffer */
+       char *ret;      /* return result */
+       ubig start;     /* where to start generating */
+       ubig stop;      /* don't generate at or above this value */
+
+       /*
+        * parse args
+        */
+       program = argv[0];
+       start = 0;
+       stop = BIG;
+       if (argc == 3) {
+               /* convert low and high args */
+               if (read_num_buf(NULL, argv[1]) == NULL) {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       exit(1);
+               }
+               if (read_num_buf(NULL, argv[2]) == NULL) {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       exit(1);
+               }
+               if (sscanf(argv[1], "%ld", &start) != 1) {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       exit(1);
+               }
+               if (sscanf(argv[2], "%ld", &stop) != 1) {
+                       fprintf(stderr, "%s: ouch\n", program);

                        exit(1);
                }
                        exit(1);
                }

-       } else {                        /* get starting no. from stdin */
-               while ((i = scanf("%f", &fstart)) != 1
-                   || fstart < 0.0 || fstart >= BIG) {
-                       if (i == EOF)
+
+       } else if (argc == 2) {
+               /* convert low arg */
+               if (read_num_buf(NULL, argv[1]) == NULL) {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       exit(1);
+               }
+               if (sscanf(argv[1], "%ld", &start) != 1) {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       exit(1);
+               }
+
+       } else {
+               /* read input until we get a good line */
+               if (read_num_buf(stdin, buf) != NULL) {
+
+                       /* convert the buffer */
+                       if (sscanf(buf, "%ld", &start) != 1) {
+                               fprintf(stderr, "%s: ouch\n", program);

                                exit(1);
                                exit(1);

-                       ouch();
+                       }
+               } else {
+                       exit(0);

                }
        }
                }
        }

-       start = (unsigned)fstart;
+       if (start > stop) {
+               fprintf(stderr, "%s: ouch\n", program);
+               exit(1);
+       }
+       primes(start, stop);
+       exit(0);
+}
+
+/*
+ * read_num_buf - read a number buffer from a stream
+ *
+ * Read a number on a line of the form:
+ *
+ *     ^[ \t]*\(+?[0-9][0-9]\)*.*$
+ *
+ * where ? is a 1-or-0 operator and the number is within \( \).
+ *
+ * If does not match the above pattern, it is ignored and a new
+ * line is read.  If the number is too large or small, we will
+ * print ouch and read a new line.
+ *
+ * We have to be very careful on how we check the magnitude of the
+ * input.  We can not use numeric checks because of the need to
+ * check values against maximum numeric values.
+ *
+ * This routine will return a line containing a ascii number between
+ * 0 and BIG, or it will return NULL.
+ *
+ * If the stream is NULL then buf will be processed as if were
+ * a single line stream.
+ *
+ * returns:
+ *     char *  pointer to leading digit or +
+ *     NULL    EOF or error
+ */
+char *
+read_num_buf(input, buf)
+       FILE *input;            /* input stream or NULL */
+       char *buf;              /* input buffer */
+{
+       static char limit[MAX_LINE+1];  /* ascii value of BIG */
+       static int limit_len;           /* digit count of limit */
+       int len;                        /* digits in input (excluding +/-) */
+       char *s;        /* line start marker */
+       char *d;        /* first digit, skip +/- */
+       char *p;        /* scan pointer */
+       char *z;        /* zero scan pointer */

 
 

+       /* form the ascii value of SEMIBIG if needed */
+       if (!isascii(limit[0]) || !isdigit(limit[0])) {
+               sprintf(limit, "%ld", SEMIBIG);
+               limit_len = strlen(limit);
+       }
+       

        /*
        /*

-        * Quick list of primes < 100
+        * the search for a good line

         */
         */

-       if (start <= 97) {
-               for (fp = pt; *fp < start; fp++)
-                       ;
-               do
-                       printf("%u\n", *fp);
-               while (++fp < &pt[sizeof(pt) / sizeof(*pt)]);
-               start = 100;
+       if (input != NULL && fgets(buf, MAX_LINE, input) == NULL) {
+               /* error or EOF */
+               return NULL;

        }
        }

-       quot = start/2;
-       start = quot * 2 + 1;
+       do {
+
+               /* ignore leading whitespace */
+               for (s=buf; *s && s < buf+MAX_LINE; ++s) {
+                       if (!isascii(*s) || !isspace(*s)) {
+                               break;
+                       }
+               }
+
+               /* object if - */
+               if (*s == '-') {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       continue;
+               }
+
+               /* skip over any leading + */
+               if (*s == '+') {
+                       d = s+1;
+               } else {
+                       d = s;
+               }
+
+               /* note leading zeros */
+               for (z=d; *z && z < buf+MAX_LINE; ++z) {
+                       if (*z != '0') {
+                               break;
+                       }
+               }
+
+               /* scan for the first non-digit/non-plus/non-minus */
+               for (p=d; *p && p < buf+MAX_LINE; ++p) {
+                       if (!isascii(*p) || !isdigit(*p)) {
+                               break;
+                       }
+               }
+
+               /* ignore empty lines */
+               if (p == d) {
+                       continue;
+               }
+               *p = '\0';
+
+               /* object if too many digits */
+               len = strlen(z);
+               len = (len<=0) ? 1 : len;
+               /* accept if digit count is below limit */
+               if (len < limit_len) {
+                       /* we have good input */
+                       return s;
+
+               /* reject very large numbers */
+               } else if (len > limit_len) {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       continue;
+
+               /* carefully check against near limit numbers */
+               } else if (strcmp(z, limit) > 0) {
+                       fprintf(stderr, "%s: ouch\n", program);
+                       continue;
+               }
+               /* number is near limit, but is under it */
+               return s;
+       } while (input != NULL && fgets(buf, MAX_LINE, input) != NULL);
+
+       /* error or EOF */
+       return NULL;
+}

 
 /*
 
 /*

- * Loop forever:
+ * primes - sieve and print primes from start up to and but not including stop

  */
  */

-    for (;;) {
+void
+primes(start, stop)
+       ubig start;     /* where to start generating */
+       ubig stop;      /* don't generate at or above this value */
+{
+       register char *q;               /* sieve spot */
+       register ubig factor;           /* index and factor */
+       register char *tab_lim;         /* the limit to sieve on the table */
+       register ubig *p;               /* prime table pointer */
+       register ubig fact_lim;         /* highest prime for current block */
+

        /*
        /*

-        * Generate primes via sieve of Eratosthenes
+        * A number of systems can not convert double values 
+        * into unsigned longs when the values are larger than
+        * the largest signed value.  Thus we take case when
+        * the double is larger than the value SEMIBIG. *sigh*

         */
         */

-       for (p = table; p < &table[TABSIZE]; p++)       /* clear sieve */
-               *p = '\0';
-       v = (unsigned)sqrt((float)(start + tabbits)); /* highest useful factor */
-       sieve(3);
-       sieve(5);
-       sieve(7);
-       factor = 11;
-       fp = &factab[1];
-       do {
-               sieve(factor);
-               factor += *fp;
-               if (++fp >= &factab[sizeof(factab) / sizeof(*factab)])
-                       fp = factab;
-       } while (factor <= v);
+       if (start < 3) {
+               start = (ubig)2;
+       }
+       if (stop < 3) {
+               stop = (ubig)2;
+       }
+       if (stop <= start) {
+               return;
+       }
+

        /*
        /*

-        * Print generated primes
+        * be sure that the values are odd, or 2

         */
         */

-       for (i = 0; i < 8*TABSIZE; i += 2) {
-               if ((table[i>>3] & bittab[i&07]) == 0)
-                       printf("%u\n", start);
-               start += 2;
+       if (start != 2 && (start&0x1) == 0) {
+               ++start;
+       }
+       if (stop != 2 && (stop&0x1) == 0) {
+               ++stop;

        }
        }

-    }
-}

 
 

-/*
- * Insert all multiples of given factor into the sieve
- */
-sieve(factor)
-unsigned factor;
-{
-       register int    i;
-       unsigned        off;
-       unsigned        quot;
-
-       quot = start / factor;
-       off = (quot * factor) - start;
-       if ((int)off < 0)
-               off += factor;
-       while (off < tabbits ) {
-               i = (int)off;
-               table[i>>3] |= bittab[i&07];
-               off += factor;
+       /*
+        * quick list of primes <= pr_limit
+        */
+       if (start <= *pr_limit) {
+               /* skip primes up to the start value */
+               for (p = &prime[0], factor = prime[0];
+                    factor < stop && p <= pr_limit; 
+                    factor = *(++p)) {
+                       if (factor >= start) {
+                               printf("%u\n", factor);
+                       }
+               }
+               /* return early if we are done */
+               if (p <= pr_limit) {
+                       return;
+               }
+               start = *pr_limit+2;

        }
        }

-}

 
 

-/*
- * Error message
- */
-ouch()
-{
-       fprintf(stderr, "Ouch.\n");
+       /*
+        * we shall sieve a bytemap window, note primes and move the window
+        * upward until we pass the stop point
+        */
+       while (start < stop) {
+               /*
+                * factor out 3, 5, 7, 11 and 13
+                */
+               /* initial pattern copy */
+               factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
+               memcpy(table, &pattern[factor], pattern_size-factor);
+               /* main block pattern copies */
+               for (fact_lim=pattern_size-factor;
+                    fact_lim+pattern_size<=TABSIZE;
+                    fact_lim+=pattern_size) {
+                       memcpy(&table[fact_lim], pattern, pattern_size);
+               }
+               /* final block pattern copy */
+               memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
+
+               /*
+                * sieve for primes 17 and higher
+                */
+               /* note highest useful factor and sieve spot */
+               if (stop-start > TABSIZE+TABSIZE) {
+                       tab_lim = &table[TABSIZE]; /* sieve it all */
+                       fact_lim = (int)sqrt(
+                                       (double)(start)+TABSIZE+TABSIZE+1.0);
+               } else {
+                       tab_lim = &table[(stop-start)/2]; /* partial sieve */
+                       fact_lim = (int)sqrt((double)(stop)+1.0);
+               }
+               /* sieve for factors >= 17 */
+               factor = 17;    /* 17 is first prime to use */
+               p = &prime[7];  /* 19 is next prime, pi(19)=7 */
+               do {
+                       /* determine the factor's initial sieve point */
+                       q = (char *)(start%factor); /* temp storage for mod */
+                       if ((int)q & 0x1) {
+                               q = &table[(factor-(int)q)/2];
+                       } else {
+                               q = &table[q ? factor-((int)q/2) : 0];
+                       }
+                       /* sive for our current factor */
+                       for ( ; q < tab_lim; q += factor) {
+                               *q = '\0'; /* sieve out a spot */
+                       }
+               } while ((factor=(ubig)(*(p++))) <= fact_lim);
+
+               /*
+                * print generated primes
+                */
+               for (q = table; q < tab_lim; ++q, start+=2) {
+                       if (*q) {
+                               printf("%u\n", start);
+                       }
+               }
+       }

 }
 }