If primes(6) will generate 32-bit primes, factor should factor 64-bit numbers.

3 files changed, 112 insertions, 13 deletions

diff --git a/games/factor/Makefile b/games/factor/Makefile
index 361bf56aa01..e5b7697eee5 100644
--- a/games/factor/Makefile
+++ b/games/factor/Makefile
@@ -1,10 +1,12 @@
-#	$OpenBSD: Makefile,v 1.2 1997/09/21 11:35:51 deraadt Exp $
+#	$OpenBSD: Makefile,v 1.3 2001/08/19 23:19:46 pjanzen Exp $
 
 PROG=	factor
-SRCS=	factor.c pr_tbl.c
+SRCS=	factor.c pattern.c pr_tbl.c
 CFLAGS+=-I${.CURDIR}/../primes
 MAN=	factor.6
 MLINKS+=factor.6 primes.6
+DPADD=  ${LIBMATH}
+LDADD=  -lm
 .PATH:	${.CURDIR}/../primes
 
 .include <bsd.prog.mk>
diff --git a/games/factor/factor.6 b/games/factor/factor.6
index 9859c9d163d..5aaaa4585ed 100644
--- a/games/factor/factor.6
+++ b/games/factor/factor.6
@@ -1,4 +1,4 @@
-.\"	$OpenBSD: factor.6,v 1.7 2001/08/19 16:55:02 pjanzen Exp $
+.\"	$OpenBSD: factor.6,v 1.8 2001/08/19 23:19:46 pjanzen Exp $
 .\"
 .\" Copyright (c) 1989, 1993
 .\"	The Regents of the University of California.  All rights reserved.
@@ -57,7 +57,7 @@
 .Sh DESCRIPTION
 The
 .Nm
-utility will factor integers between 1 and 4294967295 inclusive.
+utility will factor 64-bit positive integers.
 When a number is factored, it is printed, followed by a
 .Dq \: ,
 and the list of factors on a single line.
diff --git a/games/factor/factor.c b/games/factor/factor.c
index 2f33b48bbe3..feab8290889 100644
--- a/games/factor/factor.c
+++ b/games/factor/factor.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: factor.c,v 1.9 2001/08/19 17:12:40 pjanzen Exp $	*/
+/*	$OpenBSD: factor.c,v 1.10 2001/08/19 23:19:46 pjanzen Exp $	*/
 /*	$NetBSD: factor.c,v 1.5 1995/03/23 08:28:07 cgd Exp $	*/
 
 /*
@@ -47,7 +47,7 @@ static char copyright[] =
 #if 0
 static char sccsid[] = "@(#)factor.c	8.4 (Berkeley) 5/4/95";
 #else
-static char rcsid[] = "$OpenBSD: factor.c,v 1.9 2001/08/19 17:12:40 pjanzen Exp $";
+static char rcsid[] = "$OpenBSD: factor.c,v 1.10 2001/08/19 23:19:46 pjanzen Exp $";
 #endif
 #endif /* not lint */
 
@@ -75,8 +75,10 @@ static char rcsid[] = "$OpenBSD: factor.c,v 1.9 2001/08/19 17:12:40 pjanzen Exp
 #include <ctype.h>
 #include <errno.h>
 #include <limits.h>
+#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include <unistd.h>
 
 #include "primes.h"
@@ -89,8 +91,11 @@ static char rcsid[] = "$OpenBSD: factor.c,v 1.9 2001/08/19 17:12:40 pjanzen Exp
  */
 extern const ubig prime[];
 extern const ubig *pr_limit;		/* largest prime in the prime array */
+extern const char pattern[];
+extern const int pattern_size;
 
-void	pr_fact __P((ubig));	/* print factors of a value */
+void	pr_fact __P((u_int64_t));	/* print factors of a value */
+void	pr_bigfact __P((u_int64_t));
 void	usage __P((void));
 
 int
@@ -98,7 +103,7 @@ main(argc, argv)
 	int argc;
 	char *argv[];
 {
-	ubig val;
+	u_int64_t val;
 	int ch;
 	char *p, buf[100];		/* > max number of digits. */
 
@@ -131,7 +136,7 @@ main(argc, argv)
 			if (*p == '-')
 				errx(1, "negative numbers aren't permitted.");
 			errno = 0;
-			val = strtoul(buf, &p, 10);
+			val = strtouq(buf, &p, 10);
 			if (errno)
 				err(1, "%s", buf);
 			for (; isblank(*p); ++p);
@@ -145,7 +150,7 @@ main(argc, argv)
 			if (argv[0][0] == '-')
 				errx(1, "negative numbers aren't permitted.");
 			errno = 0;
-			val = strtoul(argv[0], &p, 10);
+			val = strtouq(argv[0], &p, 10);
 			if (errno)
 				err(1, "%s", argv[0]);
 			if (*p != '\0')
@@ -170,7 +175,7 @@ main(argc, argv)
  */
 void
 pr_fact(val)
-	ubig val;		/* Factor this value. */
+	u_int64_t val;		/* Factor this value. */
 {
 	const ubig *fact;	/* The factor found. */
 
@@ -183,7 +188,8 @@ pr_fact(val)
 	}
 
 	/* Factor value. */
-	(void)printf("%lu:", (unsigned long) val);
+	(void)printf("%llu:", val);
+	fflush(stdout);
 	for (fact = &prime[0]; val > 1; ++fact) {
 		/* Look for the smallest factor. */
 		do {
@@ -193,7 +199,10 @@ pr_fact(val)
 
 		/* Watch for primes larger than the table. */
 		if (fact > pr_limit) {
-			(void)printf(" %lu", (unsigned long) val);
+			if (val > BIG)
+				pr_bigfact(val);
+			else
+				(void)printf(" %llu", val);
 			break;
 		}
 
@@ -209,6 +218,94 @@ pr_fact(val)
 	(void)putchar('\n');
 }
 
+
+/* At this point, our number may have factors greater than those in primes[];
+ * however, we can generate primes up to 32 bits (see primes(6)), which is
+ * sufficient to factor a 64-bit quad.
+ */
+void
+pr_bigfact(val)
+	u_int64_t val;		/* Factor this value. */
+{
+	ubig start, stop, factor;
+	char *q;
+	const ubig *p;
+	ubig fact_lim, mod;
+	char *tab_lim;
+	char table[TABSIZE];	/* Eratosthenes sieve of odd numbers */
+
+	start = *pr_limit + 2;
+	stop  = (ubig)sqrt((double)val);
+	if ((stop & 0x1) == 0)
+		stop++;
+	/*
+	 * Following code barely modified from that in primes(6)
+	 *
+	 * we shall sieve a bytemap window, note primes and move the window
+	 * upward until we pass the stop point
+	 */
+	while (start < stop) {
+		/*
+		 * factor out 3, 5, 7, 11 and 13
+		 */
+		/* initial pattern copy */
+		factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */
+		memcpy(table, &pattern[factor], pattern_size-factor);
+		/* main block pattern copies */
+		for (fact_lim = pattern_size - factor;
+		    fact_lim + pattern_size <= TABSIZE; fact_lim += pattern_size) {
+			memcpy(&table[fact_lim], pattern, pattern_size);
+		}
+		/* final block pattern copy */
+		memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim);
+
+		if (stop-start > TABSIZE+TABSIZE) {
+			tab_lim = &table[TABSIZE]; /* sieve it all */
+			fact_lim = (int)sqrt(
+					(double)(start)+TABSIZE+TABSIZE+1.0);
+		} else {
+			tab_lim = &table[(stop - start)/2]; /* partial sieve */
+			fact_lim = (int)sqrt((double)(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 */
+			mod = start % factor;
+			if (mod & 0x1)
+				q = &table[(factor-mod)/2];
+			else
+				q = &table[mod ? factor-(mod/2) : 0];
+			/* sieve for our current factor */
+			for ( ; q < tab_lim; q += factor) {
+				*q = '\0'; /* sieve out a spot */
+			}
+		} while ((factor=(ubig)(*(p++))) <= fact_lim);
+
+		/*
+		 * use generated primes
+		 */
+		for (q = table; q < tab_lim; ++q, start+=2) {
+			if (*q) {
+				if (val % start == 0) {
+					do {
+						(void)printf(" %lu", (unsigned long) start);
+						val /= start;
+					} while ((val % start) == 0);
+					(void)fflush(stdout);
+					stop  = (ubig)sqrt((double)val);
+					if ((stop & 0x1) == 0)
+						stop++;
+				}
+			}
+		}
+	}
+	if (val > 1)
+		printf(" %llu", val);
+}
+
+
 void
 usage()
 {