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.

Definition approx_inv_thresh := (1/70368744177664)%R.

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.
    
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.