From dd2278ec1161079f2d28079da92398cb4743972b Mon Sep 17 00:00:00 2001
From: Guillaume Melquiond
Date: Tue, 23 Jul 2019 17:48:29 +0200
Subject: Transpose the C code of uint63_div21 to the OCaml implementations.

---
 kernel/uint63_i386_31.ml | 56 ++++++++++++------------------------------------
 1 file changed, 14 insertions(+), 42 deletions(-)

(limited to 'kernel/uint63_i386_31.ml')

diff --git a/kernel/uint63_i386_31.ml b/kernel/uint63_i386_31.ml
index c3279779e1..91a1e397a8 100644
--- a/kernel/uint63_i386_31.ml
+++ b/kernel/uint63_i386_31.ml
@@ -86,53 +86,25 @@ let rem x y =
 let addmuldiv p x y =
   l_or (l_sl x p) (l_sr y Int64.(sub (of_int uint_size) p))
 
-(* A few helper functions on 128 bits *)
-let lt128 xh xl yh yl =
-  lt xh yh || (xh = yh && lt xl yl)
-
-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
-  (* 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;
-    (* 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
-  (* 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 *)
-      (* 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;
-    maskl := l_or (l_sr !maskl one) (l_sl !maskh (of_int (uint_size - 1)));
-    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 *)
+  let nh = ref (Int64.rem xh y) in
+  let nl = ref xl in
+  let q = ref 0L in
+  for _i = 0 to 62 do
+    (* invariants: 0 <= nh < y, nl = (xl*2^i) % 2^64,
+       (q*y + nh) * 2^(63-i) + (xl % 2^(63-i)) = (xh%y) * 2^63 + xl *)
+    nl := Int64.shift_left !nl 1;
+    nh := Int64.logor (Int64.shift_left !nh 1) (Int64.shift_right_logical !nl 63);
+    q := Int64.shift_left !q 1;
+    (* TODO: use "Int64.unsigned_compare !nh y >= 0",
+       once OCaml 4.08 becomes the minimal required version *)
+    if Int64.compare !nh 0L < 0 || Int64.compare !nh y >= 0 then
+      begin q := Int64.logor !q 1L; nh := Int64.sub !nh y; end
   done;
-  !quotientl, !reml
+  !q, !nh
 
 let div21 xh xl y = if Int64.equal y zero then zero, zero else div21 xh xl y
 
-- 
cgit v1.2.3


From 654254ca2e6ac94d3e0c62a127f92caff4b5938c Mon Sep 17 00:00:00 2001
From: Guillaume Melquiond
Date: Mon, 29 Jul 2019 21:24:11 +0200
Subject: Use the precondition of diveucl_21 to avoid massaging the dividend.

All the implementations now return (0, 0) when the dividend is so large
that the quotient would overflow.
---
 kernel/uint63_i386_31.ml | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

(limited to 'kernel/uint63_i386_31.ml')

diff --git a/kernel/uint63_i386_31.ml b/kernel/uint63_i386_31.ml
index 91a1e397a8..2306427b03 100644
--- a/kernel/uint63_i386_31.ml
+++ b/kernel/uint63_i386_31.ml
@@ -87,10 +87,10 @@ let addmuldiv p x y =
   l_or (l_sl x p) (l_sr y Int64.(sub (of_int uint_size) p))
 
     (* division of two numbers by one *)
-(* precondition: y <> 0 *)
-(* outputs: q % 2^63, r s.t. x = q * y + r, r < y *)
+(* precondition: xh < y *)
+(* outputs: q, r s.t. x = q * y + r, r < y *)
 let div21 xh xl y =
-  let nh = ref (Int64.rem xh y) in
+  let nh = ref xh in
   let nl = ref xl in
   let q = ref 0L in
   for _i = 0 to 62 do
@@ -106,7 +106,8 @@ let div21 xh xl y =
   done;
   !q, !nh
 
-let div21 xh xl y = if Int64.equal y zero then zero, zero else div21 xh xl y
+let div21 xh xl y =
+  if Int64.compare y xh <= 0 then zero, zero else div21 xh xl y
 
      (* exact multiplication *)
 let mulc x y =
-- 
cgit v1.2.3