From 62b9314c61c9071ddbc1ae28aea86872b82251e0 Mon Sep 17 00:00:00 2001 From: tb Date: Mon, 11 Jul 2016 18:30:21 +0000 Subject: [PATCH] Instead of using the floating point square root, use an integer version of the Newton method from ping.c. Fixes a rounding issue that caused failure to factor numbers close to 2^64, e.g. 18446744030759878681. While there, fix an off by one error that caused 4295360521 to be reported as a prime. Issues reported by Paul Stoeber and Michael Bozon. ok tedu, deraadt --- games/factor/Makefile | 4 +-- games/factor/factor.c | 57 ++++++++++++++++++++++++++++--------------- 2 files changed, 39 insertions(+), 22 deletions(-) diff --git a/games/factor/Makefile b/games/factor/Makefile index b38a839511d..261e7a45b19 100644 --- a/games/factor/Makefile +++ b/games/factor/Makefile @@ -1,11 +1,9 @@ -# $OpenBSD: Makefile,v 1.5 2016/03/30 06:38:40 jmc Exp $ +# $OpenBSD: Makefile,v 1.6 2016/07/11 18:30:21 tb Exp $ PROG= factor SRCS= factor.c pattern.c pr_tbl.c CFLAGS+=-I${.CURDIR}/../primes MAN= factor.6 -DPADD= ${LIBM} -LDADD= -lm .PATH: ${.CURDIR}/../primes .include diff --git a/games/factor/factor.c b/games/factor/factor.c index fc86828ca25..4b00a043207 100644 --- a/games/factor/factor.c +++ b/games/factor/factor.c @@ -1,4 +1,4 @@ -/* $OpenBSD: factor.c,v 1.27 2016/03/07 12:07:56 mestre Exp $ */ +/* $OpenBSD: factor.c,v 1.28 2016/07/11 18:30:21 tb Exp $ */ /* $NetBSD: factor.c,v 1.5 1995/03/23 08:28:07 cgd Exp $ */ /* @@ -55,7 +55,6 @@ #include #include #include -#include #include #include #include @@ -66,7 +65,7 @@ /* * prime[i] is the (i+1)th prime. * - * We are able to sieve 2^32-1 because this byte table yields all primes + * 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[]; @@ -74,9 +73,10 @@ extern const ubig *pr_limit; /* largest prime in the prime array */ extern const char pattern[]; extern const int pattern_size; -void pr_fact(u_int64_t); /* print factors of a value */ -void pr_bigfact(u_int64_t); -__dead void usage(void); +static void pr_fact(u_int64_t); /* print factors of a value */ +static void pr_bigfact(u_int64_t); +static u_int32_t usqrt(u_int64_t); +static void __dead usage(void); int main(int argc, char *argv[]) @@ -152,7 +152,7 @@ main(int argc, char *argv[]) * * Prime factors are printed with leading spaces. */ -void +static void pr_fact(u_int64_t val) /* Factor this value. */ { const ubig *fact; /* The factor found. */ @@ -196,15 +196,16 @@ pr_fact(u_int64_t val) /* Factor this value. */ (void)putchar('\n'); } - -/* At this point, our number may have factors greater than those in primes[]; +/* + * 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 +static void pr_bigfact(u_int64_t val) /* Factor this value. */ { - ubig start, stop, factor; + u_int64_t start, stop; + ubig factor; char *q; const ubig *p; ubig fact_lim, mod; @@ -212,7 +213,7 @@ pr_bigfact(u_int64_t val) /* Factor this value. */ char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ start = *pr_limit + 2; - stop = (ubig)sqrt((double)val); + stop = usqrt(val) + 1; if ((stop & 0x1) == 0) stop++; /* @@ -230,7 +231,8 @@ pr_bigfact(u_int64_t val) /* Factor this value. */ 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) { + fact_lim + pattern_size <= TABSIZE; + fact_lim += pattern_size) { memcpy(&table[fact_lim], pattern, pattern_size); } /* final block pattern copy */ @@ -238,11 +240,10 @@ pr_bigfact(u_int64_t val) /* Factor this value. */ if (stop-start > TABSIZE+TABSIZE) { tab_lim = &table[TABSIZE]; /* sieve it all */ - fact_lim = (int)sqrt( - (double)(start)+TABSIZE+TABSIZE+1.0); + fact_lim = usqrt(start + TABSIZE + TABSIZE + 1); } else { tab_lim = &table[(stop - start)/2]; /* partial sieve */ - fact_lim = (int)sqrt((double)(stop) + 1.0); + fact_lim = usqrt(stop + 1); } /* sieve for factors >= 17 */ factor = 17; /* 17 is first prime to use */ @@ -267,11 +268,11 @@ pr_bigfact(u_int64_t val) /* Factor this value. */ if (*q) { if (val % start == 0) { do { - (void)printf(" %lu", (unsigned long) start); + printf(" %llu", start); val /= start; } while ((val % start) == 0); (void)fflush(stdout); - stop = (ubig)sqrt((double)val); + stop = usqrt(val) + 1; if ((stop & 0x1) == 0) stop++; } @@ -282,8 +283,26 @@ pr_bigfact(u_int64_t val) /* Factor this value. */ printf(" %llu", val); } +/* Code taken from ping.c */ +static u_int32_t +usqrt(u_int64_t n) +{ + u_int64_t y, x = 1; + + if (n == 0 || n == 1) + return n; + + do { /* newton was a stinker */ + y = x; + x = n / x; + x += y; + x /= 2; + } while (((y < x) && (x - y) > 1) || (y - x) > 1); + + return (u_int32_t)x; +} -void +static void __dead usage(void) { (void)fprintf(stderr, "usage: %s [number ...]\n", getprogname()); -- 2.20.1