Merge PR #10046: [primitive integers] Make div21 implems consistent with its specification

Ack-by: Zimmi48
Ack-by: herbelin
Ack-by: maximedenes
Ack-by: proux01
Reviewed-by: vbgl

10 files changed, 160 insertions, 111 deletions

diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index 2293ae9dfd..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();
@@ -1374,40 +1361,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;
       }
@@ -1616,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 d431dc1e5c..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 */
@@ -109,37 +110,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.