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/uint63_amd64.ml | 26 +++++++++++++++++++-------
 1 file changed, 19 insertions(+), 7 deletions(-)

(limited to 'kernel/uint63_amd64.ml')

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 *)
-- 
cgit v1.2.3