factor/factor.c

   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