factor/factor.c

   1 /*      $NetBSD: factor.c,v 1.17 2007/12/15 19:44:40 perry 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\n\
  38         The Regents of the University of California.  All rights reserved.\n");
  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.17 2007/12/15 19:44:40 perry Exp $");
  46 #endif
  47 #endif /* not lint */
  48
  49 /*
  50  * factor - factor a number into primes
  51  *
  52  * By: Landon Curt Noll   chongo@toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
  53  *
  54  *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
  55  *
  56  * usage:
  57  *      factor [number] ...
  58  *
  59  * The form of the output is:
  60  *
  61  *      number: factor1 factor1 factor2 factor3 factor3 factor3 ...
  62  *
  63  * where factor1 <= factor2 <= factor3 <= ...
  64  *
  65  * If no args are given, the list of numbers are read from stdin.
  66  */
  67
  68 #include <ctype.h>
  69 #include <err.h>
  70 #include <errno.h>
  71 #include <limits.h>
  72 #include <stdio.h>
  73 #include <stdlib.h>
  74 #include <unistd.h>
  75
  76 #ifdef HAVE_OPENSSL
  77 #include <openssl/bn.h>
  78 #else
  79 typedef long    BIGNUM;
  80 typedef u_long  BN_ULONG;
  81 int     BN_dec2bn(BIGNUM **a, const char *str);
  82 #define BN_new()                ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
  83 #define BN_is_zero(v)           (*(v) == 0)
  84 #define BN_is_one(v)            (*(v) == 1)
  85 #define BN_new()                ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
  86 #define BN_is_zero(v)           (*(v) == 0)
  87 #define BN_is_one(v)            (*(v) == 1)
  88 #define BN_mod_word(a, b)       (*(a) % (b))
  89 #endif
  90
  91 #include "primes.h"
  92
  93 /*
  94  * prime[i] is the (i-1)th prime.
  95  *
  96  * We are able to sieve 2^32-1 because this byte table yields all primes
  97  * up to 65537 and 65537^2 > 2^32-1.
  98  */
  99 extern const ubig prime[];
 100 extern const ubig *pr_limit;            /* largest prime in the prime array */
 101
 102 #define PRIME_CHECKS    5
 103
 104 #ifdef HAVE_OPENSSL
 105 BN_CTX *ctx;                            /* just use a global context */
 106 #endif
 107
 108 int     main(int, char *[]);
 109 void    pr_fact(BIGNUM *);              /* print factors of a value */
 110 void    BN_print_dec_fp(FILE *, const BIGNUM *);
 111 void    usage(void) __dead;
 112 #ifdef HAVE_OPENSSL
 113 void    pollard_pminus1(BIGNUM *);      /* print factors for big numbers */
 114 #else
 115 char    *BN_bn2dec(const BIGNUM *);
 116 BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
 117 #endif
 118
 119
 120 #ifndef HAVE_OPENSSL
 121 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;
 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         while ((ch = getopt(argc, argv, "")) != -1)
 149                 switch (ch) {
 150                 case '?':
 151                 default:
 152                         usage();
 153                 }
 154         argc -= optind;
 155         argv += optind;
 156
 157         /* No args supplied, read numbers from stdin. */
 158         if (argc == 0)
 159                 for (;;) {
 160                         if (fgets(buf, sizeof(buf), stdin) == NULL) {
 161                                 if (ferror(stdin))
 162                                         err(1, "stdin");
 163                                 exit (0);
 164                         }
 165                         for (p = buf; isblank(*p); ++p);
 166                         if (*p == '\n' || *p == '\0')
 167                                 continue;
 168                         if (*p == '-')
 169                                 errx(1, "negative numbers aren't permitted.");
 170                         if (BN_dec2bn(&val, buf) == 0)
 171                                 errx(1, "%s: illegal numeric format.", argv[0]);
 172                         pr_fact(val);
 173                 }
 174         /* Factor the arguments. */
 175         else
 176                 for (; *argv != NULL; ++argv) {
 177                         if (argv[0][0] == '-')
 178                                 errx(1, "negative numbers aren't permitted.");
 179                         if (BN_dec2bn(&val, argv[0]) == 0)
 180                                 errx(1, "%s: illegal numeric format.", argv[0]);
 181                         pr_fact(val);
 182                 }
 183         exit(0);
 184 }
 185
 186 /*
 187  * pr_fact - print the factors of a number
 188  *
 189  * If the number is 0 or 1, then print the number and return.
 190  * If the number is < 0, print -1, negate the number and continue
 191  * processing.
 192  *
 193  * Print the factors of the number, from the lowest to the highest.
 194  * A factor will be printed numtiple times if it divides the value
 195  * multiple times.
 196  *
 197  * Factors are printed with leading tabs.
 198  */
 199 void
 200 pr_fact(BIGNUM *val)
 201 {
 202         const ubig *fact;               /* The factor found. */
 203
 204         /* Firewall - catch 0 and 1. */
 205         if (BN_is_zero(val))    /* Historical practice; 0 just exits. */
 206                 exit(0);
 207         if (BN_is_one(val)) {
 208                 printf("1: 1\n");
 209                 return;
 210         }
 211
 212         /* Factor value. */
 213
 214         BN_print_dec_fp(stdout, val);
 215         putchar(':');
 216         for (fact = &prime[0]; !BN_is_one(val); ++fact) {
 217                 /* Look for the smallest factor. */
 218                 do {
 219                         if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
 220                                 break;
 221                 } while (++fact <= pr_limit);
 222
 223                 /* Watch for primes larger than the table. */
 224                 if (fact > pr_limit) {
 225 #ifdef HAVE_OPENSSL
 226                         BIGNUM *bnfact;
 227
 228                         bnfact = BN_new();
 229                         BN_set_word(bnfact, *(fact - 1));
 230                         BN_sqr(bnfact, bnfact, ctx);
 231                         if (BN_cmp(bnfact, val) > 0
 232                             || BN_is_prime(val, PRIME_CHECKS, NULL, NULL,
 233                                            NULL) == 1) {
 234                                 putchar(' ');
 235                                 BN_print_dec_fp(stdout, val);
 236                         } else
 237                                 pollard_pminus1(val);
 238 #else
 239                         printf(" %s", BN_bn2dec(val));
 240 #endif
 241                         break;
 242                 }
 243
 244                 /* Divide factor out until none are left. */
 245                 do {
 246                         printf(" %lu", *fact);
 247                         BN_div_word(val, (BN_ULONG)*fact);
 248                 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
 249
 250                 /* Let the user know we're doing something. */
 251                 fflush(stdout);
 252         }
 253         putchar('\n');
 254 }
 255
 256 /*
 257  * Sigh..  No _decimal_ output to file functions in BN.
 258  */
 259 void
 260 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
 261 {
 262         char *buf;
 263
 264         buf = BN_bn2dec(num);
 265         if (buf == NULL)
 266                 return; /* XXX do anything here? */
 267         fprintf(fp, buf);
 268         free(buf);
 269 }
 270
 271 void
 272 usage(void)
 273 {
 274         fprintf(stderr, "usage: factor [value ...]\n");
 275         exit (0);
 276 }
 277
 278
 279
 280
 281 #ifdef HAVE_OPENSSL
 282 /* pollard p-1, algorithm from Jim Gillogly, May 2000 */
 283
 284 void
 285 pollard_pminus1(BIGNUM *val)
 286 {
 287         BIGNUM *base, *rbase, *num, *i, *x;
 288
 289         base = BN_new();
 290         rbase = BN_new();
 291         num = BN_new();
 292         i = BN_new();
 293         x = BN_new();
 294
 295         BN_set_word(rbase, 1);
 296  newbase:
 297         BN_add_word(rbase, 1);
 298         BN_set_word(i, 2);
 299         BN_copy(base, rbase);
 300
 301         for (;;) {
 302                 BN_mod_exp(base, base, i, val, ctx);
 303                 if (BN_is_one(base))
 304                         goto newbase;
 305
 306                 BN_copy(x, base);
 307                 BN_sub_word(x, 1);
 308                 BN_gcd(x, x, val, ctx);
 309
 310                 if (!BN_is_one(x)) {
 311                         if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
 312                             NULL) == 1) {
 313                                 putchar(' ');
 314                                 BN_print_dec_fp(stdout, x);
 315                         } else
 316                                 pollard_pminus1(x);
 317                         fflush(stdout);
 318
 319                         BN_div(num, NULL, val, x, ctx);
 320                         if (BN_is_one(num))
 321                                 return;
 322                         if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
 323                             NULL) == 1) {
 324                                 putchar(' ');
 325                                 BN_print_dec_fp(stdout, num);
 326                                 fflush(stdout);
 327                                 return;
 328                         }
 329                         BN_copy(val, num);
 330                 }
 331                 BN_add_word(i, 1);
 332         }
 333 }
 334 #else
 335 char *
 336 BN_bn2dec(const BIGNUM *val)
 337 {
 338         char *buf;
 339
 340         buf = malloc(100);
 341         if (!buf)
 342                 return buf;
 343         snprintf(buf, 100, "%ld", (long)*val);
 344         return buf;
 345 }
 346
 347 BN_ULONG
 348 BN_div_word(BIGNUM *a, BN_ULONG b)
 349 {
 350         BN_ULONG mod;
 351
 352         mod = *a % b;
 353         *a /= b;
 354         return mod;
 355 }
 356 #endif