1/*
2 * Copyright (c) 1989, 1993
3 *	The Regents of the University of California.  All rights reserved.
4 *
5 * This code is derived from software contributed to Berkeley by
6 * Landon Curt Noll.
7 *
8 * Redistribution and use in source and binary forms, with or without
9 * modification, are permitted provided that the following conditions
10 * are met:
11 * 1. Redistributions of source code must retain the above copyright
12 *    notice, this list of conditions and the following disclaimer.
13 * 2. Redistributions in binary form must reproduce the above copyright
14 *    notice, this list of conditions and the following disclaimer in the
15 *    documentation and/or other materials provided with the distribution.
16 * 3. Neither the name of the University nor the names of its contributors
17 *    may be used to endorse or promote products derived from this software
18 *    without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
21 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
24 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
25 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
26 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
27 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
28 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
29 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
30 * SUCH DAMAGE.
31 */
32
33#ifndef lint
34#include <sys/cdefs.h>
35#ifdef __COPYRIGHT
36__COPYRIGHT("@(#) Copyright (c) 1989, 1993\
37	The Regents of the University of California.  All rights reserved.");
38#endif
39#ifdef __SCCSID
40__SCCSID("@(#)factor.c	8.4 (Berkeley) 5/4/95");
41#endif
42#ifdef __RCSID
43__RCSID("$NetBSD: factor.c,v 1.19 2009/08/12 05:54:31 dholland Exp $");
44#endif
45#ifdef __FBSDID
46__FBSDID("$FreeBSD$");
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 [-h] [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 <ctype.h>
70#include <err.h>
71#include <errno.h>
72#include <limits.h>
73#include <stdio.h>
74#include <stdlib.h>
75#include <unistd.h>
76
77#include "primes.h"
78
79#ifdef HAVE_OPENSSL
80
81#include <openssl/bn.h>
82
83#define	PRIME_CHECKS	5
84
85static void	pollard_pminus1(BIGNUM *); /* print factors for big numbers */
86
87#else
88
89typedef ubig	BIGNUM;
90typedef u_long	BN_ULONG;
91
92#define BN_CTX			int
93#define BN_CTX_new()		NULL
94#define BN_new()		((BIGNUM *)calloc(sizeof(BIGNUM), 1))
95#define BN_is_zero(v)		(*(v) == 0)
96#define BN_is_one(v)		(*(v) == 1)
97#define BN_mod_word(a, b)	(*(a) % (b))
98
99static int	BN_dec2bn(BIGNUM **a, const char *str);
100static int	BN_hex2bn(BIGNUM **a, const char *str);
101static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
102static void	BN_print_fp(FILE *, const BIGNUM *);
103
104#endif
105
106static void	BN_print_dec_fp(FILE *, const BIGNUM *);
107
108static void	pr_fact(BIGNUM *);	/* print factors of a value */
109static void	pr_print(BIGNUM *);	/* print a prime */
110static void	usage(void);
111
112static BN_CTX	*ctx;			/* just use a global context */
113static int	hflag;
114
115int
116main(int argc, char *argv[])
117{
118	BIGNUM *val;
119	int ch;
120	char *p, buf[LINE_MAX];		/* > max number of digits. */
121
122	ctx = BN_CTX_new();
123	val = BN_new();
124	if (val == NULL)
125		errx(1, "can't initialise bignum");
126
127	while ((ch = getopt(argc, argv, "h")) != -1)
128		switch (ch) {
129		case 'h':
130			hflag++;
131			break;
132		case '?':
133		default:
134			usage();
135		}
136	argc -= optind;
137	argv += optind;
138
139	/* No args supplied, read numbers from stdin. */
140	if (argc == 0)
141		for (;;) {
142			if (fgets(buf, sizeof(buf), stdin) == NULL) {
143				if (ferror(stdin))
144					err(1, "stdin");
145				exit (0);
146			}
147			for (p = buf; isblank(*p); ++p);
148			if (*p == '\n' || *p == '\0')
149				continue;
150			if (*p == '-')
151				errx(1, "negative numbers aren't permitted.");
152			if (BN_dec2bn(&val, buf) == 0 &&
153			    BN_hex2bn(&val, buf) == 0)
154				errx(1, "%s: illegal numeric format.", buf);
155			pr_fact(val);
156		}
157	/* Factor the arguments. */
158	else
159		for (; *argv != NULL; ++argv) {
160			if (argv[0][0] == '-')
161				errx(1, "negative numbers aren't permitted.");
162			if (BN_dec2bn(&val, argv[0]) == 0 &&
163			    BN_hex2bn(&val, argv[0]) == 0)
164				errx(1, "%s: illegal numeric format.", argv[0]);
165			pr_fact(val);
166		}
167	exit(0);
168}
169
170/*
171 * pr_fact - print the factors of a number
172 *
173 * Print the factors of the number, from the lowest to the highest.
174 * A factor will be printed multiple times if it divides the value
175 * multiple times.
176 *
177 * Factors are printed with leading tabs.
178 */
179static void
180pr_fact(BIGNUM *val)
181{
182	const ubig *fact;	/* The factor found. */
183
184	/* Firewall - catch 0 and 1. */
185	if (BN_is_zero(val))	/* Historical practice; 0 just exits. */
186		exit(0);
187	if (BN_is_one(val)) {
188		printf("1: 1\n");
189		return;
190	}
191
192	/* Factor value. */
193
194	if (hflag) {
195		fputs("0x", stdout);
196		BN_print_fp(stdout, val);
197	} else
198		BN_print_dec_fp(stdout, val);
199	putchar(':');
200	for (fact = &prime[0]; !BN_is_one(val); ++fact) {
201		/* Look for the smallest factor. */
202		do {
203			if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
204				break;
205		} while (++fact <= pr_limit);
206
207		/* Watch for primes larger than the table. */
208		if (fact > pr_limit) {
209#ifdef HAVE_OPENSSL
210			BIGNUM *bnfact;
211
212			bnfact = BN_new();
213			BN_set_word(bnfact, *(fact - 1));
214			if (!BN_sqr(bnfact, bnfact, ctx))
215				errx(1, "error in BN_sqr()");
216			if (BN_cmp(bnfact, val) > 0 ||
217			    BN_is_prime(val, PRIME_CHECKS,
218					NULL, NULL, NULL) == 1)
219				pr_print(val);
220			else
221				pollard_pminus1(val);
222#else
223			pr_print(val);
224#endif
225			break;
226		}
227
228		/* Divide factor out until none are left. */
229		do {
230			printf(hflag ? " 0x%lx" : " %lu", *fact);
231			BN_div_word(val, (BN_ULONG)*fact);
232		} while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
233
234		/* Let the user know we're doing something. */
235		fflush(stdout);
236	}
237	putchar('\n');
238}
239
240static void
241pr_print(BIGNUM *val)
242{
243	if (hflag) {
244		fputs(" 0x", stdout);
245		BN_print_fp(stdout, val);
246	} else {
247		putchar(' ');
248		BN_print_dec_fp(stdout, val);
249	}
250}
251
252static void
253usage(void)
254{
255	fprintf(stderr, "usage: factor [-h] [value ...]\n");
256	exit(1);
257}
258
259#ifdef HAVE_OPENSSL
260
261/* pollard p-1, algorithm from Jim Gillogly, May 2000 */
262static void
263pollard_pminus1(BIGNUM *val)
264{
265	BIGNUM *base, *rbase, *num, *i, *x;
266
267	base = BN_new();
268	rbase = BN_new();
269	num = BN_new();
270	i = BN_new();
271	x = BN_new();
272
273	BN_set_word(rbase, 1);
274newbase:
275	if (!BN_add_word(rbase, 1))
276		errx(1, "error in BN_add_word()");
277	BN_set_word(i, 2);
278	BN_copy(base, rbase);
279
280	for (;;) {
281		BN_mod_exp(base, base, i, val, ctx);
282		if (BN_is_one(base))
283			goto newbase;
284
285		BN_copy(x, base);
286		BN_sub_word(x, 1);
287		if (!BN_gcd(x, x, val, ctx))
288			errx(1, "error in BN_gcd()");
289
290		if (!BN_is_one(x)) {
291			if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
292			    NULL) == 1)
293				pr_print(x);
294			else
295				pollard_pminus1(x);
296			fflush(stdout);
297
298			BN_div(num, NULL, val, x, ctx);
299			if (BN_is_one(num))
300				return;
301			if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
302			    NULL) == 1) {
303				pr_print(num);
304				fflush(stdout);
305				return;
306			}
307			BN_copy(val, num);
308		}
309		if (!BN_add_word(i, 1))
310			errx(1, "error in BN_add_word()");
311	}
312}
313
314/*
315 * Sigh..  No _decimal_ output to file functions in BN.
316 */
317static void
318BN_print_dec_fp(FILE *fp, const BIGNUM *num)
319{
320	char *buf;
321
322	buf = BN_bn2dec(num);
323	if (buf == NULL)
324		return;	/* XXX do anything here? */
325	fprintf(fp, "%s", buf);
326	free(buf);
327}
328
329#else
330
331static void
332BN_print_fp(FILE *fp, const BIGNUM *num)
333{
334	fprintf(fp, "%lx", (unsigned long)*num);
335}
336
337static void
338BN_print_dec_fp(FILE *fp, const BIGNUM *num)
339{
340	fprintf(fp, "%lu", (unsigned long)*num);
341}
342
343static int
344BN_dec2bn(BIGNUM **a, const char *str)
345{
346	char *p;
347
348	errno = 0;
349	**a = strtoul(str, &p, 10);
350	return (errno == 0 && (*p == '\n' || *p == '\0'));
351}
352
353static int
354BN_hex2bn(BIGNUM **a, const char *str)
355{
356	char *p;
357
358	errno = 0;
359	**a = strtoul(str, &p, 16);
360	return (errno == 0 && (*p == '\n' || *p == '\0'));
361}
362
363static BN_ULONG
364BN_div_word(BIGNUM *a, BN_ULONG b)
365{
366	BN_ULONG mod;
367
368	mod = *a % b;
369	*a /= b;
370	return mod;
371}
372
373#endif
374