Division_u16.v 3.11 KB
Newer Older
1
Require Import Reals Psatz.
2
Require Import Fcore Gappa.Gappa_tactic.
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128

Open Scope R_scope.

Lemma Rdiv_le_l :
  forall a b c,
  0 < a -> b <= c / a -> a * b <= c.
Proof.
intros a b c Ha Hb.
apply Rmult_le_reg_r with (/a).
now apply Rinv_0_lt_compat.
replace (a * b * /a) with b.
exact Hb.
field.
now apply Rgt_not_eq.
Qed.

Lemma FIX_format_Z2R :
  forall beta x, FIX_format beta 0 (Z2R x).
Proof.
intros beta x.
exists (Float beta x 0).
split.
apply sym_eq, Rmult_1_r.
apply eq_refl.
Qed.

Lemma Rmult_le_lt_reg_l :
  forall a b c d, 0 < a ->
  a * b <= a * c < a * d ->
  b <= c < d.
Proof.
intros a b c d Ha H.
split.
now apply Rmult_le_reg_l with (1 := Ha).
now apply Rmult_lt_reg_l with (1 := Ha).
Qed.

Lemma Rplus_le_lt_reg_r :
  forall a b c d,
  b + a <= c + a < d + a ->
  b <= c < d.
Proof.
intros a b c d H.
split.
now apply Rplus_le_reg_r with a.
now apply Rplus_lt_reg_r with a.
Qed.

Coercion Z2R : Z >-> R.

Lemma Z2R_le_le :
  forall a b c,
  (a <= b <= c)%Z ->
  a <= b <= c.
Proof.
intros a b c.
now split ; apply Z2R_le.
Qed.

Notation pow2 := (bpow radix2).
Notation register_fmt := (FLT_exp (-65597) 64).
Notation fma x y z := (round radix2 register_fmt rndNE (x * y + z)).
Notation fnma x y z := (round radix2 register_fmt rndNE (z - x * y)).

Axiom frcpa : R -> R.
Axiom frcpa_spec : forall x : R,
  1 <= Rabs x <= 65536 ->
  generic_format radix2 (FLT_exp (-65597) 11) (frcpa x) /\
  Rabs (frcpa x - 1 / x) <= 4433 * pow2 (-21) * Rabs (1 / x).

Definition div_u16 a b :=
  let y0 := frcpa b in
  let q0 := fma a y0 0 in
  let e0 := fnma b y0 (1 + pow2 (-17)) in
  let q1 := fma e0 q0 q0 in
  Zfloor q1.

Lemma div_u16_spec : forall a b,
  (1 <= a <= 65535)%Z ->
  (1 <= b <= 65535)%Z ->
  div_u16 a b = (a / b)%Z.
Proof.
intros a b Ba Bb.
unfold div_u16.
set (y0 := frcpa b).
set (q0 := fma a y0 0).
set (e0 := fnma b y0 (1 + pow2 (-17))).
set (q1 := fma e0 q0 q0).
apply Zfloor_imp.
rewrite Z2R_plus.
simpl.
apply Rmult_le_lt_reg_l with b.
  apply (Z2R_lt 0) ; lia.
apply Rplus_le_lt_reg_r with (-a).
replace (b * (a / b)%Z + - a) with (-(a - (a / b)%Z * b)) by ring.
replace (b * ((a / b)%Z + 1) + - a) with (b - (a - (a / b)%Z * b)) by ring.
rewrite <- Z2R_mult, <- Z2R_minus.
rewrite <- Zmod_eq_full by lia.
cut (0 <= b * q1 - a < 1).
  cut (0 <= Zmod a b <= b - 1).
    lra.
  change 1 with (Z2R 1).
  rewrite <- Z2R_minus.
  apply (Z2R_le_le 0).
  generalize (Z_mod_lt a b).
  lia.
assert (Ba': 1 <= a <= 65535).
  change (1%Z <= a <= 65535%Z).
  now split ; apply Z2R_le.
assert (Bb': 1 <= b <= 65535).
  change (1%Z <= b <= 65535%Z).
  now split ; apply Z2R_le.
refine (let '(conj H1 H2) := _ in conj H1 (Rnot_le_lt _ _ H2)).
set (err := (q1 - a / b) / (a / b)).
replace (b * q1 - a) with (a * err) by abstract (unfold err ; field ; lra).
set (Mq0 := a * y0 + 0).
set (Me0 := 1 + pow2 (-17) - b * y0).
set (Mq1 := Me0 * Mq0 + Mq0).
set (eps0 := (y0 - 1 / b) / (1 / b)).
assert (H: (Mq1 - a / b) / (a / b) = -(eps0 * eps0) + (1 + eps0) * pow2 (-17))
  by abstract (unfold Mq1, Me0, Mq0, eps0 ; field ; lra).
revert H.
unfold Mq1, Me0, Mq0, eps0, err, q1, e0, q0, y0.
generalize (frcpa_spec b) (FIX_format_Z2R radix2 a) (FIX_format_Z2R radix2 b).
gappa.
Qed.