From dd60b4a292b870e08c23ddcb363630cbb2ed1227 Mon Sep 17 00:00:00 2001
From: Pierre Roux
Date: Thu, 2 May 2019 08:16:37 +0200
Subject: [primitive integers] Make div21 implems consistent with its
 specification

There are three implementations of this primitive:
* one in OCaml on 63 bits integer in kernel/uint63_amd64.ml
* one in OCaml on Int64 in kernel/uint63_x86.ml
* one in C on unsigned 64 bit integers in kernel/byterun/coq_uint63_native.h
Its specification is the axiom `diveucl_21_spec` in
theories/Numbers/Cyclic/Int63/Int63.v

* comment the implementations with loop invariants to enable an easy
  pen&paper proof of correctness (note to reviewers: the one in
  uint63_amd64.ml might be the easiest to read)
* make sure the three implementations are equivalent
* fix the specification in Int63.v
  (only the lowest part of the result is actually returned)
* make a little optimisation in div21 enabled by the proof of correctness
  (cmp is computed at the end of the first loop rather than at the beginning,
   potentially saving one loop iteration while remaining correct)
* update the proofs in Int63.v and Cyclic63.v to take into account the
  new specifiation of div21
* add a test
---
 kernel/byterun/coq_interp.c        | 39 ++++------------------------
 kernel/byterun/coq_uint63_native.h | 53 ++++++++++++++++++++++++++------------
 kernel/uint63.mli                  |  4 +++
 kernel/uint63_amd64.ml             | 26 ++++++++++++++-----
 kernel/uint63_x86.ml               | 25 +++++++++++++-----
 5 files changed, 82 insertions(+), 65 deletions(-)

(limited to 'kernel')

diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index 2293ae9dfd..e38d458b36 100644
--- a/kernel/byterun/coq_interp.c
+++ b/kernel/byterun/coq_interp.c
@@ -1374,40 +1374,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;
       }
diff --git a/kernel/byterun/coq_uint63_native.h b/kernel/byterun/coq_uint63_native.h
index d431dc1e5c..8dee0d69d3 100644
--- a/kernel/byterun/coq_uint63_native.h
+++ b/kernel/byterun/coq_uint63_native.h
@@ -109,37 +109,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))
diff --git a/kernel/uint63.mli b/kernel/uint63.mli
index b5f40ca804..f25f24512d 100644
--- a/kernel/uint63.mli
+++ b/kernel/uint63.mli
@@ -40,6 +40,10 @@ val rem     : t -> t -> t
       (* Specific arithmetic operations *)
 val mulc    : t -> t -> t * t
 val addmuldiv : t -> t -> t -> t
+
+(** [div21 xh xl y] returns [q % 2^63, r]
+    s.t. [xh * 2^63 + xl = q * y + r] and [r < y].
+    When [y] is [0], returns [0, 0]. *)
 val div21   : t -> t -> t -> t * t
 
       (* comparison *)
diff --git a/kernel/uint63_amd64.ml b/kernel/uint63_amd64.ml
index 010b594de8..2d4d685775 100644
--- a/kernel/uint63_amd64.ml
+++ b/kernel/uint63_amd64.ml
@@ -102,26 +102,35 @@ let le128 xh xl yh yl =
   lt xh yh || (xh = yh && le xl yl)
 
     (* division of two numbers by one *)
+(* precondition: y <> 0 *)
+(* outputs: q % 2^63, r s.t. x = q * y + r, r < y *)
 let div21 xh xl y =
   let maskh = ref 0 in
   let maskl = ref 1 in
   let dh = ref 0 in
   let dl = ref y in
   let cmp = ref true in
-  while !dh >= 0 && !cmp do
-    cmp := lt128 !dh !dl xh xl;
+  (* n = ref 0 *)
+  (* loop invariant: mask = 2^n, d = mask * y, (2 * d <= x -> cmp), n >= 0 *)
+  while !dh >= 0 && !cmp do (* dh >= 0 tests that dh highest bit is zero *)
     (* We don't use addmuldiv below to avoid checks on 1 *)
     dh := (!dh lsl 1) lor (!dl lsr (uint_size - 1));
     dl := !dl lsl 1;
     maskh := (!maskh lsl 1) lor (!maskl lsr (uint_size - 1));
-    maskl := !maskl lsl 1
-  done; (* mask = 2^N, d = 2^N * d, d >= x *)
+    maskl := !maskl lsl 1;
+    (* incr n *)
+    cmp := lt128 !dh !dl xh xl;
+  done; (* mask = 2^n, d = 2^n * y, 2 * d > x *)
   let remh = ref xh in
   let reml = ref xl in
-  let quotient = ref 0 in
+  (* quotienth = ref 0 *)
+  let quotientl = ref 0 in
+  (* loop invariant: x = quotient * y + rem, y * 2^(n+1) > r,
+     mask = floor(2^n), d = mask * y, n >= -1 *)
   while !maskh lor !maskl <> 0 do
     if le128 !dh !dl !remh !reml then begin (* if rem >= d, add one bit and subtract d *)
-      quotient := !quotient lor !maskl;
+      (* quotienth := !quotienth lor !maskh *)
+      quotientl := !quotientl lor !maskl;
       remh := if lt !reml !dl then !remh - !dh - 1 else !remh - !dh;
       reml := !reml - !dl;
     end;
@@ -129,8 +138,11 @@ let div21 xh xl y =
     maskh := !maskh lsr 1;
     dl := (!dl lsr 1) lor (!dh lsl (uint_size - 1));
     dh := !dh lsr 1;
+    (* decr n *)
   done;
-  !quotient, !reml
+  !quotientl, !reml
+
+let div21 xh xl y = if y = 0 then 0, 0 else div21 xh xl y
 
      (* exact multiplication *)
 (* TODO: check that none of these additions could be a logical or *)
diff --git a/kernel/uint63_x86.ml b/kernel/uint63_x86.ml
index 461184c432..fa45c90241 100644
--- a/kernel/uint63_x86.ml
+++ b/kernel/uint63_x86.ml
@@ -94,26 +94,35 @@ let le128 xh xl yh yl =
   lt xh yh || (xh = yh && le xl yl)
 
     (* division of two numbers by one *)
+(* precondition: y <> 0 *)
+(* outputs: q % 2^63, r s.t. x = q * y + r, r < y *)
 let div21 xh xl y =
   let maskh = ref zero in
   let maskl = ref one in
   let dh = ref zero in
   let dl = ref y in
   let cmp = ref true in
-  while le zero !dh && !cmp do
-    cmp := lt128 !dh !dl xh xl;
+  (* n = ref 0 *)
+  (* loop invariant: mask = 2^n, d = mask * y, (2 * d <= x -> cmp), n >= 0 *)
+  while Int64.equal (l_sr !dh (of_int (uint_size - 1))) zero && !cmp do
     (* We don't use addmuldiv below to avoid checks on 1 *)
     dh := l_or (l_sl !dh one) (l_sr !dl (of_int (uint_size - 1)));
     dl := l_sl !dl one;
     maskh := l_or (l_sl !maskh one) (l_sr !maskl (of_int (uint_size - 1)));
-    maskl := l_sl !maskl one
-  done; (* mask = 2^N, d = 2^N * d, d >= x *)
+    maskl := l_sl !maskl one;
+    (* incr n *)
+    cmp := lt128 !dh !dl xh xl;
+  done; (* mask = 2^n, d = 2^n * d, 2 * d > x *)
   let remh = ref xh in
   let reml = ref xl in
-  let quotient = ref zero in
+  (* quotienth = ref 0 *)
+  let quotientl = ref zero in
+  (* loop invariant: x = quotient * y + rem, y * 2^(n+1) > r,
+     mask = floor(2^n), d = mask * y, n >= -1 *)
   while not (Int64.equal (l_or !maskh !maskl) zero) do
     if le128 !dh !dl !remh !reml then begin (* if rem >= d, add one bit and subtract d *)
-      quotient := l_or !quotient !maskl;
+      (* quotienth := !quotienth lor !maskh *)
+      quotientl := l_or !quotientl !maskl;
       remh := if lt !reml !dl then sub (sub !remh !dh) one else sub !remh !dh;
       reml := sub !reml !dl
     end;
@@ -121,9 +130,11 @@ let div21 xh xl y =
     maskh := l_sr !maskh one;
     dl := l_or (l_sr !dl one) (l_sl !dh (of_int (uint_size - 1)));
     dh := l_sr !dh one
+    (* decr n *)
   done;
-  !quotient, !reml
+  !quotientl, !reml
 
+let div21 xh xl y = if Int64.equal y zero then zero, zero else div21 xh xl y
 
      (* exact multiplication *)
 let mulc x y =
-- 
cgit v1.2.3


From 09cdf7b1fad8761cdf7048bf38a468c8558eb0d5 Mon Sep 17 00:00:00 2001
From: Pierre Roux
Date: Thu, 2 May 2019 08:41:45 +0200
Subject: Remove now useless commented code

---
 kernel/byterun/coq_interp.c        | 15 +--------------
 kernel/byterun/coq_uint63_emul.h   |  2 ++
 kernel/byterun/coq_uint63_native.h |  1 +
 3 files changed, 4 insertions(+), 14 deletions(-)

(limited to 'kernel')

diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index e38d458b36..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();
@@ -1587,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 8dee0d69d3..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 */
-- 
cgit v1.2.3