diff options
author | David Monniaux <David.Monniaux@univ-grenoble-alpes.fr> | 2021-12-13 20:40:25 +0100 |
---|---|---|
committer | David Monniaux <David.Monniaux@univ-grenoble-alpes.fr> | 2021-12-13 20:40:25 +0100 |
commit | 9007714f50a8ba49e2e6188cddada22a9fceed11 (patch) | |
tree | 6d06089baabea55ee10c76288fecf88b906eb31e /kvx/ExtValues.v | |
parent | a90018e65444bcd7a4c85ddd815d315cd852e120 (diff) | |
download | compcert-kvx-9007714f50a8ba49e2e6188cddada22a9fceed11.tar.gz compcert-kvx-9007714f50a8ba49e2e6188cddada22a9fceed11.zip |
begin div algo
Diffstat (limited to 'kvx/ExtValues.v')
-rw-r--r-- | kvx/ExtValues.v | 65 |
1 files changed, 65 insertions, 0 deletions
diff --git a/kvx/ExtValues.v b/kvx/ExtValues.v index b4e14898..c478f70b 100644 --- a/kvx/ExtValues.v +++ b/kvx/ExtValues.v @@ -754,3 +754,68 @@ Definition fmaddfs := triple_op_single (fun f1 f2 f3 => Float32.fma f2 f3 f1). Definition fmsubf := triple_op_float (fun f1 f2 f3 => Float.fma (Float.neg f2) f3 f1). Definition fmsubfs := triple_op_single (fun f1 f2 f3 => Float32.fma (Float32.neg f2) f3 f1). + +From Flocq Require Import Core Digits Operations Round Bracket Sterbenz + Binary Round_odd. +Require Import IEEE754_extra Zdiv Psatz Floats ExtFloats. + +Definition div_approx_reals (a b : Z) (x : R) := + let q:=ZnearestE x in + let r:=a-q*b in + if r <? 0 + then q-1 + else q. + +Lemma floor_ball1: + forall x : R, forall y : Z, + (Rabs (x - IZR y) < 1)%R -> Zfloor x = (y-1)%Z \/ Zfloor x = y. +Proof. + intros x y BALL. + apply Rabs_lt_inv in BALL. + case (Rcompare_spec x (IZR y)); intro CMP. + - left. apply Zfloor_imp. + ring_simplify (y-1+1). + rewrite minus_IZR. lra. + - subst. + rewrite Zfloor_IZR. right. reflexivity. + - right. apply Zfloor_imp. + rewrite plus_IZR. lra. +Qed. + +Theorem div_approx_reals_correct: + forall a b : Z, forall x : R, + b > 0 -> + (Rabs (x - IZR a/ IZR b) < 1/2)%R -> + div_approx_reals a b x = (a/b)%Z. +Proof. + intros a b x bPOS GAP. + assert (0 < IZR b)%R by (apply IZR_lt ; lia). + unfold div_approx_reals. + pose proof (Znearest_imp2 (fun x => negb (Z.even x)) x) as NEAR. + assert (Rabs (IZR (ZnearestE x) - IZR a/ IZR b) < 1)%R as BALL. + { pose proof (Rabs_triang (IZR (ZnearestE x) - x) + (x - IZR a/ IZR b)) as TRI. + ring_simplify (IZR (ZnearestE x) - x + (x - IZR a / IZR b))%R in TRI. + lra. + } + clear GAP NEAR. + rewrite Rabs_minus_sym in BALL. + pose proof (floor_ball1 _ _ BALL) as FLOOR. + clear BALL. + rewrite Zfloor_div in FLOOR by lia. + pose proof (Z_div_mod_eq_full a b) as DIV_MOD. + assert (0 < b) as bPOS' by lia. + pose proof (Z.mod_pos_bound a b bPOS') as MOD_BOUNDS. + case Z.ltb_spec; intro; destruct FLOOR; lia. +Qed. + +Definition my_div (a b : val) := + let b_d := Val.maketotal (Val.floatofintu b) in + let invb_d := Val.floatofsingle (invfs (Val.maketotal (Val.singleofintu b))) in + let alpha := fmsubf (Vfloat ExtFloat.one) invb_d b_d in + let x := fmaddf invb_d alpha invb_d in + Val.mulf (Val.maketotal (Val.floatofintu a)) x. + +(* +Compute (my_div (Vint (Int.repr 3)) (Vint (Int.repr 5))). +*) |