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@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