1 /* $NetBSD: factor.c,v 1.29 2018/02/06 16:53:27 christos 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 43 static char sccsid[] = "@(#)factor.c 8.4 (Berkeley) 5/4/95"; 44 #else 45 __RCSID("$NetBSD: factor.c,v 1.29 2018/02/06 16:53:27 christos Exp $"); 46 #endif 47 #endif /* not lint */ 48 49 /* 50 * factor - factor a number into primes 51 * 52 * By Landon Curt Noll, http://www.isthe.com/chongo/index.html /\oo/\ 53 * 54 * usage: 55 * factor [-h] [number] ... 56 * 57 * By Default, the form of the output is: 58 * 59 * number: factor1 factor1 factor2 factor3 factor3 factor3 ... 60 * 61 * where factor1 <= factor2 <= factor3 <= ... 62 * 63 * If the -h flag is specified, the output is in "human-readable" format. 64 * Duplicate factors are printed in the form of x^n. 65 * 66 * If no number 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 #include <inttypes.h> 77 78 #ifdef HAVE_OPENSSL 79 #include <openssl/bn.h> 80 #else 81 typedef long BIGNUM; 82 typedef u_long BN_ULONG; 83 static int BN_dec2bn(BIGNUM **a, const char *str); 84 #define BN_new() ((BIGNUM *)calloc(sizeof(BIGNUM), 1)) 85 #define BN_is_zero(v) (*(v) == 0) 86 #define BN_is_one(v) (*(v) == 1) 87 #define BN_mod_word(a, b) (*(a) % (b)) 88 #endif 89 90 #include "primes.h" 91 92 /* 93 * prime[i] is the (i-1)th prime. 94 * 95 * We are able to sieve 2^32-1 because this byte table yields all primes 96 * up to 65537 and 65537^2 > 2^32-1. 97 */ 98 99 #if 0 /* debugging: limit table use to stress the "pollard" code */ 100 #define pr_limit &prime[0] 101 #endif 102 103 #define PRIME_CHECKS 5 104 105 #ifdef HAVE_OPENSSL 106 static BN_CTX *ctx; /* just use a global context */ 107 #endif 108 109 static void pr_fact(BIGNUM *, int); /* print factors of a value */ 110 static void BN_print_dec_fp(FILE *, const BIGNUM *); 111 static void usage(void) __dead; 112 #ifdef HAVE_OPENSSL 113 static void pollard_rho(BIGNUM *); /* print factors for big numbers */ 114 #else 115 static char *BN_bn2dec(const BIGNUM *); 116 static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG); 117 #endif 118 119 120 #ifndef HAVE_OPENSSL 121 static int 122 BN_dec2bn(BIGNUM **a, const char *str) 123 { 124 char *p; 125 126 errno = 0; 127 **a = strtoul(str, &p, 10); 128 if (errno) 129 err(1, "%s", str); 130 return (*p == '\n' || *p == '\0'); 131 } 132 #endif 133 134 int 135 main(int argc, char *argv[]) 136 { 137 BIGNUM *val; 138 int ch, hflag; 139 char *p, buf[LINE_MAX]; /* > max number of digits. */ 140 141 #ifdef HAVE_OPENSSL 142 ctx = BN_CTX_new(); 143 #endif 144 val = BN_new(); 145 if (val == NULL) 146 errx(1, "can't initialise bignum"); 147 148 hflag = 0; 149 while ((ch = getopt(argc, argv, "h")) != -1) 150 switch (ch) { 151 case 'h': 152 hflag = 1; 153 break; 154 case '?': 155 default: 156 usage(); 157 } 158 argc -= optind; 159 argv += optind; 160 161 /* No args supplied, read numbers from stdin. */ 162 if (argc == 0) 163 for (;;) { 164 if (fgets(buf, sizeof(buf), stdin) == NULL) { 165 if (ferror(stdin)) 166 err(1, "stdin"); 167 exit (0); 168 } 169 for (p = buf; isblank((unsigned char)*p); ++p); 170 if (*p == '\n' || *p == '\0') 171 continue; 172 if (*p == '-') 173 errx(1, "negative numbers aren't permitted."); 174 if (BN_dec2bn(&val, buf) == 0) 175 errx(1, "%s: illegal numeric format.", argv[0]); 176 pr_fact(val, hflag); 177 } 178 /* Factor the arguments. */ 179 else 180 for (; *argv != NULL; ++argv) { 181 if (argv[0][0] == '-') 182 errx(1, "numbers <= 1 aren't permitted."); 183 if (BN_dec2bn(&val, argv[0]) == 0) 184 errx(1, "%s: illegal numeric format.", argv[0]); 185 pr_fact(val, hflag); 186 } 187 exit(0); 188 } 189 190 /* 191 * pr_fact - print the factors of a number 192 * 193 * If the number is 0 or 1, then print the number and return. 194 * If the number is < 0, print -1, negate the number and continue 195 * processing. 196 * 197 * Print the factors of the number, from the lowest to the highest. 198 * By default, a factor will be printed numtiple times if it divides 199 * the value multiple times. 200 * 201 * If hflag is specified, duplicate factors are printed in "human- 202 * readable" form of x^n. 203 * 204 * Factors are printed with leading tabs. 205 */ 206 static void 207 pr_fact(BIGNUM *val, int hflag) 208 { 209 const uint64_t *fact; /* The factor found. */ 210 int i; 211 212 /* Firewall - catch 0 and 1. */ 213 if (BN_is_zero(val) || BN_is_one(val)) 214 errx(1, "numbers <= 1 aren't permitted."); 215 216 /* Factor value. */ 217 218 BN_print_dec_fp(stdout, val); 219 putchar(':'); 220 for (fact = &prime[0]; !BN_is_one(val); ++fact) { 221 /* Look for the smallest factor. */ 222 while (fact <= pr_limit) { 223 if (BN_mod_word(val, (BN_ULONG)*fact) == 0) 224 break; 225 fact++; 226 } 227 228 /* Watch for primes larger than the table. */ 229 if (fact > pr_limit) { 230 #ifdef HAVE_OPENSSL 231 BIGNUM *bnfact; 232 233 bnfact = BN_new(); 234 BN_set_word(bnfact, (BN_ULONG)*(fact - 1)); 235 BN_sqr(bnfact, bnfact, ctx); 236 if (BN_cmp(bnfact, val) > 0 237 || BN_is_prime_ex(val, PRIME_CHECKS, NULL, NULL) 238 == 1) { 239 putchar(' '); 240 BN_print_dec_fp(stdout, val); 241 } else 242 pollard_rho(val); 243 #else 244 printf(" %s", BN_bn2dec(val)); 245 #endif 246 break; 247 } 248 249 /* Divide factor out until none are left. */ 250 i = 0; 251 do { 252 i++; 253 if (!hflag) 254 printf(" %" PRIu64, *fact); 255 BN_div_word(val, (BN_ULONG)*fact); 256 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0); 257 258 if (hflag) { 259 printf(" %" PRIu64, *fact); 260 if (i > 1) 261 printf("^%d", i); 262 } 263 264 /* Let the user know we're doing something. */ 265 fflush(stdout); 266 } 267 putchar('\n'); 268 } 269 270 /* 271 * Sigh.. No _decimal_ output to file functions in BN. 272 */ 273 static void 274 BN_print_dec_fp(FILE *fp, const BIGNUM *num) 275 { 276 char *buf; 277 278 buf = BN_bn2dec(num); 279 if (buf == NULL) 280 return; /* XXX do anything here? */ 281 fprintf(fp, "%s", buf); 282 free(buf); 283 } 284 285 void 286 usage(void) 287 { 288 fprintf(stderr, "usage: factor [-h] [value ...]\n"); 289 exit (0); 290 } 291 292 293 294 295 #ifdef HAVE_OPENSSL 296 static void 297 pollard_rho(BIGNUM *val) 298 { 299 BIGNUM *x, *y, *tmp, *num; 300 BN_ULONG a; 301 unsigned int steps_taken, steps_limit; 302 303 x = BN_new(); 304 y = BN_new(); 305 tmp = BN_new(); 306 num = BN_new(); 307 a = 1; 308 restart: 309 steps_taken = 0; 310 steps_limit = 2; 311 BN_set_word(x, 1); 312 BN_copy(y, x); 313 314 for (;;) { 315 BN_sqr(tmp, x, ctx); 316 BN_add_word(tmp, a); 317 BN_mod(x, tmp, val, ctx); 318 BN_sub(tmp, x, y); 319 if (BN_is_zero(tmp)) { 320 #ifdef DEBUG 321 printf(" (loop)"); 322 #endif 323 a++; 324 goto restart; 325 } 326 BN_gcd(tmp, tmp, val, ctx); 327 328 if (!BN_is_one(tmp)) { 329 if (BN_is_prime_ex(tmp, PRIME_CHECKS, NULL, NULL) == 1) 330 { 331 putchar(' '); 332 BN_print_dec_fp(stdout, tmp); 333 } else { 334 #ifdef DEBUG 335 printf(" (recurse for "); 336 BN_print_dec_fp(stdout, tmp); 337 putchar(')'); 338 #endif 339 pollard_rho(BN_dup(tmp)); 340 #ifdef DEBUG 341 printf(" (back)"); 342 #endif 343 } 344 fflush(stdout); 345 346 BN_div(num, NULL, val, tmp, ctx); 347 if (BN_is_one(num)) 348 return; 349 if (BN_is_prime_ex(num, PRIME_CHECKS, NULL, NULL) == 1) 350 { 351 putchar(' '); 352 BN_print_dec_fp(stdout, num); 353 fflush(stdout); 354 return; 355 } 356 BN_copy(val, num); 357 goto restart; 358 } 359 steps_taken++; 360 if (steps_taken == steps_limit) { 361 BN_copy(y, x); /* teleport the turtle */ 362 steps_taken = 0; 363 steps_limit *= 2; 364 if (steps_limit == 0) { 365 #ifdef DEBUG 366 printf(" (overflow)"); 367 #endif 368 a++; 369 goto restart; 370 } 371 } 372 } 373 } 374 #else 375 char * 376 BN_bn2dec(const BIGNUM *val) 377 { 378 char *buf; 379 380 buf = malloc(100); 381 if (!buf) 382 return buf; 383 snprintf(buf, 100, "%ld", (long)*val); 384 return buf; 385 } 386 387 static BN_ULONG 388 BN_div_word(BIGNUM *a, BN_ULONG b) 389 { 390 BN_ULONG mod; 391 392 mod = *a % b; 393 *a /= b; 394 return mod; 395 } 396 #endif 397