factor/factor.c

   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