-/* $NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $ */
+/* $NetBSD: factor.c,v 1.24 2010/05/15 21:22:39 joerg Exp $ */
/*
* Copyright (c) 1989, 1993
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- * must display the following acknowledgement:
- * This product includes software developed by the University of
- * California, Berkeley and its contributors.
- * 4. Neither the name of the University nor the names of its contributors
+ * 3. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
#include <sys/cdefs.h>
#ifndef lint
-__COPYRIGHT("@(#) Copyright (c) 1989, 1993\n\
- The Regents of the University of California. All rights reserved.\n");
+__COPYRIGHT("@(#) Copyright (c) 1989, 1993\
+ The Regents of the University of California. All rights reserved.");
#endif /* not lint */
#ifndef lint
#if 0
static char sccsid[] = "@(#)factor.c 8.4 (Berkeley) 5/4/95";
#else
-__RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $");
+__RCSID("$NetBSD: factor.c,v 1.24 2010/05/15 21:22:39 joerg Exp $");
#endif
#endif /* not lint */
*
* number: factor1 factor1 factor2 factor3 factor3 factor3 ...
*
- * where factor1 < factor2 < factor3 < ...
+ * where factor1 <= factor2 <= factor3 <= ...
*
* If no args are given, the list of numbers are read from stdin.
*/
#else
typedef long BIGNUM;
typedef u_long BN_ULONG;
-int BN_dec2bn(BIGNUM **a, const char *str);
+static int BN_dec2bn(BIGNUM **a, const char *str);
#define BN_new() ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
#define BN_is_zero(v) (*(v) == 0)
#define BN_is_one(v) (*(v) == 1)
*/
extern const ubig prime[];
extern const ubig *pr_limit; /* largest prime in the prime array */
+#if 0 /* debugging: limit table use to stress the "pollard" code */
+#define pr_limit &prime[0]
+#endif
#define PRIME_CHECKS 5
#ifdef HAVE_OPENSSL
-BN_CTX *ctx; /* just use a global context */
+static BN_CTX *ctx; /* just use a global context */
#endif
-int main(int, char *[]);
-void pr_fact(BIGNUM *); /* print factors of a value */
-void BN_print_dec_fp(FILE *, const BIGNUM *);
-void usage(void) __attribute__((__noreturn__));
+static void pr_fact(BIGNUM *); /* print factors of a value */
+static void BN_print_dec_fp(FILE *, const BIGNUM *);
+static void usage(void) __dead;
#ifdef HAVE_OPENSSL
-void pollard_pminus1(BIGNUM *); /* print factors for big numbers */
+static void pollard_rho(BIGNUM *); /* print factors for big numbers */
#else
-char *BN_bn2dec(const BIGNUM *);
-BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
+static char *BN_bn2dec(const BIGNUM *);
+static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
#endif
#ifndef HAVE_OPENSSL
-int
+static int
BN_dec2bn(BIGNUM **a, const char *str)
{
char *p;
err(1, "stdin");
exit (0);
}
- for (p = buf; isblank(*p); ++p);
+ for (p = buf; isblank((unsigned char)*p); ++p);
if (*p == '\n' || *p == '\0')
continue;
if (*p == '-')
else
for (; *argv != NULL; ++argv) {
if (argv[0][0] == '-')
- errx(1, "negative numbers aren't permitted.");
+ errx(1, "numbers <= 1 aren't permitted.");
if (BN_dec2bn(&val, argv[0]) == 0)
errx(1, "%s: illegal numeric format.", argv[0]);
pr_fact(val);
*
* Factors are printed with leading tabs.
*/
-void
+static void
pr_fact(BIGNUM *val)
{
const ubig *fact; /* The factor found. */
/* Firewall - catch 0 and 1. */
- if (BN_is_zero(val)) /* Historical practice; 0 just exits. */
- exit(0);
- if (BN_is_one(val)) {
- printf("1: 1\n");
- return;
- }
+ if (BN_is_zero(val) || BN_is_one(val))
+ errx(1, "numbers <= 1 aren't permitted.");
/* Factor value. */
putchar(':');
for (fact = &prime[0]; !BN_is_one(val); ++fact) {
/* Look for the smallest factor. */
- do {
+ while (fact <= pr_limit) {
if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
break;
- } while (++fact <= pr_limit);
+ fact++;
+ }
/* Watch for primes larger than the table. */
if (fact > pr_limit) {
BIGNUM *bnfact;
bnfact = BN_new();
- BN_set_word(bnfact, *(fact - 1));
+ BN_set_word(bnfact, (BN_ULONG)*(fact - 1));
BN_sqr(bnfact, bnfact, ctx);
- if (BN_cmp(bnfact, val) > 0) {
+ if (BN_cmp(bnfact, val) > 0
+ || BN_is_prime(val, PRIME_CHECKS, NULL, NULL,
+ NULL) == 1) {
putchar(' ');
BN_print_dec_fp(stdout, val);
} else
- pollard_pminus1(val);
+ pollard_rho(val);
#else
printf(" %s", BN_bn2dec(val));
#endif
/*
* Sigh.. No _decimal_ output to file functions in BN.
*/
-void
+static void
BN_print_dec_fp(FILE *fp, const BIGNUM *num)
{
char *buf;
#ifdef HAVE_OPENSSL
-/* pollard rho, algorithm from Jim Gillogly, May 2000 */
-
-void
-pollard_pminus1(BIGNUM *val)
+static void
+pollard_rho(BIGNUM *val)
{
- BIGNUM *base, *num, *i, *x;
+ BIGNUM *x, *y, *tmp, *num;
+ BN_ULONG a;
+ unsigned int steps_taken, steps_limit;
- base = BN_new();
- num = BN_new();
- i = BN_new();
x = BN_new();
-
- BN_set_word(i, 2);
- BN_set_word(base, 2);
+ y = BN_new();
+ tmp = BN_new();
+ num = BN_new();
+ a = 1;
+restart:
+ steps_taken = 0;
+ steps_limit = 2;
+ BN_set_word(x, 1);
+ BN_copy(y, x);
for (;;) {
- BN_mod_exp(base, base, i, val, ctx);
-
- BN_copy(x, base);
- BN_sub_word(x, 1);
- BN_gcd(x, x, val, ctx);
+ BN_sqr(tmp, x, ctx);
+ BN_add_word(tmp, a);
+ BN_mod(x, tmp, val, ctx);
+ BN_sub(tmp, x, y);
+ if (BN_is_zero(tmp)) {
+#ifdef DEBUG
+ printf(" (loop)");
+#endif
+ a++;
+ goto restart;
+ }
+ BN_gcd(tmp, tmp, val, ctx);
- if (!BN_is_one(x)) {
- if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
+ if (!BN_is_one(tmp)) {
+ if (BN_is_prime(tmp, PRIME_CHECKS, NULL, NULL,
NULL) == 1) {
putchar(' ');
- BN_print_dec_fp(stdout, x);
- } else
- pollard_pminus1(x);
+ BN_print_dec_fp(stdout, tmp);
+ } else {
+#ifdef DEBUG
+ printf(" (recurse for ");
+ BN_print_dec_fp(stdout, tmp);
+ putchar(')');
+#endif
+ pollard_rho(BN_dup(tmp));
+#ifdef DEBUG
+ printf(" (back)");
+#endif
+ }
fflush(stdout);
- BN_div(num, NULL, val, x, ctx);
+ BN_div(num, NULL, val, tmp, ctx);
if (BN_is_one(num))
return;
if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
return;
}
BN_copy(val, num);
+ goto restart;
+ }
+ steps_taken++;
+ if (steps_taken == steps_limit) {
+ BN_copy(y, x); /* teleport the turtle */
+ steps_taken = 0;
+ steps_limit *= 2;
+ if (steps_limit == 0) {
+#ifdef DEBUG
+ printf(" (overflow)");
+#endif
+ a++;
+ goto restart;
+ }
}
- BN_add_word(i, 1);
}
}
#else
return buf;
}
-BN_ULONG
+static BN_ULONG
BN_div_word(BIGNUM *a, BN_ULONG b)
{
BN_ULONG mod;