Home | History | Annotate | Line # | Download | only in factor
factor.c revision 1.30
      1 /*	$NetBSD: factor.c,v 1.30 2020/10/03 22:27:00 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.30 2020/10/03 22:27:00 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  * If no args are given, the list of numbers are read from stdin.
     74  */
     75 
     76 #include <ctype.h>
     77 #include <err.h>
     78 #include <errno.h>
     79 #include <inttypes.h>
     80 #include <limits.h>
     81 #include <stdbool.h>
     82 #include <stdio.h>
     83 #include <stdlib.h>
     84 #include <unistd.h>
     85 
     86 #include "primes.h"
     87 
     88 #ifdef HAVE_OPENSSL
     89 
     90 #include <openssl/bn.h>
     91 
     92 #define	PRIME_CHECKS	5
     93 
     94 static void	pollard_pminus1(BIGNUM *, int); /* print factors for big numbers */
     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()		((BIGNUM *)calloc(sizeof(BIGNUM), 1))
    104 #define BN_is_zero(v)		(*(v) == 0)
    105 #define BN_is_one(v)		(*(v) == 1)
    106 #define BN_mod_word(a, b)	(*(a) % (b))
    107 
    108 static int	BN_dec2bn(BIGNUM **, const char *);
    109 static int	BN_hex2bn(BIGNUM **, const char *);
    110 static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
    111 static void	BN_print_fp(FILE *, const BIGNUM *);
    112 
    113 #endif
    114 
    115 static void	BN_print_dec_fp(FILE *, const BIGNUM *);
    116 static void	convert_str2bn(BIGNUM **, char *);
    117 static bool	is_hex_str(char *);
    118 static void	pr_fact(BIGNUM *, int, int);	/* print factors of a value */
    119 static void	pr_print(BIGNUM *, int);	/* print a prime */
    120 static void	usage(void);
    121 
    122 static BN_CTX	*ctx;			/* just use a global context */
    123 
    124 int
    125 main(int argc, char *argv[])
    126 {
    127 	BIGNUM *val;
    128 	int ch, hflag, xflag;
    129 	char *p, buf[LINE_MAX];		/* > max number of digits. */
    130 
    131 	hflag = xflag = 0;
    132 	ctx = BN_CTX_new();
    133 	val = BN_new();
    134 	if (val == NULL)
    135 		errx(1, "can't initialise bignum");
    136 
    137 	while ((ch = getopt(argc, argv, "hx")) != -1)
    138 		switch (ch) {
    139 		case 'h':
    140 			hflag++;
    141 			break;
    142 		case 'x':
    143 			xflag++;
    144 			break;
    145 		case '?':
    146 		default:
    147 			usage();
    148 		}
    149 	argc -= optind;
    150 	argv += optind;
    151 
    152 	/* No args supplied, read numbers from stdin. */
    153 	if (argc == 0)
    154 		for (;;) {
    155 			if (fgets(buf, sizeof(buf), stdin) == NULL) {
    156 				if (ferror(stdin))
    157 					err(1, "stdin");
    158 				exit (0);
    159 			}
    160 			for (p = buf; isblank((unsigned char)*p); ++p);
    161 			if (*p == '\n' || *p == '\0')
    162 				continue;
    163 			convert_str2bn(&val, p);
    164 			pr_fact(val, hflag, xflag);
    165 		}
    166 	/* Factor the arguments. */
    167 	else
    168 		for (p = *argv; p != NULL; p = *++argv) {
    169 			convert_str2bn(&val, p);
    170 			pr_fact(val, hflag, xflag);
    171 		}
    172 	exit(0);
    173 }
    174 
    175 /*
    176  * pr_fact - print the factors of a number
    177  *
    178  * Print the factors of the number, from the lowest to the highest.
    179  * A factor will be printed multiple times if it divides the value
    180  * multiple times.
    181  *
    182  * If hflag is specified, duplicate factors are printed in "human-
    183  * readable" form of x^n.
    184  *
    185  * If the xflag is specified, numbers will be printed in hex.
    186  *
    187  * Factors are printed with leading tabs.
    188  */
    189 static void
    190 pr_fact(BIGNUM *val, int hflag, int xflag)
    191 {
    192 	int i;
    193 	const uint64_t *fact;	/* The factor found. */
    194 
    195 	/* Firewall - catch 0 and 1. */
    196 	if (BN_is_zero(val))	/* Historical practice; 0 just exits. */
    197 		exit(0);
    198 	if (BN_is_one(val)) {
    199 		printf("1: 1\n");
    200 		return;
    201 	}
    202 
    203 	/* Factor value. */
    204 
    205 	if (xflag) {
    206 		fputs("0x", stdout);
    207 		BN_print_fp(stdout, val);
    208 	} else
    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 		do {
    214 			if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
    215 				break;
    216 		} while (++fact <= pr_limit);
    217 
    218 		/* Watch for primes larger than the table. */
    219 		if (fact > pr_limit) {
    220 #ifdef HAVE_OPENSSL
    221 			BIGNUM *bnfact;
    222 
    223 			bnfact = BN_new();
    224 			BN_set_word(bnfact, *(fact - 1));
    225 			if (!BN_sqr(bnfact, bnfact, ctx))
    226 				errx(1, "error in BN_sqr()");
    227 			if (BN_cmp(bnfact, val) > 0 ||
    228 			    BN_is_prime_ex(val, PRIME_CHECKS, NULL, NULL) == 1)
    229 				pr_print(val, xflag);
    230 			else
    231 				pollard_pminus1(val, xflag);
    232 #else
    233 			pr_print(val, xflag);
    234 #endif
    235 			break;
    236 		}
    237 
    238 		/* Divide factor out until none are left. */
    239 		i = 0;
    240 		do {
    241 			i++;
    242 			if (!hflag)
    243 				printf(xflag ? " 0x%" PRIx64 : " %" PRIu64,
    244 				    *fact);
    245 			BN_div_word(val, (BN_ULONG)*fact);
    246 		} while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
    247 
    248 		if (hflag) {
    249 			printf(xflag ? " 0x%" PRIx64 : " %" PRIu64, *fact);
    250 			if (i > 1)
    251 				printf(xflag ? "^0x%x" : "^%d", i);
    252 		}
    253 
    254 		/* Let the user know we're doing something. */
    255 		fflush(stdout);
    256 	}
    257 	putchar('\n');
    258 }
    259 
    260 static void
    261 pr_print(BIGNUM *val, int xflag)
    262 {
    263 	if (xflag) {
    264 		fputs(" 0x", stdout);
    265 		BN_print_fp(stdout, val);
    266 	} else {
    267 		putchar(' ');
    268 		BN_print_dec_fp(stdout, val);
    269 	}
    270 }
    271 
    272 static void
    273 usage(void)
    274 {
    275 	fprintf(stderr, "usage: factor [-h] [value ...]\n");
    276 	exit(1);
    277 }
    278 
    279 #ifdef HAVE_OPENSSL
    280 
    281 /* pollard p-1, algorithm from Jim Gillogly, May 2000 */
    282 static void
    283 pollard_pminus1(BIGNUM *val, int xflag)
    284 {
    285 	BIGNUM *base, *rbase, *num, *i, *x;
    286 
    287 	base = BN_new();
    288 	rbase = BN_new();
    289 	num = BN_new();
    290 	i = BN_new();
    291 	x = BN_new();
    292 
    293 	BN_set_word(rbase, 1);
    294 newbase:
    295 	if (!BN_add_word(rbase, 1))
    296 		errx(1, "error in BN_add_word()");
    297 	BN_set_word(i, 2);
    298 	BN_copy(base, rbase);
    299 
    300 	for (;;) {
    301 		BN_mod_exp(base, base, i, val, ctx);
    302 		if (BN_is_one(base))
    303 			goto newbase;
    304 
    305 		BN_copy(x, base);
    306 		BN_sub_word(x, 1);
    307 		if (!BN_gcd(x, x, val, ctx))
    308 			errx(1, "error in BN_gcd()");
    309 
    310 		if (!BN_is_one(x)) {
    311 			if (BN_is_prime_ex(x, PRIME_CHECKS, NULL, NULL) == 1)
    312 				pr_print(x, xflag);
    313 			else
    314 				pollard_pminus1(x, xflag);
    315 			fflush(stdout);
    316 
    317 			BN_div(num, NULL, val, x, ctx);
    318 			if (BN_is_one(num))
    319 				return;
    320 			if (BN_is_prime_ex(num, PRIME_CHECKS, NULL,
    321 			    NULL) == 1) {
    322 				pr_print(num, xflag);
    323 				fflush(stdout);
    324 				return;
    325 			}
    326 			BN_copy(val, num);
    327 		}
    328 		if (!BN_add_word(i, 1))
    329 			errx(1, "error in BN_add_word()");
    330 	}
    331 }
    332 
    333 /*
    334  * Sigh..  No _decimal_ output to file functions in BN.
    335  */
    336 static void
    337 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
    338 {
    339 	char *buf;
    340 
    341 	buf = BN_bn2dec(num);
    342 	if (buf == NULL)
    343 		return;	/* XXX do anything here? */
    344 	fprintf(fp, "%s", buf);
    345 	free(buf);
    346 }
    347 
    348 #else
    349 
    350 static void
    351 BN_print_fp(FILE *fp, const BIGNUM *num)
    352 {
    353 	fprintf(fp, "%lx", (unsigned long)*num);
    354 }
    355 
    356 static void
    357 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
    358 {
    359 	fprintf(fp, "%lu", (unsigned long)*num);
    360 }
    361 
    362 static int
    363 BN_dec2bn(BIGNUM **a, const char *str)
    364 {
    365 	char *p;
    366 
    367 	errno = 0;
    368 	**a = strtoul(str, &p, 10);
    369 	return (errno == 0 ? 1 : 0);	/* OpenSSL returns 0 on error! */
    370 }
    371 
    372 static int
    373 BN_hex2bn(BIGNUM **a, const char *str)
    374 {
    375 	char *p;
    376 
    377 	errno = 0;
    378 	**a = strtoul(str, &p, 16);
    379 	return (errno == 0 ? 1 : 0);	/* OpenSSL returns 0 on error! */
    380 }
    381 
    382 static BN_ULONG
    383 BN_div_word(BIGNUM *a, BN_ULONG b)
    384 {
    385 	BN_ULONG mod;
    386 
    387 	mod = *a % b;
    388 	*a /= b;
    389 	return mod;
    390 }
    391 
    392 #endif
    393 
    394 /*
    395  * Scan the string from left-to-right to see if the longest substring
    396  * is a valid hexadecimal number.
    397  */
    398 static bool
    399 is_hex_str(char *str)
    400 {
    401 	char c, *p;
    402 	bool saw_hex = false;
    403 
    404 	for (p = str; *p; p++) {
    405 		if (isdigit((unsigned char)*p))
    406 			continue;
    407 		c = tolower((unsigned char)*p);
    408 		if (c >= 'a' && c <= 'f') {
    409 			saw_hex = true;
    410 			continue;
    411 		}
    412 		break;	/* Not a hexadecimal digit. */
    413 	}
    414 	return saw_hex;
    415 }
    416 
    417 /* Convert string pointed to by *str to a bignum.  */
    418 static void
    419 convert_str2bn(BIGNUM **val, char *p)
    420 {
    421 	int n = 0;
    422 
    423 	if (*p == '+') p++;
    424 	if (*p == '-')
    425 		errx(1, "negative numbers aren't permitted.");
    426 	if (*p == '0') {
    427 		p++;
    428 		if (*p == 'x' || *p == 'X')
    429 			n = BN_hex2bn(val, ++p);
    430 	} else {
    431 		n = is_hex_str(p) ? BN_hex2bn(val, p) : BN_dec2bn(val, p);
    432 	}
    433 	if (n == 0)
    434 		errx(1, "%s: illegal numeric format.", p);
    435 }
    436