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