Merge PR #10046: [primitive integers] Make div21 implems consistent with its specification

Ack-by: Zimmi48
Ack-by: herbelin
Ack-by: maximedenes
Ack-by: proux01
Reviewed-by: vbgl

3 files changed, 45 insertions, 65 deletions

diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index 2293ae9dfd..1b348ae777 100644
--- a/kernel/byterun/coq_interp.c
+++ b/kernel/byterun/coq_interp.c
@@ -29,13 +29,6 @@
 #include "coq_uint63_emul.h"
 #endif
 
-/* spiwack: I append here a few macros for value/number manipulation */
-#define uint32_of_value(val) (((uint32_t)(val)) >> 1)
-#define value_of_uint32(i)   ((value)((((uint32_t)(i)) << 1) | 1))
-#define UI64_of_uint32(lo) ((uint64_t)((uint32_t)(lo)))
-#define UI64_of_value(val) (UI64_of_uint32(uint32_of_value(val)))
-/* /spiwack */
-
 
 
 /* Registers for the abstract machine:
@@ -1298,12 +1291,6 @@ value coq_interprete
         /*returns the multiplication on a pair */
         print_instr("MULCINT63");
         CheckInt2();
-        /*accu = 2v+1, *sp=2w+1 ==> p = 2v*w */
-        /* TODO: implement
-        p = I64_mul (UI64_of_value (accu), UI64_of_uint32 ((*sp++)^1));
-        AllocPair(); */
-        /* Field(accu, 0) = (value)(I64_lsr(p,31)|1) ; */ /*higher part*/
-        /* Field(accu, 1) = (value)(I64_to_int32(p)|1); */ /*lower part*/
         Uint63_mulc(accu, *sp, sp);
         *--sp = accu;
         AllocPair();
@@ -1374,40 +1361,11 @@ value coq_interprete
       Instruct (CHECKDIV21INT63) {
         print_instr("DIV21INT63");
         CheckInt3();
-        /* spiwack: takes three int31 (the two first ones represent an
-                    int62) and performs the euclidian division of the
-                    int62 by the int31 */
-        /* TODO: implement this
-        bigint = UI64_of_value(accu);
-        bigint = I64_or(I64_lsl(bigint, 31),UI64_of_value(*sp++));
-        uint64 divisor;
-        divisor = UI64_of_value(*sp++);
-        Alloc_small(accu, 2, 1); */ /* ( _ , arity, tag ) */
-        /* if (I64_is_zero (divisor)) {
-           Field(accu, 0) = 1; */ /* 2*0+1 */
-        /* Field(accu, 1) = 1; */ /* 2*0+1 */
-        /* }
-        else {
-          uint64 quo, mod;
-          I64_udivmod(bigint, divisor, &quo, &mod);
-          Field(accu, 0) = value_of_uint32(I64_to_int32(quo));
-          Field(accu, 1) = value_of_uint32(I64_to_int32(mod));
-        } */
-        int b;
-        Uint63_eq0(b, sp[1]);
-        if (b) {
-          AllocPair();
-          Field(accu, 0) = sp[1];
-          Field(accu, 1) = sp[1];
-	}
-        else {
-          Uint63_div21(accu, sp[0], sp[1], sp);
-          sp[1] = sp[0];
-          Swap_accu_sp;
-          AllocPair();
-          Field(accu, 0) = sp[1];
-          Field(accu, 1) = sp[0];
-	}
+        Uint63_div21(accu, sp[0], sp[1], &(sp[1]));
+        Swap_accu_sp;
+        AllocPair();
+        Field(accu, 0) = sp[1];
+        Field(accu, 1) = sp[0];
         sp += 2;
         Next;
       }
@@ -1616,7 +1574,7 @@ value coq_push_vstack(value stk, value max_stack_size) {
   print_instr("push_vstack");print_int(len);
    for(i = 0; i < len; i++) coq_sp[i] = Field(stk,i);
   sp = coq_sp;
-  CHECK_STACK(uint32_of_value(max_stack_size));
+  CHECK_STACK(uint_of_value(max_stack_size));
   return Val_unit;
 }
 
diff --git a/kernel/byterun/coq_uint63_emul.h b/kernel/byterun/coq_uint63_emul.h
index d982f67566..528cc6fc1f 100644
--- a/kernel/byterun/coq_uint63_emul.h
+++ b/kernel/byterun/coq_uint63_emul.h
@@ -6,6 +6,8 @@
 
 #define Is_uint63(v) (Tag_val(v) == Custom_tag)
 
+#define uint_of_value(val) (((uint32_t)(val)) >> 1)
+
 # define DECLARE_NULLOP(name) \
 value uint63_##name() { \
   static value* cb = 0; \
diff --git a/kernel/byterun/coq_uint63_native.h b/kernel/byterun/coq_uint63_native.h
index d431dc1e5c..1fdafc9d8f 100644
--- a/kernel/byterun/coq_uint63_native.h
+++ b/kernel/byterun/coq_uint63_native.h
@@ -1,5 +1,6 @@
 #define Is_uint63(v) (Is_long(v))
 
+#define uint_of_value(val) (((uint64_t)(val)) >> 1)
 #define uint63_of_value(val) ((uint64_t)(val) >> 1)
 
 /* 2^63 * y + x as a value */
@@ -109,37 +110,56 @@ value uint63_mulc(value x, value y, value* h) {
 #define lt128(xh,xl,yh,yl) (uint63_lt(xh,yh) || (uint63_eq(xh,yh) && uint63_lt(xl,yl)))
 #define le128(xh,xl,yh,yl) (uint63_lt(xh,yh) || (uint63_eq(xh,yh) && uint63_leq(xl,yl)))
 
-value uint63_div21(value xh, value xl, value y, value* q) {
-  xh = (uint64_t)xh >> 1;
-  xl = ((uint64_t)xl >> 1) | ((uint64_t)xh << 63);
-  xh = (uint64_t)xh >> 1;
+#define maxuint63 ((uint64_t)0x7FFFFFFFFFFFFFFF)
+/* precondition: y <> 0 */
+/* outputs r and sets ql to q % 2^63 s.t. x = q * y + r, r < y */
+static value uint63_div21_aux(value xh, value xl, value y, value* ql) {
+  xh = uint63_of_value(xh);
+  xl = uint63_of_value(xl);
+  y = uint63_of_value(y);
   uint64_t maskh = 0;
   uint64_t maskl = 1;
   uint64_t dh = 0;
-  uint64_t dl = (uint64_t)y >> 1;
+  uint64_t dl = y;
   int cmp = 1;
-  while (dh >= 0 && cmp) {
+  /* int n = 0 */
+  /* loop invariant: mask = 2^n, d = mask * y, (2 * d <= x -> cmp), n >= 0, d < 2^(2*63) */
+  while (!(dh >> (63 - 1)) && cmp) {
+    dh = (dh << 1) | (dl >> (63 - 1));
+    dl = (dl << 1) & maxuint63;
+    maskh = (maskh << 1) | (maskl >> (63 - 1));
+    maskl = (maskl << 1) & maxuint63;
+    /* ++n */
     cmp = lt128(dh,dl,xh,xl);
-    dh = (dh << 1) | (dl >> 63);
-    dl = dl << 1;
-    maskh = (maskh << 1) | (maskl >> 63);
-    maskl = maskl << 1;
   }
   uint64_t remh = xh;
   uint64_t reml = xl;
-  uint64_t quotient = 0;
+  /* uint64_t quotienth = 0; */
+  uint64_t quotientl = 0;
+  /* loop invariant: x = quotient * y + rem, y * 2^(n+1) > r,
+     mask = floor(2^n), d = mask * y, n >= -1 */
   while (maskh | maskl) {
-    if (le128(dh,dl,remh,reml)) {
-      quotient = quotient | maskl;
-      if (uint63_lt(reml,dl)) {remh = remh - dh - 1;} else {remh = remh - dh;}
+    if (le128(dh,dl,remh,reml)) { /* if rem >= d, add one bit and subtract d */
+      /* quotienth = quotienth | maskh */
+      quotientl = quotientl | maskl;
+      remh = (uint63_lt(reml,dl)) ? (remh - dh - 1) : (remh - dh);
       reml = reml - dl;
     }
-    maskl = (maskl >> 1) | (maskh << 63);
+    maskl = (maskl >> 1) | ((maskh << (63 - 1)) & maxuint63);
     maskh = maskh >> 1;
-    dl = (dl >> 1) | (dh << 63);
+    dl = (dl >> 1) | ((dh << (63 - 1)) & maxuint63);
     dh = dh >> 1;
+    /* decr n */
   }
-  *q = Val_int(quotient);
+  *ql = Val_int(quotientl);
   return Val_int(reml);
 }
+value uint63_div21(value xh, value xl, value y, value* ql) {
+  if (uint63_of_value(y) == 0) {
+    *ql = Val_int(0);
+    return Val_int(0);
+  } else {
+    return uint63_div21_aux(xh, xl, y, ql);
+  }
+}
 #define Uint63_div21(xh, xl, y, q) (accu = uint63_div21(xh, xl, y, q))