Merge PR #10553: Use a more traditional definition for uint63_div21 (fixes #10551).

Reviewed-by: maximedenes
Reviewed-by: proux01

5 files changed, 66 insertions, 154 deletions

diff --git a/kernel/byterun/coq_uint63_native.h b/kernel/byterun/coq_uint63_native.h
index 1fdafc9d8f..9fbd3f83d8 100644
--- a/kernel/byterun/coq_uint63_native.h
+++ b/kernel/byterun/coq_uint63_native.h
@@ -111,51 +111,26 @@ value uint63_mulc(value x, value y, value* h) {
 #define le128(xh,xl,yh,yl) (uint63_lt(xh,yh) || (uint63_eq(xh,yh) && uint63_leq(xl,yl)))
 
 #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 */
+/* precondition: xh < y */
+/* outputs r and sets ql to q 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);
+  uint64_t nh = uint63_of_value(xh);
+  uint64_t nl = uint63_of_value(xl);
   y = uint63_of_value(y);
-  uint64_t maskh = 0;
-  uint64_t maskl = 1;
-  uint64_t dh = 0;
-  uint64_t dl = y;
-  int cmp = 1;
-  /* 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);
+  uint64_t q = 0;
+  for (int i = 0; i < 63; ++i) {
+    // 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 <<= 1;
+    nh = (nh << 1) | (nl >> 63);
+    q <<= 1;
+    if (nh >= y) { q |= 1; nh -= y; }
   }
-  uint64_t remh = xh;
-  uint64_t reml = xl;
-  /* 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)) { /* 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 - 1)) & maxuint63);
-    maskh = maskh >> 1;
-    dl = (dl >> 1) | ((dh << (63 - 1)) & maxuint63);
-    dh = dh >> 1;
-    /* decr n */
-  }
-  *ql = Val_int(quotientl);
-  return Val_int(reml);
+  *ql = Val_int(q);
+  return Val_int(nh);
 }
 value uint63_div21(value xh, value xl, value y, value* ql) {
-  if (uint63_of_value(y) == 0) {
+  if (uint63_leq(y, xh)) {
     *ql = Val_int(0);
     return Val_int(0);
   } else {
diff --git a/kernel/uint63_amd64_63.ml b/kernel/uint63_amd64_63.ml
index d6b077a9f5..5c4028e1c8 100644
--- a/kernel/uint63_amd64_63.ml
+++ b/kernel/uint63_amd64_63.ml
@@ -96,55 +96,32 @@ let le (x : int) (y : int) =
   (x lxor 0x4000000000000000) <= (y lxor 0x4000000000000000)
 [@@ocaml.inline always]
 
-(* 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 *)
+(* precondition: xh < y *)
+(* outputs: q, 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
-  (* 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;
-    (* 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
-  (* 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 *)
-      (* quotienth := !quotienth lor !maskh *)
-      quotientl := !quotientl lor !maskl;
-      remh := if lt !reml !dl then !remh - !dh - 1 else !remh - !dh;
-      reml := !reml - !dl;
-    end;
-    maskl := (!maskl lsr 1) lor (!maskh lsl (uint_size - 1));
-    maskh := !maskh lsr 1;
-    dl := (!dl lsr 1) lor (!dh lsl (uint_size - 1));
-    dh := !dh lsr 1;
-    (* decr n *)
+  (* nh might temporarily grow as large as 2*y - 1 in the loop body,
+     so we store it as a 64-bit unsigned integer *)
+  let nh = ref xh in
+  let nl = ref xl in
+  let q = ref 0 in
+  for _i = 0 to 62 do
+    (* invariants: 0 <= nh < y, nl = (xl*2^i) % 2^63,
+       (q*y + nh) * 2^(63-i) + (xl % 2^(63-i)) = (xh%y) * 2^63 + xl *)
+    nh := Int64.logor (Int64.shift_left !nh 1) (Int64.of_int (!nl lsr 62));
+    nl := !nl lsl 1;
+    q := !q lsl 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 := !q lor 1; nh := Int64.sub !nh y; end
   done;
-  !quotientl, !reml
+  !q, Int64.to_int !nh
 
-let div21 xh xl y = if y = 0 then 0, 0 else div21 xh xl y
+let div21 xh xl y =
+  let xh = to_uint64 xh in
+  let y = to_uint64 y in
+  if Int64.compare y xh <= 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_i386_31.ml b/kernel/uint63_i386_31.ml
index 2a3fc75ec1..b8eccd19fb 100644
--- a/kernel/uint63_i386_31.ml
+++ b/kernel/uint63_i386_31.ml
@@ -88,55 +88,28 @@ let diveucl x y = (div x y, 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 *)
+(* precondition: xh < y *)
+(* outputs: q, 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 xh 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
+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 =
diff --git a/test-suite/arithmetic/diveucl_21.v b/test-suite/arithmetic/diveucl_21.v
index b888c97be3..b12dba429c 100644
--- a/test-suite/arithmetic/diveucl_21.v
+++ b/test-suite/arithmetic/diveucl_21.v
@@ -10,11 +10,11 @@ Check (eq_refl (4611686018427387904,1) <<: diveucl_21 1 1 2 = (46116860184273879
 Definition compute1 := Eval compute in diveucl_21 1 1 2.
 Check (eq_refl compute1 : (4611686018427387904,1) = (4611686018427387904,1)).
 
-Check (eq_refl : diveucl_21 3 1 2 = (4611686018427387904, 1)).
-Check (eq_refl (4611686018427387904, 1) <: diveucl_21 3 1 2 = (4611686018427387904, 1)).
-Check (eq_refl (4611686018427387904, 1) <<: diveucl_21 3 1 2 = (4611686018427387904, 1)).
+Check (eq_refl : diveucl_21 3 1 2 = (0, 0)).
+Check (eq_refl (0, 0) <: diveucl_21 3 1 2 = (0, 0)).
+Check (eq_refl (0, 0) <<: diveucl_21 3 1 2 = (0, 0)).
 Definition compute2 := Eval compute in diveucl_21 3 1 2.
-Check (eq_refl compute2 : (4611686018427387904, 1) = (4611686018427387904, 1)).
+Check (eq_refl compute2 : (0, 0) = (0, 0)).
 
 Check (eq_refl : diveucl_21 1 1 0 = (0,0)).
 Check (eq_refl (0,0) <: diveucl_21 1 1 0 = (0,0)).
@@ -23,3 +23,7 @@ Check (eq_refl (0,0) <<: diveucl_21 1 1 0 = (0,0)).
 Check (eq_refl : diveucl_21 9223372036854775807 0 1 = (0,0)).
 Check (eq_refl (0,0) <: diveucl_21 9223372036854775807 0 1 = (0,0)).
 Check (eq_refl (0,0) <<: diveucl_21 9223372036854775807 0 1 = (0,0)).
+
+Check (eq_refl : diveucl_21 9305446873517 1793572051078448654 4930380657631323783 = (17407905077428, 3068214991893055266)).
+Check (eq_refl (17407905077428, 3068214991893055266) <: diveucl_21 9305446873517 1793572051078448654 4930380657631323783 = (17407905077428, 3068214991893055266)).
+Check (eq_refl (17407905077428, 3068214991893055266) <<: diveucl_21 9305446873517 1793572051078448654 4930380657631323783 = (17407905077428, 3068214991893055266)).
diff --git a/theories/Numbers/Cyclic/Int63/Int63.v b/theories/Numbers/Cyclic/Int63/Int63.v
index c81ba02230..9e9481341f 100644
--- a/theories/Numbers/Cyclic/Int63/Int63.v
+++ b/theories/Numbers/Cyclic/Int63/Int63.v
@@ -388,7 +388,7 @@ Axiom diveucl_def_spec : forall x y, diveucl x y = diveucl_def x y.
 Axiom diveucl_21_spec :  forall a1 a2 b,
    let (q,r) := diveucl_21 a1 a2 b in
    let (q',r') := Z.div_eucl ([|a1|] * wB + [|a2|]) [|b|] in
-   [|q|] = Z.modulo q' wB /\ [|r|] = r'.
+   [|a1|] < [|b|] -> [|q|] = q' /\ [|r|] = r'.
 
 Axiom addmuldiv_def_spec : forall p x y,
   addmuldiv p x y = addmuldiv_def p x y.
@@ -1421,26 +1421,9 @@ Proof.
  generalize (Z_div_mod ([|a1|]*wB+[|a2|]) [|b|] H).
  revert W.
  destruct (diveucl_21 a1 a2 b); destruct (Z.div_eucl ([|a1|]*wB+[|a2|]) [|b|]).
- intros (H', H''); rewrite H', H''; clear H' H''.
+ intros (H', H''); auto; rewrite H', H''; clear H' H''.
  intros (H', H''); split; [ |exact H''].
- rewrite H', Zmult_comm, Z.mod_small; [reflexivity| ].
- split.
- { revert H'; case z; [now simpl..|intros p H'].
-   exfalso; apply (Z.lt_irrefl 0), (Z.le_lt_trans _ ([|a1|] * wB + [|a2|])).
-   { now apply Z.add_nonneg_nonneg; [apply Z.mul_nonneg_nonneg| ]. }
-   rewrite H'; apply (Zplus_lt_reg_r _ _ (- z0)); ring_simplify.
-   apply (Z.le_lt_trans _ (- [|b|])); [ |now auto with zarith].
-   rewrite Z.opp_eq_mul_m1; apply Zmult_le_compat_l; [ |now apply Wb].
-   rewrite <-!Pos2Z.opp_pos, <-Z.opp_le_mono.
-   now change 1 with (Z.succ 0); apply Zlt_le_succ. }
- rewrite <-Z.nle_gt; intro Hz; revert H2; apply Zle_not_lt.
- rewrite (Z.div_unique_pos (wB * [|a1|] + [|a2|]) wB [|a1|] [|a2|]);
-   [ |now simpl..].
- rewrite Z.mul_comm, H'.
- rewrite (Z.div_unique_pos (wB * [|b|] + z0) wB [|b|] z0) at 1;
-   [ |split; [ |apply (Z.lt_trans _ [|b|])]; now simpl|reflexivity].
- apply Z_div_le; [now simpl| ]; rewrite Z.mul_comm; apply Zplus_le_compat_r.
- now apply Zmult_le_compat_l.
+ now rewrite H', Zmult_comm.
 Qed.
 
 Lemma div2_phi ih il j: (2^62 <= [|j|] -> [|ih|] < [|j|] ->