Last fall there was a thread containing some examples of composite numbers that were not factored properly by factor(6):
https://marc.info/?t=144155842300002&r=1&w=2 Some suggestions were made and some (incomplete) diffs were posted, but no action was taken. Below is a diff that uses Newton's method from ping.c, as suggested by Otto, and which fixes both bugs that were apparent in the thread: 1. The bigger problem is rounding: stop = (ubig)sqrt((double)val); is zero if val is close to 2^64. The reason is that the square root will be rounded up to 2^32 so that the cast with ubig (u_int32_t) results in zero. 2. The other problem is an off-by-one in the loop condition in pr_bigfact(): for val = 4295360521 = 65539 * 65539, both start and stop are equal to 65539 and hence the loop while (start < stop) { ... } is never entered so that 4295360521 is falsely reported as a prime. Both problems are only relevant if val has no prime factors smaller than 65539. Finally, as mentioned by naddy, another option would be to think about the FreeBSD way and make use of libcrypto's BN functions. Index: Makefile =================================================================== RCS file: /var/cvs/src/games/factor/Makefile,v retrieving revision 1.5 diff -u -p -r1.5 Makefile --- Makefile 30 Mar 2016 06:38:40 -0000 1.5 +++ Makefile 11 Jul 2016 12:04:27 -0000 @@ -4,8 +4,6 @@ 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 <bsd.prog.mk> Index: factor.c =================================================================== RCS file: /var/cvs/src/games/factor/factor.c,v retrieving revision 1.27 diff -u -p -r1.27 factor.c --- factor.c 7 Mar 2016 12:07:56 -0000 1.27 +++ factor.c 11 Jul 2016 15:43:35 -0000 @@ -55,7 +55,6 @@ #include <ctype.h> #include <err.h> #include <errno.h> -#include <math.h> #include <stdio.h> #include <stdlib.h> #include <string.h> @@ -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 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 v (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 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 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 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 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 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());
