[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

9 files changed, 156 insertions, 97 deletions

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 =
diff --git a/test-suite/arithmetic/diveucl_21.v b/test-suite/arithmetic/diveucl_21.v
index 7e12a08906..b888c97be3 100644
--- a/test-suite/arithmetic/diveucl_21.v
+++ b/test-suite/arithmetic/diveucl_21.v
@@ -15,3 +15,11 @@ Check (eq_refl (4611686018427387904, 1) <: diveucl_21 3 1 2 = (46116860184273879
 Check (eq_refl (4611686018427387904, 1) <<: diveucl_21 3 1 2 = (4611686018427387904, 1)).
 Definition compute2 := Eval compute in diveucl_21 3 1 2.
 Check (eq_refl compute2 : (4611686018427387904, 1) = (4611686018427387904, 1)).
+
+Check (eq_refl : diveucl_21 1 1 0 = (0,0)).
+Check (eq_refl (0,0) <: diveucl_21 1 1 0 = (0,0)).
+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)).
diff --git a/test-suite/bugs/closed/bug_10031.v b/test-suite/bugs/closed/bug_10031.v
new file mode 100644
index 0000000000..15b53de00d
--- /dev/null
+++ b/test-suite/bugs/closed/bug_10031.v
@@ -0,0 +1,9 @@
+Require Import Int63 ZArith.
+
+Open Scope int63_scope.
+
+Goal False.
+cut (let (q, r) := (0, 0) in ([|q|], [|r|]) = (9223372036854775808%Z, 0%Z));
+  [discriminate| ].
+Fail (change (0, 0) with (diveucl_21 1 0 1); apply diveucl_21_spec).
+Abort.
diff --git a/theories/Numbers/Cyclic/Int63/Cyclic63.v b/theories/Numbers/Cyclic/Int63/Cyclic63.v
index 3b431d5b47..c03e6615cb 100644
--- a/theories/Numbers/Cyclic/Int63/Cyclic63.v
+++ b/theories/Numbers/Cyclic/Int63/Cyclic63.v
@@ -177,21 +177,6 @@ Proof.
  inversion W;rewrite Zmult_comm;trivial.
 Qed.
 
-Lemma diveucl_21_spec_aux : forall a1 a2 b,
-      wB/2 <= [|b|] ->
-      [|a1|] < [|b|] ->
-      let (q,r) := diveucl_21 a1 a2 b in
-      [|a1|] *wB+ [|a2|] = [|q|] * [|b|] + [|r|] /\
-      0 <= [|r|] < [|b|].
-Proof.
- intros a1 a2 b H1 H2;assert (W:= diveucl_21_spec a1 a2 b).
- assert (W1:= to_Z_bounded a1).
- assert ([|b|]>0) by (auto with zarith).
- generalize (Z_div_mod ([|a1|]*wB+[|a2|]) [|b|] H).
- destruct (diveucl_21 a1 a2 b);destruct (Z.div_eucl ([|a1|]*wB+[|a2|]) [|b|]).
- inversion W;rewrite (Zmult_comm [|b|]);trivial.
-Qed.
-
 Lemma shift_unshift_mod_2 : forall n p a, 0 <= p <= n ->
    ((a * 2 ^ (n - p)) mod (2^n) / 2 ^ (n - p)) mod (2^n) =
    a mod 2 ^ p.
diff --git a/theories/Numbers/Cyclic/Int63/Int63.v b/theories/Numbers/Cyclic/Int63/Int63.v
index eac26add03..3c96130bf3 100644
--- a/theories/Numbers/Cyclic/Int63/Int63.v
+++ b/theories/Numbers/Cyclic/Int63/Int63.v
@@ -387,7 +387,8 @@ 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
-   ([|q|],[|r|]) = Z.div_eucl ([|a1|] * wB + [|a2|]) [|b|].
+   let (q',r') := Z.div_eucl ([|a1|] * wB + [|a2|]) [|b|] in
+   [|q|] = Z.modulo q' wB /\ [|r|] = r'.
 
 Axiom addmuldiv_def_spec : forall p x y,
   addmuldiv p x y = addmuldiv_def p x y.
@@ -1413,12 +1414,51 @@ Proof.
  apply Z.le_trans with ([|ih|] * wB)%Z;try rewrite Z.pow_2_r; auto with zarith.
 Qed.
 
-Lemma div2_phi ih il j:
-  [|fst (diveucl_21 ih il j)|] = [|| WW ih il||] /[|j|].
-Proof.
- generalize (diveucl_21_spec ih il j).
- case diveucl_21; intros q r Heq.
- simpl zn2z_to_Z;unfold Z.div;rewrite <- Heq;trivial.
+Lemma diveucl_21_spec_aux : forall a1 a2 b,
+      wB/2 <= [|b|] ->
+      [|a1|] < [|b|] ->
+      let (q,r) := diveucl_21 a1 a2 b in
+      [|a1|] *wB+ [|a2|] = [|q|] * [|b|] + [|r|] /\
+      0 <= [|r|] < [|b|].
+Proof.
+ intros a1 a2 b H1 H2;assert (W:= diveucl_21_spec a1 a2 b).
+ assert (W1:= to_Z_bounded a1).
+ assert (W2:= to_Z_bounded a2).
+ assert (Wb:= to_Z_bounded b).
+ assert ([|b|]>0) by (auto with zarith).
+ 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''); 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.
+Qed.
+
+Lemma div2_phi ih il j: (2^62 <= [|j|] -> [|ih|] < [|j|] ->
+  [|fst (diveucl_21 ih il j)|] = [|| WW ih il||] /[|j|])%Z.
+Proof.
+ intros Hj Hj1.
+ generalize (diveucl_21_spec_aux ih il j Hj Hj1).
+ case diveucl_21; intros q r (Hq, Hr).
+ apply Zdiv_unique with [|r|]; auto with zarith.
+ simpl @fst; apply eq_trans with (1 := Hq); ring.
 Qed.
 
 Lemma sqrt2_step_correct rec ih il j:
@@ -1436,9 +1476,9 @@ Proof.
  case (to_Z_bounded il); intros Hil1 _.
  case (to_Z_bounded j); intros _ Hj1.
  assert (Hp3: (0 < [||WW ih il||])).
-  simpl zn2z_to_Z;apply Z.lt_le_trans with ([|ih|] * wB)%Z; auto with zarith.
+ {simpl zn2z_to_Z;apply Z.lt_le_trans with ([|ih|] * wB)%Z; auto with zarith.
   apply Zmult_lt_0_compat; auto with zarith.
-  refine (Z.lt_le_trans _ _ _ _ Hih); auto with zarith.
+  refine (Z.lt_le_trans _ _ _ _ Hih); auto with zarith. }
  cbv zeta.
  case_eq (ih < j)%int63;intros Heq.
  rewrite -> ltb_spec in Heq.
@@ -1450,28 +1490,28 @@ Proof.
  2: rewrite Zmult_comm, Z_div_plus_full_l; unfold base; auto with zarith.
  case (Zle_or_lt (2^(Z_of_nat size -1)) [|j|]); intros Hjj.
  case_eq (fst (diveucl_21 ih il j) < j)%int63;intros Heq0.
- 2: rewrite <-not_true_iff_false, ltb_spec, div2_phi in Heq0.
+ 2: rewrite <-not_true_iff_false, ltb_spec, (div2_phi _ _ _ Hjj Heq) in Heq0.
  2: split; auto; apply sqrt_test_true; auto with zarith.
- rewrite -> ltb_spec, div2_phi in Heq0.
+ rewrite -> ltb_spec, (div2_phi _ _ _ Hjj Heq) in Heq0.
  match goal with |- context[rec _ _ ?X] =>
   set (u := X)
  end.
  assert (H: [|u|] = ([|j|] + ([||WW ih il||])/([|j|]))/2).
-  unfold u; generalize (addc_spec j (fst (diveucl_21 ih il j)));
-  case addc;unfold interp_carry;rewrite div2_phi;simpl zn2z_to_Z.
-  intros i H;rewrite lsr_spec, H;trivial.
+ { unfold u; generalize (addc_spec j (fst (diveucl_21 ih il j)));
+  case addc;unfold interp_carry;rewrite (div2_phi _ _ _ Hjj Heq);simpl zn2z_to_Z.
+  { intros i H;rewrite lsr_spec, H;trivial. }
   intros i H;rewrite <- H.
   case (to_Z_bounded i); intros H1i H2i.
   rewrite -> add_spec, Zmod_small, lsr_spec.
-  change (1 * wB) with ([|(1 << (digits -1))|] * 2)%Z.
-  rewrite Z_div_plus_full_l; auto with zarith.
+  { change (1 * wB) with ([|(1 << (digits -1))|] * 2)%Z.
+    rewrite Z_div_plus_full_l; auto with zarith. }
   change wB with (2 * (wB/2))%Z; auto.
   replace [|(1 << (digits - 1))|] with (wB/2); auto.
   rewrite lsr_spec; auto.
   replace (2^[|1|]) with 2%Z; auto.
   split; auto with zarith.
   assert ([|i|]/2 < wB/2); auto with zarith.
-  apply Zdiv_lt_upper_bound; auto with zarith.
+  apply Zdiv_lt_upper_bound; auto with zarith. }
  apply Hrec; rewrite H; clear u H.
  assert (Hf1: 0 <= [||WW ih il||]/ [|j|]) by (apply Z_div_pos; auto with zarith).
  case (Zle_lt_or_eq 1 ([|j|])); auto with zarith; intros Hf2.