Home | History | Annotate | Line # | Download | only in factor 
      1 /*	$NetBSD: factor.c,v 1.38 2020/10/12 13:54:51 christos Exp $	*/
      2 /*
      3  * Copyright (c) 1989, 1993
      4  *	The Regents of the University of California.  All rights reserved.
      5  *
      6  * This code is derived from software contributed to Berkeley by
      7  * Landon Curt Noll.
      8  *
      9  * Redistribution and use in source and binary forms, with or without
     10  * modification, are permitted provided that the following conditions
     11  * are met:
     12  * 1. Redistributions of source code must retain the above copyright
     13  *    notice, this list of conditions and the following disclaimer.
     14  * 2. Redistributions in binary form must reproduce the above copyright
     15  *    notice, this list of conditions and the following disclaimer in the
     16  *    documentation and/or other materials provided with the distribution.
     17  * 3. Neither the name of the University nor the names of its contributors
     18  *    may be used to endorse or promote products derived from this software
     19  *    without specific prior written permission.
     20  *
     21  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
     22  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
     25  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     26  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     27  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     28  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     29  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     30  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     31  * SUCH DAMAGE.
     32  */
     33 
     34 #ifndef lint
     35 #include <sys/cdefs.h>
     36 #ifdef __COPYRIGHT
     37 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
     38 	The Regents of the University of California.  All rights reserved.");
     39 #endif
     40 #ifdef __SCCSID
     41 __SCCSID("@(#)factor.c	8.4 (Berkeley) 5/4/95");
     42 #endif
     43 #ifdef __RCSID
     44 __RCSID("$NetBSD: factor.c,v 1.38 2020/10/12 13:54:51 christos Exp $");
     45 #endif
     46 #ifdef __FBSDID
     47 __FBSDID("$FreeBSD: head/usr.bin/factor/factor.c 356666 2020-01-12 20:25:11Z gad $");
     48 #endif
     49 #endif /* not lint */
     50 
     51 /*
     52  * factor - factor a number into primes
     53  *
     54  * By: Landon Curt Noll   chongo (at) toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
     55  *
     56  *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
     57  *
     58  * usage:
     59  *	factor [-hx] [number] ...
     60  *
     61  * The form of the output is:
     62  *
     63  *	number: factor1 factor1 factor2 factor3 factor3 factor3 ...
     64  *
     65  * where factor1 <= factor2 <= factor3 <= ...
     66  *
     67  * If the -h flag is specified, the output is in "human-readable" format.
     68  * Duplicate factors are printed in the form of x^n.
     69  *
     70  * If the -x flag is specified numbers are printed in hex.
     71  *
     72  * If no number args are given, the list of numbers are read from stdin.
     73  */
     74 
     75 #include <ctype.h>
     76 #include <err.h>
     77 #include <errno.h>
     78 #include <inttypes.h>
     79 #include <limits.h>
     80 #include <stdbool.h>
     81 #include <stdio.h>
     82 #include <stdlib.h>
     83 #include <unistd.h>
     84 
     85 #include "primes.h"
     86 
     87 #ifdef HAVE_OPENSSL
     88 
     89 #include <openssl/bn.h>
     90 
     91 #define	PRIME_CHECKS	5
     92 
     93 /* print factors for big numbers */
     94 static void	pollard_pminus1(BIGNUM *, int, int);
     95 
     96 #else
     97 
     98 typedef long	BIGNUM;
     99 typedef u_long	BN_ULONG;
    100 
    101 #define BN_CTX			int
    102 #define BN_CTX_new()		NULL
    103 #define BN_new()		calloc(sizeof(BIGNUM), 1)
    104 #define BN_free(a)		free(a)
    105 
    106 #define BN_copy(a, b)		*(a) = *(b)
    107 #define BN_cmp(a, b)		(*(a) - *(b))
    108 #define BN_is_zero(v)		(*(v) == 0)
    109 #define BN_is_one(v)		(*(v) == 1)
    110 #define BN_mod_word(a, b)	(*(a) % (b))
    111 
    112 static BIGNUM  *BN_dup(const BIGNUM *);
    113 static int	BN_dec2bn(BIGNUM **, const char *);
    114 static int	BN_hex2bn(BIGNUM **, const char *);
    115 static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
    116 static void	BN_print_fp(FILE *, const BIGNUM *);
    117 
    118 #endif
    119 
    120 static void	BN_print_dec_fp(FILE *, const BIGNUM *);
    121 static void	convert_str2bn(BIGNUM **, char *);
    122 static void	pr_fact(BIGNUM *, int, int);	/* print factors of a value */
    123 static void	pr_print(BIGNUM *, int, int);	/* print a prime */
    124 static void	usage(void) __dead;
    125 
    126 static BN_CTX	*ctx;			/* just use a global context */
    127 
    128 int
    129 main(int argc, char *argv[])
    130 {
    131 	BIGNUM *val;
    132 	int ch, hflag, xflag;
    133 	char *p, buf[LINE_MAX];		/* > max number of digits. */
    134 
    135 	hflag = xflag = 0;
    136 	ctx = BN_CTX_new();
    137 	val = BN_new();
    138 	if (val == NULL)
    139 		errx(1, "can't initialise bignum");
    140 
    141 	while ((ch = getopt(argc, argv, "hx")) != -1)
    142 		switch (ch) {
    143 		case 'h':
    144 			hflag++;
    145 			break;
    146 		case 'x':
    147 			xflag++;
    148 			break;
    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 			convert_str2bn(&val, p);
    168 			pr_fact(val, hflag, xflag);
    169 		}
    170 	/* Factor the arguments. */
    171 	else
    172 		for (p = *argv; p != NULL; p = *++argv) {
    173 			convert_str2bn(&val, p);
    174 			pr_fact(val, hflag, xflag);
    175 		}
    176 	exit(0);
    177 }
    178 
    179 static void
    180 pr_exp(int i, int xflag)
    181 {
    182 	printf(xflag && i > 9 ? "^0x%x" : "^%d", i);
    183 }
    184 
    185 static void
    186 pr_uint64(uint64_t i, int xflag)
    187 {
    188 	printf(xflag ? " 0x%" PRIx64 : " %" PRIu64, i);
    189 }
    190 
    191 /*
    192  * pr_fact - print the factors of a number
    193  *
    194  * Print the factors of the number, from the lowest to the highest.
    195  * A factor will be printed multiple times if it divides the value
    196  * multiple times.
    197  *
    198  * If hflag is specified, duplicate factors are printed in "human-
    199  * readable" form of x^n.
    200  *
    201  * If the xflag is specified, numbers will be printed in hex.
    202  *
    203  * Factors are printed with leading tabs.
    204  */
    205 static void
    206 pr_fact(BIGNUM *val, int hflag, int xflag)
    207 {
    208 	int i;
    209 	const uint64_t *fact;	/* The factor found. */
    210 
    211 	/* Firewall - catch 0 and 1. */
    212 	if (BN_is_zero(val))	/* Historical practice; 0 just exits. */
    213 		exit(0);
    214 	if (BN_is_one(val)) {
    215 		printf("1: 1\n");
    216 		return;
    217 	}
    218 
    219 	/* Factor value. */
    220 
    221 	if (xflag) {
    222 		fputs("0x", stdout);
    223 		BN_print_fp(stdout, val);
    224 	} else
    225 		BN_print_dec_fp(stdout, val);
    226 	putchar(':');
    227 	fflush(stdout);
    228 	for (fact = &prime[0]; !BN_is_one(val); ++fact) {
    229 		/* Look for the smallest factor. */
    230 		do {
    231 			if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
    232 				break;
    233 		} while (++fact <= pr_limit);
    234 
    235 		/* Watch for primes larger than the table. */
    236 		if (fact > pr_limit) {
    237 #ifdef HAVE_OPENSSL
    238 			BIGNUM *bnfact;
    239 
    240 			bnfact = BN_new();
    241 			BN_set_word(bnfact, *(fact - 1));
    242 			if (!BN_sqr(bnfact, bnfact, ctx))
    243 				errx(1, "error in BN_sqr()");
    244 			if (BN_cmp(bnfact, val) > 0 ||
    245 			    BN_is_prime_ex(val, PRIME_CHECKS, NULL, NULL) == 1)
    246 				pr_print(val, hflag, xflag);
    247 			else
    248 				pollard_pminus1(val, hflag, xflag);
    249 #else
    250 			pr_print(val, hflag, xflag);
    251 #endif
    252 			pr_print(NULL, hflag, xflag);
    253 			break;
    254 		}
    255 
    256 		/* Divide factor out until none are left. */
    257 		i = 0;
    258 		do {
    259 			i++;
    260 			if (!hflag)
    261 				pr_uint64(*fact, xflag);
    262 			BN_div_word(val, (BN_ULONG)*fact);
    263 		} while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
    264 
    265 		if (hflag) {
    266 			pr_uint64(*fact, xflag);
    267 			if (i > 1)
    268 				pr_exp(i, xflag);
    269 		}
    270 
    271 		/* Let the user know we're doing something. */
    272 		fflush(stdout);
    273 	}
    274 	putchar('\n');
    275 }
    276 
    277 static void
    278 pr_print(BIGNUM *val, int hflag, int xflag)
    279 {
    280 	static BIGNUM *sval;
    281 	static int ex = 1;
    282 	BIGNUM *pval;
    283 
    284 	if (hflag) {
    285 		if (sval == NULL) {
    286 			sval = BN_dup(val);
    287 			return;
    288 		}
    289 
    290 		if (val != NULL && BN_cmp(val, sval) == 0) {
    291 			ex++;
    292 			return;
    293 		}
    294 		pval = sval;
    295 	} else if (val == NULL) {
    296 		return;
    297 	} else {
    298 		pval = val;
    299 	}
    300 
    301 	if (xflag) {
    302 		fputs(" 0x", stdout);
    303 		BN_print_fp(stdout, pval);
    304 	} else {
    305 		putchar(' ');
    306 		BN_print_dec_fp(stdout, pval);
    307 	}
    308 
    309 	if (hflag) {
    310 		if (ex > 1)
    311 			pr_exp(ex, xflag);
    312 
    313 		if (val != NULL) {
    314 			BN_copy(sval, val);
    315 		} else {
    316 			BN_free(sval);
    317 			sval = NULL;
    318 		}
    319 		ex = 1;
    320 	}
    321 }
    322 
    323 static void
    324 usage(void)
    325 {
    326 	fprintf(stderr, "Usage: %s [-hx] [value ...]\n", getprogname());
    327 	exit(1);
    328 }
    329 
    330 #ifdef HAVE_OPENSSL
    331 
    332 /* pollard p-1, algorithm from Jim Gillogly, May 2000 */
    333 static void
    334 pollard_pminus1(BIGNUM *val, int hflag, int xflag)
    335 {
    336 	BIGNUM *base, *rbase, *num, *i, *x;
    337 
    338 	base = BN_new();
    339 	rbase = BN_new();
    340 	num = BN_new();
    341 	i = BN_new();
    342 	x = BN_new();
    343 
    344 	BN_set_word(rbase, 1);
    345 newbase:
    346 	if (!BN_add_word(rbase, 1))
    347 		errx(1, "error in BN_add_word()");
    348 	BN_set_word(i, 2);
    349 	BN_copy(base, rbase);
    350 
    351 	for (;;) {
    352 		BN_mod_exp(base, base, i, val, ctx);
    353 		if (BN_is_one(base))
    354 			goto newbase;
    355 
    356 		BN_copy(x, base);
    357 		BN_sub_word(x, 1);
    358 		if (!BN_gcd(x, x, val, ctx))
    359 			errx(1, "error in BN_gcd()");
    360 
    361 		if (!BN_is_one(x)) {
    362 			if (BN_is_prime_ex(x, PRIME_CHECKS, NULL, NULL) == 1)
    363 				pr_print(x, hflag, xflag);
    364 			else
    365 				pollard_pminus1(x, hflag, xflag);
    366 			fflush(stdout);
    367 
    368 			BN_div(num, NULL, val, x, ctx);
    369 			if (BN_is_one(num))
    370 				return;
    371 			if (BN_is_prime_ex(num, PRIME_CHECKS, NULL,
    372 			    NULL) == 1) {
    373 				pr_print(num, hflag, xflag);
    374 				fflush(stdout);
    375 				return;
    376 			}
    377 			BN_copy(val, num);
    378 		}
    379 		if (!BN_add_word(i, 1))
    380 			errx(1, "error in BN_add_word()");
    381 	}
    382 }
    383 
    384 /*
    385  * Sigh..  No _decimal_ output to file functions in BN.
    386  */
    387 static void
    388 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
    389 {
    390 	char *buf;
    391 
    392 	buf = BN_bn2dec(num);
    393 	if (buf == NULL)
    394 		return;	/* XXX do anything here? */
    395 	fprintf(fp, "%s", buf);
    396 	free(buf);
    397 }
    398 
    399 #else
    400 
    401 static void
    402 BN_print_fp(FILE *fp, const BIGNUM *num)
    403 {
    404 	fprintf(fp, "%lx", (unsigned long)*num);
    405 }
    406 
    407 static void
    408 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
    409 {
    410 	fprintf(fp, "%lu", (unsigned long)*num);
    411 }
    412 
    413 static int
    414 BN_dec2bn(BIGNUM **a, const char *str)
    415 {
    416 	char *p;
    417 
    418 	errno = 0;
    419 	**a = strtoul(str, &p, 10);
    420 	return (errno == 0 ? 1 : 0);	/* OpenSSL returns 0 on error! */
    421 }
    422 
    423 static int
    424 BN_hex2bn(BIGNUM **a, const char *str)
    425 {
    426 	char *p;
    427 
    428 	errno = 0;
    429 	**a = strtoul(str, &p, 16);
    430 	return (errno == 0 ? 1 : 0);	/* OpenSSL returns 0 on error! */
    431 }
    432 
    433 static BN_ULONG
    434 BN_div_word(BIGNUM *a, BN_ULONG b)
    435 {
    436 	BN_ULONG mod;
    437 
    438 	mod = *a % b;
    439 	*a /= b;
    440 	return mod;
    441 }
    442 
    443 static BIGNUM *
    444 BN_dup(const BIGNUM *a)
    445 {
    446 	BIGNUM *b = BN_new();
    447 	BN_copy(b, a);
    448 	return b;
    449 }
    450 
    451 #endif
    452 
    453 /* Convert string pointed to by *str to a bignum.  */
    454 static void
    455 convert_str2bn(BIGNUM **val, char *p)
    456 {
    457 	int n = 0;
    458 
    459 	if (*p == '+') p++;
    460 	if (*p == '-')
    461 		errx(1, "negative numbers aren't permitted.");
    462 	if (*p == '0' && (p[1] == 'x' || p[1] == 'X')) {
    463 		n = BN_hex2bn(val, p + 2);
    464 	} else {
    465 		n = BN_dec2bn(val, p);
    466 	}
    467 	if (n == 0)
    468 		errx(1, "%s: illegal numeric format.", p);
    469 }
    470