home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / games / primes / primes.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-04-08  |  10.9 KB  |  411 lines
  1. /*
  2.  * Copyright (c) 1989 The Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * This code is derived from software contributed to Berkeley by
  6.  * Landon Curt Noll.
  7.  *
  8.  * Redistribution and use in source and binary forms, with or without
  9.  * modification, are permitted provided that the following conditions
  10.  * are met:
  11.  * 1. Redistributions of source code must retain the above copyright
  12.  *    notice, this list of conditions and the following disclaimer.
  13.  * 2. Redistributions in binary form must reproduce the above copyright
  14.  *    notice, this list of conditions and the following disclaimer in the
  15.  *    documentation and/or other materials provided with the distribution.
  16.  * 3. All advertising materials mentioning features or use of this software
  17.  *    must display the following acknowledgement:
  18.  *    This product includes software developed by the University of
  19.  *    California, Berkeley and its contributors.
  20.  * 4. Neither the name of the University nor the names of its contributors
  21.  *    may be used to endorse or promote products derived from this software
  22.  *    without specific prior written permission.
  23.  *
  24.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  25.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  26.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  27.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  28.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  29.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  30.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  31.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  32.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  33.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  34.  * SUCH DAMAGE.
  35.  */
  36.  
  37. #ifndef lint
  38. char copyright[] =
  39. "@(#) Copyright (c) 1989 The Regents of the University of California.\n\
  40.  All rights reserved.\n";
  41. #endif /* not lint */
  42.  
  43. #ifndef lint
  44. static char sccsid[] = "@(#)primes.c    5.4 (Berkeley) 6/1/90";
  45. #endif /* not lint */
  46.  
  47. /*
  48.  * primes - generate a table of primes between two values
  49.  *
  50.  * By: Landon Curt Noll   chongo@toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
  51.  *
  52.  *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
  53.  *
  54.  * usage:
  55.  *    primes [start [stop]]
  56.  *
  57.  *    Print primes >= start and < stop.  If stop is omitted,
  58.  *    the value 4294967295 (2^32-1) is assumed.  If start is
  59.  *    omitted, start is read from standard input.
  60.  *
  61.  *    Prints "ouch" if start or stop is bogus.
  62.  *
  63.  * validation check: there are 664579 primes between 0 and 10^7
  64.  */
  65.  
  66. #include <stdio.h>
  67. #include <math.h>
  68. #include <memory.h>
  69. #include <ctype.h>
  70. #include "primes.h"
  71.  
  72. /*
  73.  * Eratosthenes sieve table
  74.  *
  75.  * We only sieve the odd numbers.  The base of our sieve windows are always
  76.  * odd.  If the base of table is 1, table[i] represents 2*i-1.  After the
  77.  * sieve, table[i] == 1 if and only iff 2*i-1 is prime.
  78.  *
  79.  * We make TABSIZE large to reduce the overhead of inner loop setup.
  80.  */
  81. char table[TABSIZE];     /* Eratosthenes sieve of odd numbers */
  82.  
  83. /*
  84.  * prime[i] is the (i-1)th prime.
  85.  *
  86.  * We are able to sieve 2^32-1 because this byte table yields all primes 
  87.  * up to 65537 and 65537^2 > 2^32-1.
  88.  */
  89. extern ubig prime[];
  90. extern ubig *pr_limit;        /* largest prime in the prime array */
  91.  
  92. /*
  93.  * To avoid excessive sieves for small factors, we use the table below to 
  94.  * setup our sieve blocks.  Each element represents a odd number starting 
  95.  * with 1.  All non-zero elements are factors of 3, 5, 7, 11 and 13.
  96.  */
  97. extern char pattern[];
  98. extern int pattern_size;    /* length of pattern array */
  99.  
  100. #define MAX_LINE 255    /* max line allowed on stdin */
  101.  
  102. char *read_num_buf();     /* read a number buffer */
  103. void primes();         /* print the primes in range */
  104. char *program;         /* our name */
  105.  
  106. main(argc, argv)
  107.     int argc;    /* arg count */
  108.     char *argv[];    /* args */
  109. {
  110.     char buf[MAX_LINE+1];   /* input buffer */
  111.     char *ret;    /* return result */
  112.     ubig start;    /* where to start generating */
  113.     ubig stop;    /* don't generate at or above this value */
  114.  
  115.     /*
  116.      * parse args
  117.      */
  118.     program = argv[0];
  119.     start = 0;
  120.     stop = BIG;
  121.     if (argc == 3) {
  122.         /* convert low and high args */
  123.         if (read_num_buf(NULL, argv[1]) == NULL) {
  124.             fprintf(stderr, "%s: ouch\n", program);
  125.             exit(1);
  126.         }
  127.         if (read_num_buf(NULL, argv[2]) == NULL) {
  128.             fprintf(stderr, "%s: ouch\n", program);
  129.             exit(1);
  130.         }
  131.         if (sscanf(argv[1], "%ld", &start) != 1) {
  132.             fprintf(stderr, "%s: ouch\n", program);
  133.             exit(1);
  134.         }
  135.         if (sscanf(argv[2], "%ld", &stop) != 1) {
  136.             fprintf(stderr, "%s: ouch\n", program);
  137.             exit(1);
  138.         }
  139.  
  140.     } else if (argc == 2) {
  141.         /* convert low arg */
  142.         if (read_num_buf(NULL, argv[1]) == NULL) {
  143.             fprintf(stderr, "%s: ouch\n", program);
  144.             exit(1);
  145.         }
  146.         if (sscanf(argv[1], "%ld", &start) != 1) {
  147.             fprintf(stderr, "%s: ouch\n", program);
  148.             exit(1);
  149.         }
  150.  
  151.     } else {
  152.         /* read input until we get a good line */
  153.         if (read_num_buf(stdin, buf) != NULL) {
  154.  
  155.             /* convert the buffer */
  156.             if (sscanf(buf, "%ld", &start) != 1) {
  157.                 fprintf(stderr, "%s: ouch\n", program);
  158.                 exit(1);
  159.             }
  160.         } else {
  161.             exit(0);
  162.         }
  163.     }
  164.     if (start > stop) {
  165.         fprintf(stderr, "%s: ouch\n", program);
  166.         exit(1);
  167.     }
  168.     primes(start, stop);
  169.     exit(0);
  170. }
  171.  
  172. /*
  173.  * read_num_buf - read a number buffer from a stream
  174.  *
  175.  * Read a number on a line of the form:
  176.  *
  177.  *    ^[ \t]*\(+?[0-9][0-9]\)*.*$
  178.  *
  179.  * where ? is a 1-or-0 operator and the number is within \( \).
  180.  *
  181.  * If does not match the above pattern, it is ignored and a new
  182.  * line is read.  If the number is too large or small, we will
  183.  * print ouch and read a new line.
  184.  *
  185.  * We have to be very careful on how we check the magnitude of the
  186.  * input.  We can not use numeric checks because of the need to
  187.  * check values against maximum numeric values.
  188.  *
  189.  * This routine will return a line containing a ascii number between
  190.  * 0 and BIG, or it will return NULL.
  191.  *
  192.  * If the stream is NULL then buf will be processed as if were
  193.  * a single line stream.
  194.  *
  195.  * returns:
  196.  *    char *    pointer to leading digit or +
  197.  *    NULL    EOF or error
  198.  */
  199. char *
  200. read_num_buf(input, buf)
  201.     FILE *input;        /* input stream or NULL */
  202.     char *buf;        /* input buffer */
  203. {
  204.     static char limit[MAX_LINE+1];    /* ascii value of BIG */
  205.     static int limit_len;        /* digit count of limit */
  206.     int len;            /* digits in input (excluding +/-) */
  207.     char *s;    /* line start marker */
  208.     char *d;    /* first digit, skip +/- */
  209.     char *p;    /* scan pointer */
  210.     char *z;    /* zero scan pointer */
  211.  
  212.     /* form the ascii value of SEMIBIG if needed */
  213.     if (!isascii(limit[0]) || !isdigit(limit[0])) {
  214.         sprintf(limit, "%ld", SEMIBIG);
  215.         limit_len = strlen(limit);
  216.     }
  217.     
  218.     /*
  219.      * the search for a good line
  220.      */
  221.     if (input != NULL && fgets(buf, MAX_LINE, input) == NULL) {
  222.         /* error or EOF */
  223.         return NULL;
  224.     }
  225.     do {
  226.  
  227.         /* ignore leading whitespace */
  228.         for (s=buf; *s && s < buf+MAX_LINE; ++s) {
  229.             if (!isascii(*s) || !isspace(*s)) {
  230.                 break;
  231.             }
  232.         }
  233.  
  234.         /* object if - */
  235.         if (*s == '-') {
  236.             fprintf(stderr, "%s: ouch\n", program);
  237.             continue;
  238.         }
  239.  
  240.         /* skip over any leading + */
  241.         if (*s == '+') {
  242.             d = s+1;
  243.         } else {
  244.             d = s;
  245.         }
  246.  
  247.         /* note leading zeros */
  248.         for (z=d; *z && z < buf+MAX_LINE; ++z) {
  249.             if (*z != '0') {
  250.                 break;
  251.             }
  252.         }
  253.  
  254.         /* scan for the first non-digit/non-plus/non-minus */
  255.         for (p=d; *p && p < buf+MAX_LINE; ++p) {
  256.             if (!isascii(*p) || !isdigit(*p)) {
  257.                 break;
  258.             }
  259.         }
  260.  
  261.         /* ignore empty lines */
  262.         if (p == d) {
  263.             continue;
  264.         }
  265.         *p = '\0';
  266.  
  267.         /* object if too many digits */
  268.         len = strlen(z);
  269.         len = (len<=0) ? 1 : len;
  270.         /* accept if digit count is below limit */
  271.         if (len < limit_len) {
  272.             /* we have good input */
  273.             return s;
  274.  
  275.         /* reject very large numbers */
  276.         } else if (len > limit_len) {
  277.             fprintf(stderr, "%s: ouch\n", program);
  278.             continue;
  279.  
  280.         /* carefully check against near limit numbers */
  281.         } else if (strcmp(z, limit) > 0) {
  282.             fprintf(stderr, "%s: ouch\n", program);
  283.             continue;
  284.         }
  285.         /* number is near limit, but is under it */
  286.         return s;
  287.     } while (input != NULL && fgets(buf, MAX_LINE, input) != NULL);
  288.  
  289.     /* error or EOF */
  290.     return NULL;
  291. }
  292.  
  293. /*
  294.  * primes - sieve and print primes from start up to and but not including stop
  295.  */
  296. void
  297. primes(start, stop)
  298.     ubig start;    /* where to start generating */
  299.     ubig stop;    /* don't generate at or above this value */
  300. {
  301.     register char *q;        /* sieve spot */
  302.     register ubig factor;        /* index and factor */
  303.     register char *tab_lim;        /* the limit to sieve on the table */
  304.     register ubig *p;        /* prime table pointer */
  305.     register ubig fact_lim;        /* highest prime for current block */
  306.  
  307.     /*
  308.      * A number of systems can not convert double values 
  309.      * into unsigned longs when the values are larger than
  310.      * the largest signed value.  Thus we take case when
  311.      * the double is larger than the value SEMIBIG. *sigh*
  312.      */
  313.     if (start < 3) {
  314.         start = (ubig)2;
  315.     }
  316.     if (stop < 3) {
  317.         stop = (ubig)2;
  318.     }
  319.     if (stop <= start) {
  320.         return;
  321.     }
  322.  
  323.     /*
  324.      * be sure that the values are odd, or 2
  325.      */
  326.     if (start != 2 && (start&0x1) == 0) {
  327.         ++start;
  328.     }
  329.     if (stop != 2 && (stop&0x1) == 0) {
  330.         ++stop;
  331.     }
  332.  
  333.     /*
  334.      * quick list of primes <= pr_limit
  335.      */
  336.     if (start <= *pr_limit) {
  337.         /* skip primes up to the start value */
  338.         for (p = &prime[0], factor = prime[0];
  339.              factor < stop && p <= pr_limit; 
  340.              factor = *(++p)) {
  341.             if (factor >= start) {
  342.                 printf("%u\n", factor);
  343.             }
  344.         }
  345.         /* return early if we are done */
  346.         if (p <= pr_limit) {
  347.             return;
  348.         }
  349.         start = *pr_limit+2;
  350.     }
  351.  
  352.     /*
  353.      * we shall sieve a bytemap window, note primes and move the window
  354.      * upward until we pass the stop point
  355.      */
  356.     while (start < stop) {
  357.         /*
  358.          * factor out 3, 5, 7, 11 and 13
  359.          */
  360.         /* initial pattern copy */
  361.         factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
  362.         memcpy(table, &pattern[factor], pattern_size-factor);
  363.         /* main block pattern copies */
  364.         for (fact_lim=pattern_size-factor;
  365.              fact_lim+pattern_size<=TABSIZE;
  366.              fact_lim+=pattern_size) {
  367.             memcpy(&table[fact_lim], pattern, pattern_size);
  368.         }
  369.         /* final block pattern copy */
  370.         memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
  371.  
  372.         /*
  373.          * sieve for primes 17 and higher
  374.          */
  375.         /* note highest useful factor and sieve spot */
  376.         if (stop-start > TABSIZE+TABSIZE) {
  377.             tab_lim = &table[TABSIZE]; /* sieve it all */
  378.             fact_lim = (int)sqrt(
  379.                     (double)(start)+TABSIZE+TABSIZE+1.0);
  380.         } else {
  381.             tab_lim = &table[(stop-start)/2]; /* partial sieve */
  382.             fact_lim = (int)sqrt((double)(stop)+1.0);
  383.         }
  384.         /* sieve for factors >= 17 */
  385.         factor = 17;    /* 17 is first prime to use */
  386.         p = &prime[7];    /* 19 is next prime, pi(19)=7 */
  387.         do {
  388.             /* determine the factor's initial sieve point */
  389.             q = (char *)(start%factor); /* temp storage for mod */
  390.             if ((int)q & 0x1) {
  391.                 q = &table[(factor-(int)q)/2];
  392.             } else {
  393.                 q = &table[q ? factor-((int)q/2) : 0];
  394.             }
  395.             /* sive for our current factor */
  396.             for ( ; q < tab_lim; q += factor) {
  397.                 *q = '\0'; /* sieve out a spot */
  398.             }
  399.         } while ((factor=(ubig)(*(p++))) <= fact_lim);
  400.  
  401.         /*
  402.          * print generated primes
  403.          */
  404.         for (q = table; q < tab_lim; ++q, start+=2) {
  405.             if (*q) {
  406.                 printf("%u\n", start);
  407.             }
  408.         }
  409.     }
  410. }
  411.