factor/factor.c

   1 /*      $NetBSD: factor.c,v 1.12 2002/06/17 15:43:52 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.12 2002/06/17 15:43:52 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 #define BN_new()                ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
  86 #define BN_dec2bn(pp, str)      (**(pp) = atol(str))
  87 #define BN_is_zero(v)           (*(v) == 0)
  88 #define BN_is_one(v)            (*(v) == 1)
  89 #define BN_mod_word(a, b)       (*(a) % (b))
  90 #endif
  91
  92 #include "primes.h"
  93
  94 /*
  95  * prime[i] is the (i-1)th prime.
  96  *
  97  * We are able to sieve 2^32-1 because this byte table yields all primes
  98  * up to 65537 and 65537^2 > 2^32-1.
  99  */
 100 extern const ubig prime[];
 101 extern const ubig *pr_limit;            /* largest prime in the prime array */
 102
 103 #define PRIME_CHECKS    5
 104
 105 #ifdef HAVE_OPENSSL
 106 BN_CTX *ctx;                            /* just use a global context */
 107 #endif
 108
 109 int     main(int, char *[]);
 110 void    pr_fact(BIGNUM *);              /* print factors of a value */
 111 void    BN_print_dec_fp(FILE *, const BIGNUM *);
 112 void    usage(void) __attribute__((__noreturn__));
 113 #ifdef HAVE_OPENSSL
 114 void    pollard_pminus1(BIGNUM *);      /* print factors for big numbers */
 115 #else
 116 char    *BN_bn2dec(const BIGNUM *);
 117 BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
 118 #endif
 119
 120 int
 121 main(int argc, char *argv[])
 122 {
 123         BIGNUM *val;
 124         int ch;
 125         char *p, buf[LINE_MAX];         /* > max number of digits. */
 126
 127 #ifdef HAVE_OPENSSL
 128         ctx = BN_CTX_new();
 129 #endif
 130         val = BN_new();
 131         if (val == NULL)
 132                 errx(1, "can't initialise bignum");
 133
 134         while ((ch = getopt(argc, argv, "")) != -1)
 135                 switch (ch) {
 136                 case '?':
 137                 default:
 138                         usage();
 139                 }
 140         argc -= optind;
 141         argv += optind;
 142
 143         /* No args supplied, read numbers from stdin. */
 144         if (argc == 0)
 145                 for (;;) {
 146                         if (fgets(buf, sizeof(buf), stdin) == NULL) {
 147                                 if (ferror(stdin))
 148                                         err(1, "stdin");
 149                                 exit (0);
 150                         }
 151                         for (p = buf; isblank(*p); ++p);
 152                         if (*p == '\n' || *p == '\0')
 153                                 continue;
 154                         if (*p == '-')
 155                                 errx(1, "negative numbers aren't permitted.");
 156                         if (BN_dec2bn(&val, buf) == 0)
 157                                 errx(1, "%s: illegal numeric format.", argv[0]);
 158                         pr_fact(val);
 159                 }
 160         /* Factor the arguments. */
 161         else
 162                 for (; *argv != NULL; ++argv) {
 163                         if (argv[0][0] == '-')
 164                                 errx(1, "negative numbers aren't permitted.");
 165                         if (BN_dec2bn(&val, argv[0]) == 0)
 166                                 errx(1, "%s: illegal numeric format.", argv[0]);
 167                         pr_fact(val);
 168                 }
 169         exit(0);
 170 }
 171
 172 /*
 173  * pr_fact - print the factors of a number
 174  *
 175  * If the number is 0 or 1, then print the number and return.
 176  * If the number is < 0, print -1, negate the number and continue
 177  * processing.
 178  *
 179  * Print the factors of the number, from the lowest to the highest.
 180  * A factor will be printed numtiple times if it divides the value
 181  * multiple times.
 182  *
 183  * Factors are printed with leading tabs.
 184  */
 185 void
 186 pr_fact(BIGNUM *val)
 187 {
 188         const ubig *fact;               /* The factor found. */
 189
 190         /* Firewall - catch 0 and 1. */
 191         if (BN_is_zero(val))    /* Historical practice; 0 just exits. */
 192                 exit(0);
 193         if (BN_is_one(val)) {
 194                 printf("1: 1\n");
 195                 return;
 196         }
 197
 198         /* Factor value. */
 199
 200         BN_print_dec_fp(stdout, val);
 201         putchar(':');
 202         for (fact = &prime[0]; !BN_is_one(val); ++fact) {
 203                 /* Look for the smallest factor. */
 204                 do {
 205                         if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
 206                                 break;
 207                 } while (++fact <= pr_limit);
 208
 209                 /* Watch for primes larger than the table. */
 210                 if (fact > pr_limit) {
 211 #ifdef HAVE_OPENSSL
 212                         BIGNUM *bnfact;
 213
 214                         bnfact = BN_new();
 215                         BN_set_word(bnfact, *(fact - 1));
 216                         BN_sqr(bnfact, bnfact, ctx);
 217                         if (BN_cmp(bnfact, val) > 0) {
 218                                 putchar(' ');
 219                                 BN_print_dec_fp(stdout, val);
 220                         } else
 221                                 pollard_pminus1(val);
 222 #else
 223                         printf(" %s", BN_bn2dec(val));
 224 #endif
 225                         break;
 226                 }
 227
 228                 /* Divide factor out until none are left. */
 229                 do {
 230                         printf(" %lu", *fact);
 231                         BN_div_word(val, (BN_ULONG)*fact);
 232                 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
 233
 234                 /* Let the user know we're doing something. */
 235                 fflush(stdout);
 236         }
 237         putchar('\n');
 238 }
 239
 240 /*
 241  * Sigh..  No _decimal_ output to file functions in BN.
 242  */
 243 void
 244 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
 245 {
 246         char *buf;
 247
 248         buf = BN_bn2dec(num);
 249         if (buf == NULL)
 250                 return; /* XXX do anything here? */
 251         fprintf(fp, buf);
 252         free(buf);
 253 }
 254
 255 void
 256 usage(void)
 257 {
 258         fprintf(stderr, "usage: factor [value ...]\n");
 259         exit (0);
 260 }
 261
 262
 263
 264
 265 #ifdef HAVE_OPENSSL
 266 /* pollard rho, algorithm from Jim Gillogly, May 2000 */
 267
 268 void
 269 pollard_pminus1(BIGNUM *val)
 270 {
 271         BIGNUM *base, *num, *i, *x;
 272
 273         base = BN_new();
 274         num = BN_new();
 275         i = BN_new();
 276         x = BN_new();
 277
 278         BN_set_word(i, 2);
 279         BN_set_word(base, 2);
 280
 281         for (;;) {
 282                 BN_mod_exp(base, base, i, val, ctx);
 283
 284                 BN_copy(x, base);
 285                 BN_sub_word(x, 1);
 286                 BN_gcd(x, x, val, ctx);
 287
 288                 if (!BN_is_one(x)) {
 289                         if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
 290                             NULL) == 1) {
 291                                 putchar(' ');
 292                                 BN_print_dec_fp(stdout, x);
 293                         } else
 294                                 pollard_pminus1(x);
 295                         fflush(stdout);
 296
 297                         BN_div(num, NULL, val, x, ctx);
 298                         if (BN_is_one(num))
 299                                 return;
 300                         if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
 301                             NULL) == 1) {
 302                                 putchar(' ');
 303                                 BN_print_dec_fp(stdout, num);
 304                                 fflush(stdout);
 305                                 return;
 306                         }
 307                         BN_copy(val, num);
 308                 }
 309                 BN_add_word(i, 1);
 310         }
 311 }
 312 #else
 313 char *
 314 BN_bn2dec(const BIGNUM *val)
 315 {
 316         char *buf;
 317
 318         buf = malloc(100);
 319         if (!buf)
 320                 return buf;
 321         snprintf(buf, 100, "%ld", (long)*val);
 322         return buf;
 323 }
 324
 325 BN_ULONG
 326 BN_div_word(BIGNUM *a, BN_ULONG b)
 327 {
 328         BN_ULONG mod;
 329
 330         mod = *a % b;
 331         *a /= b;
 332         return mod;
 333 }
 334 #endif