primes.c revision 104728
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 3853920Sbillfstatic const 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 4453920Sbillf#if 0 4523726Speterstatic char sccsid[] = "@(#)primes.c 8.5 (Berkeley) 5/10/95"; 4653920Sbillf#endif 4753920Sbillfstatic const char rcsid[] = 4853920Sbillf "$FreeBSD: head/games/primes/primes.c 104728 2002-10-09 20:42:40Z fanf $"; 492490Sjkh#endif /* not lint */ 502490Sjkh 512490Sjkh/* 522490Sjkh * primes - generate a table of primes between two values 532490Sjkh * 542490Sjkh * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo 552490Sjkh * 562490Sjkh * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\ 572490Sjkh * 582490Sjkh * usage: 59104720Sfanf * primes [-h] [start [stop]] 602490Sjkh * 612490Sjkh * Print primes >= start and < stop. If stop is omitted, 622490Sjkh * the value 4294967295 (2^32-1) is assumed. If start is 632490Sjkh * omitted, start is read from standard input. 642490Sjkh * 652490Sjkh * validation check: there are 664579 primes between 0 and 10^7 662490Sjkh */ 672490Sjkh 682490Sjkh#include <ctype.h> 692490Sjkh#include <err.h> 702490Sjkh#include <errno.h> 712490Sjkh#include <limits.h> 722490Sjkh#include <math.h> 732490Sjkh#include <stdio.h> 742490Sjkh#include <stdlib.h> 75104728Sfanf#include <string.h> 7623726Speter#include <unistd.h> 772490Sjkh 782490Sjkh#include "primes.h" 792490Sjkh 802490Sjkh/* 812490Sjkh * Eratosthenes sieve table 822490Sjkh * 832490Sjkh * We only sieve the odd numbers. The base of our sieve windows are always 842490Sjkh * odd. If the base of table is 1, table[i] represents 2*i-1. After the 8588530Sroam * sieve, table[i] == 1 if and only if 2*i-1 is prime. 862490Sjkh * 872490Sjkh * We make TABSIZE large to reduce the overhead of inner loop setup. 882490Sjkh */ 8991027Sbillfstatic char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ 902490Sjkh 9191027Sbillfstatic int hflag; 9242338Simp 9391027Sbillfstatic void primes(ubig, ubig); 9491027Sbillfstatic ubig read_num_buf(void); 9591027Sbillfstatic void usage(void); 962490Sjkh 972490Sjkhint 9891027Sbillfmain(int argc, char *argv[]) 992490Sjkh{ 1002490Sjkh ubig start; /* where to start generating */ 1012490Sjkh ubig stop; /* don't generate at or above this value */ 1022490Sjkh int ch; 1032490Sjkh char *p; 1042490Sjkh 10542338Simp while ((ch = getopt(argc, argv, "h")) != -1) 1062490Sjkh switch (ch) { 10742338Simp case 'h': 10842338Simp hflag++; 10942338Simp break; 1102490Sjkh case '?': 1112490Sjkh default: 1122490Sjkh usage(); 1132490Sjkh } 1142490Sjkh argc -= optind; 1152490Sjkh argv += optind; 1162490Sjkh 1172490Sjkh start = 0; 1182490Sjkh stop = BIG; 1192490Sjkh 1202490Sjkh /* 1212490Sjkh * Convert low and high args. Strtoul(3) sets errno to 1222490Sjkh * ERANGE if the number is too large, but, if there's 1232490Sjkh * a leading minus sign it returns the negation of the 1242490Sjkh * result of the conversion, which we'd rather disallow. 1252490Sjkh */ 1262490Sjkh switch (argc) { 1272490Sjkh case 2: 1282490Sjkh /* Start and stop supplied on the command line. */ 1292490Sjkh if (argv[0][0] == '-' || argv[1][0] == '-') 1302490Sjkh errx(1, "negative numbers aren't permitted."); 1312490Sjkh 1322490Sjkh errno = 0; 13342338Simp start = strtoul(argv[0], &p, 0); 1342490Sjkh if (errno) 1352490Sjkh err(1, "%s", argv[0]); 1362490Sjkh if (*p != '\0') 1372490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1382490Sjkh 1392490Sjkh errno = 0; 14042338Simp stop = strtoul(argv[1], &p, 0); 1412490Sjkh if (errno) 1422490Sjkh err(1, "%s", argv[1]); 1432490Sjkh if (*p != '\0') 1442490Sjkh errx(1, "%s: illegal numeric format.", argv[1]); 1452490Sjkh break; 1462490Sjkh case 1: 1472490Sjkh /* Start on the command line. */ 1482490Sjkh if (argv[0][0] == '-') 1492490Sjkh errx(1, "negative numbers aren't permitted."); 1502490Sjkh 1512490Sjkh errno = 0; 15242338Simp start = strtoul(argv[0], &p, 0); 1532490Sjkh if (errno) 1542490Sjkh err(1, "%s", argv[0]); 1552490Sjkh if (*p != '\0') 1562490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1572490Sjkh break; 1582490Sjkh case 0: 15938199Sphk start = read_num_buf(); 1602490Sjkh break; 1612490Sjkh default: 1622490Sjkh usage(); 1632490Sjkh } 1642490Sjkh 1652490Sjkh if (start > stop) 1662490Sjkh errx(1, "start value must be less than stop value."); 16738199Sphk primes(start, stop); 16891027Sbillf return (0); 1692490Sjkh} 1702490Sjkh 1712490Sjkh/* 1722490Sjkh * read_num_buf -- 1732490Sjkh * This routine returns a number n, where 0 <= n && n <= BIG. 1742490Sjkh */ 17591027Sbillfstatic ubig 17691027Sbillfread_num_buf(void) 1772490Sjkh{ 1782490Sjkh ubig val; 179104720Sfanf char *p, buf[LINE_MAX]; /* > max number of digits. */ 1802490Sjkh 1812490Sjkh for (;;) { 1822490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1832490Sjkh if (ferror(stdin)) 1842490Sjkh err(1, "stdin"); 1852490Sjkh exit(0); 1862490Sjkh } 1872490Sjkh for (p = buf; isblank(*p); ++p); 1882490Sjkh if (*p == '\n' || *p == '\0') 1892490Sjkh continue; 1902490Sjkh if (*p == '-') 1912490Sjkh errx(1, "negative numbers aren't permitted."); 1922490Sjkh errno = 0; 19342338Simp val = strtoul(buf, &p, 0); 1942490Sjkh if (errno) 1952490Sjkh err(1, "%s", buf); 1962490Sjkh if (*p != '\n') 1972490Sjkh errx(1, "%s: illegal numeric format.", buf); 1982490Sjkh return (val); 1992490Sjkh } 2002490Sjkh} 2012490Sjkh 2022490Sjkh/* 2032490Sjkh * primes - sieve and print primes from start up to and but not including stop 2042490Sjkh */ 20591027Sbillfstatic void 20691027Sbillfprimes(ubig start, ubig stop) 2072490Sjkh{ 20853210Sbillf char *q; /* sieve spot */ 20953210Sbillf ubig factor; /* index and factor */ 21053210Sbillf char *tab_lim; /* the limit to sieve on the table */ 211104720Sfanf const ubig *p; /* prime table pointer */ 21253210Sbillf ubig fact_lim; /* highest prime for current block */ 213104720Sfanf ubig mod; /* temp storage for mod */ 2142490Sjkh 2152490Sjkh /* 2162490Sjkh * A number of systems can not convert double values into unsigned 2172490Sjkh * longs when the values are larger than the largest signed value. 2182490Sjkh * We don't have this problem, so we can go all the way to BIG. 2192490Sjkh */ 2202490Sjkh if (start < 3) { 2212490Sjkh start = (ubig)2; 2222490Sjkh } 2232490Sjkh if (stop < 3) { 2242490Sjkh stop = (ubig)2; 2252490Sjkh } 2262490Sjkh if (stop <= start) { 2272490Sjkh return; 2282490Sjkh } 2292490Sjkh 2302490Sjkh /* 2312490Sjkh * be sure that the values are odd, or 2 2322490Sjkh */ 2332490Sjkh if (start != 2 && (start&0x1) == 0) { 2342490Sjkh ++start; 2352490Sjkh } 2362490Sjkh if (stop != 2 && (stop&0x1) == 0) { 2372490Sjkh ++stop; 2382490Sjkh } 2392490Sjkh 2402490Sjkh /* 2412490Sjkh * quick list of primes <= pr_limit 2422490Sjkh */ 2432490Sjkh if (start <= *pr_limit) { 2442490Sjkh /* skip primes up to the start value */ 2452490Sjkh for (p = &prime[0], factor = prime[0]; 2462490Sjkh factor < stop && p <= pr_limit; factor = *(++p)) { 2472490Sjkh if (factor >= start) { 24842357Simp printf(hflag ? "0x%lx\n" : "%lu\n", factor); 2492490Sjkh } 2502490Sjkh } 2512490Sjkh /* return early if we are done */ 2522490Sjkh if (p <= pr_limit) { 2532490Sjkh return; 2542490Sjkh } 2552490Sjkh start = *pr_limit+2; 2562490Sjkh } 2572490Sjkh 2582490Sjkh /* 2592490Sjkh * we shall sieve a bytemap window, note primes and move the window 2602490Sjkh * upward until we pass the stop point 2612490Sjkh */ 2622490Sjkh while (start < stop) { 2632490Sjkh /* 2642490Sjkh * factor out 3, 5, 7, 11 and 13 2652490Sjkh */ 2662490Sjkh /* initial pattern copy */ 2672490Sjkh factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 2682490Sjkh memcpy(table, &pattern[factor], pattern_size-factor); 2692490Sjkh /* main block pattern copies */ 2702490Sjkh for (fact_lim=pattern_size-factor; 2712490Sjkh fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { 2722490Sjkh memcpy(&table[fact_lim], pattern, pattern_size); 2732490Sjkh } 2742490Sjkh /* final block pattern copy */ 2752490Sjkh memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 2762490Sjkh 2772490Sjkh /* 2782490Sjkh * sieve for primes 17 and higher 2792490Sjkh */ 2802490Sjkh /* note highest useful factor and sieve spot */ 2812490Sjkh if (stop-start > TABSIZE+TABSIZE) { 2822490Sjkh tab_lim = &table[TABSIZE]; /* sieve it all */ 283104720Sfanf fact_lim = sqrt(start+1.0+TABSIZE+TABSIZE); 2842490Sjkh } else { 2852490Sjkh tab_lim = &table[(stop-start)/2]; /* partial sieve */ 286104720Sfanf fact_lim = sqrt(stop+1.0); 2872490Sjkh } 2882490Sjkh /* sieve for factors >= 17 */ 2892490Sjkh factor = 17; /* 17 is first prime to use */ 2902490Sjkh p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 2912490Sjkh do { 2922490Sjkh /* determine the factor's initial sieve point */ 293104720Sfanf mod = start%factor; 294104720Sfanf if (mod & 0x1) { 295104720Sfanf q = &table[(factor-mod)/2]; 2962490Sjkh } else { 297104720Sfanf q = &table[mod ? factor-(mod/2) : 0]; 2982490Sjkh } 2992490Sjkh /* sive for our current factor */ 3002490Sjkh for ( ; q < tab_lim; q += factor) { 3012490Sjkh *q = '\0'; /* sieve out a spot */ 3022490Sjkh } 303104720Sfanf factor = *p++; 304104720Sfanf } while (factor <= fact_lim); 3052490Sjkh 3062490Sjkh /* 3072490Sjkh * print generated primes 3082490Sjkh */ 3092490Sjkh for (q = table; q < tab_lim; ++q, start+=2) { 3102490Sjkh if (*q) { 31142357Simp printf(hflag ? "0x%lx\n" : "%lu\n", start); 3122490Sjkh } 3132490Sjkh } 3142490Sjkh } 3152490Sjkh} 3162490Sjkh 31791027Sbillfstatic void 31891027Sbillfusage(void) 3192490Sjkh{ 320104720Sfanf fprintf(stderr, "usage: primes [-h] [start [stop]]\n"); 3212490Sjkh exit(1); 3222490Sjkh} 323