factor.c revision 1.11 
1/*	$OpenBSD: factor.c,v 1.11 2001/10/24 14:32:56 deraadt Exp $	*/
2/*	$NetBSD: factor.c,v 1.5 1995/03/23 08:28:07 cgd Exp $	*/
3
4/*
5 * Copyright (c) 1989, 1993
6 *	The Regents of the University of California.  All rights reserved.
7 *
8 * This code is derived from software contributed to Berkeley by
9 * Landon Curt Noll.
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted provided that the following conditions
13 * are met:
14 * 1. Redistributions of source code must retain the above copyright
15 *    notice, this list of conditions and the following disclaimer.
16 * 2. Redistributions in binary form must reproduce the above copyright
17 *    notice, this list of conditions and the following disclaimer in the
18 *    documentation and/or other materials provided with the distribution.
19 * 3. All advertising materials mentioning features or use of this software
20 *    must display the following acknowledgement:
21 *	This product includes software developed by the University of
22 *	California, Berkeley and its contributors.
23 * 4. Neither the name of the University nor the names of its contributors
24 *    may be used to endorse or promote products derived from this software
25 *    without specific prior written permission.
26 *
27 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
28 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
30 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
31 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
32 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
33 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
34 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
35 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
36 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
37 * SUCH DAMAGE.
38 */
39
40#ifndef lint
41static char copyright[] =
42"@(#) Copyright (c) 1989, 1993\n\
43	The Regents of the University of California.  All rights reserved.\n";
44#endif /* not lint */
45
46#ifndef lint
47#if 0
48static char sccsid[] = "@(#)factor.c	8.4 (Berkeley) 5/4/95";
49#else
50static char rcsid[] = "$OpenBSD: factor.c,v 1.11 2001/10/24 14:32:56 deraadt Exp $";
51#endif
52#endif /* not lint */
53
54/*
55 * factor - factor a number into primes
56 *
57 * By: Landon Curt Noll   chongo@toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
58 *
59 *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
60 *
61 * usage:
62 *	factor [number] ...
63 *
64 * The form of the output is:
65 *
66 *	number: factor1 factor1 factor2 factor3 factor3 factor3 ...
67 *
68 * where factor1 < factor2 < factor3 < ...
69 *
70 * If no args are given, the list of numbers are read from stdin.
71 */
72
73#include <sys/types.h>
74#include <err.h>
75#include <ctype.h>
76#include <errno.h>
77#include <limits.h>
78#include <math.h>
79#include <stdio.h>
80#include <stdlib.h>
81#include <string.h>
82#include <unistd.h>
83
84#include "primes.h"
85
86/*
87 * prime[i] is the (i+1)th prime.
88 *
89 * We are able to sieve 2^32-1 because this byte table yields all primes
90 * up to 65537 and 65537^2 > 2^32-1.
91 */
92extern const ubig prime[];
93extern const ubig *pr_limit;		/* largest prime in the prime array */
94extern const char pattern[];
95extern const int pattern_size;
96
97void	pr_fact __P((u_int64_t));	/* print factors of a value */
98void	pr_bigfact __P((u_int64_t));
99void	usage __P((void));
100
101int
102main(argc, argv)
103	int argc;
104	char *argv[];
105{
106	u_int64_t val;
107	int ch;
108	char *p, buf[100];		/* > max number of digits. */
109
110	/* revoke privs */
111	setegid(getgid());
112	setgid(getgid());
113
114	while ((ch = getopt(argc, argv, "")) != -1)
115		switch (ch) {
116		case '?':
117		default:
118			usage();
119		}
120	argc -= optind;
121	argv += optind;
122
123	/* No args supplied, read numbers from stdin. */
124	if (argc == 0)
125		for (;;) {
126			if (fgets(buf, sizeof(buf), stdin) == NULL) {
127				if (ferror(stdin))
128					err(1, "stdin");
129				exit (0);
130			}
131			if (*(p = buf + strlen(buf) - 1) == '\n')
132				*p = '\0';
133			for (p = buf; isblank(*p); ++p);
134			if (*p == '\0')
135				continue;
136			if (*p == '-')
137				errx(1, "negative numbers aren't permitted.");
138			errno = 0;
139			val = strtouq(buf, &p, 10);
140			if (errno)
141				err(1, "%s", buf);
142			for (; isblank(*p); ++p);
143			if (*p != '\0')
144				errx(1, "%s: illegal numeric format.", buf);
145			pr_fact(val);
146		}
147	/* Factor the arguments. */
148	else
149		for (; *argv != NULL; ++argv) {
150			if (argv[0][0] == '-')
151				errx(1, "negative numbers aren't permitted.");
152			errno = 0;
153			val = strtouq(argv[0], &p, 10);
154			if (errno)
155				err(1, "%s", argv[0]);
156			if (*p != '\0')
157				errx(1, "%s: illegal numeric format.", argv[0]);
158			pr_fact(val);
159		}
160	exit(0);
161}
162
163/*
164 * pr_fact - print the factors of a number
165 *
166 * If the number is 0 or 1, then print the number and return.
167 * If the number is < 0, print -1, negate the number and continue
168 * processing.
169 *
170 * Print the factors of the number, from the lowest to the highest.
171 * A factor will be printed multiple times if it divides the value
172 * multiple times.
173 *
174 * Factors are printed with leading tabs.
175 */
176void
177pr_fact(val)
178	u_int64_t val;		/* Factor this value. */
179{
180	const ubig *fact;	/* The factor found. */
181
182	/* Firewall - catch 0 and 1. */
183	if (val == 0)		/* Historical practice; 0 just exits. */
184		exit(0);
185	if (val == 1) {
186		(void)printf("1: 1\n");
187		return;
188	}
189
190	/* Factor value. */
191	(void)printf("%llu:", val);
192	fflush(stdout);
193	for (fact = &prime[0]; val > 1; ++fact) {
194		/* Look for the smallest factor. */
195		do {
196			if (val % (long)*fact == 0)
197				break;
198		} while (++fact <= pr_limit);
199
200		/* Watch for primes larger than the table. */
201		if (fact > pr_limit) {
202			if (val > BIG)
203				pr_bigfact(val);
204			else
205				(void)printf(" %llu", val);
206			break;
207		}
208
209		/* Divide factor out until none are left. */
210		do {
211			(void)printf(" %lu", (unsigned long) *fact);
212			val /= (long)*fact;
213		} while ((val % (long)*fact) == 0);
214
215		/* Let the user know we're doing something. */
216		(void)fflush(stdout);
217	}
218	(void)putchar('\n');
219}
220
221
222/* At this point, our number may have factors greater than those in primes[];
223 * however, we can generate primes up to 32 bits (see primes(6)), which is
224 * sufficient to factor a 64-bit quad.
225 */
226void
227pr_bigfact(val)
228	u_int64_t val;		/* Factor this value. */
229{
230	ubig start, stop, factor;
231	char *q;
232	const ubig *p;
233	ubig fact_lim, mod;
234	char *tab_lim;
235	char table[TABSIZE];	/* Eratosthenes sieve of odd numbers */
236
237	start = *pr_limit + 2;
238	stop  = (ubig)sqrt((double)val);
239	if ((stop & 0x1) == 0)
240		stop++;
241	/*
242	 * Following code barely modified from that in primes(6)
243	 *
244	 * we shall sieve a bytemap window, note primes and move the window
245	 * upward until we pass the stop point
246	 */
247	while (start < stop) {
248		/*
249		 * factor out 3, 5, 7, 11 and 13
250		 */
251		/* initial pattern copy */
252		factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
253		memcpy(table, &pattern[factor], pattern_size-factor);
254		/* main block pattern copies */
255		for (fact_lim = pattern_size - factor;
256		    fact_lim + pattern_size <= TABSIZE; fact_lim += pattern_size) {
257			memcpy(&table[fact_lim], pattern, pattern_size);
258		}
259		/* final block pattern copy */
260		memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
261
262		if (stop-start > TABSIZE+TABSIZE) {
263			tab_lim = &table[TABSIZE]; /* sieve it all */
264			fact_lim = (int)sqrt(
265					(double)(start)+TABSIZE+TABSIZE+1.0);
266		} else {
267			tab_lim = &table[(stop - start)/2]; /* partial sieve */
268			fact_lim = (int)sqrt((double)(stop) + 1.0);
269		}
270		/* sieve for factors >= 17 */
271		factor = 17;	/* 17 is first prime to use */
272		p = &prime[7];	/* 19 is next prime, pi(19)=7 */
273		do {
274			/* determine the factor's initial sieve point */
275			mod = start % factor;
276			if (mod & 0x1)
277				q = &table[(factor-mod)/2];
278			else
279				q = &table[mod ? factor-(mod/2) : 0];
280			/* sieve for our current factor */
281			for ( ; q < tab_lim; q += factor) {
282				*q = '\0'; /* sieve out a spot */
283			}
284		} while ((factor=(ubig)(*(p++))) <= fact_lim);
285
286		/*
287		 * use generated primes
288		 */
289		for (q = table; q < tab_lim; ++q, start+=2) {
290			if (*q) {
291				if (val % start == 0) {
292					do {
293						(void)printf(" %lu", (unsigned long) start);
294						val /= start;
295					} while ((val % start) == 0);
296					(void)fflush(stdout);
297					stop  = (ubig)sqrt((double)val);
298					if ((stop & 0x1) == 0)
299						stop++;
300				}
301			}
302		}
303	}
304	if (val > 1)
305		printf(" %llu", val);
306}
307
308
309void
310usage()
311{
312	(void)fprintf(stderr, "usage: factor [value ...]\n");
313	exit (0);
314}
315