factor/factor.c

   1 /*      $NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb 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. All advertising materials mentioning features or use of this software
  19  *    must display the following acknowledgement:
  20  *      This product includes software developed by the University of
  21  *      California, Berkeley and its contributors.
  22  * 4. Neither the name of the University nor the names of its contributors
  23  *    may be used to endorse or promote products derived from this software
  24  *    without specific prior written permission.
  25  *
  26  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  27  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  28  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  29  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  30  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  31  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  32  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  33  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  34  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  35  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  36  * SUCH DAMAGE.
  37  */
  38
  39 #include <sys/cdefs.h>
  40 #ifndef lint
  41 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\n\
  42         The Regents of the University of California.  All rights reserved.\n");
  43 #endif /* not lint */
  44
  45 #ifndef lint
  46 #if 0
  47 static char sccsid[] = "@(#)factor.c    8.4 (Berkeley) 5/4/95";
  48 #else
  49 __RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $");
  50 #endif
  51 #endif /* not lint */
  52
  53 /*
  54  * factor - factor a number into primes
  55  *
  56  * By: Landon Curt Noll   chongo@toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
  57  *
  58  *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
  59  *
  60  * usage:
  61  *      factor [number] ...
  62  *
  63  * The form of the output is:
  64  *
  65  *      number: factor1 factor1 factor2 factor3 factor3 factor3 ...
  66  *
  67  * where factor1 < factor2 < factor3 < ...
  68  *
  69  * If no args are given, the list of numbers are read from stdin.
  70  */
  71
  72 #include <ctype.h>
  73 #include <err.h>
  74 #include <errno.h>
  75 #include <limits.h>
  76 #include <stdio.h>
  77 #include <stdlib.h>
  78 #include <unistd.h>
  79
  80 #ifdef HAVE_OPENSSL
  81 #include <openssl/bn.h>
  82 #else
  83 typedef long    BIGNUM;
  84 typedef u_long  BN_ULONG;
  85 int     BN_dec2bn(BIGNUM **a, const char *str);
  86 #define BN_new()                ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
  87 #define BN_is_zero(v)           (*(v) == 0)
  88 #define BN_is_one(v)            (*(v) == 1)
  89 #define BN_new()                ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
  90 #define BN_is_zero(v)           (*(v) == 0)
  91 #define BN_is_one(v)            (*(v) == 1)
  92 #define BN_mod_word(a, b)       (*(a) % (b))
  93 #endif
  94
  95 #include "primes.h"
  96
  97 /*
  98  * prime[i] is the (i-1)th prime.
  99  *
 100  * We are able to sieve 2^32-1 because this byte table yields all primes
 101  * up to 65537 and 65537^2 > 2^32-1.
 102  */
 103 extern const ubig prime[];
 104 extern const ubig *pr_limit;            /* largest prime in the prime array */
 105
 106 #define PRIME_CHECKS    5
 107
 108 #ifdef HAVE_OPENSSL
 109 BN_CTX *ctx;                            /* just use a global context */
 110 #endif
 111
 112 int     main(int, char *[]);
 113 void    pr_fact(BIGNUM *);              /* print factors of a value */
 114 void    BN_print_dec_fp(FILE *, const BIGNUM *);
 115 void    usage(void) __attribute__((__noreturn__));
 116 #ifdef HAVE_OPENSSL
 117 void    pollard_pminus1(BIGNUM *);      /* print factors for big numbers */
 118 #else
 119 char    *BN_bn2dec(const BIGNUM *);
 120 BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
 121 #endif
 122
 123
 124 #ifndef HAVE_OPENSSL
 125 int
 126 BN_dec2bn(BIGNUM **a, const char *str)
 127 {
 128         char *p;
 129
 130         errno = 0;
 131         **a = strtoul(str, &p, 10);
 132         if (errno)
 133                 err(1, "%s", str);
 134         return (*p == '\n' || *p == '\0');
 135 }
 136 #endif
 137
 138 int
 139 main(int argc, char *argv[])
 140 {
 141         BIGNUM *val;
 142         int ch;
 143         char *p, buf[LINE_MAX];         /* > max number of digits. */
 144
 145 #ifdef HAVE_OPENSSL
 146         ctx = BN_CTX_new();
 147 #endif
 148         val = BN_new();
 149         if (val == NULL)
 150                 errx(1, "can't initialise bignum");
 151
 152         while ((ch = getopt(argc, argv, "")) != -1)
 153                 switch (ch) {
 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(*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);
 177                 }
 178         /* Factor the arguments. */
 179         else
 180                 for (; *argv != NULL; ++argv) {
 181                         if (argv[0][0] == '-')
 182                                 errx(1, "negative numbers aren't permitted.");
 183                         if (BN_dec2bn(&val, argv[0]) == 0)
 184                                 errx(1, "%s: illegal numeric format.", argv[0]);
 185                         pr_fact(val);
 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  * A factor will be printed numtiple times if it divides the value
 199  * multiple times.
 200  *
 201  * Factors are printed with leading tabs.
 202  */
 203 void
 204 pr_fact(BIGNUM *val)
 205 {
 206         const ubig *fact;               /* The factor found. */
 207
 208         /* Firewall - catch 0 and 1. */
 209         if (BN_is_zero(val))    /* Historical practice; 0 just exits. */
 210                 exit(0);
 211         if (BN_is_one(val)) {
 212                 printf("1: 1\n");
 213                 return;
 214         }
 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                 do {
 223                         if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
 224                                 break;
 225                 } while (++fact <= pr_limit);
 226
 227                 /* Watch for primes larger than the table. */
 228                 if (fact > pr_limit) {
 229 #ifdef HAVE_OPENSSL
 230                         BIGNUM *bnfact;
 231
 232                         bnfact = BN_new();
 233                         BN_set_word(bnfact, *(fact - 1));
 234                         BN_sqr(bnfact, bnfact, ctx);
 235                         if (BN_cmp(bnfact, val) > 0) {
 236                                 putchar(' ');
 237                                 BN_print_dec_fp(stdout, val);
 238                         } else
 239                                 pollard_pminus1(val);
 240 #else
 241                         printf(" %s", BN_bn2dec(val));
 242 #endif
 243                         break;
 244                 }
 245
 246                 /* Divide factor out until none are left. */
 247                 do {
 248                         printf(" %lu", *fact);
 249                         BN_div_word(val, (BN_ULONG)*fact);
 250                 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
 251
 252                 /* Let the user know we're doing something. */
 253                 fflush(stdout);
 254         }
 255         putchar('\n');
 256 }
 257
 258 /*
 259  * Sigh..  No _decimal_ output to file functions in BN.
 260  */
 261 void
 262 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
 263 {
 264         char *buf;
 265
 266         buf = BN_bn2dec(num);
 267         if (buf == NULL)
 268                 return; /* XXX do anything here? */
 269         fprintf(fp, buf);
 270         free(buf);
 271 }
 272
 273 void
 274 usage(void)
 275 {
 276         fprintf(stderr, "usage: factor [value ...]\n");
 277         exit (0);
 278 }
 279
 280
 281
 282
 283 #ifdef HAVE_OPENSSL
 284 /* pollard rho, algorithm from Jim Gillogly, May 2000 */
 285
 286 void
 287 pollard_pminus1(BIGNUM *val)
 288 {
 289         BIGNUM *base, *num, *i, *x;
 290
 291         base = BN_new();
 292         num = BN_new();
 293         i = BN_new();
 294         x = BN_new();
 295
 296         BN_set_word(i, 2);
 297         BN_set_word(base, 2);
 298
 299         for (;;) {
 300                 BN_mod_exp(base, base, i, val, ctx);
 301
 302                 BN_copy(x, base);
 303                 BN_sub_word(x, 1);
 304                 BN_gcd(x, x, val, ctx);
 305
 306                 if (!BN_is_one(x)) {
 307                         if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
 308                             NULL) == 1) {
 309                                 putchar(' ');
 310                                 BN_print_dec_fp(stdout, x);
 311                         } else
 312                                 pollard_pminus1(x);
 313                         fflush(stdout);
 314
 315                         BN_div(num, NULL, val, x, ctx);
 316                         if (BN_is_one(num))
 317                                 return;
 318                         if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
 319                             NULL) == 1) {
 320                                 putchar(' ');
 321                                 BN_print_dec_fp(stdout, num);
 322                                 fflush(stdout);
 323                                 return;
 324                         }
 325                         BN_copy(val, num);
 326                 }
 327                 BN_add_word(i, 1);
 328         }
 329 }
 330 #else
 331 char *
 332 BN_bn2dec(const BIGNUM *val)
 333 {
 334         char *buf;
 335
 336         buf = malloc(100);
 337         if (!buf)
 338                 return buf;
 339         snprintf(buf, 100, "%ld", (long)*val);
 340         return buf;
 341 }
 342
 343 BN_ULONG
 344 BN_div_word(BIGNUM *a, BN_ULONG b)
 345 {
 346         BN_ULONG mod;
 347
 348         mod = *a % b;
 349         *a /= b;
 350         return mod;
 351 }
 352 #endif