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> 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)); 214216598Suqs if (!BN_sqr(bnfact, bnfact, ctx)) 215216598Suqs errx(1, "error in BN_sqr()"); 216199815Sfanf if (BN_cmp(bnfact, val) > 0 || 217199815Sfanf BN_is_prime(val, PRIME_CHECKS, 218199815Sfanf NULL, NULL, NULL) == 1) 219104722Sfanf pr_print(val); 220104722Sfanf else 221104722Sfanf pollard_pminus1(val); 222104722Sfanf#else 223104722Sfanf pr_print(val); 224104722Sfanf#endif 2252490Sjkh break; 2262490Sjkh } 2272490Sjkh 2282490Sjkh /* Divide factor out until none are left. */ 2292490Sjkh do { 230104720Sfanf printf(hflag ? " 0x%lx" : " %lu", *fact); 231104722Sfanf BN_div_word(val, (BN_ULONG)*fact); 232104722Sfanf } while (BN_mod_word(val, (BN_ULONG)*fact) == 0); 2332490Sjkh 2342490Sjkh /* Let the user know we're doing something. */ 235104720Sfanf fflush(stdout); 2362490Sjkh } 237104720Sfanf putchar('\n'); 2382490Sjkh} 2392490Sjkh 240104720Sfanfstatic void 241104722Sfanfpr_print(BIGNUM *val) 242104722Sfanf{ 243104722Sfanf if (hflag) { 244104722Sfanf fputs(" 0x", stdout); 245104722Sfanf BN_print_fp(stdout, val); 246104722Sfanf } else { 247104722Sfanf putchar(' '); 248104722Sfanf BN_print_dec_fp(stdout, val); 249104722Sfanf } 250104722Sfanf} 251104722Sfanf 252104722Sfanfstatic void 253104720Sfanfusage(void) 2542490Sjkh{ 255104720Sfanf fprintf(stderr, "usage: factor [-h] [value ...]\n"); 256104720Sfanf exit(1); 2572490Sjkh} 258104722Sfanf 259104722Sfanf#ifdef HAVE_OPENSSL 260104722Sfanf 261199815Sfanf/* pollard p-1, algorithm from Jim Gillogly, May 2000 */ 262104722Sfanfstatic void 263104722Sfanfpollard_pminus1(BIGNUM *val) 264104722Sfanf{ 265199815Sfanf BIGNUM *base, *rbase, *num, *i, *x; 266104722Sfanf 267104722Sfanf base = BN_new(); 268199815Sfanf rbase = BN_new(); 269104722Sfanf num = BN_new(); 270104722Sfanf i = BN_new(); 271104722Sfanf x = BN_new(); 272104722Sfanf 273199815Sfanf BN_set_word(rbase, 1); 274199815Sfanfnewbase: 275216598Suqs if (!BN_add_word(rbase, 1)) 276216598Suqs errx(1, "error in BN_add_word()"); 277104722Sfanf BN_set_word(i, 2); 278199815Sfanf BN_copy(base, rbase); 279104722Sfanf 280104722Sfanf for (;;) { 281104722Sfanf BN_mod_exp(base, base, i, val, ctx); 282199815Sfanf if (BN_is_one(base)) 283199815Sfanf goto newbase; 284104722Sfanf 285104722Sfanf BN_copy(x, base); 286104722Sfanf BN_sub_word(x, 1); 287216598Suqs if (!BN_gcd(x, x, val, ctx)) 288216598Suqs errx(1, "error in BN_gcd()"); 289104722Sfanf 290104722Sfanf if (!BN_is_one(x)) { 291104722Sfanf if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL, 292104722Sfanf NULL) == 1) 293104722Sfanf pr_print(x); 294104722Sfanf else 295104722Sfanf pollard_pminus1(x); 296104722Sfanf fflush(stdout); 297104722Sfanf 298104722Sfanf BN_div(num, NULL, val, x, ctx); 299104722Sfanf if (BN_is_one(num)) 300104722Sfanf return; 301104722Sfanf if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL, 302104722Sfanf NULL) == 1) { 303104722Sfanf pr_print(num); 304104722Sfanf fflush(stdout); 305104722Sfanf return; 306104722Sfanf } 307104722Sfanf BN_copy(val, num); 308104722Sfanf } 309216598Suqs if (!BN_add_word(i, 1)) 310216598Suqs errx(1, "error in BN_add_word()"); 311104722Sfanf } 312104722Sfanf} 313104722Sfanf 314104722Sfanf/* 315104722Sfanf * Sigh.. No _decimal_ output to file functions in BN. 316104722Sfanf */ 317104722Sfanfstatic void 318104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 319104722Sfanf{ 320104722Sfanf char *buf; 321104722Sfanf 322104722Sfanf buf = BN_bn2dec(num); 323104722Sfanf if (buf == NULL) 324104722Sfanf return; /* XXX do anything here? */ 325228596Sdim fprintf(fp, "%s", buf); 326104722Sfanf free(buf); 327104722Sfanf} 328104722Sfanf 329104722Sfanf#else 330104722Sfanf 331104722Sfanfstatic void 332104722SfanfBN_print_fp(FILE *fp, const BIGNUM *num) 333104722Sfanf{ 334104722Sfanf fprintf(fp, "%lx", (unsigned long)*num); 335104722Sfanf} 336104722Sfanf 337104722Sfanfstatic void 338104722SfanfBN_print_dec_fp(FILE *fp, const BIGNUM *num) 339104722Sfanf{ 340104722Sfanf fprintf(fp, "%lu", (unsigned long)*num); 341104722Sfanf} 342104722Sfanf 343104722Sfanfstatic int 344104722SfanfBN_dec2bn(BIGNUM **a, const char *str) 345104722Sfanf{ 346104722Sfanf char *p; 347104722Sfanf 348104722Sfanf errno = 0; 349104722Sfanf **a = strtoul(str, &p, 10); 350104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 351104722Sfanf} 352104722Sfanf 353104722Sfanfstatic int 354104722SfanfBN_hex2bn(BIGNUM **a, const char *str) 355104722Sfanf{ 356104722Sfanf char *p; 357104722Sfanf 358104722Sfanf errno = 0; 359104722Sfanf **a = strtoul(str, &p, 16); 360104722Sfanf return (errno == 0 && (*p == '\n' || *p == '\0')); 361104722Sfanf} 362104722Sfanf 363104722Sfanfstatic BN_ULONG 364104722SfanfBN_div_word(BIGNUM *a, BN_ULONG b) 365104722Sfanf{ 366104722Sfanf BN_ULONG mod; 367104722Sfanf 368104722Sfanf mod = *a % b; 369104722Sfanf *a /= b; 370104722Sfanf return mod; 371104722Sfanf} 372104722Sfanf 373104722Sfanf#endif 374