factor.c revision 199815
12490Sjkh/* 22490Sjkh * Copyright (c) 1989, 1993 32490Sjkh * The Regents of the University of California. All rights reserved. 42490Sjkh * 52490Sjkh * This code is derived from software contributed to Berkeley by 62490Sjkh * Landon Curt Noll. 72490Sjkh * 82490Sjkh * Redistribution and use in source and binary forms, with or without 92490Sjkh * modification, are permitted provided that the following conditions 102490Sjkh * are met: 112490Sjkh * 1. Redistributions of source code must retain the above copyright 122490Sjkh * notice, this list of conditions and the following disclaimer. 132490Sjkh * 2. Redistributions in binary form must reproduce the above copyright 142490Sjkh * notice, this list of conditions and the following disclaimer in the 152490Sjkh * documentation and/or other materials provided with the distribution. 16199815Sfanf * 3. Neither the name of the University nor the names of its contributors 172490Sjkh * may be used to endorse or promote products derived from this software 182490Sjkh * without specific prior written permission. 192490Sjkh * 202490Sjkh * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 212490Sjkh * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 222490Sjkh * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 232490Sjkh * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 242490Sjkh * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 252490Sjkh * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 262490Sjkh * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 272490Sjkh * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 282490Sjkh * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 292490Sjkh * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 302490Sjkh * SUCH DAMAGE. 312490Sjkh */ 322490Sjkh 332490Sjkh#ifndef lint 34199815Sfanf#include <sys/cdefs.h> 35199815Sfanf#ifdef __COPYRIGHT 36199815Sfanf__COPYRIGHT("@(#) Copyright (c) 1989, 1993\ 37199815Sfanf The Regents of the University of California. All rights reserved."); 3853920Sbillf#endif 39199815Sfanf#ifdef __SCCSID 40199815Sfanf__SCCSID("@(#)factor.c 8.4 (Berkeley) 5/4/95"); 41199815Sfanf#endif 42199815Sfanf#ifdef __RCSID 43199815Sfanf__RCSID("$NetBSD: factor.c,v 1.19 2009/08/12 05:54:31 dholland Exp $"); 44199815Sfanf#endif 45199815Sfanf#ifdef __FBSDID 46199815Sfanf__FBSDID("$FreeBSD: head/games/factor/factor.c 199815 2009-11-26 00:38:13Z fanf $"); 47199815Sfanf#endif 482490Sjkh#endif /* not lint */ 492490Sjkh 502490Sjkh/* 512490Sjkh * factor - factor a number into primes 522490Sjkh * 532490Sjkh * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo 542490Sjkh * 552490Sjkh * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\ 562490Sjkh * 572490Sjkh * usage: 58104720Sfanf * factor [-h] [number] ... 592490Sjkh * 602490Sjkh * The form of the output is: 612490Sjkh * 622490Sjkh * number: factor1 factor1 factor2 factor3 factor3 factor3 ... 632490Sjkh * 64199815Sfanf * where factor1 <= factor2 <= factor3 <= ... 652490Sjkh * 662490Sjkh * If no args are given, the list of numbers are read from stdin. 672490Sjkh */ 682490Sjkh 69104720Sfanf#include <ctype.h> 702490Sjkh#include <err.h> 712490Sjkh#include <errno.h> 722490Sjkh#include <limits.h> 732490Sjkh#include <stdio.h> 742490Sjkh#include <stdlib.h> 7523726Speter#include <unistd.h> 762490Sjkh 772490Sjkh#include "primes.h" 782490Sjkh 79104722Sfanf#ifdef HAVE_OPENSSL 802490Sjkh 81104722Sfanf#include <openssl/bn.h> 82104722Sfanf 83104722Sfanf#define PRIME_CHECKS 5 84104722Sfanf 85104722Sfanfstatic void pollard_pminus1(BIGNUM *); /* print factors for big numbers */ 86104722Sfanf 87104722Sfanf#else 88104722Sfanf 89104722Sfanftypedef ubig BIGNUM; 90104722Sfanftypedef u_long BN_ULONG; 91104722Sfanf 92104722Sfanf#define BN_CTX int 93104722Sfanf#define BN_CTX_new() NULL 94104722Sfanf#define BN_new() ((BIGNUM *)calloc(sizeof(BIGNUM), 1)) 95104722Sfanf#define BN_is_zero(v) (*(v) == 0) 96104722Sfanf#define BN_is_one(v) (*(v) == 1) 97104722Sfanf#define BN_mod_word(a, b) (*(a) % (b)) 98104722Sfanf 99104722Sfanfstatic int BN_dec2bn(BIGNUM **a, const char *str); 100104722Sfanfstatic int BN_hex2bn(BIGNUM **a, const char *str); 101104722Sfanfstatic BN_ULONG BN_div_word(BIGNUM *, BN_ULONG); 102104722Sfanfstatic void BN_print_fp(FILE *, const BIGNUM *); 103104722Sfanf 104104722Sfanf#endif 105104722Sfanf 106104722Sfanfstatic void BN_print_dec_fp(FILE *, const BIGNUM *); 107104722Sfanf 108104722Sfanfstatic void pr_fact(BIGNUM *); /* print factors of a value */ 109104722Sfanfstatic void pr_print(BIGNUM *); /* print a prime */ 110104720Sfanfstatic void usage(void); 11142338Simp 112104722Sfanfstatic BN_CTX *ctx; /* just use a global context */ 113104722Sfanfstatic int hflag; 114104722Sfanf 1152490Sjkhint 116104720Sfanfmain(int argc, char *argv[]) 1172490Sjkh{ 118104722Sfanf BIGNUM *val; 1192490Sjkh int ch; 120104720Sfanf char *p, buf[LINE_MAX]; /* > max number of digits. */ 1212490Sjkh 122104722Sfanf ctx = BN_CTX_new(); 123104722Sfanf val = BN_new(); 124104722Sfanf if (val == NULL) 125104722Sfanf errx(1, "can't initialise bignum"); 126104722Sfanf 12742338Simp while ((ch = getopt(argc, argv, "h")) != -1) 1282490Sjkh switch (ch) { 12942338Simp case 'h': 13042338Simp hflag++; 13142338Simp break; 1322490Sjkh case '?': 1332490Sjkh default: 1342490Sjkh usage(); 1352490Sjkh } 1362490Sjkh argc -= optind; 1372490Sjkh argv += optind; 1382490Sjkh 1392490Sjkh /* No args supplied, read numbers from stdin. */ 1402490Sjkh if (argc == 0) 1412490Sjkh for (;;) { 1422490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1432490Sjkh if (ferror(stdin)) 1442490Sjkh err(1, "stdin"); 1452490Sjkh exit (0); 1462490Sjkh } 1472490Sjkh for (p = buf; isblank(*p); ++p); 1482490Sjkh if (*p == '\n' || *p == '\0') 1492490Sjkh continue; 1502490Sjkh if (*p == '-') 1512490Sjkh errx(1, "negative numbers aren't permitted."); 152104722Sfanf if (BN_dec2bn(&val, buf) == 0 && 153104722Sfanf BN_hex2bn(&val, buf) == 0) 1542490Sjkh errx(1, "%s: illegal numeric format.", buf); 1552490Sjkh pr_fact(val); 1562490Sjkh } 1572490Sjkh /* Factor the arguments. */ 1582490Sjkh else 1592490Sjkh for (; *argv != NULL; ++argv) { 1602490Sjkh if (argv[0][0] == '-') 1612490Sjkh errx(1, "negative numbers aren't permitted."); 162104722Sfanf if (BN_dec2bn(&val, argv[0]) == 0 && 163104722Sfanf BN_hex2bn(&val, argv[0]) == 0) 1642490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1652490Sjkh pr_fact(val); 1662490Sjkh } 1672490Sjkh exit(0); 1682490Sjkh} 1692490Sjkh 1702490Sjkh/* 1712490Sjkh * pr_fact - print the factors of a number 1722490Sjkh * 1732490Sjkh * Print the factors of the number, from the lowest to the highest. 17485408Sroam * A factor will be printed multiple times if it divides the value 1752490Sjkh * multiple times. 1762490Sjkh * 1772490Sjkh * Factors are printed with leading tabs. 1782490Sjkh */ 179104720Sfanfstatic void 180104722Sfanfpr_fact(BIGNUM *val) 1812490Sjkh{ 182104720Sfanf const ubig *fact; /* The factor found. */ 1832490Sjkh 1842490Sjkh /* Firewall - catch 0 and 1. */ 185104722Sfanf if (BN_is_zero(val)) /* Historical practice; 0 just exits. */ 1862490Sjkh exit(0); 187104722Sfanf if (BN_is_one(val)) { 188104720Sfanf printf("1: 1\n"); 1892490Sjkh return; 1902490Sjkh } 1912490Sjkh 1922490Sjkh /* Factor value. */ 193104722Sfanf 194104722Sfanf if (hflag) { 195104722Sfanf fputs("0x", stdout); 196104722Sfanf BN_print_fp(stdout, val); 197104722Sfanf } else 198104722Sfanf BN_print_dec_fp(stdout, val); 199104722Sfanf putchar(':'); 200104722Sfanf for (fact = &prime[0]; !BN_is_one(val); ++fact) { 2012490Sjkh /* Look for the smallest factor. */ 2022490Sjkh do { 203104722Sfanf if (BN_mod_word(val, (BN_ULONG)*fact) == 0) 2042490Sjkh break; 2052490Sjkh } while (++fact <= pr_limit); 2062490Sjkh 2072490Sjkh /* Watch for primes larger than the table. */ 2082490Sjkh if (fact > pr_limit) { 209104722Sfanf#ifdef HAVE_OPENSSL 210104722Sfanf BIGNUM *bnfact; 211104722Sfanf 212104722Sfanf bnfact = BN_new(); 213104722Sfanf BN_set_word(bnfact, *(fact - 1)); 214104722Sfanf BN_sqr(bnfact, bnfact, ctx); 215199815Sfanf if (BN_cmp(bnfact, val) > 0 || 216199815Sfanf BN_is_prime(val, PRIME_CHECKS, 217199815Sfanf NULL, NULL, NULL) == 1) 218104722Sfanf pr_print(val); 219104722Sfanf else 220104722Sfanf pollard_pminus1(val); 221104722Sfanf#else 222104722Sfanf pr_print(val); 223104722Sfanf#endif 2242490Sjkh break; 2252490Sjkh } 2262490Sjkh 2272490Sjkh /* Divide factor out until none are left. */ 2282490Sjkh do { 229104720Sfanf printf(hflag ? " 0x%lx" : " %lu", *fact); 230104722Sfanf BN_div_word(val, (BN_ULONG)*fact); 231104722Sfanf } while (BN_mod_word(val, (BN_ULONG)*fact) == 0); 2322490Sjkh 2332490Sjkh /* Let the user know we're doing something. */ 234104720Sfanf fflush(stdout); 2352490Sjkh } 236104720Sfanf putchar('\n'); 2372490Sjkh} 2382490Sjkh 239104720Sfanfstatic void 240104722Sfanfpr_print(BIGNUM *val) 241104722Sfanf{ 242104722Sfanf if (hflag) { 243104722Sfanf fputs(" 0x", stdout); 244104722Sfanf BN_print_fp(stdout, val); 245104722Sfanf } else { 246104722Sfanf putchar(' '); 247104722Sfanf BN_print_dec_fp(stdout, val); 248104722Sfanf } 249104722Sfanf} 250104722Sfanf 251104722Sfanfstatic void 252104720Sfanfusage(void) 2532490Sjkh{ 254104720Sfanf fprintf(stderr, "usage: factor [-h] [value ...]\n"); 255104720Sfanf exit(1); 2562490Sjkh} 257104722Sfanf 258104722Sfanf#ifdef HAVE_OPENSSL 259104722Sfanf 260199815Sfanf/* pollard p-1, algorithm from Jim Gillogly, May 2000 */ 261104722Sfanfstatic void 262104722Sfanfpollard_pminus1(BIGNUM *val) 263104722Sfanf{ 264199815Sfanf BIGNUM *base, *rbase, *num, *i, *x; 265104722Sfanf 266104722Sfanf base = BN_new(); 267199815Sfanf rbase = BN_new(); 268104722Sfanf num = BN_new(); 269104722Sfanf i = BN_new(); 270104722Sfanf x = BN_new(); 271104722Sfanf 272199815Sfanf BN_set_word(rbase, 1); 273199815Sfanfnewbase: 274199815Sfanf BN_add_word(rbase, 1); 275104722Sfanf BN_set_word(i, 2); 276199815Sfanf BN_copy(base, rbase); 277104722Sfanf 278104722Sfanf for (;;) { 279104722Sfanf BN_mod_exp(base, base, i, val, ctx); 280199815Sfanf if (BN_is_one(base)) 281199815Sfanf goto newbase; 282104722Sfanf 283104722Sfanf BN_copy(x, base); 284104722Sfanf BN_sub_word(x, 1); 285104722Sfanf BN_gcd(x, x, val, ctx); 286104722Sfanf 287104722Sfanf if (!BN_is_one(x)) { 288104722Sfanf if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL, 289104722Sfanf NULL) == 1) 290104722Sfanf pr_print(x); 291104722Sfanf else 292104722Sfanf pollard_pminus1(x); 293104722Sfanf fflush(stdout); 294104722Sfanf 295104722Sfanf BN_div(num, NULL, val, x, ctx); 296104722Sfanf if (BN_is_one(num)) 297104722Sfanf return; 298104722Sfanf if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL, 299104722Sfanf NULL) == 1) { 300104722Sfanf pr_print(num); 301104722Sfanf fflush(stdout); 302104722Sfanf return; 303104722Sfanf } 304104722Sfanf BN_copy(val, num); 305104722Sfanf } 306104722Sfanf BN_add_word(i, 1); 307104722Sfanf } 308104722Sfanf} 309104722Sfanf 310104722Sfanf/* 311104722Sfanf * Sigh.. No _decimal_ output to file functions in BN. 312104722Sfanf */ 313104722Sfanfstatic void 314104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 315104722Sfanf{ 316104722Sfanf char *buf; 317104722Sfanf 318104722Sfanf buf = BN_bn2dec(num); 319104722Sfanf if (buf == NULL) 320104722Sfanf return; /* XXX do anything here? */ 321104722Sfanf fprintf(fp, buf); 322104722Sfanf free(buf); 323104722Sfanf} 324104722Sfanf 325104722Sfanf#else 326104722Sfanf 327104722Sfanfstatic void 328104722SfanfBN_print_fp(FILE *fp, const BIGNUM *num) 329104722Sfanf{ 330104722Sfanf fprintf(fp, "%lx", (unsigned long)*num); 331104722Sfanf} 332104722Sfanf 333104722Sfanfstatic void 334104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 335104722Sfanf{ 336104722Sfanf fprintf(fp, "%lu", (unsigned long)*num); 337104722Sfanf} 338104722Sfanf 339104722Sfanfstatic int 340104722SfanfBN_dec2bn(BIGNUM **a, const char *str) 341104722Sfanf{ 342104722Sfanf char *p; 343104722Sfanf 344104722Sfanf errno = 0; 345104722Sfanf **a = strtoul(str, &p, 10); 346104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 347104722Sfanf} 348104722Sfanf 349104722Sfanfstatic int 350104722SfanfBN_hex2bn(BIGNUM **a, const char *str) 351104722Sfanf{ 352104722Sfanf char *p; 353104722Sfanf 354104722Sfanf errno = 0; 355104722Sfanf **a = strtoul(str, &p, 16); 356104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 357104722Sfanf} 358104722Sfanf 359104722Sfanfstatic BN_ULONG 360104722SfanfBN_div_word(BIGNUM *a, BN_ULONG b) 361104722Sfanf{ 362104722Sfanf BN_ULONG mod; 363104722Sfanf 364104722Sfanf mod = *a % b; 365104722Sfanf *a /= b; 366104722Sfanf return mod; 367104722Sfanf} 368104722Sfanf 369104722Sfanf#endif 370