1: /*      $NetBSD: factor.c,v 1.15 2004/02/08 11:47:36 jsm 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.15 2004/02/08 11:47:36 jsm 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) __attribute__((__noreturn__));
  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:         /* Revoke setgid privileges */
  149:         setregid(getgid(), getgid());
  150: 
  151:         while ((ch = getopt(argc, argv, "")) != -1)
  152:                 switch (ch) {
  153:                 case '?':
  154:                 default:
  155:                         usage();
  156:                 }
  157:         argc -= optind;
  158:         argv += optind;
  159: 
  160:         /* No args supplied, read numbers from stdin. */
  161:         if (argc == 0)
  162:                 for (;;) {
  163:                         if (fgets(buf, sizeof(buf), stdin) == NULL) {
  164:                                 if (ferror(stdin))
  165:                                         err(1, "stdin");
  166:                                 exit (0);
  167:                         }
  168:                         for (p = buf; isblank(*p); ++p);
  169:                         if (*p == '\n' || *p == '\0')
  170:                                 continue;
  171:                         if (*p == '-')
  172:                                 errx(1, "negative numbers aren't permitted.");
  173:                         if (BN_dec2bn(&val, buf) == 0)
  174:                                 errx(1, "%s: illegal numeric format.", argv[0]);
  175:                         pr_fact(val);
  176:                 }
  177:         /* Factor the arguments. */
  178:         else
  179:                 for (; *argv != NULL; ++argv) {
  180:                         if (argv[0][0] == '-')
  181:                                 errx(1, "negative numbers aren't permitted.");
  182:                         if (BN_dec2bn(&val, argv[0]) == 0)
  183:                                 errx(1, "%s: illegal numeric format.", argv[0]);
  184:                         pr_fact(val);
  185:                 }
  186:         exit(0);
  187: }
  188: 
  189: /*
  190:  * pr_fact - print the factors of a number
  191:  *
  192:  * If the number is 0 or 1, then print the number and return.
  193:  * If the number is < 0, print -1, negate the number and continue
  194:  * processing.
  195:  *
  196:  * Print the factors of the number, from the lowest to the highest.
  197:  * A factor will be printed numtiple times if it divides the value
  198:  * multiple times.
  199:  *
  200:  * Factors are printed with leading tabs.
  201:  */
  202: void
  203: pr_fact(BIGNUM *val)
  204: {
  205:         const ubig *fact;              /* The factor found. */
  206: 
  207:         /* Firewall - catch 0 and 1. */
  208:         if (BN_is_zero(val))   /* Historical practice; 0 just exits. */
  209:                 exit(0);
  210:         if (BN_is_one(val)) {
  211:                 printf("1: 1\n");
  212:                 return;
  213:         }
  214: 
  215:         /* Factor value. */
  216: 
  217:         BN_print_dec_fp(stdout, val);
  218:         putchar(':');
  219:         for (fact = &prime[0]; !BN_is_one(val); ++fact) {
  220:                 /* Look for the smallest factor. */
  221:                 do {
  222:                         if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
  223:                                 break;
  224:                 } while (++fact <= pr_limit);
  225: 
  226:                 /* Watch for primes larger than the table. */
  227:                 if (fact > pr_limit) {
  228: #ifdef HAVE_OPENSSL
  229:                         BIGNUM *bnfact;
  230: 
  231:                         bnfact = BN_new();
  232:                         BN_set_word(bnfact, *(fact - 1));
  233:                         BN_sqr(bnfact, bnfact, ctx);
  234:                         if (BN_cmp(bnfact, val) > 0
  235:                             || BN_is_prime(val, PRIME_CHECKS, NULL, NULL,
  236:                                            NULL) == 1) {
  237:                                 putchar(' ');
  238:                                 BN_print_dec_fp(stdout, val);
  239:                         } else
  240:                                 pollard_pminus1(val);
  241: #else
  242:                         printf(" %s", BN_bn2dec(val));
  243: #endif
  244:                         break;
  245:                 }
  246: 
  247:                 /* Divide factor out until none are left. */
  248:                 do {
  249:                         printf(" %lu", *fact);
  250:                         BN_div_word(val, (BN_ULONG)*fact);
  251:                 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
  252: 
  253:                 /* Let the user know we're doing something. */
  254:                 fflush(stdout);
  255:         }
  256:         putchar('\n');
  257: }
  258: 
  259: /*
  260:  * Sigh..  No _decimal_ output to file functions in BN.
  261:  */
  262: void
  263: BN_print_dec_fp(FILE *fp, const BIGNUM *num)
  264: {
  265:         char *buf;
  266:         
  267:         buf = BN_bn2dec(num);
  268:         if (buf == NULL)
  269:                 return;       /* XXX do anything here? */
  270:         fprintf(fp, buf);
  271:         free(buf);
  272: }
  273: 
  274: void
  275: usage(void)
  276: {
  277:         fprintf(stderr, "usage: factor [value ...]\n");
  278:         exit (0);
  279: }
  280: 
  281: 
  282: 
  283: 
  284: #ifdef HAVE_OPENSSL
  285: /* pollard p-1, algorithm from Jim Gillogly, May 2000 */
  286: 
  287: void
  288: pollard_pminus1(BIGNUM *val)
  289: {
  290:         BIGNUM *base, *rbase, *num, *i, *x;
  291: 
  292:         base = BN_new();
  293:         rbase = BN_new();
  294:         num = BN_new();
  295:         i = BN_new();
  296:         x = BN_new();
  297: 
  298:         BN_set_word(rbase, 1);
  299:  newbase:
  300:         BN_add_word(rbase, 1);
  301:         BN_set_word(i, 2);
  302:         BN_copy(base, rbase);
  303: 
  304:         for (;;) {
  305:                 BN_mod_exp(base, base, i, val, ctx);
  306:                 if (BN_is_one(base))
  307:                         goto newbase;
  308: 
  309:                 BN_copy(x, base);
  310:                 BN_sub_word(x, 1);
  311:                 BN_gcd(x, x, val, ctx);
  312: 
  313:                 if (!BN_is_one(x)) {
  314:                         if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
  315:                             NULL) == 1) {
  316:                                 putchar(' ');
  317:                                 BN_print_dec_fp(stdout, x);
  318:                         } else
  319:                                 pollard_pminus1(x);
  320:                         fflush(stdout);
  321: 
  322:                         BN_div(num, NULL, val, x, ctx);
  323:                         if (BN_is_one(num))
  324:                                 return;
  325:                         if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
  326:                             NULL) == 1) {
  327:                                 putchar(' ');
  328:                                 BN_print_dec_fp(stdout, num);
  329:                                 fflush(stdout);
  330:                                 return;
  331:                         }
  332:                         BN_copy(val, num);
  333:                 }
  334:                 BN_add_word(i, 1);
  335:         }
  336: }
  337: #else
  338: char *
  339: BN_bn2dec(const BIGNUM *val)
  340: {
  341:         char *buf;
  342: 
  343:         buf = malloc(100);
  344:         if (!buf)
  345:                 return buf;
  346:         snprintf(buf, 100, "%ld", (long)*val);
  347:         return buf;
  348: }
  349: 
  350: BN_ULONG
  351: BN_div_word(BIGNUM *a, BN_ULONG b)
  352: {
  353:         BN_ULONG mod;
  354: 
  355:         mod = *a % b;
  356:         *a /= b;
  357:         return mod;
  358: }
  359: #endif