Below is my proposed patch to primes(6) and factor(6) which I plan to commit in one go since the changes are somewhat inter-dependent. Feedback is welcomed. I'm in the process of fixing the manual.
Merge changes from NetBSD and perform some cleaning up. primes: const-correctness and removal of wacky casting (from NetBSD); move declarations of external tables into primes.h instead of repeating them in primes.c and factor.c. Still limited to ULONG_MAX though. factor: style fixes (#include ordering, ANSI functions, const-correctness, staticisation -- all but the latter from NetBSD); remove bogus comment; fix usage() exit value; and the biggie: if OpenSSL is available, use bignums using clever code from NetBSD. I have cleaned the latter up somewhat so that it supports FreeBSD's -h feature and doesn't introduce regressions for N >= 2^31 in the non-OpenSSL case. Tony. -- f.a.n.finch <[EMAIL PROTECTED]> http://dotat.at/ PORTLAND: EAST OR SOUTHEAST BACKING NORTHEAST 5 OR 6, OCCASIONALLY 7. OCCASIONAL RAIN. GOOD. --- factor/Makefile 26 Mar 2001 14:20:55 -0000 1.4 +++ factor/Makefile 8 Oct 2002 19:31:03 -0000 @@ -4,6 +4,13 @@ PROG= factor SRCS= factor.c pr_tbl.c CFLAGS+=-I${.CURDIR}/../primes + +.if !defined(NO_OPENSSL) +CFLAGS+=-DHAVE_OPENSSL +LDADD+= -lcrypto +DPADD+= ${LIBCRYPTO} +.endif + MAN= factor.6 MLINKS+=factor.6 primes.6 .PATH: ${.CURDIR}/../primes --- factor/factor.c 18 Feb 2002 05:15:15 -0000 1.11 +++ factor/factor.c 9 Oct 2002 15:28:21 -0000 @@ -43,6 +43,7 @@ #ifndef lint #if 0 static char sccsid[] = "@(#)factor.c 8.4 (Berkeley) 5/4/95"; +__RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $"); #endif static const char rcsid[] = "$FreeBSD: src/games/factor/factor.c,v 1.11 2002/02/18 05:15:15 imp Exp $"; @@ -67,8 +68,8 @@ * If no args are given, the list of numbers are read from stdin. */ -#include <err.h> #include <ctype.h> +#include <err.h> #include <errno.h> #include <limits.h> #include <stdio.h> @@ -77,28 +78,53 @@ #include "primes.h" -/* - * prime[i] is the (i-1)th prime. - * - * We are able to sieve 2^32-1 because this byte table yields all primes - * up to 65537 and 65537^2 > 2^32-1. - */ -extern ubig prime[]; -extern ubig *pr_limit; /* largest prime in the prime array */ +#ifdef HAVE_OPENSSL + +#include <openssl/bn.h> + +#define PRIME_CHECKS 5 + +static void pollard_pminus1(BIGNUM *); /* print factors for big numbers */ + +#else + +typedef ubig BIGNUM; +typedef u_long BN_ULONG; + +#define BN_CTX int +#define BN_CTX_new() NULL +#define BN_new() ((BIGNUM *)calloc(sizeof(BIGNUM), 1)) +#define BN_is_zero(v) (*(v) == 0) +#define BN_is_one(v) (*(v) == 1) +#define BN_mod_word(a, b) (*(a) % (b)) + +static int BN_dec2bn(BIGNUM **a, const char *str); +static int BN_hex2bn(BIGNUM **a, const char *str); +static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG); +static void BN_print_fp(FILE *, const BIGNUM *); + +#endif + +static void BN_print_dec_fp(FILE *, const BIGNUM *); -int hflag; +static void pr_fact(BIGNUM *); /* print factors of a value */ +static void pr_print(BIGNUM *); /* print a prime */ +static void usage(void); -void pr_fact(ubig); /* print factors of a value */ -void usage(void); +static BN_CTX *ctx; /* just use a global context */ +static int hflag; int -main(argc, argv) - int argc; - char *argv[]; +main(int argc, char *argv[]) { - ubig val; + BIGNUM *val; int ch; - char *p, buf[100]; /* > max number of digits. */ + char *p, buf[LINE_MAX]; /* > max number of digits. */ + + ctx = BN_CTX_new(); + val = BN_new(); + if (val == NULL) + errx(1, "can't initialise bignum"); while ((ch = getopt(argc, argv, "h")) != -1) switch (ch) { @@ -125,12 +151,9 @@ continue; if (*p == '-') errx(1, "negative numbers aren't permitted."); - errno = 0; - val = strtoul(buf, &p, 0); - if (errno) - err(1, "%s", buf); - if (*p != '\n') - errx(1, "%s: illegal numeric format.", buf); + if (BN_dec2bn(&val, buf) == 0 && + BN_hex2bn(&val, buf) == 0) + errx(1, "%s: illegal numeric format.", argv[0]); pr_fact(val); } /* Factor the arguments. */ @@ -138,11 +161,8 @@ for (; *argv != NULL; ++argv) { if (argv[0][0] == '-') errx(1, "negative numbers aren't permitted."); - errno = 0; - val = strtoul(argv[0], &p, 0); - if (errno) - err(1, "%s", argv[0]); - if (*p != '\0') + if (BN_dec2bn(&val, argv[0]) == 0 && + BN_hex2bn(&val, argv[0]) == 0) errx(1, "%s: illegal numeric format.", argv[0]); pr_fact(val); } @@ -152,60 +172,192 @@ /* * pr_fact - print the factors of a number * - * If the number is 0 or 1, then print the number and return. - * If the number is < 0, print -1, negate the number and continue - * processing. - * * Print the factors of the number, from the lowest to the highest. * A factor will be printed multiple times if it divides the value * multiple times. * * Factors are printed with leading tabs. */ -void -pr_fact(val) - ubig val; /* Factor this value. */ +static void +pr_fact(BIGNUM *val) { - ubig *fact; /* The factor found. */ + const ubig *fact; /* The factor found. */ /* Firewall - catch 0 and 1. */ - if (val == 0) /* Historical practice; 0 just exits. */ + if (BN_is_zero(val)) /* Historical practice; 0 just exits. */ exit(0); - if (val == 1) { - (void)printf("1: 1\n"); + if (BN_is_one(val)) { + printf("1: 1\n"); return; } /* Factor value. */ - (void)printf(hflag ? "0x%lx:" : "%lu:", val); - for (fact = &prime[0]; val > 1; ++fact) { + + if (hflag) { + fputs("0x", stdout); + BN_print_fp(stdout, val); + } else + BN_print_dec_fp(stdout, val); + putchar(':'); + for (fact = &prime[0]; !BN_is_one(val); ++fact) { /* Look for the smallest factor. */ do { - if (val % (long)*fact == 0) + if (BN_mod_word(val, (BN_ULONG)*fact) == 0) break; } while (++fact <= pr_limit); /* Watch for primes larger than the table. */ if (fact > pr_limit) { - (void)printf(hflag ? " 0x%lx" : " %lu", val); +#ifdef HAVE_OPENSSL + BIGNUM *bnfact; + + bnfact = BN_new(); + BN_set_word(bnfact, *(fact - 1)); + BN_sqr(bnfact, bnfact, ctx); + if (BN_cmp(bnfact, val) > 0) + pr_print(val); + else + pollard_pminus1(val); +#else + pr_print(val); +#endif break; } /* Divide factor out until none are left. */ do { - (void)printf(hflag ? " 0x%lx" : " %lu", *fact); - val /= *fact; - } while ((val % *fact) == 0); + printf(hflag ? " 0x%lx" : " %lu", *fact); + BN_div_word(val, (BN_ULONG)*fact); + } while (BN_mod_word(val, (BN_ULONG)*fact) == 0); /* Let the user know we're doing something. */ - (void)fflush(stdout); + fflush(stdout); + } + putchar('\n'); +} + +static void +pr_print(BIGNUM *val) +{ + if (hflag) { + fputs(" 0x", stdout); + BN_print_fp(stdout, val); + } else { + putchar(' '); + BN_print_dec_fp(stdout, val); + } +} + +static void +usage(void) +{ + fprintf(stderr, "usage: factor [-h] [value ...]\n"); + exit(1); +} + +#ifdef HAVE_OPENSSL + +/* pollard rho, algorithm from Jim Gillogly, May 2000 */ +static void +pollard_pminus1(BIGNUM *val) +{ + BIGNUM *base, *num, *i, *x; + + base = BN_new(); + num = BN_new(); + i = BN_new(); + x = BN_new(); + + BN_set_word(i, 2); + BN_set_word(base, 2); + + for (;;) { + BN_mod_exp(base, base, i, val, ctx); + + BN_copy(x, base); + BN_sub_word(x, 1); + BN_gcd(x, x, val, ctx); + + if (!BN_is_one(x)) { + if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL, + NULL) == 1) + pr_print(x); + else + pollard_pminus1(x); + fflush(stdout); + + BN_div(num, NULL, val, x, ctx); + if (BN_is_one(num)) + return; + if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL, + NULL) == 1) { + pr_print(num); + fflush(stdout); + return; + } + BN_copy(val, num); + } + BN_add_word(i, 1); } - (void)putchar('\n'); } -void -usage() +/* + * Sigh.. No _decimal_ output to file functions in BN. + */ +static void +BN_print_dec_fp(FILE *fp, const BIGNUM *num) +{ + char *buf; + + buf = BN_bn2dec(num); + if (buf == NULL) + return; /* XXX do anything here? */ + fprintf(fp, buf); + free(buf); +} + +#else + +static void +BN_print_fp(FILE *fp, const BIGNUM *num) +{ + fprintf(fp, "%lx", (unsigned long)*num); +} + +static void +BN_print_dec_fp(FILE *fp, const BIGNUM *num) +{ + fprintf(fp, "%lu", (unsigned long)*num); +} + +static int +BN_dec2bn(BIGNUM **a, const char *str) { - (void)fprintf(stderr, "usage: factor -h [value ...]\n"); - exit (0); + char *p; + + errno = 0; + **a = strtoul(str, &p, 10); + return (errno == 0 && (*p == '\n' || *p == '\0')); } + +static int +BN_hex2bn(BIGNUM **a, const char *str) +{ + char *p; + + errno = 0; + **a = strtoul(str, &p, 16); + return (errno == 0 && (*p == '\n' || *p == '\0')); +} + +static BN_ULONG +BN_div_word(BIGNUM *a, BN_ULONG b) +{ + BN_ULONG mod; + + mod = *a % b; + *a /= b; + return mod; +} + +#endif --- primes/pattern.c 21 Feb 2002 18:13:31 -0000 1.4 +++ primes/pattern.c 8 Oct 2002 20:14:47 -0000 @@ -56,7 +56,9 @@ #include <stddef.h> -char pattern[] = { +#include "primes.h" + +const char pattern[] = { 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,0,1,1,0,0,1,0,1,1,0,0, 1,0,1,0,0,1,0,0,0,1,0,1,1,0,1,1,0,1,0,0,0,0,0,0,1,0,1,0,0,1,1,0,0,0,0,1,1,0,0, 1,0,0,1,0,1,0,0,1,0,0,1,1,0,0,0,0,1,1,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,1, @@ -443,4 +445,4 @@ 0,0,1,1,0,0,0,0,1,1,0,0,1,0,1,0,0,0,0,0,0,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0,1,0,1, 0,0,1,1,0,1,0,0,1,1,0,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,0,0,0,0,0,0,1 }; -size_t pattern_size = (sizeof(pattern)/sizeof(pattern[0])); +const size_t pattern_size = (sizeof(pattern)/sizeof(pattern[0])); --- primes/pr_tbl.c 21 Feb 2002 17:33:56 -0000 1.4 +++ primes/pr_tbl.c 8 Oct 2002 20:14:58 -0000 @@ -53,9 +53,11 @@ * and 65537^2 > 2^32-1. */ +#include <stddef.h> + #include "primes.h" -ubig prime[] = { +const ubig prime[] = { 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103, 107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199, 211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313, @@ -547,4 +549,4 @@ }; /* pr_limit - largest prime in the prime table */ -ubig *pr_limit = &prime[(sizeof(prime)/sizeof(prime[0]))-1]; +const ubig *const pr_limit = &prime[(sizeof(prime)/sizeof(prime[0]))-1]; --- primes/primes.c 21 Feb 2002 18:13:31 -0000 1.18 +++ primes/primes.c 9 Oct 2002 15:02:05 -0000 @@ -56,7 +56,7 @@ * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\ * * usage: - * primes [start [stop]] + * primes [-h] [start [stop]] * * Print primes >= start and < stop. If stop is omitted, * the value 4294967295 (2^32-1) is assumed. If start is @@ -88,23 +88,6 @@ */ static char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ -/* - * prime[i] is the (i-1)th prime. - * - * We are able to sieve 2^32-1 because this byte table yields all primes - * up to 65537 and 65537^2 > 2^32-1. - */ -extern ubig prime[]; -extern ubig *pr_limit; /* largest prime in the prime array */ - -/* - * To avoid excessive sieves for small factors, we use the table below to - * setup our sieve blocks. Each element represents a odd number starting - * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13. - */ -extern char pattern[]; -extern size_t pattern_size; /* length of pattern array */ - static int hflag; static void primes(ubig, ubig); @@ -225,8 +208,9 @@ char *q; /* sieve spot */ ubig factor; /* index and factor */ char *tab_lim; /* the limit to sieve on the table */ - ubig *p; /* prime table pointer */ + const ubig *p; /* prime table pointer */ ubig fact_lim; /* highest prime for current block */ + ubig mod; /* temp storage for mod */ /* * A number of systems can not convert double values into unsigned @@ -296,22 +280,21 @@ /* note highest useful factor and sieve spot */ if (stop-start > TABSIZE+TABSIZE) { tab_lim = &table[TABSIZE]; /* sieve it all */ - fact_lim = (int)sqrt( - (double)(start)+TABSIZE+TABSIZE+1.0); + fact_lim = sqrt(start+1.0+TABSIZE+TABSIZE); } else { tab_lim = &table[(stop-start)/2]; /* partial sieve */ - fact_lim = (int)sqrt((double)(stop)+1.0); + fact_lim = sqrt(stop+1.0); } /* sieve for factors >= 17 */ factor = 17; /* 17 is first prime to use */ p = &prime[7]; /* 19 is next prime, pi(19)=7 */ do { /* determine the factor's initial sieve point */ - q = (char *)(start%factor); /* temp storage for mod */ - if ((long)q & 0x1) { - q = &table[(factor-(long)q)/2]; + mod = start%factor; + if (mod & 0x1) { + q = &table[(factor-mod)/2]; } else { - q = &table[q ? factor-((long)q/2) : 0]; + q = &table[mod ? factor-(mod/2) : 0]; } /* sive for our current factor */ for ( ; q < tab_lim; q += factor) { @@ -333,6 +316,6 @@ static void usage(void) { - (void)fprintf(stderr, "usage: primes [-h] [start [stop]]\n"); + fprintf(stderr, "usage: primes [-h] [start [stop]]\n"); exit(1); } --- primes/primes.h 4 Sep 1994 04:03:08 -0000 1.1.1.1 +++ primes/primes.h 8 Oct 2002 20:11:29 -0000 @@ -50,3 +50,20 @@ /* bytes in sieve table (must be > 3*5*7*11) */ #define TABSIZE 256*1024 + +/* + * prime[i] is the (i-1)th prime. + * + * We are able to sieve 2^32-1 because this byte table yields all primes + * up to 65537 and 65537^2 > 2^32-1. + */ +extern const ubig prime[]; +extern const ubig *const pr_limit; /* largest prime in the prime array */ + +/* + * To avoid excessive sieves for small factors, we use the table below to + * setup our sieve blocks. Each element represents a odd number starting + * with 1. All non-zero elements are factors of 3, 5, 7, 11 and 13. + */ +extern const char pattern[]; +extern const size_t pattern_size; /* length of pattern array */ To Unsubscribe: send mail to [EMAIL PROTECTED] with "unsubscribe freebsd-current" in the body of the message
