Add __divdi3(), remove __udivdi3() kernel dependency

Up until now no SPL consumer attempted to perform signed 64-bit
division so there was no need to support this.  That has now
changed so I adding 64-bit division support for 32-bit platforms.
The signed implementation is based on the unsigned version.

Since the have been several bug reports in the past concerning
correct 64-bit division on 32-bit platforms I added some long
over due regression tests.  Much to my surprise the unsigned
64-bit division regression tests failed.

This was surprising because __udivdi3() was implemented by simply
calling div64_u64() which is provided by the kernel.  This meant
that the linux kernels 64-bit division algorithm on 32-bit platforms
was flawed.  After some investigation this turned out to be exactly
the case.

Because of this I was forced to abandon the kernel helper and
instead to fully implement 64-bit division in the spl.  There are
several published implementation out there on how to do this
properly and I settled on one proposed in the book Hacker's Delight.
Their proposed algoritm is freely available without restriction
and I have just modified it to be linux kernel friendly.

The update implementation now passed all the unsigned and signed
regression tests.  This should be functional, but not fast, which is
good enough for out purposes.  If you want fast too I'd strongly
suggest you upgrade to a 64-bit platform.  I have also reported the
kernel bug and we'll see if we can't get it fixed up stream.

1 files changed, 96 insertions, 23 deletions

diff --git a/module/spl/spl-generic.c b/module/spl/spl-generic.c
index b9ab7109a..6a3f49ae0 100644
--- a/module/spl/spl-generic.c
+++ b/module/spl/spl-generic.c
@@ -94,39 +94,111 @@ highbit(unsigned long i)
 }
 EXPORT_SYMBOL(highbit);
 
+#if BITS_PER_LONG == 32
 /*
- * Implementation of 64 bit division for 32-bit machines.
+ * Support 64/64 => 64 division on a 32-bit platform.  While the kernel
+ * provides a div64_u64() function for this we do not use it because the
+ * implementation is flawed.  There are cases which return incorrect
+ * results as late as linux-2.6.35.  Until this is fixed upstream the
+ * spl must provide its own implementation.
+ *
+ * This implementation is a slightly modified version of the algorithm
+ * proposed by the book 'Hacker's Delight'.  The original source can be
+ * found here and is available for use without restriction.
+ *
+ * http://www.hackersdelight.org/HDcode/newCode/divDouble.c
  */
-#if BITS_PER_LONG == 32
-uint64_t
-__udivdi3(uint64_t dividend, uint64_t divisor)
-{
-#if defined(HAVE_DIV64_64) /* 2.6.22 - 2.6.25 API */
-	return div64_64(dividend, divisor);
-#elif defined(HAVE_DIV64_U64) /* 2.6.26 - 2.6.x API */
-	return div64_u64(dividend, divisor);
-#else
-	/* Implementation from 2.6.30 kernel */
-	uint32_t high, d;
 
-	high = divisor >> 32;
-	if (high) {
-		unsigned int shift = fls(high);
+/*
+ * Calculate number of leading of zeros for a 64-bit value.
+ */
+static int
+nlz64(uint64_t x) {
+	register int n = 0;
 
-		d = divisor >> shift;
-		dividend >>= shift;
-	} else
-		d = divisor;
+	if (x == 0)
+		return 64;
 
-	do_div(dividend, d);
+	if (x <= 0x00000000FFFFFFFFULL) {n = n + 32; x = x << 32;}
+	if (x <= 0x0000FFFFFFFFFFFFULL) {n = n + 16; x = x << 16;}
+	if (x <= 0x00FFFFFFFFFFFFFFULL) {n = n +  8; x = x <<  8;}
+	if (x <= 0x0FFFFFFFFFFFFFFFULL) {n = n +  4; x = x <<  4;}
+	if (x <= 0x3FFFFFFFFFFFFFFFULL) {n = n +  2; x = x <<  2;}
+	if (x <= 0x7FFFFFFFFFFFFFFFULL) {n = n +  1;}
+
+	return n;
+}
 
-	return dividend;
-#endif /* HAVE_DIV64_64, HAVE_DIV64_U64 */
+/*
+ * Newer kernels have a div_u64() function but we define our own
+ * to simplify portibility between kernel versions.
+ */
+static inline uint64_t
+__div_u64(uint64_t u, uint32_t v)
+{
+	(void) do_div(u, v);
+	return u;
+}
+
+/*
+ * Implementation of 64-bit unsigned division for 32-bit machines.
+ *
+ * First the procedure takes care of the case in which the divisor is a
+ * 32-bit quantity. There are two subcases: (1) If the left half of the
+ * dividend is less than the divisor, one execution of do_div() is all that
+ * is required (overflow is not possible). (2) Otherwise it does two
+ * divisions, using the grade school method.
+ */
+uint64_t
+__udivdi3(uint64_t u, uint64_t v)
+{
+	uint64_t u0, u1, v1, q0, q1, k;
+	int n;
+
+	if (v >> 32 == 0) {			// If v < 2**32:
+		if (u >> 32 < v) {		// If u/v cannot overflow,
+			return __div_u64(u, v);	// just do one division.
+		} else {			// If u/v would overflow:
+			u1 = u >> 32;		// Break u into two halves.
+			u0 = u & 0xFFFFFFFF;
+			q1 = __div_u64(u1, v);	// First quotient digit.
+			k  = u1 - q1 * v;	// First remainder, < v.
+			u0 += (k << 32);
+			q0 = __div_u64(u0, v);	// Seconds quotient digit.
+			return (q1 << 32) + q0;
+		}
+	} else {				// If v >= 2**32:
+		n = nlz64(v);			// 0 <= n <= 31.
+		v1 = (v << n) >> 32;		// Normalize divisor, MSB is 1.
+		u1 = u >> 1;			// To ensure no overflow.
+		q1 = __div_u64(u1, v1);		// Get quotient from
+		q0 = (q1 << n) >> 31;		// Undo normalization and
+						// division of u by 2.
+		if (q0 != 0)			// Make q0 correct or
+			q0 = q0 - 1;		// too small by 1.
+		if ((u - q0 * v) >= v)
+			q0 = q0 + 1;		// Now q0 is correct.
+	
+		return q0;
+	}
 }
 EXPORT_SYMBOL(__udivdi3);
 
 /*
- * Implementation of 64 bit modulo for 32-bit machines.
+ * Implementation of 64-bit signed division for 32-bit machines.
+ */
+int64_t
+__divdi3(int64_t u, int64_t v)
+{
+	int64_t q, t;
+	q = __udivdi3(abs64(u), abs64(v));
+	t = (u ^ v) >> 63;	// If u, v have different
+	return (q ^ t) - t;	// signs, negate q.
+}
+EXPORT_SYMBOL(__divdi3);
+
+/*
+ * Implementation of 64-bit unsigned modulo for 32-bit machines.
  */
 uint64_t
 __umoddi3(uint64_t dividend, uint64_t divisor)
@@ -134,6 +206,7 @@ __umoddi3(uint64_t dividend, uint64_t divisor)
 	return (dividend - (divisor * __udivdi3(dividend, divisor)));
 }
 EXPORT_SYMBOL(__umoddi3);
+
 #endif /* BITS_PER_LONG */
 
 /* NOTE: The strtoxx behavior is solely based on my reading of the Solaris