Return to primes.c CVS log

Up to [CSRG BSD Unix] / 43BSDReno / games / primes

Annotation of 43BSDReno/games/primes/primes.c, revision 1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.