ARM: Implement long division.

* lib/mes/div.c (__mesabi_log2i): New procedure.
(__mesabi_uldiv1): New procedure.
(__mesabi_uldiv): Use it.

Co-Authored-By: Nathalie Kopaczewski <natkopa@gmail.com>

lib/mes/div.c

@@ -2,6 +2,7 @@
  * GNU Mes --- Maxwell Equations of Software
  * Copyright © 2016,2017,2018,2019 Jan (janneke) Nieuwenhuizen <janneke@gnu.org>
  * Copyright © 2019 Danny Milosavljevic <dannym@scratchpost.org>
+ * Copyright © 2020 Nathalie Kopaczewski <natkopa@gmail.com>
  *
  * This file is part of GNU Mes.
  *
@@ -24,77 +25,211 @@
 #include <limits.h>
 #include <signal.h>
 
-typedef struct
+struct ldiv_t
 {
   long quot;
   long rem;
-} ldiv_t;
+};
 
-int __raise(int);
+int __raise (int);
 
 void
 __mesabi_div0 (void)
 {
-  if (__raise(SIGFPE) < 0) { /* could not raise SIGFPE */
+  if (__raise (SIGFPE) < 0) /* could not raise SIGFPE */
-    /* Fail in any way possible */
+    {
-    unsigned char* x = (unsigned char*) 0;
+      /* Fail in any way possible */
-    *x = 2;
+      unsigned char *x = (unsigned char *) 0;
-  }
+      *x = 2;
+    }
 }
 
+#define ULONG_HIGHBITMASK LONG_MIN
+#define ULONG_BITCOUNT (sizeof (unsigned long)*8)
+
+/** Compute the logarithm of base 2 of D.  The result is rounded down.
+    That is equal to the highest-index set bit in D.
+
+    The idea is to shift D to the right in order to find the index i of the first most-significant digit > 0.
+
+    The computation is done by bisection, for speed.
+
+    Recurse:
+      Two halves are determined of the remaining slice.
+      The first half checked is the higher-significant half.
+        If that higher-significant half is not zero, recurse on that one.
+        Otherwise, recurse on the lower-significant half.
+
+    Precondition: D > 0 */
+static unsigned int
+__mesabi_log2i (unsigned long D)
+{
+  unsigned int n = ULONG_BITCOUNT;
+  unsigned int i = 0U;
+  unsigned long D1;
+  while (n >= 2U)
+    {                           /* while still two halves possible */
+      n >>= 1U;
+      /* D1: higher-significant half of D */
+      D1 = D >> n;
+      if (D1 > 0UL)
+        {
+          /* We know that the resulting index has to be in the higher-significant half.
+             In that case, lower-significant half of D is superfluous for determination of i,
+             therefore scroll and continue with higher-significant half. */
+          D = D1;
+          i += n;
+        }
+    }
+  return i;
+}
+
+#if 0
+static void
+test_log2i (void)
+{
+  assert (log2i (1) == 0);
+  assert (log2i (1) == 0);
+  assert (log2i (2) == 1);
+  assert (log2i (3) == 1);
+  assert (log2i (4) == 2);
+  assert (log2i (5) == 2);
+  assert (log2i (6) == 2);
+  assert (log2i (7) == 2);
+  assert (log2i (8) == 3);
+  assert (log2i (9) == 3);
+  assert (log2i (10) == 3);
+  assert (log2i (11) == 3);
+  assert (log2i (12) == 3);
+  assert (log2i (13) == 3);
+  assert (log2i (71) == 6);
+  assert (log2i (72) == 6);
+  assert (log2i (73) == 6);
+  assert (log2i (74) == 6);
+  assert (log2i (75) == 6);
+  assert (log2i (99) == 6);
+  assert (log2i (2147483648) == 31);
+  assert (log2i (3221225471) == 31);
+  assert (log2i (4294967294) == 31);
+  assert (log2i (4294967295) == 31);
+}
+#endif
+
+/** Perform unsigned integer division of N by D; store the remainder
+    into *REMAINDER; return the quotient.
+
+    This is currently implemented as long division.
+
+    R is the remainder.  R >= 0.  R starts at N.
+
+    QUOTIENT is built up bit by bit starting at the most significant bit [*].
+
+    Values D', starting at D << ULONG_BITCOUNT [*], going down to 1,
+    divided by 2 each time, are iterated over, doing: If R >= D',
+    subtract D' from R, and append new LSB 1 to the QUOTIENT.
+    Otherwise, subtract 0 from R (implicit), and append new LSB 0 to the
+    QUOTIENT (0 is the implicit default).
+
+    [*] As a special consideration for C throwing away bits when
+        left-shifting, D' starts at the highest value that will not lose
+        bits in this way instead.  (ULONG_BITCOUNT - log2i(D) - 1) is
+        the number of leading zeroes in D in binary radix.
+
+    Precondition: D > 0 */
+static unsigned long
+__mesabi_uldiv1 (unsigned long N, unsigned long D, unsigned long *remainder)
+{
+  // Note: __mesabi_log2i(D) < ULONG_BITCOUNT
+  unsigned int j = ULONG_BITCOUNT - __mesabi_log2i (D); /* Note: Or j = __mesabi_log2i(N) + 1 - __mesabi_log2i(D) */
+  // Note: assert(j - 1 == __builtin_clzl(D));  on GCC
+  unsigned long quotient = 0UL;
+  unsigned long R = N;
+  for (D <<= (j - 1); j > 0U; --j, D >>= 1U)
+    {
+      quotient <<= 1U;
+      if (R >= D)
+        {
+          R -= D;
+          quotient |= 1UL;
+        }
+    }
+  *remainder = R;
+  return quotient;
+}
+
+#if 0
+static void
+assert_uldiv (unsigned long N, unsigned long D, unsigned long expected_quotient,
+              unsigned long expected_remainder)
+{
+  unsigned long remainder;
+  unsigned long quotient;
+  quotient = uldiv (N, D, &remainder);
+  printf ("%lu/%lu = %lu;%lu\n", N, D, quotient, remainder);
+  assert (quotient == expected_quotient);
+  assert (remainder == expected_remainder);
+}
+
+static void
+test_uldiv (void)
+{
+  //assert_uldiv(0, 0, 0, 0);
+  assert_uldiv (0, 1, 0, 0);
+  assert_uldiv (1, 1, 1, 0);
+  assert_uldiv (72, 5, 14, 2);
+  assert_uldiv (0xffffffff, 1, 0xffffffff, 0);
+  assert_uldiv (0xffffffff, 2, 0x7fffffff, 1);
+}
+#endif
+
 /* Compare gcc: __udivmoddi4 */
 unsigned long
-__mesabi_uldiv (unsigned long a, unsigned long b, unsigned long* remainder)
+__mesabi_uldiv (unsigned long a, unsigned long b, unsigned long *remainder)
 {
   unsigned long tmp;
   if (!remainder)
     remainder = &tmp;
   *remainder = 0;
-  switch (b) {
+  switch (b)
-  case 64UL:
-    *remainder = a & 63UL;
-    return a >> 6UL;
-  case 32UL:
-    *remainder = a & 31UL;
-    return a >> 5UL;
-  case 16UL:
-    *remainder = a & 15UL;
-    return a >> 4UL;
-  case 8UL:
-    *remainder = a & 7UL;
-    return a >> 3UL;
-  case 4UL:
-    *remainder = a & 3UL;
-    return a >> 2UL;
-  case 2UL:
-    *remainder = a & 1UL;
-    return a >> 1UL;
-  case 1UL:
-    *remainder = 0;
-    return a;
-  case 0UL:
-    __mesabi_div0();
-    return 0UL;
-  default:
     {
-      unsigned long x;
+    case 64UL:
-      for (x = 0; a >= b; a -= b)
+      *remainder = a & 63UL;
-        ++x;
+      return a >> 6UL;
-      *remainder = a;
+    case 32UL:
-      return x;
+      *remainder = a & 31UL;
+      return a >> 5UL;
+    case 16UL:
+      *remainder = a & 15UL;
+      return a >> 4UL;
+    case 8UL:
+      *remainder = a & 7UL;
+      return a >> 3UL;
+    case 4UL:
+      *remainder = a & 3UL;
+      return a >> 2UL;
+    case 2UL:
+      *remainder = a & 1UL;
+      return a >> 1UL;
+    case 1UL:
+      *remainder = 0;
+      return a;
+    case 0UL:
+      __mesabi_div0 ();
+      return 0UL;
+    default:
+      return __mesabi_uldiv1 (a, b, remainder);
     }
-  }
 }
 
 /* Note: Rounds towards zero.
    Maintainer: Be careful to satisfy quot * b + rem == a.
                That means that rem can be negative. */
 void
-__mesabi_ldiv(long a, long b, ldiv_t* result)
+__mesabi_ldiv (long a, long b, struct ldiv_t *result)
 {
   int negate_result = (a < 0) ^ (b < 0);
   if (b == LONG_MIN)
-    __mesabi_div0();
+    __mesabi_div0 ();
   if (a != LONG_MIN)
     {
       int negative_a = (a < 0);
@@ -102,7 +237,7 @@ __mesabi_ldiv(long a, long b, ldiv_t* result)
         a = -a;
       if (b < 0)
         b = -b;
-      result->quot = __mesabi_uldiv(a, b, &result->rem);
+      result->quot = __mesabi_uldiv (a, b, &result->rem);
       if (negate_result)
         result->quot = -result->quot;
       if (negative_a)
@@ -120,13 +255,13 @@ __mesabi_ldiv(long a, long b, ldiv_t* result)
           negate_result = !negate_result;
         }
       else if (b == 0)
-        __mesabi_div0();
+        __mesabi_div0 ();
       else
         {
           long x;
           for (x = 0; a <= -b; a += b)
             ++x;
-          result->rem = a; /* negative */
+          result->rem = a;      /* negative */
           result->quot = x;
         }
       if (negate_result)
@@ -137,15 +272,15 @@ __mesabi_ldiv(long a, long b, ldiv_t* result)
 long
 __mesabi_imod (long a, long b)
 {
-  ldiv_t result;
+  struct ldiv_t result;
-  __mesabi_ldiv(a, b, &result);
+  __mesabi_ldiv (a, b, &result);
   return result.rem;
 }
 
 int
 __mesabi_idiv (int a, int b)
 {
-  ldiv_t result;
+  struct ldiv_t result;
-  __mesabi_ldiv(a, b, &result);
+  __mesabi_ldiv (a, b, &result);
   return result.quot;
 }