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. 16203932Simp * 3. Neither the name of the University nor the names of its contributors 172490Sjkh * may be used to endorse or promote products derived from this software 182490Sjkh * without specific prior written permission. 192490Sjkh * 202490Sjkh * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 212490Sjkh * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 222490Sjkh * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 232490Sjkh * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 242490Sjkh * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 252490Sjkh * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 262490Sjkh * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 272490Sjkh * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 282490Sjkh * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 292490Sjkh * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 302490Sjkh * SUCH DAMAGE. 312490Sjkh */ 322490Sjkh 332490Sjkh#ifndef lint 3453920Sbillfstatic const char copyright[] = 352490Sjkh"@(#) Copyright (c) 1989, 1993\n\ 362490Sjkh The Regents of the University of California. All rights reserved.\n"; 372490Sjkh#endif /* not lint */ 382490Sjkh 392490Sjkh#ifndef lint 4053920Sbillf#if 0 4123726Speterstatic char sccsid[] = "@(#)primes.c 8.5 (Berkeley) 5/10/95"; 4253920Sbillf#endif 4353920Sbillfstatic const char rcsid[] = 4453920Sbillf "$FreeBSD: releng/10.3/games/primes/primes.c 203932 2010-02-15 18:46:02Z imp $"; 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: 55104720Sfanf * primes [-h] [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 <stdio.h> 702490Sjkh#include <stdlib.h> 71104728Sfanf#include <string.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 8188530Sroam * sieve, table[i] == 1 if and only if 2*i-1 is prime. 822490Sjkh * 832490Sjkh * We make TABSIZE large to reduce the overhead of inner loop setup. 842490Sjkh */ 8591027Sbillfstatic char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ 862490Sjkh 8791027Sbillfstatic int hflag; 8842338Simp 8991027Sbillfstatic void primes(ubig, ubig); 9091027Sbillfstatic ubig read_num_buf(void); 9191027Sbillfstatic void usage(void); 922490Sjkh 932490Sjkhint 9491027Sbillfmain(int argc, char *argv[]) 952490Sjkh{ 962490Sjkh ubig start; /* where to start generating */ 972490Sjkh ubig stop; /* don't generate at or above this value */ 982490Sjkh int ch; 992490Sjkh char *p; 1002490Sjkh 10142338Simp while ((ch = getopt(argc, argv, "h")) != -1) 1022490Sjkh switch (ch) { 10342338Simp case 'h': 10442338Simp hflag++; 10542338Simp break; 1062490Sjkh case '?': 1072490Sjkh default: 1082490Sjkh usage(); 1092490Sjkh } 1102490Sjkh argc -= optind; 1112490Sjkh argv += optind; 1122490Sjkh 1132490Sjkh start = 0; 1142490Sjkh stop = BIG; 1152490Sjkh 1162490Sjkh /* 1172490Sjkh * Convert low and high args. Strtoul(3) sets errno to 1182490Sjkh * ERANGE if the number is too large, but, if there's 1192490Sjkh * a leading minus sign it returns the negation of the 1202490Sjkh * result of the conversion, which we'd rather disallow. 1212490Sjkh */ 1222490Sjkh switch (argc) { 1232490Sjkh case 2: 1242490Sjkh /* Start and stop supplied on the command line. */ 1252490Sjkh if (argv[0][0] == '-' || argv[1][0] == '-') 1262490Sjkh errx(1, "negative numbers aren't permitted."); 1272490Sjkh 1282490Sjkh errno = 0; 12942338Simp start = strtoul(argv[0], &p, 0); 1302490Sjkh if (errno) 1312490Sjkh err(1, "%s", argv[0]); 1322490Sjkh if (*p != '\0') 1332490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1342490Sjkh 1352490Sjkh errno = 0; 13642338Simp stop = strtoul(argv[1], &p, 0); 1372490Sjkh if (errno) 1382490Sjkh err(1, "%s", argv[1]); 1392490Sjkh if (*p != '\0') 1402490Sjkh errx(1, "%s: illegal numeric format.", argv[1]); 1412490Sjkh break; 1422490Sjkh case 1: 1432490Sjkh /* Start on the command line. */ 1442490Sjkh if (argv[0][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 break; 1542490Sjkh case 0: 15538199Sphk start = read_num_buf(); 1562490Sjkh break; 1572490Sjkh default: 1582490Sjkh usage(); 1592490Sjkh } 1602490Sjkh 1612490Sjkh if (start > stop) 1622490Sjkh errx(1, "start value must be less than stop value."); 16338199Sphk primes(start, stop); 16491027Sbillf return (0); 1652490Sjkh} 1662490Sjkh 1672490Sjkh/* 1682490Sjkh * read_num_buf -- 1692490Sjkh * This routine returns a number n, where 0 <= n && n <= BIG. 1702490Sjkh */ 17191027Sbillfstatic ubig 17291027Sbillfread_num_buf(void) 1732490Sjkh{ 1742490Sjkh ubig val; 175104720Sfanf char *p, buf[LINE_MAX]; /* > max number of digits. */ 1762490Sjkh 1772490Sjkh for (;;) { 1782490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1792490Sjkh if (ferror(stdin)) 1802490Sjkh err(1, "stdin"); 1812490Sjkh exit(0); 1822490Sjkh } 1832490Sjkh for (p = buf; isblank(*p); ++p); 1842490Sjkh if (*p == '\n' || *p == '\0') 1852490Sjkh continue; 1862490Sjkh if (*p == '-') 1872490Sjkh errx(1, "negative numbers aren't permitted."); 1882490Sjkh errno = 0; 18942338Simp val = strtoul(buf, &p, 0); 1902490Sjkh if (errno) 1912490Sjkh err(1, "%s", buf); 1922490Sjkh if (*p != '\n') 1932490Sjkh errx(1, "%s: illegal numeric format.", buf); 1942490Sjkh return (val); 1952490Sjkh } 1962490Sjkh} 1972490Sjkh 1982490Sjkh/* 1992490Sjkh * primes - sieve and print primes from start up to and but not including stop 2002490Sjkh */ 20191027Sbillfstatic void 20291027Sbillfprimes(ubig start, ubig stop) 2032490Sjkh{ 20453210Sbillf char *q; /* sieve spot */ 20553210Sbillf ubig factor; /* index and factor */ 20653210Sbillf char *tab_lim; /* the limit to sieve on the table */ 207104720Sfanf const ubig *p; /* prime table pointer */ 20853210Sbillf ubig fact_lim; /* highest prime for current block */ 209104720Sfanf ubig mod; /* temp storage for mod */ 2102490Sjkh 2112490Sjkh /* 2122490Sjkh * A number of systems can not convert double values into unsigned 2132490Sjkh * longs when the values are larger than the largest signed value. 2142490Sjkh * We don't have this problem, so we can go all the way to BIG. 2152490Sjkh */ 2162490Sjkh if (start < 3) { 2172490Sjkh start = (ubig)2; 2182490Sjkh } 2192490Sjkh if (stop < 3) { 2202490Sjkh stop = (ubig)2; 2212490Sjkh } 2222490Sjkh if (stop <= start) { 2232490Sjkh return; 2242490Sjkh } 2252490Sjkh 2262490Sjkh /* 2272490Sjkh * be sure that the values are odd, or 2 2282490Sjkh */ 2292490Sjkh if (start != 2 && (start&0x1) == 0) { 2302490Sjkh ++start; 2312490Sjkh } 2322490Sjkh if (stop != 2 && (stop&0x1) == 0) { 2332490Sjkh ++stop; 2342490Sjkh } 2352490Sjkh 2362490Sjkh /* 2372490Sjkh * quick list of primes <= pr_limit 2382490Sjkh */ 2392490Sjkh if (start <= *pr_limit) { 2402490Sjkh /* skip primes up to the start value */ 2412490Sjkh for (p = &prime[0], factor = prime[0]; 2422490Sjkh factor < stop && p <= pr_limit; factor = *(++p)) { 2432490Sjkh if (factor >= start) { 24442357Simp printf(hflag ? "0x%lx\n" : "%lu\n", factor); 2452490Sjkh } 2462490Sjkh } 2472490Sjkh /* return early if we are done */ 2482490Sjkh if (p <= pr_limit) { 2492490Sjkh return; 2502490Sjkh } 2512490Sjkh start = *pr_limit+2; 2522490Sjkh } 2532490Sjkh 2542490Sjkh /* 2552490Sjkh * we shall sieve a bytemap window, note primes and move the window 2562490Sjkh * upward until we pass the stop point 2572490Sjkh */ 2582490Sjkh while (start < stop) { 2592490Sjkh /* 2602490Sjkh * factor out 3, 5, 7, 11 and 13 2612490Sjkh */ 2622490Sjkh /* initial pattern copy */ 2632490Sjkh factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 2642490Sjkh memcpy(table, &pattern[factor], pattern_size-factor); 2652490Sjkh /* main block pattern copies */ 2662490Sjkh for (fact_lim=pattern_size-factor; 2672490Sjkh fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { 2682490Sjkh memcpy(&table[fact_lim], pattern, pattern_size); 2692490Sjkh } 2702490Sjkh /* final block pattern copy */ 2712490Sjkh memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 2722490Sjkh 2732490Sjkh /* 2742490Sjkh * sieve for primes 17 and higher 2752490Sjkh */ 2762490Sjkh /* note highest useful factor and sieve spot */ 2772490Sjkh if (stop-start > TABSIZE+TABSIZE) { 2782490Sjkh tab_lim = &table[TABSIZE]; /* sieve it all */ 279104720Sfanf fact_lim = sqrt(start+1.0+TABSIZE+TABSIZE); 2802490Sjkh } else { 2812490Sjkh tab_lim = &table[(stop-start)/2]; /* partial sieve */ 282104720Sfanf fact_lim = sqrt(stop+1.0); 2832490Sjkh } 2842490Sjkh /* sieve for factors >= 17 */ 2852490Sjkh factor = 17; /* 17 is first prime to use */ 2862490Sjkh p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 2872490Sjkh do { 2882490Sjkh /* determine the factor's initial sieve point */ 289104720Sfanf mod = start%factor; 290104720Sfanf if (mod & 0x1) { 291104720Sfanf q = &table[(factor-mod)/2]; 2922490Sjkh } else { 293104720Sfanf q = &table[mod ? factor-(mod/2) : 0]; 2942490Sjkh } 2952490Sjkh /* sive for our current factor */ 2962490Sjkh for ( ; q < tab_lim; q += factor) { 2972490Sjkh *q = '\0'; /* sieve out a spot */ 2982490Sjkh } 299104720Sfanf factor = *p++; 300104720Sfanf } while (factor <= fact_lim); 3012490Sjkh 3022490Sjkh /* 3032490Sjkh * print generated primes 3042490Sjkh */ 3052490Sjkh for (q = table; q < tab_lim; ++q, start+=2) { 3062490Sjkh if (*q) { 30742357Simp printf(hflag ? "0x%lx\n" : "%lu\n", start); 3082490Sjkh } 3092490Sjkh } 3102490Sjkh } 3112490Sjkh} 3122490Sjkh 31391027Sbillfstatic void 31491027Sbillfusage(void) 3152490Sjkh{ 316104720Sfanf fprintf(stderr, "usage: primes [-h] [start [stop]]\n"); 3172490Sjkh exit(1); 3182490Sjkh} 319