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: stable/11/usr.bin/primes/primes.c 320218 2017-06-22 05:26:08Z cperciva $"; 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> 67272207Scperciva#include <inttypes.h> 682490Sjkh#include <limits.h> 692490Sjkh#include <math.h> 702490Sjkh#include <stdio.h> 712490Sjkh#include <stdlib.h> 72104728Sfanf#include <string.h> 7323726Speter#include <unistd.h> 742490Sjkh 752490Sjkh#include "primes.h" 762490Sjkh 772490Sjkh/* 782490Sjkh * Eratosthenes sieve table 792490Sjkh * 802490Sjkh * We only sieve the odd numbers. The base of our sieve windows are always 812490Sjkh * odd. If the base of table is 1, table[i] represents 2*i-1. After the 8288530Sroam * sieve, table[i] == 1 if and only if 2*i-1 is prime. 832490Sjkh * 842490Sjkh * We make TABSIZE large to reduce the overhead of inner loop setup. 852490Sjkh */ 8691027Sbillfstatic char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ 872490Sjkh 8891027Sbillfstatic int hflag; 8942338Simp 9091027Sbillfstatic void primes(ubig, ubig); 9191027Sbillfstatic ubig read_num_buf(void); 9291027Sbillfstatic void usage(void); 932490Sjkh 942490Sjkhint 9591027Sbillfmain(int argc, char *argv[]) 962490Sjkh{ 972490Sjkh ubig start; /* where to start generating */ 982490Sjkh ubig stop; /* don't generate at or above this value */ 992490Sjkh int ch; 1002490Sjkh char *p; 1012490Sjkh 10242338Simp while ((ch = getopt(argc, argv, "h")) != -1) 1032490Sjkh switch (ch) { 10442338Simp case 'h': 10542338Simp hflag++; 10642338Simp break; 1072490Sjkh case '?': 1082490Sjkh default: 1092490Sjkh usage(); 1102490Sjkh } 1112490Sjkh argc -= optind; 1122490Sjkh argv += optind; 1132490Sjkh 1142490Sjkh start = 0; 115320218Scperciva stop = (uint64_t)(-1); 1162490Sjkh 1172490Sjkh /* 118272207Scperciva * Convert low and high args. Strtoumax(3) sets errno to 1192490Sjkh * ERANGE if the number is too large, but, if there's 1202490Sjkh * a leading minus sign it returns the negation of the 1212490Sjkh * result of the conversion, which we'd rather disallow. 1222490Sjkh */ 1232490Sjkh switch (argc) { 1242490Sjkh case 2: 1252490Sjkh /* Start and stop supplied on the command line. */ 1262490Sjkh if (argv[0][0] == '-' || argv[1][0] == '-') 1272490Sjkh errx(1, "negative numbers aren't permitted."); 1282490Sjkh 1292490Sjkh errno = 0; 130272207Scperciva start = strtoumax(argv[0], &p, 0); 1312490Sjkh if (errno) 1322490Sjkh err(1, "%s", argv[0]); 1332490Sjkh if (*p != '\0') 1342490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1352490Sjkh 1362490Sjkh errno = 0; 137272207Scperciva stop = strtoumax(argv[1], &p, 0); 1382490Sjkh if (errno) 1392490Sjkh err(1, "%s", argv[1]); 1402490Sjkh if (*p != '\0') 1412490Sjkh errx(1, "%s: illegal numeric format.", argv[1]); 1422490Sjkh break; 1432490Sjkh case 1: 1442490Sjkh /* Start on the command line. */ 1452490Sjkh if (argv[0][0] == '-') 1462490Sjkh errx(1, "negative numbers aren't permitted."); 1472490Sjkh 1482490Sjkh errno = 0; 149272207Scperciva start = strtoumax(argv[0], &p, 0); 1502490Sjkh if (errno) 1512490Sjkh err(1, "%s", argv[0]); 1522490Sjkh if (*p != '\0') 1532490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1542490Sjkh break; 1552490Sjkh case 0: 15638199Sphk start = read_num_buf(); 1572490Sjkh break; 1582490Sjkh default: 1592490Sjkh usage(); 1602490Sjkh } 1612490Sjkh 1622490Sjkh if (start > stop) 1632490Sjkh errx(1, "start value must be less than stop value."); 16438199Sphk primes(start, stop); 16591027Sbillf return (0); 1662490Sjkh} 1672490Sjkh 1682490Sjkh/* 1692490Sjkh * read_num_buf -- 1702490Sjkh * This routine returns a number n, where 0 <= n && n <= BIG. 1712490Sjkh */ 17291027Sbillfstatic ubig 17391027Sbillfread_num_buf(void) 1742490Sjkh{ 1752490Sjkh ubig val; 176104720Sfanf char *p, buf[LINE_MAX]; /* > max number of digits. */ 1772490Sjkh 1782490Sjkh for (;;) { 1792490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1802490Sjkh if (ferror(stdin)) 1812490Sjkh err(1, "stdin"); 1822490Sjkh exit(0); 1832490Sjkh } 1842490Sjkh for (p = buf; isblank(*p); ++p); 1852490Sjkh if (*p == '\n' || *p == '\0') 1862490Sjkh continue; 1872490Sjkh if (*p == '-') 1882490Sjkh errx(1, "negative numbers aren't permitted."); 1892490Sjkh errno = 0; 190272207Scperciva val = strtoumax(buf, &p, 0); 1912490Sjkh if (errno) 1922490Sjkh err(1, "%s", buf); 1932490Sjkh if (*p != '\n') 1942490Sjkh errx(1, "%s: illegal numeric format.", buf); 1952490Sjkh return (val); 1962490Sjkh } 1972490Sjkh} 1982490Sjkh 1992490Sjkh/* 2002490Sjkh * primes - sieve and print primes from start up to and but not including stop 2012490Sjkh */ 20291027Sbillfstatic void 20391027Sbillfprimes(ubig start, ubig stop) 2042490Sjkh{ 20553210Sbillf char *q; /* sieve spot */ 20653210Sbillf ubig factor; /* index and factor */ 20753210Sbillf char *tab_lim; /* the limit to sieve on the table */ 208104720Sfanf const ubig *p; /* prime table pointer */ 20953210Sbillf ubig fact_lim; /* highest prime for current block */ 210104720Sfanf ubig mod; /* temp storage for mod */ 2112490Sjkh 2122490Sjkh /* 2132490Sjkh * A number of systems can not convert double values into unsigned 2142490Sjkh * longs when the values are larger than the largest signed value. 2152490Sjkh * We don't have this problem, so we can go all the way to BIG. 2162490Sjkh */ 2172490Sjkh if (start < 3) { 2182490Sjkh start = (ubig)2; 2192490Sjkh } 2202490Sjkh if (stop < 3) { 2212490Sjkh stop = (ubig)2; 2222490Sjkh } 2232490Sjkh if (stop <= start) { 2242490Sjkh return; 2252490Sjkh } 2262490Sjkh 2272490Sjkh /* 2282490Sjkh * be sure that the values are odd, or 2 2292490Sjkh */ 2302490Sjkh if (start != 2 && (start&0x1) == 0) { 2312490Sjkh ++start; 2322490Sjkh } 2332490Sjkh if (stop != 2 && (stop&0x1) == 0) { 2342490Sjkh ++stop; 2352490Sjkh } 2362490Sjkh 2372490Sjkh /* 2382490Sjkh * quick list of primes <= pr_limit 2392490Sjkh */ 2402490Sjkh if (start <= *pr_limit) { 2412490Sjkh /* skip primes up to the start value */ 2422490Sjkh for (p = &prime[0], factor = prime[0]; 2432490Sjkh factor < stop && p <= pr_limit; factor = *(++p)) { 2442490Sjkh if (factor >= start) { 245272207Scperciva printf(hflag ? "%" PRIx64 "\n" : "%" PRIu64 "\n", factor); 2462490Sjkh } 2472490Sjkh } 2482490Sjkh /* return early if we are done */ 2492490Sjkh if (p <= pr_limit) { 2502490Sjkh return; 2512490Sjkh } 2522490Sjkh start = *pr_limit+2; 2532490Sjkh } 2542490Sjkh 2552490Sjkh /* 2562490Sjkh * we shall sieve a bytemap window, note primes and move the window 2572490Sjkh * upward until we pass the stop point 2582490Sjkh */ 2592490Sjkh while (start < stop) { 2602490Sjkh /* 2612490Sjkh * factor out 3, 5, 7, 11 and 13 2622490Sjkh */ 2632490Sjkh /* initial pattern copy */ 2642490Sjkh factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 2652490Sjkh memcpy(table, &pattern[factor], pattern_size-factor); 2662490Sjkh /* main block pattern copies */ 2672490Sjkh for (fact_lim=pattern_size-factor; 2682490Sjkh fact_lim+pattern_size<=TABSIZE; fact_lim+=pattern_size) { 2692490Sjkh memcpy(&table[fact_lim], pattern, pattern_size); 2702490Sjkh } 2712490Sjkh /* final block pattern copy */ 2722490Sjkh memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 2732490Sjkh 2742490Sjkh /* 2752490Sjkh * sieve for primes 17 and higher 2762490Sjkh */ 2772490Sjkh /* note highest useful factor and sieve spot */ 2782490Sjkh if (stop-start > TABSIZE+TABSIZE) { 2792490Sjkh tab_lim = &table[TABSIZE]; /* sieve it all */ 280104720Sfanf fact_lim = sqrt(start+1.0+TABSIZE+TABSIZE); 2812490Sjkh } else { 2822490Sjkh tab_lim = &table[(stop-start)/2]; /* partial sieve */ 283104720Sfanf fact_lim = sqrt(stop+1.0); 2842490Sjkh } 2852490Sjkh /* sieve for factors >= 17 */ 2862490Sjkh factor = 17; /* 17 is first prime to use */ 2872490Sjkh p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 2882490Sjkh do { 2892490Sjkh /* determine the factor's initial sieve point */ 290104720Sfanf mod = start%factor; 291104720Sfanf if (mod & 0x1) { 292104720Sfanf q = &table[(factor-mod)/2]; 2932490Sjkh } else { 294104720Sfanf q = &table[mod ? factor-(mod/2) : 0]; 2952490Sjkh } 2962490Sjkh /* sive for our current factor */ 2972490Sjkh for ( ; q < tab_lim; q += factor) { 2982490Sjkh *q = '\0'; /* sieve out a spot */ 2992490Sjkh } 300104720Sfanf factor = *p++; 301104720Sfanf } while (factor <= fact_lim); 3022490Sjkh 3032490Sjkh /* 3042490Sjkh * print generated primes 3052490Sjkh */ 3062490Sjkh for (q = table; q < tab_lim; ++q, start+=2) { 3072490Sjkh if (*q) { 308272207Scperciva if (start > SIEVEMAX) { 309272166Scperciva if (!isprime(start)) 310272166Scperciva continue; 311272166Scperciva } 312272207Scperciva printf(hflag ? "%" PRIx64 "\n" : "%" PRIu64 "\n", start); 3132490Sjkh } 3142490Sjkh } 3152490Sjkh } 3162490Sjkh} 3172490Sjkh 31891027Sbillfstatic void 31991027Sbillfusage(void) 3202490Sjkh{ 321104720Sfanf fprintf(stderr, "usage: primes [-h] [start [stop]]\n"); 3222490Sjkh exit(1); 3232490Sjkh} 324