factor.c revision 104722
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. 162490Sjkh * 3. All advertising materials mentioning features or use of this software 172490Sjkh * must display the following acknowledgement: 182490Sjkh * This product includes software developed by the University of 192490Sjkh * California, Berkeley and its contributors. 202490Sjkh * 4. Neither the name of the University nor the names of its contributors 212490Sjkh * may be used to endorse or promote products derived from this software 222490Sjkh * without specific prior written permission. 232490Sjkh * 242490Sjkh * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 252490Sjkh * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 262490Sjkh * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 272490Sjkh * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 282490Sjkh * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 292490Sjkh * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 302490Sjkh * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 312490Sjkh * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 322490Sjkh * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 332490Sjkh * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 342490Sjkh * SUCH DAMAGE. 352490Sjkh */ 362490Sjkh 372490Sjkh#ifndef lint 3853920Sbillfstatic const char copyright[] = 392490Sjkh"@(#) Copyright (c) 1989, 1993\n\ 402490Sjkh The Regents of the University of California. All rights reserved.\n"; 412490Sjkh#endif /* not lint */ 422490Sjkh 432490Sjkh#ifndef lint 4453920Sbillf#if 0 4523726Speterstatic char sccsid[] = "@(#)factor.c 8.4 (Berkeley) 5/4/95"; 46104722Sfanf__RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $"); 4753920Sbillf#endif 4853920Sbillfstatic const char rcsid[] = 4953920Sbillf "$FreeBSD: head/games/factor/factor.c 104722 2002-10-09 19:55:04Z fanf $"; 502490Sjkh#endif /* not lint */ 512490Sjkh 522490Sjkh/* 532490Sjkh * factor - factor a number into primes 542490Sjkh * 552490Sjkh * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo 562490Sjkh * 572490Sjkh * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\ 582490Sjkh * 592490Sjkh * usage: 60104720Sfanf * factor [-h] [number] ... 612490Sjkh * 622490Sjkh * The form of the output is: 632490Sjkh * 642490Sjkh * number: factor1 factor1 factor2 factor3 factor3 factor3 ... 652490Sjkh * 662490Sjkh * where factor1 < factor2 < factor3 < ... 672490Sjkh * 682490Sjkh * If no args are given, the list of numbers are read from stdin. 692490Sjkh */ 702490Sjkh 71104720Sfanf#include <ctype.h> 722490Sjkh#include <err.h> 732490Sjkh#include <errno.h> 742490Sjkh#include <limits.h> 752490Sjkh#include <stdio.h> 762490Sjkh#include <stdlib.h> 7723726Speter#include <unistd.h> 782490Sjkh 792490Sjkh#include "primes.h" 802490Sjkh 81104722Sfanf#ifdef HAVE_OPENSSL 822490Sjkh 83104722Sfanf#include <openssl/bn.h> 84104722Sfanf 85104722Sfanf#define PRIME_CHECKS 5 86104722Sfanf 87104722Sfanfstatic void pollard_pminus1(BIGNUM *); /* print factors for big numbers */ 88104722Sfanf 89104722Sfanf#else 90104722Sfanf 91104722Sfanftypedef ubig BIGNUM; 92104722Sfanftypedef u_long BN_ULONG; 93104722Sfanf 94104722Sfanf#define BN_CTX int 95104722Sfanf#define BN_CTX_new() NULL 96104722Sfanf#define BN_new() ((BIGNUM *)calloc(sizeof(BIGNUM), 1)) 97104722Sfanf#define BN_is_zero(v) (*(v) == 0) 98104722Sfanf#define BN_is_one(v) (*(v) == 1) 99104722Sfanf#define BN_mod_word(a, b) (*(a) % (b)) 100104722Sfanf 101104722Sfanfstatic int BN_dec2bn(BIGNUM **a, const char *str); 102104722Sfanfstatic int BN_hex2bn(BIGNUM **a, const char *str); 103104722Sfanfstatic BN_ULONG BN_div_word(BIGNUM *, BN_ULONG); 104104722Sfanfstatic void BN_print_fp(FILE *, const BIGNUM *); 105104722Sfanf 106104722Sfanf#endif 107104722Sfanf 108104722Sfanfstatic void BN_print_dec_fp(FILE *, const BIGNUM *); 109104722Sfanf 110104722Sfanfstatic void pr_fact(BIGNUM *); /* print factors of a value */ 111104722Sfanfstatic void pr_print(BIGNUM *); /* print a prime */ 112104720Sfanfstatic void usage(void); 11342338Simp 114104722Sfanfstatic BN_CTX *ctx; /* just use a global context */ 115104722Sfanfstatic int hflag; 116104722Sfanf 1172490Sjkhint 118104720Sfanfmain(int argc, char *argv[]) 1192490Sjkh{ 120104722Sfanf BIGNUM *val; 1212490Sjkh int ch; 122104720Sfanf char *p, buf[LINE_MAX]; /* > max number of digits. */ 1232490Sjkh 124104722Sfanf ctx = BN_CTX_new(); 125104722Sfanf val = BN_new(); 126104722Sfanf if (val == NULL) 127104722Sfanf errx(1, "can't initialise bignum"); 128104722Sfanf 12942338Simp while ((ch = getopt(argc, argv, "h")) != -1) 1302490Sjkh switch (ch) { 13142338Simp case 'h': 13242338Simp hflag++; 13342338Simp break; 1342490Sjkh case '?': 1352490Sjkh default: 1362490Sjkh usage(); 1372490Sjkh } 1382490Sjkh argc -= optind; 1392490Sjkh argv += optind; 1402490Sjkh 1412490Sjkh /* No args supplied, read numbers from stdin. */ 1422490Sjkh if (argc == 0) 1432490Sjkh for (;;) { 1442490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1452490Sjkh if (ferror(stdin)) 1462490Sjkh err(1, "stdin"); 1472490Sjkh exit (0); 1482490Sjkh } 1492490Sjkh for (p = buf; isblank(*p); ++p); 1502490Sjkh if (*p == '\n' || *p == '\0') 1512490Sjkh continue; 1522490Sjkh if (*p == '-') 1532490Sjkh errx(1, "negative numbers aren't permitted."); 154104722Sfanf if (BN_dec2bn(&val, buf) == 0 && 155104722Sfanf BN_hex2bn(&val, buf) == 0) 1562490Sjkh errx(1, "%s: illegal numeric format.", buf); 1572490Sjkh pr_fact(val); 1582490Sjkh } 1592490Sjkh /* Factor the arguments. */ 1602490Sjkh else 1612490Sjkh for (; *argv != NULL; ++argv) { 1622490Sjkh if (argv[0][0] == '-') 1632490Sjkh errx(1, "negative numbers aren't permitted."); 164104722Sfanf if (BN_dec2bn(&val, argv[0]) == 0 && 165104722Sfanf BN_hex2bn(&val, argv[0]) == 0) 1662490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1672490Sjkh pr_fact(val); 1682490Sjkh } 1692490Sjkh exit(0); 1702490Sjkh} 1712490Sjkh 1722490Sjkh/* 1732490Sjkh * pr_fact - print the factors of a number 1742490Sjkh * 1752490Sjkh * Print the factors of the number, from the lowest to the highest. 17685408Sroam * A factor will be printed multiple times if it divides the value 1772490Sjkh * multiple times. 1782490Sjkh * 1792490Sjkh * Factors are printed with leading tabs. 1802490Sjkh */ 181104720Sfanfstatic void 182104722Sfanfpr_fact(BIGNUM *val) 1832490Sjkh{ 184104720Sfanf const ubig *fact; /* The factor found. */ 1852490Sjkh 1862490Sjkh /* Firewall - catch 0 and 1. */ 187104722Sfanf if (BN_is_zero(val)) /* Historical practice; 0 just exits. */ 1882490Sjkh exit(0); 189104722Sfanf if (BN_is_one(val)) { 190104720Sfanf printf("1: 1\n"); 1912490Sjkh return; 1922490Sjkh } 1932490Sjkh 1942490Sjkh /* Factor value. */ 195104722Sfanf 196104722Sfanf if (hflag) { 197104722Sfanf fputs("0x", stdout); 198104722Sfanf BN_print_fp(stdout, val); 199104722Sfanf } else 200104722Sfanf BN_print_dec_fp(stdout, val); 201104722Sfanf putchar(':'); 202104722Sfanf for (fact = &prime[0]; !BN_is_one(val); ++fact) { 2032490Sjkh /* Look for the smallest factor. */ 2042490Sjkh do { 205104722Sfanf if (BN_mod_word(val, (BN_ULONG)*fact) == 0) 2062490Sjkh break; 2072490Sjkh } while (++fact <= pr_limit); 2082490Sjkh 2092490Sjkh /* Watch for primes larger than the table. */ 2102490Sjkh if (fact > pr_limit) { 211104722Sfanf#ifdef HAVE_OPENSSL 212104722Sfanf BIGNUM *bnfact; 213104722Sfanf 214104722Sfanf bnfact = BN_new(); 215104722Sfanf BN_set_word(bnfact, *(fact - 1)); 216104722Sfanf BN_sqr(bnfact, bnfact, ctx); 217104722Sfanf if (BN_cmp(bnfact, val) > 0) 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 260104722Sfanf/* pollard rho, algorithm from Jim Gillogly, May 2000 */ 261104722Sfanfstatic void 262104722Sfanfpollard_pminus1(BIGNUM *val) 263104722Sfanf{ 264104722Sfanf BIGNUM *base, *num, *i, *x; 265104722Sfanf 266104722Sfanf base = BN_new(); 267104722Sfanf num = BN_new(); 268104722Sfanf i = BN_new(); 269104722Sfanf x = BN_new(); 270104722Sfanf 271104722Sfanf BN_set_word(i, 2); 272104722Sfanf BN_set_word(base, 2); 273104722Sfanf 274104722Sfanf for (;;) { 275104722Sfanf BN_mod_exp(base, base, i, val, ctx); 276104722Sfanf 277104722Sfanf BN_copy(x, base); 278104722Sfanf BN_sub_word(x, 1); 279104722Sfanf BN_gcd(x, x, val, ctx); 280104722Sfanf 281104722Sfanf if (!BN_is_one(x)) { 282104722Sfanf if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL, 283104722Sfanf NULL) == 1) 284104722Sfanf pr_print(x); 285104722Sfanf else 286104722Sfanf pollard_pminus1(x); 287104722Sfanf fflush(stdout); 288104722Sfanf 289104722Sfanf BN_div(num, NULL, val, x, ctx); 290104722Sfanf if (BN_is_one(num)) 291104722Sfanf return; 292104722Sfanf if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL, 293104722Sfanf NULL) == 1) { 294104722Sfanf pr_print(num); 295104722Sfanf fflush(stdout); 296104722Sfanf return; 297104722Sfanf } 298104722Sfanf BN_copy(val, num); 299104722Sfanf } 300104722Sfanf BN_add_word(i, 1); 301104722Sfanf } 302104722Sfanf} 303104722Sfanf 304104722Sfanf/* 305104722Sfanf * Sigh.. No _decimal_ output to file functions in BN. 306104722Sfanf */ 307104722Sfanfstatic void 308104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 309104722Sfanf{ 310104722Sfanf char *buf; 311104722Sfanf 312104722Sfanf buf = BN_bn2dec(num); 313104722Sfanf if (buf == NULL) 314104722Sfanf return; /* XXX do anything here? */ 315104722Sfanf fprintf(fp, buf); 316104722Sfanf free(buf); 317104722Sfanf} 318104722Sfanf 319104722Sfanf#else 320104722Sfanf 321104722Sfanfstatic void 322104722SfanfBN_print_fp(FILE *fp, const BIGNUM *num) 323104722Sfanf{ 324104722Sfanf fprintf(fp, "%lx", (unsigned long)*num); 325104722Sfanf} 326104722Sfanf 327104722Sfanfstatic void 328104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 329104722Sfanf{ 330104722Sfanf fprintf(fp, "%lu", (unsigned long)*num); 331104722Sfanf} 332104722Sfanf 333104722Sfanfstatic int 334104722SfanfBN_dec2bn(BIGNUM **a, const char *str) 335104722Sfanf{ 336104722Sfanf char *p; 337104722Sfanf 338104722Sfanf errno = 0; 339104722Sfanf **a = strtoul(str, &p, 10); 340104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 341104722Sfanf} 342104722Sfanf 343104722Sfanfstatic int 344104722SfanfBN_hex2bn(BIGNUM **a, const char *str) 345104722Sfanf{ 346104722Sfanf char *p; 347104722Sfanf 348104722Sfanf errno = 0; 349104722Sfanf **a = strtoul(str, &p, 16); 350104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 351104722Sfanf} 352104722Sfanf 353104722Sfanfstatic BN_ULONG 354104722SfanfBN_div_word(BIGNUM *a, BN_ULONG b) 355104722Sfanf{ 356104722Sfanf BN_ULONG mod; 357104722Sfanf 358104722Sfanf mod = *a % b; 359104722Sfanf *a /= b; 360104722Sfanf return mod; 361104722Sfanf} 362104722Sfanf 363104722Sfanf#endif 364