games/factor/factor.c

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