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