factor/factor.c

   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 }