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