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$"); 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> 72272210Ssbruno#include <inttypes.h> 732490Sjkh#include <limits.h> 742490Sjkh#include <stdio.h> 752490Sjkh#include <stdlib.h> 7623726Speter#include <unistd.h> 772490Sjkh 782490Sjkh#include "primes.h" 792490Sjkh 80104722Sfanf#ifdef HAVE_OPENSSL 812490Sjkh 82104722Sfanf#include <openssl/bn.h> 83104722Sfanf 84104722Sfanf#define PRIME_CHECKS 5 85104722Sfanf 86104722Sfanfstatic void pollard_pminus1(BIGNUM *); /* print factors for big numbers */ 87104722Sfanf 88104722Sfanf#else 89104722Sfanf 90104722Sfanftypedef ubig BIGNUM; 91104722Sfanftypedef u_long BN_ULONG; 92104722Sfanf 93104722Sfanf#define BN_CTX int 94104722Sfanf#define BN_CTX_new() NULL 95104722Sfanf#define BN_new() ((BIGNUM *)calloc(sizeof(BIGNUM), 1)) 96104722Sfanf#define BN_is_zero(v) (*(v) == 0) 97104722Sfanf#define BN_is_one(v) (*(v) == 1) 98104722Sfanf#define BN_mod_word(a, b) (*(a) % (b)) 99104722Sfanf 100104722Sfanfstatic int BN_dec2bn(BIGNUM **a, const char *str); 101104722Sfanfstatic int BN_hex2bn(BIGNUM **a, const char *str); 102104722Sfanfstatic BN_ULONG BN_div_word(BIGNUM *, BN_ULONG); 103104722Sfanfstatic void BN_print_fp(FILE *, const BIGNUM *); 104104722Sfanf 105104722Sfanf#endif 106104722Sfanf 107104722Sfanfstatic void BN_print_dec_fp(FILE *, const BIGNUM *); 108104722Sfanf 109104722Sfanfstatic void pr_fact(BIGNUM *); /* print factors of a value */ 110104722Sfanfstatic void pr_print(BIGNUM *); /* print a prime */ 111104720Sfanfstatic void usage(void); 11242338Simp 113104722Sfanfstatic BN_CTX *ctx; /* just use a global context */ 114104722Sfanfstatic int hflag; 115104722Sfanf 1162490Sjkhint 117104720Sfanfmain(int argc, char *argv[]) 1182490Sjkh{ 119104722Sfanf BIGNUM *val; 1202490Sjkh int ch; 121104720Sfanf char *p, buf[LINE_MAX]; /* > max number of digits. */ 1222490Sjkh 123104722Sfanf ctx = BN_CTX_new(); 124104722Sfanf val = BN_new(); 125104722Sfanf if (val == NULL) 126104722Sfanf errx(1, "can't initialise bignum"); 127104722Sfanf 12842338Simp while ((ch = getopt(argc, argv, "h")) != -1) 1292490Sjkh switch (ch) { 13042338Simp case 'h': 13142338Simp hflag++; 13242338Simp break; 1332490Sjkh case '?': 1342490Sjkh default: 1352490Sjkh usage(); 1362490Sjkh } 1372490Sjkh argc -= optind; 1382490Sjkh argv += optind; 1392490Sjkh 1402490Sjkh /* No args supplied, read numbers from stdin. */ 1412490Sjkh if (argc == 0) 1422490Sjkh for (;;) { 1432490Sjkh if (fgets(buf, sizeof(buf), stdin) == NULL) { 1442490Sjkh if (ferror(stdin)) 1452490Sjkh err(1, "stdin"); 1462490Sjkh exit (0); 1472490Sjkh } 1482490Sjkh for (p = buf; isblank(*p); ++p); 1492490Sjkh if (*p == '\n' || *p == '\0') 1502490Sjkh continue; 1512490Sjkh if (*p == '-') 1522490Sjkh errx(1, "negative numbers aren't permitted."); 153104722Sfanf if (BN_dec2bn(&val, buf) == 0 && 154104722Sfanf BN_hex2bn(&val, buf) == 0) 1552490Sjkh errx(1, "%s: illegal numeric format.", buf); 1562490Sjkh pr_fact(val); 1572490Sjkh } 1582490Sjkh /* Factor the arguments. */ 1592490Sjkh else 1602490Sjkh for (; *argv != NULL; ++argv) { 1612490Sjkh if (argv[0][0] == '-') 1622490Sjkh errx(1, "negative numbers aren't permitted."); 163104722Sfanf if (BN_dec2bn(&val, argv[0]) == 0 && 164104722Sfanf BN_hex2bn(&val, argv[0]) == 0) 1652490Sjkh errx(1, "%s: illegal numeric format.", argv[0]); 1662490Sjkh pr_fact(val); 1672490Sjkh } 1682490Sjkh exit(0); 1692490Sjkh} 1702490Sjkh 1712490Sjkh/* 1722490Sjkh * pr_fact - print the factors of a number 1732490Sjkh * 1742490Sjkh * Print the factors of the number, from the lowest to the highest. 17585408Sroam * A factor will be printed multiple times if it divides the value 1762490Sjkh * multiple times. 1772490Sjkh * 1782490Sjkh * Factors are printed with leading tabs. 1792490Sjkh */ 180104720Sfanfstatic void 181104722Sfanfpr_fact(BIGNUM *val) 1822490Sjkh{ 183104720Sfanf const ubig *fact; /* The factor found. */ 1842490Sjkh 1852490Sjkh /* Firewall - catch 0 and 1. */ 186104722Sfanf if (BN_is_zero(val)) /* Historical practice; 0 just exits. */ 1872490Sjkh exit(0); 188104722Sfanf if (BN_is_one(val)) { 189104720Sfanf printf("1: 1\n"); 1902490Sjkh return; 1912490Sjkh } 1922490Sjkh 1932490Sjkh /* Factor value. */ 194104722Sfanf 195104722Sfanf if (hflag) { 196104722Sfanf fputs("0x", stdout); 197104722Sfanf BN_print_fp(stdout, val); 198104722Sfanf } else 199104722Sfanf BN_print_dec_fp(stdout, val); 200104722Sfanf putchar(':'); 201104722Sfanf for (fact = &prime[0]; !BN_is_one(val); ++fact) { 2022490Sjkh /* Look for the smallest factor. */ 2032490Sjkh do { 204104722Sfanf if (BN_mod_word(val, (BN_ULONG)*fact) == 0) 2052490Sjkh break; 2062490Sjkh } while (++fact <= pr_limit); 2072490Sjkh 2082490Sjkh /* Watch for primes larger than the table. */ 2092490Sjkh if (fact > pr_limit) { 210104722Sfanf#ifdef HAVE_OPENSSL 211104722Sfanf BIGNUM *bnfact; 212104722Sfanf 213104722Sfanf bnfact = BN_new(); 214104722Sfanf BN_set_word(bnfact, *(fact - 1)); 215216598Suqs if (!BN_sqr(bnfact, bnfact, ctx)) 216216598Suqs errx(1, "error in BN_sqr()"); 217199815Sfanf if (BN_cmp(bnfact, val) > 0 || 218199815Sfanf BN_is_prime(val, PRIME_CHECKS, 219199815Sfanf NULL, NULL, NULL) == 1) 220104722Sfanf pr_print(val); 221104722Sfanf else 222104722Sfanf pollard_pminus1(val); 223104722Sfanf#else 224104722Sfanf pr_print(val); 225104722Sfanf#endif 2262490Sjkh break; 2272490Sjkh } 2282490Sjkh 2292490Sjkh /* Divide factor out until none are left. */ 2302490Sjkh do { 231272210Ssbruno printf(hflag ? " 0x%" PRIx64 "" : " %" PRIu64 "", *fact); 232104722Sfanf BN_div_word(val, (BN_ULONG)*fact); 233104722Sfanf } while (BN_mod_word(val, (BN_ULONG)*fact) == 0); 2342490Sjkh 2352490Sjkh /* Let the user know we're doing something. */ 236104720Sfanf fflush(stdout); 2372490Sjkh } 238104720Sfanf putchar('\n'); 2392490Sjkh} 2402490Sjkh 241104720Sfanfstatic void 242104722Sfanfpr_print(BIGNUM *val) 243104722Sfanf{ 244104722Sfanf if (hflag) { 245104722Sfanf fputs(" 0x", stdout); 246104722Sfanf BN_print_fp(stdout, val); 247104722Sfanf } else { 248104722Sfanf putchar(' '); 249104722Sfanf BN_print_dec_fp(stdout, val); 250104722Sfanf } 251104722Sfanf} 252104722Sfanf 253104722Sfanfstatic void 254104720Sfanfusage(void) 2552490Sjkh{ 256104720Sfanf fprintf(stderr, "usage: factor [-h] [value ...]\n"); 257104720Sfanf exit(1); 2582490Sjkh} 259104722Sfanf 260104722Sfanf#ifdef HAVE_OPENSSL 261104722Sfanf 262199815Sfanf/* pollard p-1, algorithm from Jim Gillogly, May 2000 */ 263104722Sfanfstatic void 264104722Sfanfpollard_pminus1(BIGNUM *val) 265104722Sfanf{ 266199815Sfanf BIGNUM *base, *rbase, *num, *i, *x; 267104722Sfanf 268104722Sfanf base = BN_new(); 269199815Sfanf rbase = BN_new(); 270104722Sfanf num = BN_new(); 271104722Sfanf i = BN_new(); 272104722Sfanf x = BN_new(); 273104722Sfanf 274199815Sfanf BN_set_word(rbase, 1); 275199815Sfanfnewbase: 276216598Suqs if (!BN_add_word(rbase, 1)) 277216598Suqs errx(1, "error in BN_add_word()"); 278104722Sfanf BN_set_word(i, 2); 279199815Sfanf BN_copy(base, rbase); 280104722Sfanf 281104722Sfanf for (;;) { 282104722Sfanf BN_mod_exp(base, base, i, val, ctx); 283199815Sfanf if (BN_is_one(base)) 284199815Sfanf goto newbase; 285104722Sfanf 286104722Sfanf BN_copy(x, base); 287104722Sfanf BN_sub_word(x, 1); 288216598Suqs if (!BN_gcd(x, x, val, ctx)) 289216598Suqs errx(1, "error in BN_gcd()"); 290104722Sfanf 291104722Sfanf if (!BN_is_one(x)) { 292104722Sfanf if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL, 293104722Sfanf NULL) == 1) 294104722Sfanf pr_print(x); 295104722Sfanf else 296104722Sfanf pollard_pminus1(x); 297104722Sfanf fflush(stdout); 298104722Sfanf 299104722Sfanf BN_div(num, NULL, val, x, ctx); 300104722Sfanf if (BN_is_one(num)) 301104722Sfanf return; 302104722Sfanf if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL, 303104722Sfanf NULL) == 1) { 304104722Sfanf pr_print(num); 305104722Sfanf fflush(stdout); 306104722Sfanf return; 307104722Sfanf } 308104722Sfanf BN_copy(val, num); 309104722Sfanf } 310216598Suqs if (!BN_add_word(i, 1)) 311216598Suqs errx(1, "error in BN_add_word()"); 312104722Sfanf } 313104722Sfanf} 314104722Sfanf 315104722Sfanf/* 316104722Sfanf * Sigh.. No _decimal_ output to file functions in BN. 317104722Sfanf */ 318104722Sfanfstatic void 319104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 320104722Sfanf{ 321104722Sfanf char *buf; 322104722Sfanf 323104722Sfanf buf = BN_bn2dec(num); 324104722Sfanf if (buf == NULL) 325104722Sfanf return; /* XXX do anything here? */ 326228596Sdim fprintf(fp, "%s", buf); 327104722Sfanf free(buf); 328104722Sfanf} 329104722Sfanf 330104722Sfanf#else 331104722Sfanf 332104722Sfanfstatic void 333104722SfanfBN_print_fp(FILE *fp, const BIGNUM *num) 334104722Sfanf{ 335104722Sfanf fprintf(fp, "%lx", (unsigned long)*num); 336104722Sfanf} 337104722Sfanf 338104722Sfanfstatic void 339104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 340104722Sfanf{ 341104722Sfanf fprintf(fp, "%lu", (unsigned long)*num); 342104722Sfanf} 343104722Sfanf 344104722Sfanfstatic int 345104722SfanfBN_dec2bn(BIGNUM **a, const char *str) 346104722Sfanf{ 347104722Sfanf char *p; 348104722Sfanf 349104722Sfanf errno = 0; 350104722Sfanf **a = strtoul(str, &p, 10); 351104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 352104722Sfanf} 353104722Sfanf 354104722Sfanfstatic int 355104722SfanfBN_hex2bn(BIGNUM **a, const char *str) 356104722Sfanf{ 357104722Sfanf char *p; 358104722Sfanf 359104722Sfanf errno = 0; 360104722Sfanf **a = strtoul(str, &p, 16); 361104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 362104722Sfanf} 363104722Sfanf 364104722Sfanfstatic BN_ULONG 365104722SfanfBN_div_word(BIGNUM *a, BN_ULONG b) 366104722Sfanf{ 367104722Sfanf BN_ULONG mod; 368104722Sfanf 369104722Sfanf mod = *a % b; 370104722Sfanf *a /= b; 371104722Sfanf return mod; 372104722Sfanf} 373104722Sfanf 374104722Sfanf#endif 375