primes.c revision 38199
12490Sjkh/* 22490Sjkh * Copyright (c) 1989, 1993 32490Sjkh * The Regents of the University of California. All rights reserved. 42490Sjkh * 52490Sjkh * This code is derived from software contributed to Berkeley by 62490Sjkh * Landon Curt Noll. 72490Sjkh * 82490Sjkh * Redistribution and use in source and binary forms, with or without 92490Sjkh * modification, are permitted provided that the following conditions 102490Sjkh * are met: 112490Sjkh * 1. Redistributions of source code must retain the above copyright 122490Sjkh * notice, this list of conditions and the following disclaimer. 132490Sjkh * 2. Redistributions in binary form must reproduce the above copyright 142490Sjkh * notice, this list of conditions and the following disclaimer in the 152490Sjkh * documentation and/or other materials provided with the distribution. 162490Sjkh * 3. All advertising materials mentioning features or use of this software 172490Sjkh * must display the following acknowledgement: 182490Sjkh * This product includes software developed by the University of 192490Sjkh * California, Berkeley and its contributors. 202490Sjkh * 4. Neither the name of the University nor the names of its contributors 212490Sjkh * may be used to endorse or promote products derived from this software 222490Sjkh * without specific prior written permission. 232490Sjkh * 242490Sjkh * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 252490Sjkh * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 262490Sjkh * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 272490Sjkh * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 282490Sjkh * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 292490Sjkh * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 302490Sjkh * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 312490Sjkh * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 322490Sjkh * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 332490Sjkh * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 342490Sjkh * SUCH DAMAGE. 352490Sjkh */ 362490Sjkh 372490Sjkh#ifndef lint 382490Sjkhstatic char copyright[] = 392490Sjkh"@(#) Copyright (c) 1989, 1993\n\ 402490Sjkh The Regents of the University of California. All rights reserved.\n"; 412490Sjkh#endif /* not lint */ 422490Sjkh 432490Sjkh#ifndef lint 4423726Speterstatic char sccsid[] = "@(#)primes.c 8.5 (Berkeley) 5/10/95"; 452490Sjkh#endif /* not lint */ 462490Sjkh 472490Sjkh/* 482490Sjkh * primes - generate a table of primes between two values 492490Sjkh * 502490Sjkh * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo 512490Sjkh * 522490Sjkh * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\ 532490Sjkh * 542490Sjkh * usage: 552490Sjkh * primes [start [stop]] 562490Sjkh * 572490Sjkh * Print primes >= start and < stop. If stop is omitted, 582490Sjkh * the value 4294967295 (2^32-1) is assumed. If start is 592490Sjkh * omitted, start is read from standard input. 602490Sjkh * 612490Sjkh * validation check: there are 664579 primes between 0 and 10^7 622490Sjkh */ 632490Sjkh 642490Sjkh#include <ctype.h> 652490Sjkh#include <err.h> 662490Sjkh#include <errno.h> 672490Sjkh#include <limits.h> 682490Sjkh#include <math.h> 692490Sjkh#include <memory.h> 702490Sjkh#include <stdio.h> 712490Sjkh#include <stdlib.h> 7223726Speter#include <unistd.h> 732490Sjkh 742490Sjkh#include "primes.h" 752490Sjkh 762490Sjkh/* 772490Sjkh * Eratosthenes sieve table 782490Sjkh * 792490Sjkh * We only sieve the odd numbers. The base of our sieve windows are always 802490Sjkh * odd. If the base of table is 1, table[i] represents 2*i-1. After the 812490Sjkh * sieve, table[i] == 1 if and only iff 2*i-1 is prime. 822490Sjkh * 832490Sjkh * We make TABSIZE large to reduce the overhead of inner loop setup. 842490Sjkh */ 852490Sjkhchar table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ 862490Sjkh 872490Sjkh/* 882490Sjkh * prime[i] is the (i-1)th prime. 892490Sjkh * 908856Srgrimes * We are able to sieve 2^32-1 because this byte table yields all primes 912490Sjkh * up to 65537 and 65537^2 > 2^32-1. 922490Sjkh */ 932490Sjkhextern ubig prime[]; 942490Sjkhextern ubig *pr_limit; /* largest prime in the prime array */ 952490Sjkh 962490Sjkh/* 978856Srgrimes * To avoid excessive sieves for small factors, we use the table below to 988856Srgrimes * setup our sieve blocks. Each element represents a odd number starting 992490Sjkh * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13. 1002490Sjkh */ 1012490Sjkhextern char pattern[]; 1022490Sjkhextern int pattern_size; /* length of pattern array */ 1032490Sjkh 10438199Sphkvoid primes __P((ubig, ubig)); 10538199Sphkubig read_num_buf __P((void)); 1062490Sjkhvoid usage __P((void)); 1072490Sjkh 1082490Sjkhint 1092490Sjkhmain(argc, argv) 1102490Sjkh int argc; 1112490Sjkh char *argv[]; 1122490Sjkh{ 1132490Sjkh ubig start; /* where to start generating */ 1142490Sjkh ubig stop; /* don't generate at or above this value */ 1152490Sjkh int ch; 1162490Sjkh char *p; 1172490Sjkh 11838199Sphk while ((ch = getopt(argc, argv, "")) != -1) 1192490Sjkh switch (ch) { 1202490Sjkh case '?': 1212490Sjkh default: 1222490Sjkh usage(); 1232490Sjkh } 1242490Sjkh argc -= optind; 1252490Sjkh argv += optind; 1262490Sjkh 1272490Sjkh start = 0; 1282490Sjkh stop = BIG; 1292490Sjkh 1302490Sjkh /* 1312490Sjkh * Convert low and high args. Strtoul(3) sets errno to 1322490Sjkh * ERANGE if the number is too large, but, if there's 1332490Sjkh * a leading minus sign it returns the negation of the 1342490Sjkh * result of the conversion, which we'd rather disallow. 1352490Sjkh */ 1362490Sjkh switch (argc) { 1372490Sjkh case 2: 1382490Sjkh /* Start and stop supplied on the command line. */ 1392490Sjkh if (argv[0][0] == '-' || argv[1][0] == '-') 1402490Sjkh errx(1, "negative numbers aren't permitted."); 1412490Sjkh 1422490Sjkh errno = 0; 14338199Sphk start = strtoul(argv[0], &p, 10); 1442490Sjkh if (errno) 1452490Sjkh err(1, "%s", argv[0]); 1462490Sjkh if (*p != '\0') 1472490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1482490Sjkh 1492490Sjkh errno = 0; 15038199Sphk stop = strtoul(argv[1], &p, 10); 1512490Sjkh if (errno) 1522490Sjkh err(1, "%s", argv[1]); 1532490Sjkh if (*p != '\0') 1542490Sjkh errx(1, "%s: illegal numeric format.", argv[1]); 1552490Sjkh break; 1562490Sjkh case 1: 1572490Sjkh /* Start on the command line. */ 1582490Sjkh if (argv[0][0] == '-') 1592490Sjkh errx(1, "negative numbers aren't permitted."); 1602490Sjkh 1612490Sjkh errno = 0; 16238199Sphk start = strtoul(argv[0], &p, 10); 1632490Sjkh if (errno) 1642490Sjkh err(1, "%s", argv[0]); 1652490Sjkh if (*p != '\0') 1662490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1672490Sjkh break; 1682490Sjkh case 0: 16938199Sphk start = read_num_buf(); 1702490Sjkh break; 1712490Sjkh default: 1722490Sjkh usage(); 1732490Sjkh } 1742490Sjkh 1752490Sjkh if (start > stop) 1762490Sjkh errx(1, "start value must be less than stop value."); 17738199Sphk primes(start, stop); 1782490Sjkh exit(0); 1792490Sjkh} 1802490Sjkh 1812490Sjkh/* 1822490Sjkh * read_num_buf -- 1832490Sjkh * This routine returns a number n, where 0 <= n && n <= BIG. 1842490Sjkh */ 1852490Sjkhubig 18638199Sphkread_num_buf() 1872490Sjkh{ 1882490Sjkh ubig val; 1892490Sjkh char *p, buf[100]; /* > max number of digits. */ 1902490Sjkh 1912490Sjkh for (;;) { 1922490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1932490Sjkh if (ferror(stdin)) 1942490Sjkh err(1, "stdin"); 1952490Sjkh exit(0); 1962490Sjkh } 1972490Sjkh for (p = buf; isblank(*p); ++p); 1982490Sjkh if (*p == '\n' || *p == '\0') 1992490Sjkh continue; 2002490Sjkh if (*p == '-') 2012490Sjkh errx(1, "negative numbers aren't permitted."); 2022490Sjkh errno = 0; 20338199Sphk val = strtoul(buf, &p, 10); 2042490Sjkh if (errno) 2052490Sjkh err(1, "%s", buf); 2062490Sjkh if (*p != '\n') 2072490Sjkh errx(1, "%s: illegal numeric format.", buf); 2082490Sjkh return (val); 2092490Sjkh } 2102490Sjkh} 2112490Sjkh 2122490Sjkh/* 2132490Sjkh * primes - sieve and print primes from start up to and but not including stop 2142490Sjkh */ 2152490Sjkhvoid 21638199Sphkprimes(start, stop) 2172490Sjkh ubig start; /* where to start generating */ 2182490Sjkh ubig stop; /* don't generate at or above this value */ 2192490Sjkh{ 2202490Sjkh register char *q; /* sieve spot */ 2212490Sjkh register ubig factor; /* index and factor */ 2222490Sjkh register char *tab_lim; /* the limit to sieve on the table */ 2232490Sjkh register ubig *p; /* prime table pointer */ 2242490Sjkh register ubig fact_lim; /* highest prime for current block */ 2252490Sjkh 2262490Sjkh /* 2272490Sjkh * A number of systems can not convert double values into unsigned 2282490Sjkh * longs when the values are larger than the largest signed value. 2292490Sjkh * We don't have this problem, so we can go all the way to BIG. 2302490Sjkh */ 2312490Sjkh if (start < 3) { 2322490Sjkh start = (ubig)2; 2332490Sjkh } 2342490Sjkh if (stop < 3) { 2352490Sjkh stop = (ubig)2; 2362490Sjkh } 2372490Sjkh if (stop <= start) { 2382490Sjkh return; 2392490Sjkh } 2402490Sjkh 2412490Sjkh /* 2422490Sjkh * be sure that the values are odd, or 2 2432490Sjkh */ 2442490Sjkh if (start != 2 && (start&0x1) == 0) { 2452490Sjkh ++start; 2462490Sjkh } 2472490Sjkh if (stop != 2 && (stop&0x1) == 0) { 2482490Sjkh ++stop; 2492490Sjkh } 2502490Sjkh 2512490Sjkh /* 2522490Sjkh * quick list of primes <= pr_limit 2532490Sjkh */ 2542490Sjkh if (start <= *pr_limit) { 2552490Sjkh /* skip primes up to the start value */ 2562490Sjkh for (p = &prime[0], factor = prime[0]; 2572490Sjkh factor < stop && p <= pr_limit; factor = *(++p)) { 2582490Sjkh if (factor >= start) { 25938199Sphk printf("%lu\n", factor); 2602490Sjkh } 2612490Sjkh } 2622490Sjkh /* return early if we are done */ 2632490Sjkh if (p <= pr_limit) { 2642490Sjkh return; 2652490Sjkh } 2662490Sjkh start = *pr_limit+2; 2672490Sjkh } 2682490Sjkh 2692490Sjkh /* 2702490Sjkh * we shall sieve a bytemap window, note primes and move the window 2712490Sjkh * upward until we pass the stop point 2722490Sjkh */ 2732490Sjkh while (start < stop) { 2742490Sjkh /* 2752490Sjkh * factor out 3, 5, 7, 11 and 13 2762490Sjkh */ 2772490Sjkh /* initial pattern copy */ 2782490Sjkh factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 2792490Sjkh memcpy(table, &pattern[factor], pattern_size-factor); 2802490Sjkh /* main block pattern copies */ 2812490Sjkh for (fact_lim=pattern_size-factor; 2822490Sjkh fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { 2832490Sjkh memcpy(&table[fact_lim], pattern, pattern_size); 2842490Sjkh } 2852490Sjkh /* final block pattern copy */ 2862490Sjkh memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 2872490Sjkh 2882490Sjkh /* 2892490Sjkh * sieve for primes 17 and higher 2902490Sjkh */ 2912490Sjkh /* note highest useful factor and sieve spot */ 2922490Sjkh if (stop-start > TABSIZE+TABSIZE) { 2932490Sjkh tab_lim = &table[TABSIZE]; /* sieve it all */ 2942490Sjkh fact_lim = (int)sqrt( 2952490Sjkh (double)(start)+TABSIZE+TABSIZE+1.0); 2962490Sjkh } else { 2972490Sjkh tab_lim = &table[(stop-start)/2]; /* partial sieve */ 2982490Sjkh fact_lim = (int)sqrt((double)(stop)+1.0); 2992490Sjkh } 3002490Sjkh /* sieve for factors >= 17 */ 3012490Sjkh factor = 17; /* 17 is first prime to use */ 3022490Sjkh p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 3032490Sjkh do { 3042490Sjkh /* determine the factor's initial sieve point */ 3052490Sjkh q = (char *)(start%factor); /* temp storage for mod */ 30635892Sjb if ((long)q & 0x1) { 30735892Sjb q = &table[(factor-(long)q)/2]; 3082490Sjkh } else { 30935892Sjb q = &table[q ? factor-((long)q/2) : 0]; 3102490Sjkh } 3112490Sjkh /* sive for our current factor */ 3122490Sjkh for ( ; q < tab_lim; q += factor) { 3132490Sjkh *q = '\0'; /* sieve out a spot */ 3142490Sjkh } 3152490Sjkh } while ((factor=(ubig)(*(p++))) <= fact_lim); 3162490Sjkh 3172490Sjkh /* 3182490Sjkh * print generated primes 3192490Sjkh */ 3202490Sjkh for (q = table; q < tab_lim; ++q, start+=2) { 3212490Sjkh if (*q) { 32238199Sphk printf("%lu\n", start); 3232490Sjkh } 3242490Sjkh } 3252490Sjkh } 3262490Sjkh} 3272490Sjkh 3282490Sjkhvoid 3292490Sjkhusage() 3302490Sjkh{ 3312490Sjkh (void)fprintf(stderr, "usage: primes [start [stop]]\n"); 3322490Sjkh exit(1); 3332490Sjkh} 334