On 09/08/15 08:03, Otto Moerbeek wrote:
> My bet would be Newton iteration, which should converge faster than
> binary search. ping(8) has an implementation of integer sqrt.
>
> BTW, long double won't fix the problem on al platforms since it is not
> guaranteed to be longer than double.
Newton is in general faster. Still I'd prefer the binary search as it
is simpler and avoids issues with accuracy (LSB). Anyway, the
factorization itself is the much slower part for products of big primes.
Regards,
Michael
--- /usr/src/games/factor/factor.c Wed Oct 28 00:59:24 2009
+++ factor.c Tue Sep 8 20:06:44 2015
@@ -192,6 +192,19 @@ pr_fact(u_int64_t val) /* Factor this value. */
(void)putchar('\n');
}
+static u_int32_t
+int_sqrt(u_int64_t y)
+{
+ int i;
+ u_int32_t m_i;
+ u_int32_t m = 0;
+
+ for (i = 32; i >= 0; i--) {
+ m_i = m | (1U << i);
+ if ((u_int32_t)m_i * m_i <= y) m = m_i;
+ }
+ return m;
+}
/* 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
@@ -208,7 +221,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 = int_sqrt(val);
if ((stop & 0x1) == 0)
stop++;
/*
