primes.c revision 42357
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 10442338Simpint hflag; 10542338Simp 10638199Sphkvoid primes __P((ubig, ubig)); 10738199Sphkubig read_num_buf __P((void)); 1082490Sjkhvoid usage __P((void)); 1092490Sjkh 1102490Sjkhint 1112490Sjkhmain(argc, argv) 1122490Sjkh int argc; 1132490Sjkh char *argv[]; 1142490Sjkh{ 1152490Sjkh ubig start; /* where to start generating */ 1162490Sjkh ubig stop; /* don't generate at or above this value */ 1172490Sjkh int ch; 1182490Sjkh char *p; 1192490Sjkh 12042338Simp while ((ch = getopt(argc, argv, "h")) != -1) 1212490Sjkh switch (ch) { 12242338Simp case 'h': 12342338Simp hflag++; 12442338Simp break; 1252490Sjkh case '?': 1262490Sjkh default: 1272490Sjkh usage(); 1282490Sjkh } 1292490Sjkh argc -= optind; 1302490Sjkh argv += optind; 1312490Sjkh 1322490Sjkh start = 0; 1332490Sjkh stop = BIG; 1342490Sjkh 1352490Sjkh /* 1362490Sjkh * Convert low and high args. Strtoul(3) sets errno to 1372490Sjkh * ERANGE if the number is too large, but, if there's 1382490Sjkh * a leading minus sign it returns the negation of the 1392490Sjkh * result of the conversion, which we'd rather disallow. 1402490Sjkh */ 1412490Sjkh switch (argc) { 1422490Sjkh case 2: 1432490Sjkh /* Start and stop supplied on the command line. */ 1442490Sjkh if (argv[0][0] == '-' || argv[1][0] == '-') 1452490Sjkh errx(1, "negative numbers aren't permitted."); 1462490Sjkh 1472490Sjkh errno = 0; 14842338Simp start = strtoul(argv[0], &p, 0); 1492490Sjkh if (errno) 1502490Sjkh err(1, "%s", argv[0]); 1512490Sjkh if (*p != '\0') 1522490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1532490Sjkh 1542490Sjkh errno = 0; 15542338Simp stop = strtoul(argv[1], &p, 0); 1562490Sjkh if (errno) 1572490Sjkh err(1, "%s", argv[1]); 1582490Sjkh if (*p != '\0') 1592490Sjkh errx(1, "%s: illegal numeric format.", argv[1]); 1602490Sjkh break; 1612490Sjkh case 1: 1622490Sjkh /* Start on the command line. */ 1632490Sjkh if (argv[0][0] == '-') 1642490Sjkh errx(1, "negative numbers aren't permitted."); 1652490Sjkh 1662490Sjkh errno = 0; 16742338Simp start = strtoul(argv[0], &p, 0); 1682490Sjkh if (errno) 1692490Sjkh err(1, "%s", argv[0]); 1702490Sjkh if (*p != '\0') 1712490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1722490Sjkh break; 1732490Sjkh case 0: 17438199Sphk start = read_num_buf(); 1752490Sjkh break; 1762490Sjkh default: 1772490Sjkh usage(); 1782490Sjkh } 1792490Sjkh 1802490Sjkh if (start > stop) 1812490Sjkh errx(1, "start value must be less than stop value."); 18238199Sphk primes(start, stop); 1832490Sjkh exit(0); 1842490Sjkh} 1852490Sjkh 1862490Sjkh/* 1872490Sjkh * read_num_buf -- 1882490Sjkh * This routine returns a number n, where 0 <= n && n <= BIG. 1892490Sjkh */ 1902490Sjkhubig 19138199Sphkread_num_buf() 1922490Sjkh{ 1932490Sjkh ubig val; 1942490Sjkh char *p, buf[100]; /* > max number of digits. */ 1952490Sjkh 1962490Sjkh for (;;) { 1972490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1982490Sjkh if (ferror(stdin)) 1992490Sjkh err(1, "stdin"); 2002490Sjkh exit(0); 2012490Sjkh } 2022490Sjkh for (p = buf; isblank(*p); ++p); 2032490Sjkh if (*p == '\n' || *p == '\0') 2042490Sjkh continue; 2052490Sjkh if (*p == '-') 2062490Sjkh errx(1, "negative numbers aren't permitted."); 2072490Sjkh errno = 0; 20842338Simp val = strtoul(buf, &p, 0); 2092490Sjkh if (errno) 2102490Sjkh err(1, "%s", buf); 2112490Sjkh if (*p != '\n') 2122490Sjkh errx(1, "%s: illegal numeric format.", buf); 2132490Sjkh return (val); 2142490Sjkh } 2152490Sjkh} 2162490Sjkh 2172490Sjkh/* 2182490Sjkh * primes - sieve and print primes from start up to and but not including stop 2192490Sjkh */ 2202490Sjkhvoid 22138199Sphkprimes(start, stop) 2222490Sjkh ubig start; /* where to start generating */ 2232490Sjkh ubig stop; /* don't generate at or above this value */ 2242490Sjkh{ 2252490Sjkh register char *q; /* sieve spot */ 2262490Sjkh register ubig factor; /* index and factor */ 2272490Sjkh register char *tab_lim; /* the limit to sieve on the table */ 2282490Sjkh register ubig *p; /* prime table pointer */ 2292490Sjkh register ubig fact_lim; /* highest prime for current block */ 2302490Sjkh 2312490Sjkh /* 2322490Sjkh * A number of systems can not convert double values into unsigned 2332490Sjkh * longs when the values are larger than the largest signed value. 2342490Sjkh * We don't have this problem, so we can go all the way to BIG. 2352490Sjkh */ 2362490Sjkh if (start < 3) { 2372490Sjkh start = (ubig)2; 2382490Sjkh } 2392490Sjkh if (stop < 3) { 2402490Sjkh stop = (ubig)2; 2412490Sjkh } 2422490Sjkh if (stop <= start) { 2432490Sjkh return; 2442490Sjkh } 2452490Sjkh 2462490Sjkh /* 2472490Sjkh * be sure that the values are odd, or 2 2482490Sjkh */ 2492490Sjkh if (start != 2 && (start&0x1) == 0) { 2502490Sjkh ++start; 2512490Sjkh } 2522490Sjkh if (stop != 2 && (stop&0x1) == 0) { 2532490Sjkh ++stop; 2542490Sjkh } 2552490Sjkh 2562490Sjkh /* 2572490Sjkh * quick list of primes <= pr_limit 2582490Sjkh */ 2592490Sjkh if (start <= *pr_limit) { 2602490Sjkh /* skip primes up to the start value */ 2612490Sjkh for (p = &prime[0], factor = prime[0]; 2622490Sjkh factor < stop && p <= pr_limit; factor = *(++p)) { 2632490Sjkh if (factor >= start) { 26442357Simp printf(hflag ? "0x%lx\n" : "%lu\n", factor); 2652490Sjkh } 2662490Sjkh } 2672490Sjkh /* return early if we are done */ 2682490Sjkh if (p <= pr_limit) { 2692490Sjkh return; 2702490Sjkh } 2712490Sjkh start = *pr_limit+2; 2722490Sjkh } 2732490Sjkh 2742490Sjkh /* 2752490Sjkh * we shall sieve a bytemap window, note primes and move the window 2762490Sjkh * upward until we pass the stop point 2772490Sjkh */ 2782490Sjkh while (start < stop) { 2792490Sjkh /* 2802490Sjkh * factor out 3, 5, 7, 11 and 13 2812490Sjkh */ 2822490Sjkh /* initial pattern copy */ 2832490Sjkh factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 2842490Sjkh memcpy(table, &pattern[factor], pattern_size-factor); 2852490Sjkh /* main block pattern copies */ 2862490Sjkh for (fact_lim=pattern_size-factor; 2872490Sjkh fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { 2882490Sjkh memcpy(&table[fact_lim], pattern, pattern_size); 2892490Sjkh } 2902490Sjkh /* final block pattern copy */ 2912490Sjkh memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 2922490Sjkh 2932490Sjkh /* 2942490Sjkh * sieve for primes 17 and higher 2952490Sjkh */ 2962490Sjkh /* note highest useful factor and sieve spot */ 2972490Sjkh if (stop-start > TABSIZE+TABSIZE) { 2982490Sjkh tab_lim = &table[TABSIZE]; /* sieve it all */ 2992490Sjkh fact_lim = (int)sqrt( 3002490Sjkh (double)(start)+TABSIZE+TABSIZE+1.0); 3012490Sjkh } else { 3022490Sjkh tab_lim = &table[(stop-start)/2]; /* partial sieve */ 3032490Sjkh fact_lim = (int)sqrt((double)(stop)+1.0); 3042490Sjkh } 3052490Sjkh /* sieve for factors >= 17 */ 3062490Sjkh factor = 17; /* 17 is first prime to use */ 3072490Sjkh p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 3082490Sjkh do { 3092490Sjkh /* determine the factor's initial sieve point */ 3102490Sjkh q = (char *)(start%factor); /* temp storage for mod */ 31135892Sjb if ((long)q & 0x1) { 31235892Sjb q = &table[(factor-(long)q)/2]; 3132490Sjkh } else { 31435892Sjb q = &table[q ? factor-((long)q/2) : 0]; 3152490Sjkh } 3162490Sjkh /* sive for our current factor */ 3172490Sjkh for ( ; q < tab_lim; q += factor) { 3182490Sjkh *q = '\0'; /* sieve out a spot */ 3192490Sjkh } 3202490Sjkh } while ((factor=(ubig)(*(p++))) <= fact_lim); 3212490Sjkh 3222490Sjkh /* 3232490Sjkh * print generated primes 3242490Sjkh */ 3252490Sjkh for (q = table; q < tab_lim; ++q, start+=2) { 3262490Sjkh if (*q) { 32742357Simp printf(hflag ? "0x%lx\n" : "%lu\n", start); 3282490Sjkh } 3292490Sjkh } 3302490Sjkh } 3312490Sjkh} 3322490Sjkh 3332490Sjkhvoid 3342490Sjkhusage() 3352490Sjkh{ 33642338Simp (void)fprintf(stderr, "usage: primes [-h] [start [stop]]\n"); 3372490Sjkh exit(1); 3382490Sjkh} 339