Turn 64-bit integer division and modulus by constants into multiply-high

This trick was already implemented for 32-bit integer division and modulus.  Here we extend it to the 64-bit case.

For 32-bit target processors, the runtime library must implement 64-bit multiply-high (signed and unsigned).  Tentative implementations are provided for IA32 and PowerPC, but need testing.

6 files changed, 337 insertions, 17 deletions

diff --git a/backend/SelectDiv.vp b/backend/SelectDiv.vp
index 1fc0b689..5cc66322 100644
--- a/backend/SelectDiv.vp
+++ b/backend/SelectDiv.vp
@@ -206,14 +206,23 @@ Context {hf: helper_functions}.
 Definition modl_from_divl (equo: expr) (n: int64) :=
   subl (Eletvar O) (mullimm n equo).
 
+Definition divlu_mull (p: Z) (m: Z) :=
+  shrluimm (mullhu (Eletvar O) (Int64.repr m)) (Int.repr p).
+
 Definition divlu (e1 e2: expr) :=
   match is_longconst e2, is_longconst e1 with
   | Some n2, Some n1 => longconst (Int64.divu n1 n2)
   | Some n2, _ =>
       match Int64.is_power2' n2 with
       | Some l => shrluimm e1 l
-      | None   => divlu_base e1 e2
-      end (* TODO: multiply-high *)
+      | None   => if optim_for_size tt then
+                    divlu_base e1 e2
+                  else
+                    match divlu_mul_params (Int64.unsigned n2) with
+                    | None => divlu_base e1 e2
+                    | Some(p, m) => Elet e1 (divlu_mull p m)
+                    end
+      end
   | _, _ => divlu_base e1 e2
   end.
 
@@ -223,19 +232,41 @@ Definition modlu (e1 e2: expr) :=
   | Some n2, _ =>
       match Int64.is_power2 n2 with
       | Some l => andl e1 (longconst (Int64.sub n2 Int64.one))
-      | None   => modlu_base e1 e2
+      | None   => if optim_for_size tt then
+                    modlu_base e1 e2
+                  else
+                    match divlu_mul_params (Int64.unsigned n2) with
+                    | None => modlu_base e1 e2
+                    | Some(p, m) => Elet e1 (modl_from_divl (divlu_mull p m) n2)
+                    end
       end
   | _, _ => modlu_base e1 e2
   end.
 
+Definition divls_mull (p: Z) (m: Z) :=
+  let e2 := 
+    mullhs (Eletvar O) (Int64.repr m) in
+  let e3 :=
+    if zlt m Int64.half_modulus then e2 else addl e2 (Eletvar O) in
+  addl (shrlimm e3 (Int.repr p))
+       (shrluimm (Eletvar O) (Int.repr (Int64.zwordsize - 1))).
+
 Definition divls (e1 e2: expr) :=
   match is_longconst e2, is_longconst e1 with
   | Some n2, Some n1 => longconst (Int64.divs n1 n2)
   | Some n2, _ =>
       match Int64.is_power2' n2 with
-      | Some l => if Int.ltu l (Int.repr 63) then shrxlimm e1 l else divls_base e1 e2
-      | None   => divls_base e1 e2
-      end (* TODO: multiply-high *)
+      | Some l => if Int.ltu l (Int.repr 63)
+                  then shrxlimm e1 l
+                  else divls_base e1 e2
+      | None   => if optim_for_size tt then
+                    divls_base e1 e2
+                  else
+                    match divls_mul_params (Int64.signed n2) with
+                    | None => divls_base e1 e2
+                    | Some(p, m) => Elet e1 (divls_mull p m)
+                    end
+      end
   | _, _ => divls_base e1 e2
   end.
 
@@ -247,7 +278,13 @@ Definition modls (e1 e2: expr) :=
       | Some l => if Int.ltu l (Int.repr 63)
                   then Elet e1 (modl_from_divl (shrxlimm (Eletvar O) l) n2)
                   else modls_base e1 e2
-      | None   => modls_base e1 e2
+      | None   => if optim_for_size tt then
+                    modls_base e1 e2
+                  else
+                    match divls_mul_params (Int64.signed n2) with
+                    | None => modls_base e1 e2
+                    | Some(p, m) => Elet e1 (modl_from_divl (divls_mull p m) n2)
+                    end
       end
   | _, _ => modls_base e1 e2
   end.
diff --git a/backend/SelectDivproof.v b/backend/SelectDivproof.v
index 441f69b1..41db3c70 100644
--- a/backend/SelectDivproof.v
+++ b/backend/SelectDivproof.v
@@ -317,6 +317,165 @@ Proof.
   assert (32 < Int.max_unsigned) by (compute; auto). omega.
 Qed.
 
+(** Same, for 64-bit integers *)
+
+Lemma divls_mul_params_sound:
+  forall d m p,
+  divls_mul_params d = Some(p, m) ->
+  0 <= m < Int64.modulus /\ 0 <= p < 64 /\
+  forall n,
+  Int64.min_signed <= n <= Int64.max_signed ->
+  Z.quot n d = Zdiv (m * n) (two_p (64 + p)) + (if zlt n 0 then 1 else 0).
+Proof with (try discriminate).
+  unfold divls_mul_params; intros d m' p'.
+  destruct (find_div_mul_params Int64.wordsize
+               (Int64.half_modulus - Int64.half_modulus mod d - 1) d 64)
+  as [[p m] | ]...
+  generalize (p - 64). intro p1.  
+  destruct (zlt 0 d)...
+  destruct (zlt (two_p (64 + p1)) (m * d))...
+  destruct (zle (m * d) (two_p (64 + p1) + two_p (p1 + 1)))...
+  destruct (zle 0 m)...
+  destruct (zlt m Int64.modulus)...
+  destruct (zle 0 p1)...
+  destruct (zlt p1 64)...
+  intros EQ; inv EQ. 
+  split. auto. split. auto. intros.
+  replace (64 + p') with (63 + (p' + 1)) by omega.
+  apply Zquot_mul; try omega.
+  replace (63 + (p' + 1)) with (64 + p') by omega. omega.
+  change (Int64.min_signed <= n < Int64.half_modulus).
+  unfold Int64.max_signed in H. omega.
+Qed.
+
+Lemma divlu_mul_params_sound:
+  forall d m p,
+  divlu_mul_params d = Some(p, m) ->
+  0 <= m < Int64.modulus /\ 0 <= p < 64 /\
+  forall n,
+  0 <= n < Int64.modulus ->
+  Zdiv n d = Zdiv (m * n) (two_p (64 + p)).
+Proof with (try discriminate).
+  unfold divlu_mul_params; intros d m' p'.
+  destruct (find_div_mul_params Int64.wordsize
+               (Int64.modulus - Int64.modulus mod d - 1) d 64)
+  as [[p m] | ]...
+  generalize (p - 64); intro p1.
+  destruct (zlt 0 d)...
+  destruct (zle (two_p (64 + p1)) (m * d))...
+  destruct (zle (m * d) (two_p (64 + p1) + two_p p1))...
+  destruct (zle 0 m)...
+  destruct (zlt m Int64.modulus)...
+  destruct (zle 0 p1)...
+  destruct (zlt p1 64)...
+  intros EQ; inv EQ.
+  split. auto. split. auto. intros.
+  apply Zdiv_mul_pos; try omega. assumption.
+Qed.
+
+Remark int64_shr'_div_two_p:
+  forall x y, Int64.shr' x y = Int64.repr (Int64.signed x / two_p (Int.unsigned y)).
+Proof.
+  intros; unfold Int64.shr'. rewrite Int64.Zshiftr_div_two_p; auto. generalize (Int.unsigned_range y); omega.
+Qed.
+
+Lemma divls_mul_shift_gen:
+  forall x y m p,
+  divls_mul_params (Int64.signed y) = Some(p, m) ->
+  0 <= m < Int64.modulus /\ 0 <= p < 64 /\
+  Int64.divs x y = Int64.add (Int64.shr' (Int64.repr ((Int64.signed x * m) / Int64.modulus)) (Int.repr p))
+                             (Int64.shru x (Int64.repr 63)).
+Proof.
+  intros. set (n := Int64.signed x). set (d := Int64.signed y) in *.
+  exploit divls_mul_params_sound; eauto. intros (A & B & C).
+  split. auto. split. auto.
+  unfold Int64.divs. fold n; fold d. rewrite C by (apply Int64.signed_range).
+  rewrite two_p_is_exp by omega. rewrite <- Zdiv_Zdiv.
+  rewrite Int64.shru_lt_zero. unfold Int64.add. apply Int64.eqm_samerepr. apply Int64.eqm_add.
+  rewrite int64_shr'_div_two_p. apply Int64.eqm_unsigned_repr_r. apply Int64.eqm_refl2.
+  rewrite Int.unsigned_repr. f_equal.
+  rewrite Int64.signed_repr. rewrite Int64.modulus_power. f_equal. ring.
+  cut (Int64.min_signed <= n * m / Int64.modulus < Int64.half_modulus).
+  unfold Int64.max_signed; omega.
+  apply Zdiv_interval_1. generalize Int64.min_signed_neg; omega. apply Int64.half_modulus_pos.
+  apply Int64.modulus_pos.
+  split. apply Zle_trans with (Int64.min_signed * m). apply Zmult_le_compat_l_neg. omega. generalize Int64.min_signed_neg; omega.
+  apply Zmult_le_compat_r. unfold n; generalize (Int64.signed_range x); tauto. tauto.
+  apply Zle_lt_trans with (Int64.half_modulus * m).
+  apply Zmult_le_compat_r. generalize (Int64.signed_range x); unfold n, Int64.max_signed; omega. tauto.
+  apply Zmult_lt_compat_l. generalize Int64.half_modulus_pos; omega. tauto.
+  assert (64 < Int.max_unsigned) by (compute; auto). omega.
+  unfold Int64.lt; fold n. rewrite Int64.signed_zero. destruct (zlt n 0); apply Int64.eqm_unsigned_repr.
+  apply two_p_gt_ZERO. omega.
+  apply two_p_gt_ZERO. omega.
+Qed.
+
+Theorem divls_mul_shift_1:
+  forall x y m p,
+  divls_mul_params (Int64.signed y) = Some(p, m) ->
+  m < Int64.half_modulus ->
+  0 <= p < 64 /\
+  Int64.divs x y = Int64.add (Int64.shr' (Int64.mulhs x (Int64.repr m)) (Int.repr p))
+                             (Int64.shru' x (Int.repr 63)).
+Proof.
+  intros. exploit divls_mul_shift_gen; eauto. instantiate (1 := x).
+  intros (A & B & C). split. auto. rewrite C.
+  unfold Int64.mulhs. rewrite Int64.signed_repr. auto.
+  generalize Int64.min_signed_neg; unfold Int64.max_signed; omega.
+Qed.
+
+Theorem divls_mul_shift_2:
+  forall x y m p,
+  divls_mul_params (Int64.signed y) = Some(p, m) ->
+  m >= Int64.half_modulus ->
+  0 <= p < 64 /\
+  Int64.divs x y = Int64.add (Int64.shr' (Int64.add (Int64.mulhs x (Int64.repr m)) x) (Int.repr p))
+                             (Int64.shru' x (Int.repr 63)).
+Proof.
+  intros. exploit divls_mul_shift_gen; eauto. instantiate (1 := x).
+  intros (A & B & C). split. auto. rewrite C. f_equal. f_equal.
+  rewrite Int64.add_signed. unfold Int64.mulhs. set (n := Int64.signed x).
+  transitivity (Int64.repr (n * (m - Int64.modulus) / Int64.modulus + n)).
+  f_equal.
+  replace (n * (m - Int64.modulus)) with (n * m +  (-n) * Int64.modulus) by ring.
+  rewrite Z_div_plus. ring. apply Int64.modulus_pos.
+  apply Int64.eqm_samerepr. apply Int64.eqm_add; auto with ints.
+  apply Int64.eqm_sym. eapply Int64.eqm_trans. apply Int64.eqm_signed_unsigned.
+  apply Int64.eqm_unsigned_repr_l. apply Int64.eqm_refl2. f_equal. f_equal.
+  rewrite Int64.signed_repr_eq. rewrite Zmod_small by assumption.
+  apply zlt_false. omega.
+Qed.
+
+Remark int64_shru'_div_two_p:
+  forall x y, Int64.shru' x y = Int64.repr (Int64.unsigned x / two_p (Int.unsigned y)).
+Proof.
+  intros; unfold Int64.shru'. rewrite Int64.Zshiftr_div_two_p; auto. generalize (Int.unsigned_range y); omega.
+Qed.
+
+Theorem divlu_mul_shift:
+  forall x y m p,
+  divlu_mul_params (Int64.unsigned y) = Some(p, m) ->
+  0 <= p < 64 /\
+  Int64.divu x y = Int64.shru' (Int64.mulhu x (Int64.repr m)) (Int.repr p).
+Proof.
+  intros. exploit divlu_mul_params_sound; eauto. intros (A & B & C).
+  split. auto.
+  rewrite int64_shru'_div_two_p. rewrite Int.unsigned_repr.
+  unfold Int64.divu, Int64.mulhu. f_equal. rewrite C by apply Int64.unsigned_range.
+  rewrite two_p_is_exp by omega. rewrite <- Zdiv_Zdiv by (apply two_p_gt_ZERO; omega).
+  f_equal. rewrite (Int64.unsigned_repr m).
+  rewrite Int64.unsigned_repr. f_equal. ring.
+  cut (0 <= Int64.unsigned x * m / Int64.modulus < Int64.modulus).
+  unfold Int64.max_unsigned; omega.
+  apply Zdiv_interval_1. omega. compute; auto. compute; auto.
+  split. simpl. apply Z.mul_nonneg_nonneg. generalize (Int64.unsigned_range x); omega. omega.
+  apply Zle_lt_trans with (Int64.modulus * m).
+  apply Zmult_le_compat_r. generalize (Int64.unsigned_range x); omega. omega.
+  apply Zmult_lt_compat_l. compute; auto. omega.
+  unfold Int64.max_unsigned; omega.
+  assert (64 < Int.max_unsigned) by (compute; auto). omega.
+Qed.
+
 (** * Correctness of the smart constructors for division and modulus *)
 
 Section CMCONSTRS.
@@ -559,6 +718,22 @@ Proof.
   simpl in B1; inv B1. simpl in B2; inv B2. exact A2.
 Qed.
 
+Lemma eval_divlu_mull:
+  forall le x y p M,
+  divlu_mul_params (Int64.unsigned y) = Some(p, M) ->
+  nth_error le O = Some (Vlong x) ->
+  eval_expr ge sp e m le (divlu_mull p M) (Vlong (Int64.divu x y)).
+Proof.
+  intros. unfold divlu_mull. exploit (divlu_mul_shift x); eauto. intros [A B].
+  assert (A0: eval_expr ge sp e m le (Eletvar O) (Vlong x)) by (constructor; auto).
+  exploit eval_mullhu. eauto. eexact A0. instantiate (1 := Int64.repr M). intros (v1 & A1 & B1).
+  exploit eval_shrluimm. eauto. eexact A1. instantiate (1 := Int.repr p). intros (v2 & A2 & B2).
+  simpl in B1; inv B1. simpl in B2. replace (Int.ltu (Int.repr p) Int64.iwordsize') with true in B2. inv B2.
+  rewrite B. assumption.
+  unfold Int.ltu. rewrite Int.unsigned_repr. rewrite zlt_true; auto. tauto.
+  assert (64 < Int.max_unsigned) by (compute; auto). omega.
+Qed.
+
 Theorem eval_divlu:
   forall le a b x y z,
   eval_expr ge sp e m le a x ->
@@ -575,7 +750,12 @@ Proof.
   econstructor; split. apply eval_longconst. constructor.
 + destruct (Int64.is_power2' n2) as [l|] eqn:POW.
 * exploit Val.divlu_pow2; eauto. intros EQ; subst z. apply eval_shrluimm; auto.
-* eapply eval_divlu_base; eauto.     
+* destruct (Compopts.optim_for_size tt). eapply eval_divlu_base; eauto.
+  destruct (divlu_mul_params (Int64.unsigned n2)) as [[p M]|] eqn:PARAMS. 
+** destruct x; simpl in H1; try discriminate.
+   destruct (Int64.eq n2 Int64.zero); inv H1.
+   econstructor; split; eauto. econstructor. eauto. eapply eval_divlu_mull; eauto.  
+** eapply eval_divlu_base; eauto.     
 - eapply eval_divlu_base; eauto.
 Qed.
 
@@ -595,10 +775,56 @@ Proof.
   econstructor; split. apply eval_longconst. constructor.
 + destruct (Int64.is_power2 n2) as [l|] eqn:POW.
 * exploit Val.modlu_pow2; eauto. intros EQ; subst z. eapply eval_andl; eauto. apply eval_longconst.
-* eapply eval_modlu_base; eauto.     
+* destruct (Compopts.optim_for_size tt). eapply eval_modlu_base; eauto.
+  destruct (divlu_mul_params (Int64.unsigned n2)) as [[p M]|] eqn:PARAMS. 
+** destruct x; simpl in H1; try discriminate.
+   destruct (Int64.eq n2 Int64.zero) eqn:Z; inv H1.
+   rewrite Int64.modu_divu. 
+    econstructor; split; eauto. econstructor. eauto.
+    eapply eval_modl_from_divl; eauto.  
+    eapply eval_divlu_mull; eauto.
+    red; intros; subst n2; discriminate Z.  
+** eapply eval_modlu_base; eauto.     
 - eapply eval_modlu_base; eauto.
 Qed.
 
+Lemma eval_divls_mull:
+  forall le x y p M,
+  divls_mul_params (Int64.signed y) = Some(p, M) ->
+  nth_error le O = Some (Vlong x) ->
+  eval_expr ge sp e m le (divls_mull p M) (Vlong (Int64.divs x y)).
+Proof.
+  intros. unfold divls_mull.
+  assert (A0: eval_expr ge sp e m le (Eletvar O) (Vlong x)).
+  { constructor; auto. }
+  exploit eval_mullhs. eauto. eexact A0. instantiate (1 := Int64.repr M).  intros (v1 & A1 & B1).
+  exploit eval_addl. eexact A1. eexact A0. intros (v2 & A2 & B2).
+  exploit eval_shrluimm. eauto. eexact A0. instantiate (1 := Int.repr 63). intros (v3 & A3 & B3).
+  set (a4 := if zlt M Int64.half_modulus
+             then mullhs (Eletvar 0) (Int64.repr M)
+             else addl (mullhs (Eletvar 0) (Int64.repr M)) (Eletvar 0)).
+  set (v4 := if zlt M Int64.half_modulus then v1 else v2).
+  assert (A4: eval_expr ge sp e m le a4 v4). 
+  { unfold a4, v4; destruct (zlt M Int64.half_modulus); auto. }
+  exploit eval_shrlimm. eauto. eexact A4. instantiate (1 := Int.repr p). intros (v5 & A5 & B5).
+  exploit eval_addl. eexact A5. eexact A3. intros (v6 & A6 & B6).
+  assert (RANGE: forall x, 0 <= x < 64 -> Int.ltu (Int.repr x) Int64.iwordsize' = true).
+  { intros. unfold Int.ltu. rewrite Int.unsigned_repr. rewrite zlt_true by tauto. auto.
+    assert (64 < Int.max_unsigned) by (compute; auto). omega. }
+  simpl in B1; inv B1.
+  simpl in B2; inv B2.
+  simpl in B3; rewrite RANGE in B3 by omega; inv B3.
+  destruct (zlt M Int64.half_modulus).
+- exploit (divls_mul_shift_1 x); eauto. intros [A B].
+  simpl in B5; rewrite RANGE in B5 by auto; inv B5.
+  simpl in B6; inv B6. 
+  rewrite B; exact A6.
+- exploit (divls_mul_shift_2 x); eauto. intros [A B].
+  simpl in B5; rewrite RANGE in B5 by auto; inv B5.
+  simpl in B6; inv B6. 
+  rewrite B; exact A6.
+Qed.
+
 Theorem eval_divls:
   forall le a b x y z,
   eval_expr ge sp e m le a x ->
@@ -618,8 +844,15 @@ Proof.
 + destruct (Int64.is_power2' n2) as [l|] eqn:POW.
 * destruct (Int.ltu l (Int.repr 63)) eqn:LT.
 ** exploit Val.divls_pow2; eauto. intros EQ. eapply eval_shrxlimm; eauto.
-** eapply eval_divls_base; eauto.  
-* eapply eval_divls_base; eauto.     
+** eapply eval_divls_base; eauto.
+* destruct (Compopts.optim_for_size tt). eapply eval_divls_base; eauto.
+  destruct (divls_mul_params (Int64.signed n2)) as [[p M]|] eqn:PARAMS. 
+** destruct x; simpl in H1; try discriminate.
+   destruct (Int64.eq n2 Int64.zero
+             || Int64.eq i (Int64.repr Int64.min_signed) && Int64.eq n2 Int64.mone); inv H1.
+   econstructor; split; eauto. econstructor. eauto.
+   eapply eval_divls_mull; eauto.
+** eapply eval_divls_base; eauto.
 - eapply eval_divls_base; eauto.
 Qed.
 
@@ -654,7 +887,16 @@ Proof.
   econstructor. eauto. eapply eval_modl_from_divl. eexact A1. reflexivity.      
   rewrite Int64.mods_divs. auto.
 **eapply eval_modls_base; eauto.
-* eapply eval_modls_base; eauto.     
+* destruct (Compopts.optim_for_size tt). eapply eval_modls_base; eauto.
+  destruct (divls_mul_params (Int64.signed n2)) as [[p M]|] eqn:PARAMS. 
+** destruct x; simpl in H1; try discriminate.
+   destruct (Int64.eq n2 Int64.zero
+             || Int64.eq i (Int64.repr Int64.min_signed) && Int64.eq n2 Int64.mone); inv H1.
+   econstructor; split; eauto. econstructor. eauto.
+   rewrite Int64.mods_divs. 
+   eapply eval_modl_from_divl; auto.
+   eapply eval_divls_mull; eauto.
+** eapply eval_modls_base; eauto.     
 - eapply eval_modls_base; eauto.
 Qed.
 
diff --git a/backend/Selection.v b/backend/Selection.v
index 5cb5d119..abda1d95 100644
--- a/backend/Selection.v
+++ b/backend/Selection.v
@@ -389,10 +389,13 @@ Definition get_helpers (defmap: PTree.t globdef) : res helper_functions :=
   do i64_shl <- lookup_helper globs "__i64_shl" sig_li_l ;
   do i64_shr <- lookup_helper globs "__i64_shr" sig_li_l ;
   do i64_sar <- lookup_helper globs "__i64_sar" sig_li_l ;
+  do i64_umulh <- lookup_helper globs "__i64_umulh" sig_ll_l ;
+  do i64_smulh <- lookup_helper globs "__i64_smulh" sig_ll_l ;
   OK (mk_helper_functions
      i64_dtos i64_dtou i64_stod i64_utod i64_stof i64_utof
      i64_sdiv i64_udiv i64_smod i64_umod
-     i64_shl i64_shr i64_sar).
+     i64_shl i64_shr i64_sar
+     i64_umulh i64_smulh).
 
 (** Conversion of programs. *)
 
diff --git a/backend/SplitLong.vp b/backend/SplitLong.vp
index 5891adef..f7eeebd0 100644
--- a/backend/SplitLong.vp
+++ b/backend/SplitLong.vp
@@ -38,7 +38,9 @@ Class helper_functions := mk_helper_functions {
   i64_umod: ident;                      (**r unsigned remainder *)
   i64_shl: ident;                       (**r shift left *)
   i64_shr: ident;                       (**r shift right unsigned *)
-  i64_sar: ident                        (**r shift right signed *)
+  i64_sar: ident;                       (**r shift right signed *)
+  i64_umulh: ident;                     (**r unsigned multiply high *)
+  i64_smulh: ident;                     (**r signed multiply high *)
 }.
 
 Definition sig_l_l := mksignature (Tlong :: nil) (Some Tlong) cc_default.
@@ -255,6 +257,11 @@ Definition mull (e1 e2: expr) :=
   | _, _ => mull_base e1 e2
   end.
 
+Definition mullhu (e1: expr) (n2: int64) :=
+  Eexternal i64_umulh sig_ll_l (e1 ::: longconst n2 ::: Enil).
+Definition mullhs (e1: expr) (n2: int64) :=
+  Eexternal i64_smulh sig_ll_l (e1 ::: longconst n2 ::: Enil).
+
 Definition shrxlimm (e: expr) (n: int) :=
   if Int.eq n Int.zero then e else
     Elet e (shrlimm (addl (Eletvar O)
diff --git a/backend/SplitLongproof.v b/backend/SplitLongproof.v
index a10ee3f7..31f5db67 100644
--- a/backend/SplitLongproof.v
+++ b/backend/SplitLongproof.v
@@ -48,7 +48,9 @@ Axiom i64_helpers_correct :
  /\ (forall x y z, Val.modlu x y = Some z -> external_implements "__i64_umod" sig_ll_l (x::y::nil) z)
  /\ (forall x y, external_implements "__i64_shl" sig_li_l (x::y::nil) (Val.shll x y))
  /\ (forall x y, external_implements "__i64_shr" sig_li_l (x::y::nil) (Val.shrlu x y))
- /\ (forall x y, external_implements "__i64_sar" sig_li_l (x::y::nil) (Val.shrl x y)).
+ /\ (forall x y, external_implements "__i64_sar" sig_li_l (x::y::nil) (Val.shrl x y))
+ /\ (forall x y, external_implements "__i64_umulh" sig_ll_l (x::y::nil) (Val.mullhu x y))
+ /\ (forall x y, external_implements "__i64_smulh" sig_ll_l (x::y::nil) (Val.mullhs x y)).
 
 Definition helper_declared {F V: Type} (p: AST.program (AST.fundef F) V) (id: ident) (name: string) (sg: signature) : Prop :=
   (prog_defmap p)!id = Some (Gfun (External (EF_runtime name sg))).
@@ -66,7 +68,9 @@ Definition helper_functions_declared {F V: Type} (p: AST.program (AST.fundef F)
   /\ helper_declared p i64_umod "__i64_umod" sig_ll_l
   /\ helper_declared p i64_shl "__i64_shl" sig_li_l
   /\ helper_declared p i64_shr "__i64_shr" sig_li_l
-  /\ helper_declared p i64_sar "__i64_sar" sig_li_l.
+  /\ helper_declared p i64_sar "__i64_sar" sig_li_l
+  /\ helper_declared p i64_umulh "__i64_umulh" sig_ll_l
+  /\ helper_declared p i64_smulh "__i64_smulh" sig_ll_l.
 
 (** * Correctness of the instruction selection functions for 64-bit operators *)
 
@@ -823,6 +827,20 @@ Proof.
 - apply eval_mull_base; auto.
 Qed.
 
+Theorem eval_mullhu: 
+  forall n, unary_constructor_sound (fun a => mullhu a n) (fun v => Val.mullhu v (Vlong n)).
+Proof.
+  unfold mullhu; intros; red; intros. econstructor; split; eauto. 
+  eapply eval_helper_2; eauto. apply eval_longconst. DeclHelper; eauto. UseHelper. 
+Qed.
+
+Theorem eval_mullhs: 
+  forall n, unary_constructor_sound (fun a => mullhs a n) (fun v => Val.mullhs v (Vlong n)).
+Proof.
+  unfold mullhs; intros; red; intros. econstructor; split; eauto. 
+  eapply eval_helper_2; eauto. apply eval_longconst. DeclHelper; eauto. UseHelper. 
+Qed.
+
 Theorem eval_shrxlimm:
   forall le a n x z,
   Archi.ptr64 = false ->
diff --git a/backend/ValueDomain.v b/backend/ValueDomain.v
index bf88a450..be8bcccc 100644
--- a/backend/ValueDomain.v
+++ b/backend/ValueDomain.v
@@ -1866,6 +1866,18 @@ Lemma mull_sound:
   forall v x w y, vmatch v x -> vmatch w y -> vmatch (Val.mull v w) (mull x y).
 Proof (binop_long_sound Int64.mul).
 
+Definition mullhs := binop_long Int64.mulhs.
+
+Lemma mullhs_sound:
+  forall v x w y, vmatch v x -> vmatch w y -> vmatch (Val.mullhs v w) (mullhs x y).
+Proof (binop_long_sound Int64.mulhs).
+
+Definition mullhu := binop_long Int64.mulhu.
+
+Lemma mullhu_sound:
+  forall v x w y, vmatch v x -> vmatch w y -> vmatch (Val.mullhu v w) (mullhu x y).
+Proof (binop_long_sound Int64.mulhu).
+
 Definition divls (v w: aval) :=
   match w, v with
   | L i2, L i1 =>
@@ -4559,7 +4571,8 @@ Hint Resolve cnot_sound symbol_address_sound
        divs_sound divu_sound mods_sound modu_sound shrx_sound
        shll_sound shrl_sound shrlu_sound
        andl_sound orl_sound xorl_sound notl_sound roll_sound rorl_sound
-       negl_sound addl_sound subl_sound mull_sound
+       negl_sound addl_sound subl_sound
+       mull_sound mullhs_sound mullhu_sound
        divls_sound divlu_sound modls_sound modlu_sound shrxl_sound
        negf_sound absf_sound
        addf_sound subf_sound mulf_sound divf_sound