path: root/kvx/FPDivision64.v

(*
This needs a special gappa script

#!/bin/sh
/home/monniaux/.opam/4.12.0+flambda/bin/gappa -Eprecision=100 "$@"

in PATH before the normal gappa
 *)

From Flocq Require Import Core Digits Operations Round Bracket Sterbenz
                          Binary Round_odd Bits.
Require Archi.
Require Import Coqlib.
Require Import Compopts.
Require Import AST.
Require Import Integers.
Require Import Floats.
Require Import Op.
Require Import CminorSel.
Require Import OpHelpers.
Require Import ExtFloats.
Require Import DecBoolOps.
Require Import Chunks.
Require Import Builtins.
Require Import Values Globalenvs.
Require Compopts.
Require Import Psatz.
Require Import IEEE754_extra.

From Gappa Require Import Gappa_tactic.

Definition approx_inv_longu b :=
  let invb_s := ExtValues.invfs (Val.maketotal (Val.singleoflongu b)) in
  let invb_d := Val.floatofsingle invb_s in
  let b_d := Val.maketotal (Val.floatoflongu b) in
  let one := Vfloat (ExtFloat.one) in
  let alpha := ExtValues.fmsubf one invb_d b_d in
  ExtValues.fmaddf invb_d alpha invb_d.

Lemma Rabs_relax:
  forall b b' (INEQ : (b < b')%R) x,
    (-b <= x <= b)%R -> (Rabs x < b')%R.
Proof.
  intros.
  apply Rabs_lt.
  lra.
Qed.

Definition approx_inv_thresh := (25/2251799813685248)%R.
(* 1.11022302462516e-14 *)
    
Theorem approx_inv_longu_correct :
  forall b,
    (0 < Int64.unsigned b)%Z ->
    exists (f : float),
      (approx_inv_longu (Vlong b)) = Vfloat f /\
      is_finite _ _ f = true /\ (Rabs((B2R _ _ f) - (1 / IZR (Int64.unsigned b))) <= approx_inv_thresh)%R.
Proof.
  intros b NONZ.
  unfold approx_inv_longu.
  cbn.
  econstructor.
  split.
  reflexivity.
  Local Transparent Float.neg Float.of_single Float32.of_longu Float32.div Float.of_longu Float32.of_int Float.of_int.
  unfold Float.fma, Float.neg, Float.of_single, Float32.of_longu, ExtFloat32.inv, Float32.div, Float.of_longu, ExtFloat32.one, Float32.of_int, ExtFloat.one, Float.of_int.
  set (re := (@eq_refl Datatypes.comparison Lt)).
  change (Int.signed (Int.repr 1)) with 1%Z.
  set (b' := Int64.unsigned b) in *.
  pose proof (Int64.unsigned_range b) as RANGE.
  change Int64.modulus with 18446744073709551616%Z in RANGE.                                              
  assert(1 <= IZR b' <= 18446744073709551616)%R as RANGE'.
  { split; apply IZR_le; lia.
  }

  assert (-16777216 <= 1 <= 16777216)%Z as SILLY by lia.
  destruct (BofZ_exact 24 128 re re 1 SILLY) as (C0R & C0F & _).
  clear SILLY.
  set (one_s := (BofZ 24 128 re re 1)) in *.
  
  pose proof (BofZ_correct 24 128 re re b') as C1.
  cbn in C1.
  rewrite Rlt_bool_true in C1; cycle 1.
  { clear C1.
    eapply (Rabs_relax (IZR 18446744073709551616)).
    lra.
    set (b'' := IZR b') in *.
    gappa.
  }
  destruct C1 as (C1R & C1F & _).
  set (b_s :=  (BofZ 24 128 re re b')) in *.

  assert(1 <= B2R 24 128 b_s <= 18446744073709551616)%R as b_s_RANGE.
  { rewrite C1R.
    gappa.
  }
  assert(B2R 24 128 b_s <> 0)%R as b_s_NONZ by lra.
  
  pose proof (Bdiv_correct 24 128 re re Float32.binop_nan mode_NE one_s b_s b_s_NONZ) as C2.
  rewrite Rlt_bool_true in C2; cycle 1.
  { clear C2.
    apply Rabs_relax with (b := 1%R).
    { cbn; lra. }
    rewrite C0R.
    set (r_b_s := B2R 24 128 b_s) in *.
    cbn.
    gappa.
  }
  
  destruct C2 as (C2R & C2F & _).
  set (invb_s := (Bdiv 24 128 re re Float32.binop_nan mode_NE one_s b_s)) in *.
  rewrite C0F in C2F.

  assert ((1/18446744073709551616 <= B2R 24 128 invb_s <= 1)%R) as invb_s_RANGE.
  { rewrite C2R.
    set (r_b_s := B2R 24 128 b_s) in *.
    rewrite C0R.
    cbn.
    gappa.
  }
  
  pose proof (Bconv_correct 24 128 53 1024 re re Float.of_single_nan mode_NE invb_s C2F) as C3.
  rewrite Rlt_bool_true in C3; cycle 1.
  { clear C3.
    set (r_invb_s := (B2R 24 128 invb_s)) in *.
    apply Rabs_relax with (b := 1%R).
    { replace 1%R with (bpow radix2 0)%R by reflexivity.
      apply bpow_lt.
      lia.
    }
    cbn.
    gappa.
  }
  
  destruct C3 as (C3R & C3F & _).
  set (invb_d :=  (Bconv 24 128 53 1024 re re Float.of_single_nan mode_NE invb_s)) in *.
  assert ((1/18446744073709551616 <= B2R 53 1024 invb_d <= 1)%R) as invb_d_RANGE.
  { 
    rewrite C3R.
    set (r_invb_s := B2R 24 128 invb_s) in *.
    cbn.
    gappa.
  }

  pose proof (is_finite_Bopp 53 1024 Float.neg_nan invb_d) as opp_finite.
  rewrite C3F in opp_finite.

  pose proof (BofZ_correct 53 1024 re re 1) as C4.
  rewrite Rlt_bool_true in C4; cycle 1.
  { clear C4.
    cbn.
    eapply (Rabs_relax (IZR 18446744073709551616)).
    lra.
    set (b'' := IZR b') in *.
    gappa.
  }
  destruct C4 as (C4R & C4F & _).
  
  pose proof (BofZ_correct 53 1024 re re b') as C5.
  cbn in C5.
  rewrite Rlt_bool_true in C5; cycle 1.
  { clear C5.
    eapply (Rabs_relax (IZR 18446744073709551616)).
    lra.
    set (b'' := IZR b') in *.
    gappa.
  }
  destruct C5 as (C5R & C5F & _).
  set (b_d :=  (BofZ 53 1024 re re b')) in *.
  
  assert(1 <= B2R 53 1024 b_d <= 18446744073709551616)%R as b_d_RANGE.
  { rewrite C5R.
    gappa.
  }

  pose proof (Bfma_correct 53 1024 re re Float.fma_nan mode_NE
          (Bopp 53 1024 Float.neg_nan invb_d) (BofZ 53 1024 re re b')
          (BofZ 53 1024 re re 1) opp_finite C5F C4F) as C6.
  rewrite Rlt_bool_true in C6; cycle 1.
  { clear C6.
    rewrite C4R.
    rewrite B2R_Bopp.
    cbn.
    eapply (Rabs_relax (IZR 18446744073709551616)).
    { lra. }
    fold invb_d.
    fold b_d.
    set (r_invb_d := B2R 53 1024 invb_d) in *.
    set (r_b_d := B2R 53 1024 b_d) in *.
    gappa.
  }
  fold b_d in C6.
  destruct C6 as (C6R & C6F & _).

  pose proof (Bfma_correct 53 1024 re re Float.fma_nan mode_NE
       (Bfma 53 1024 re re Float.fma_nan mode_NE
          (Bopp 53 1024 Float.neg_nan invb_d) b_d (BofZ 53 1024 re re 1))
       invb_d invb_d C6F C3F C3F) as C7.
  rewrite Rlt_bool_true in C7; cycle 1.
  { clear C7.
    rewrite C6R.
    rewrite B2R_Bopp.
    eapply (Rabs_relax (bpow radix2 64)).
    { apply bpow_lt. lia. }
    rewrite C4R.
    cbn.
    set (r_invb_d := B2R 53 1024 invb_d) in *.
    set (r_b_d := B2R 53 1024 b_d) in *.
    gappa.
  }
  destruct C7 as (C7R & C7F & _).

  split. assumption.
  rewrite C7R.
  rewrite C6R.
  rewrite C5R.
  rewrite C4R.
  rewrite B2R_Bopp.
  rewrite C3R.
  rewrite C2R.
  rewrite C1R.
  rewrite C0R.
  cbn.
  set(b1 := IZR b') in *.
  replace (round radix2 (FLT_exp (-1074) 53) ZnearestE 1) with 1%R by gappa.
  set (bd := round radix2 (FLT_exp (-1074) 53) ZnearestE b1).
  set (x0 := round radix2 (FLT_exp (-1074) 53) ZnearestE
         (round radix2 (FLT_exp (-149) 24) ZnearestE
                    (1 / round radix2 (FLT_exp (-149) 24) ZnearestE b1))).
  set (alpha0 := (- x0 * bd + 1)%R).
  set (y1 := (round radix2 (FLT_exp (-1074) 53) ZnearestE alpha0 * x0 + x0)%R).
  set (x1 := round radix2 (FLT_exp (-1074) 53) ZnearestE y1).
  replace (x1 - 1/b1)%R with ((y1-1/b1)+(x1-y1))%R by ring.

  assert(alpha0 = -((x0-1/bd)/(1/bd)))%R as alpha0_EQ.
  { unfold alpha0.
    field.
    unfold bd.
    gappa.
  }
  assert(y1-1/b1 = ((round radix2 (FLT_exp (-1074) 53) ZnearestE alpha0)
                   - alpha0) * x0
                   + alpha0*(x0-1/b1) - ((bd-b1)/b1) * x0)%R as y1_EQ.
  { unfold y1, alpha0.
    field.
    lra.
  }
  assert(Rabs alpha0 <= 257/2147483648)%R as alpha0_ABS.
  { rewrite alpha0_EQ.
    unfold x0, bd.
    gappa.
  }
  assert (Rabs (x0 - 1 / b1) <= 3/33554432)%R as x0_delta_ABS.
  { unfold x0.
    gappa.
  }
  set (x0_delta := (x0 - 1 / b1)%R) in *.
  assert (Rabs ((bd - b1) / b1) <= 1/9007199254740992)%R as bd_delta_ABS.
  { unfold bd.
    gappa.
  }
  set (bd_delta := ((bd - b1) / b1)%R) in *.
  set (rnd_alpha0_delta := (round radix2 (FLT_exp (-1074) 53) ZnearestE alpha0 - alpha0)%R) in *.
  assert (Rabs rnd_alpha0_delta <= 1/75557863725914323419136)%R as rnd_alpha0_delta_ABS.
  { unfold rnd_alpha0_delta.
    gappa.
  }
  assert (1/18446744073709551616 <= x0 <= 1)%R as x0_RANGE.
  { unfold x0.
    gappa.
  }
  assert (Rabs (y1 - 1 / b1) <= 49/4503599627370496)%R as y1_delta_ABS.
  { rewrite y1_EQ.
    gappa.
  }
  assert (Rabs(x1 - y1) <= 1/9007199254740992)%R as x1_delta_ABS.
  { unfold x1.
    gappa.
  }
  set (y1_delta := (y1 - 1 / b1)%R) in *.
  set (x1_delta := (x1-y1)%R) in *.
  unfold approx_inv_thresh.
  gappa.
Qed.

Definition rough_approx_inv_longu b :=
  let invb_s := ExtValues.invfs (Val.maketotal (Val.singleoflongu b)) in
  Val.floatofsingle invb_s.

Definition rough_approx_inv_thresh := (3/33554432)%R.
(* 8.94069671630859e-8 *)
    
Theorem rough_approx_inv_longu_correct :
  forall b,
    (0 < Int64.unsigned b)%Z ->
    exists (f : float),
      (rough_approx_inv_longu (Vlong b)) = Vfloat f /\
        is_finite _ _ f = true /\
        (Rabs((B2R _ _ f) - (1 / IZR (Int64.unsigned b))) <= rough_approx_inv_thresh)%R /\
        Bsign _ _ f = false.
Proof.
  intros b NONZ.
  unfold rough_approx_inv_longu.
  cbn.
  econstructor.
  split.
  reflexivity.
  Local Transparent Float.neg Float.of_single Float32.of_longu Float32.div Float.of_longu Float32.of_int Float.of_int.
  unfold Float.fma, Float.neg, Float.of_single, Float32.of_longu, ExtFloat32.inv, Float32.div, Float.of_longu, ExtFloat32.one, Float32.of_int, ExtFloat.one, Float.of_int.
  set (re := (@eq_refl Datatypes.comparison Lt)).
  change (Int.signed (Int.repr 1)) with 1%Z.
  set (b' := Int64.unsigned b) in *.
  pose proof (Int64.unsigned_range b) as RANGE.
  change Int64.modulus with 18446744073709551616%Z in RANGE.                                              
  assert(1 <= IZR b' <= 18446744073709551616)%R as RANGE'.
  { split; apply IZR_le; lia.
  }

  assert (-16777216 <= 1 <= 16777216)%Z as SILLY by lia.
  destruct (BofZ_exact 24 128 re re 1 SILLY) as (C0R & C0F & _).
  clear SILLY.
  set (one_s := (BofZ 24 128 re re 1)) in *.
  
  pose proof (BofZ_correct 24 128 re re b') as C1.
  cbn in C1.
  rewrite Rlt_bool_true in C1; cycle 1.
  { clear C1.
    eapply (Rabs_relax (IZR 18446744073709551616)).
    lra.
    set (b'' := IZR b') in *.
    gappa.
  }
  rewrite (Zlt_bool_false b' 0) in C1 by lia.
  destruct C1 as (C1R & C1F & C1S).
  set (b_s :=  (BofZ 24 128 re re b')) in *.

  assert(1 <= B2R 24 128 b_s <= 18446744073709551616)%R as b_s_RANGE.
  { rewrite C1R.
    gappa.
  }
  assert(B2R 24 128 b_s <> 0)%R as b_s_NONZ by lra.
  
  pose proof (Bdiv_correct 24 128 re re Float32.binop_nan mode_NE one_s b_s b_s_NONZ) as C2.
  rewrite Rlt_bool_true in C2; cycle 1.
  { clear C2.
    apply Rabs_relax with (b := 1%R).
    { cbn; lra. }
    rewrite C0R.
    set (r_b_s := B2R 24 128 b_s) in *.
    cbn.
    gappa.
  }
  destruct C2 as (C2R & C2F & C2Sz).
  rewrite C1S in C2Sz.
  change (xorb _ _) with false in C2Sz.
  set (invb_s := (Bdiv 24 128 re re Float32.binop_nan mode_NE one_s b_s)) in *.
  rewrite C0F in C2F.
  assert (is_nan 24 128 invb_s = false) as NAN.
  { apply is_finite_not_is_nan.
    assumption.
  }
  pose proof (C2Sz NAN) as C2S.
  clear C2Sz.
  
  assert ((1/18446744073709551616 <= B2R 24 128 invb_s <= 1)%R) as invb_s_RANGE.
  { rewrite C2R.
    set (r_b_s := B2R 24 128 b_s) in *.
    rewrite C0R.
    cbn.
    gappa.
  }
  
  pose proof (Bconv_correct 24 128 53 1024 re re Float.of_single_nan mode_NE invb_s C2F) as C3.
  rewrite Rlt_bool_true in C3; cycle 1.
  { clear C3.
    set (r_invb_s := (B2R 24 128 invb_s)) in *.
    apply Rabs_relax with (b := 1%R).
    { replace 1%R with (bpow radix2 0)%R by reflexivity.
      apply bpow_lt.
      lia.
    }
    cbn.
    gappa.
  }
  
  destruct C3 as (C3R & C3F & C3S).
  set (invb_d :=  (Bconv 24 128 53 1024 re re Float.of_single_nan mode_NE invb_s)) in *.
  assert ((1/18446744073709551616 <= B2R 53 1024 invb_d <= 1)%R) as invb_d_RANGE.
  { 
    rewrite C3R.
    set (r_invb_s := B2R 24 128 invb_s) in *.
    cbn.
    gappa.
  }
  rewrite C2S in C3S.

  split. assumption. split. 2: assumption.
  unfold rough_approx_inv_thresh.
  rewrite C3R.
  rewrite C2R.
  rewrite C1R.
  rewrite C0R.
  cbn.
  gappa.
Qed.

Definition rough_approx_div_longu a b :=
  Val.maketotal (Val.longuoffloat_ne
    (Val.mulf (Val.maketotal (Val.floatoflongu a)) (rough_approx_inv_longu b))).

Definition rough_approx_div_thresh := 1649267441663%Z.

Lemma Bsign_false_nonneg:
  forall prec emax f,
    (Bsign prec emax f) = false -> (0 <= (B2R prec emax f))%R.
Proof.
  intros until f. intro SIGN.
  destruct f.
  1, 2, 3: cbn; lra.
  cbn.
  apply F2R_ge_0.
  cbn.
  cbn in SIGN.
  rewrite SIGN.
  cbn.
  lia.
Qed.

Theorem rough_approx_div_longu_correct :
  forall a b,
    (0 < Int64.unsigned b)%Z ->
    exists (q : int64),
      (rough_approx_div_longu (Vlong a) (Vlong b)) = Vlong q /\
        (Z.abs ((Int64.unsigned a) - (Int64.unsigned q) * (Int64.unsigned b)) <= rough_approx_div_thresh)%Z.
Proof.
  intros a b b_NONZ.
  unfold rough_approx_div_longu.
  destruct (rough_approx_inv_longu_correct b b_NONZ) as (f & C1R & C1F & C1D & C1S).
  rewrite C1R.
  cbn.

  Local Transparent Float.to_longu Float.mul.
  unfold Float.to_longu_ne, Float.mul, Float.of_longu.
  set (re := (@eq_refl Datatypes.comparison Lt)).

  pose proof (Int64.unsigned_range a) as a_RANGE.
  change Int64.modulus with 18446744073709551616%Z in a_RANGE.
  set (a' := Int64.unsigned a) in *.
  assert (0 <= IZR a' <= 18446744073709551615)%R as IZR_a_RANGE.
  { split ; apply IZR_le; lia.
  }

  set (a_d := (round radix2 (FLT_exp (-1074) 53) ZnearestE (IZR a'))).
  assert (0 <= a_d <= 18446744073709551616)%R as a_d_RANGE.
  { unfold a_d.
    gappa.
  }

  pose proof (Int64.unsigned_range b) as b_RANGE.
  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
  set (b' := Int64.unsigned b) in *.
  assert (1 <= IZR b' <= 18446744073709551615)%R as IZR_b_RANGE.
  { split ; apply IZR_le; lia.
  }

  set (f_r := (B2R _ _ f)) in *.
  assert (0 <= f_r <= 1)%R as f_RANGE.
  { split.
    { apply Bsign_false_nonneg. trivial .}
    unfold rough_approx_inv_thresh in C1D.
    replace f_r with ((f_r - 1 / IZR b') + 1/ IZR b')%R by (field; lra).
    set (delta := (f_r - 1 / IZR b')%R) in *.
    destruct (Z_le_gt_dec b' 1) as [LE | GT].
    { admit. }
    assert (2 <= IZR b')%R as NOT1.
    { apply IZR_le. lia. }
    gappa.
  }
  pose proof (BofZ_correct 53 1024 re re a') as C2.
  rewrite Rlt_bool_true in C2; cycle 1.
  { clear C2.
    apply Rabs_relax with (b := bpow radix2 64).
    { apply bpow_lt. lia. }
    cbn.
    gappa.
  }
  destruct C2 as (C2R & C2F & _).

  pose proof (Bmult_correct  53 1024 re re Float.binop_nan mode_NE
                             (BofZ 53 1024 re re a') f) as C3.
  rewrite Rlt_bool_true in C3; cycle 1.
  { clear C3.
    rewrite C2R.
    replace (B2R 53 1024 f)%R with ((B2R 53 1024 f - 1 / IZR b') + 1 / IZR b')%R by ring.
    unfold rough_approx_inv_thresh in C1D.
    apply Rabs_relax with (b := bpow radix2 65).
    { apply bpow_lt ; lia. }
    cbn.
    gappa.
  }
  rewrite C1F in C3.
  rewrite C2F in C3.
  cbn in C3.
  destruct C3 as (C3R & C3F & _).

  pose proof (ZofB_ne_range_correct 53 1024
               (Bmult 53 1024 re re Float.binop_nan mode_NE
                      (BofZ 53 1024 re re a') f) 0 Int64.max_unsigned) as C4.
  rewrite C3F in C4.
  rewrite C3R in C4.
  rewrite C2R in C4.

  assert(0 <=
      (round radix2 (FLT_exp (-1074) 53) ZnearestE
              (round radix2 (FLT_exp (3 - 1024 - 53) 53) 
                 (round_mode mode_NE) (IZR a') * B2R 53 1024 f))
         <= IZR Int64.max_unsigned)%R as q_RANGE.
  { clear C4.
    cbn.
    change (IZR Int64.max_unsigned) with 18446744073709551615%R.
    fold a_d.
    gappa.
    set (y := ((IZR a') * B2R 53 1024 f)%R).
    assert ((round radix2 (FLT_exp (-1074) 53) ZnearestE (IZR a') * B2R 53 1024 f) >= 0)%R.
   
    cbn.
    replace (B2R 53 1024 f)%R with ((B2R 53 1024 f - 1 / IZR b') + 1 / IZR b')%R by ring.
    unfold rough_approx_inv_thresh in C1D.
    set (delta :=  (B2R 53 1024 f - 1 / IZR b')%R) in *.
    gappa.
  }
Admitted.