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

Local Notation "'rd'" := (round radix2 (FLT_exp (-1074) 53) ZnearestE).
Local Notation "'rs'" := (round radix2 (FLT_exp (-149) 24) ZnearestE).

Definition approx_inv_rel_thresh := (1049/72057594037927936)%R.
Theorem approx_inv_longu_correct_rel :
  forall b,
    (0 < Int64.unsigned b)%Z ->
    exists (f : float),
      (approx_inv_longu (Vlong b)) = Vfloat f /\
      is_finite _ _ f = true /\ (Rabs(IZR (Int64.unsigned b) * (B2R _ _ f) - 1) <= approx_inv_rel_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 (rd 1) with 1%R by gappa.
  replace (rd (rs (1 / rs b1))) with
    ((((rd (rs (1 / rs b1)) - (/b1))/(/b1))+1)*(/ b1))%R ; cycle 1.
  { field. lra. }
  set (er0 := ((rd (rs (1 / rs b1)) - (/b1))/(/b1))%R).
  replace (rd b1) with ((((rd b1) - b1)/b1 + 1) * b1)%R; cycle 1.
  { field. lra. }
  set (er1 := (((rd b1) - b1)/b1)%R).
  replace (- ((er0 + 1) * / b1) * ((er1 + 1) * b1) + 1)%R
    with (1 - (er0 + 1)*(er1 + 1))%R ; cycle 1.
  { field. lra. }
  set (z0 := (1 - (er0 + 1) * (er1 + 1))%R).
  assert (Rabs er0 <= 257/2147483648)%R as er0_ABS.
  { unfold er0.
    gappa.
  }
  assert (Rabs er1 <= 1/9007199254740992)%R as er1_ABS.
  { unfold er1.
    gappa.
  }
  replace (rd z0) with ((rd(z0)-z0)+z0)%R by ring.
  set (ea0 := (rd(z0)-z0)%R).
  assert (Rabs ea0 <= 1/75557863725914323419136)%R as ea0_ABS.
  { unfold ea0. unfold z0.
    gappa.
  }
  set (z1 :=  ((ea0 + z0) * ((er0 + 1) * / b1) + (er0 + 1) * / b1)%R).
  replace (rd z1) with ((((rd z1)-z1)/z1+1)*z1)%R; cycle 1.
  { field.
    unfold z1.
    unfold z0.
    gappa.
  }
  set (er2 := ((rd z1 - z1) / z1)%R).
  assert (Rabs er2 <= 1/9007199254740992)%R as er2_ABS.
  { unfold er2.
    unfold z1, z0.
    gappa.
  }
  unfold z1, z0.
  replace  (b1 *
      ((er2 + 1) *
       ((ea0 + (1 - (er0 + 1) * (er1 + 1))) * ((er0 + 1) * / b1) +
          (er0 + 1) * / b1)) - 1)%R
      with (-er0*er0*er1*er2 - er0*er0*er1 + ea0*er0*er2 - er0*er0*er2 - 2*er0*er1*er2 + ea0*er0 - er0*er0 - 2*er0*er1 + ea0*er2 - er1*er2 + ea0 - er1 + er2)%R; cycle 1.
  { field. lra. }
  unfold approx_inv_rel_thresh.
  gappa.
Qed.

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

Definition step1_real_inv_thresh := (3/33554432)%R.
(* 8.94069671630859e-8 *)

Theorem step1_real_inv_longu_correct :
  forall b,
    (0 < Int64.unsigned b)%Z ->
    exists (f : float),
      (step1_real_inv_longu (Vlong b)) = Vfloat f /\
        (B2R _ _ f) = (rd (rs (1 / rs (IZR (Int64.unsigned b))))) /\
        is_finite _ _ f = true /\
        Bsign _ _ f = false.
Proof.
  intros b NONZ.
  unfold step1_real_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.
  }
  rewrite C1R in C2.
  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.
  rewrite C2R in C3R.
  rewrite C0R in C3R.

  auto.
Qed.

Theorem step1_real_inv_longu_correct1 :
  forall b,
    (Int64.unsigned b = 1%Z) ->
    exists f,
    (step1_real_inv_longu (Vlong b)) = Vfloat f /\
        (B2R _ _ f) = 1%R /\
        is_finite _ _ f = true /\
        Bsign _ _ f = false.
Proof.
  intros b EQ1.
  assert (0 < Int64.unsigned b)%Z as b_RANGE by lia.
  destruct (step1_real_inv_longu_correct b b_RANGE) as (f & C1E & C1R & C1F & C1S).
  rewrite EQ1 in C1R.
  exists f.
  repeat split; try assumption.
  rewrite C1R.
  gappa.
Qed.

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.

Lemma Znearest_IZR_le :
  forall rnd n x, (IZR n <= x)%R -> (n <= Znearest rnd x)%Z.
Proof.
  intros until x. intro ORDER.
  pose proof (Znearest_ge_floor rnd x).
  pose proof (Zfloor_le _ _ ORDER) as KK.
  rewrite Zfloor_IZR in KK.
  lia.
Qed.

Lemma Znearest_le_IZR :
  forall rnd n x, (x <= IZR n)%R -> (Znearest rnd x <= n)%Z.
Proof.
  intros until x. intro ORDER.
  pose proof (Znearest_le_ceil rnd x).
  pose proof (Zceil_le _ _ ORDER) as KK.
  rewrite Zceil_IZR in KK.
  lia.
Qed.

Definition step1_real_div_longu a b :=
  Val.mulf (Val.maketotal (Val.floatoflongu a)) (step1_real_inv_longu b).

Definition step1_div_longu a b :=
  Val.maketotal (Val.longuoffloat_ne (step1_real_div_longu a b)).

Definition step1_real_quotient (a b : R) :=
             rd ((rd (a)) * (rd (rs (1 / rs (b))))).
  
Theorem step1_real_div_longu_correct:
  forall a b,
    (1 < Int64.unsigned b)%Z ->
    exists (q : float),
      (step1_real_div_longu (Vlong a) (Vlong b)) = Vfloat q /\
        (B2R _ _ q) = step1_real_quotient (IZR (Int64.unsigned a))
                                          (IZR (Int64.unsigned b)) /\
        is_finite _ _ q = true /\
        Bsign _ _ q = false.
Proof.
  intros a b b_NON01.
  assert (0 < Int64.unsigned b)%Z as b_NON0 by lia.
  destruct (step1_real_inv_longu_correct b b_NON0) as (f & C1E & C1R & C1F & C1S).
  unfold step1_real_div_longu.
  rewrite C1E.
  cbn.
  set (b' := Int64.unsigned b) in *.
  Local Transparent Float.mul.
  unfold Float.mul, Float.of_longu.
  econstructor.
  split. reflexivity.
  set (a' := Int64.unsigned a) in *.
  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.
  assert (0 <= IZR a' <= 18446744073709551615)%R as IZR_a_RANGE.
  { split; apply IZR_le; lia. }
  pose proof (Int64.unsigned_range b) as b_RANGE.
  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
  assert (2 <= IZR b' <= 18446744073709551615)%R as IZR_b_RANGE.
  { split; apply IZR_le; lia. }
  
  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 & C2S).
  rewrite Zlt_bool_false in C2S by lia.

  pose proof (Bmult_correct 53 1024 re re Float.binop_nan mode_NE (BofZ 53 1024 re re a') f) as C3.
  rewrite C1S in C3.
  rewrite C2S in C3.
  rewrite C1F in C3.
  rewrite C2F in C3.
  rewrite C1R in C3.
  rewrite C2R in C3.
  rewrite Rlt_bool_true in C3; cycle 1.
  { apply Rabs_relax with (b := bpow radix2 64).
    { apply bpow_lt ; lia. }
    cbn.
    gappa.
  }
  cbn in C3.
  destruct C3 as (C3R & C3F & C3Sz).
  assert (is_nan 53 1024
          (Bmult 53 1024 re re Float.binop_nan mode_NE 
                 (BofZ 53 1024 re re a') f) = false) as NAN.
  { apply is_finite_not_is_nan.
    assumption. }
  pose proof (C3Sz NAN) as C3S.
  clear NAN C3Sz.

  auto.
Qed.

Definition smallb_thresh :=       4398046511104%Z.

Definition smallb_approx_real_range := 2200000000000%R.
Lemma step1_smallb_real :
  forall a b
    (a_RANGE : (1 <= a <= 18446744073709551615)%R)
    (b_RANGE : (1 <= b <= IZR smallb_thresh)%R),
    (Rabs((step1_real_quotient a b) * b - a) <= smallb_approx_real_range)%R.
Proof.
  intros.
  unfold smallb_thresh in b_RANGE.
  unfold smallb_approx_real_range.
  unfold step1_real_quotient.
  set (q := ((rd (a)) * (rd (rs (1 / rs (b)))))%R) in *.
  replace ((rd q) *b - a)%R with
     ((rd(q)-q)/q * rd(a) * (1 + (rd (rs (1 / rs (b))) - 1/b)/(1/b)) +
  (rd (a)) * ((rd (rs (1 / rs (b))) - 1 / b) / (1/b)) +
        (rd(a) - a))%R; cycle 1.
  { unfold q.
    field.
    split. lra.
    split. gappa.
    gappa.
  }
  unfold q.
  gappa.
Qed.

Lemma step1_div_longu_a0 :
  forall b,
    (0 < Int64.unsigned b)%Z ->
    (step1_div_longu (Vlong Int64.zero) (Vlong b)) = Vlong Int64.zero.
Proof.
  intros b b_NOT0.
  unfold step1_div_longu.
  unfold step1_real_div_longu.
  destruct (step1_real_inv_longu_correct b b_NOT0) as
    (f & C1E & C1R & C1F & C1S).
  rewrite C1E.
  cbn.
  unfold Float.to_longu_ne, Float.of_longu, Float.mul.
  rewrite Int64.unsigned_zero.
  set (re :=  (@eq_refl Datatypes.comparison Lt)).
  assert (- 2 ^ 53 <= 0 <= 2 ^ 53)%Z as SILLY by lia.
  destruct (BofZ_exact 53 1024 re re 0 SILLY) as (C2R & C2F & C2S).
  
  pose proof (Bmult_correct 53 1024 re re Float.binop_nan mode_NE
                            (BofZ 53 1024 re re 0) f) as C3.
  rewrite C1F in C3.
  rewrite C2F in C3.
  rewrite C1S in C3.
  rewrite C2S in C3.
  rewrite Z.ltb_irrefl in C3.
  rewrite Rlt_bool_true in C3; cycle 1.
  { clear C3.
    apply Rabs_relax with (b := bpow radix2 64).
    { apply bpow_lt. lia. }
    cbn.
    rewrite Rmult_0_l.
    gappa.
  }
  rewrite C2R in C3.
  rewrite Rmult_0_l in C3.
  destruct C3 as (C3R & C3F & C3Sz).
  change (true && true) with true in C3F.
  change (xorb false false) with false in C3Sz.
  assert (is_nan 53 1024
           (Bmult 53 1024 re re Float.binop_nan mode_NE 
                  (BofZ 53 1024 re re 0) f) = false) as NAN.
  { apply is_finite_not_is_nan.
    assumption.
  }
  pose proof (C3Sz NAN) as C3S.
  clear NAN C3Sz.
  pose proof ((ZofB_ne_range_correct 53 1024
             (Bmult 53 1024 re re Float.binop_nan mode_NE
                    (BofZ 53 1024 re re 0) f) 0 Int64.max_unsigned)) as C4.
  rewrite C3R in C4.
  replace (round radix2 (FLT_exp (3 - 1024 - 53) 53) (round_mode mode_NE) 0)
    with 0%R in C4 by (cbn ; gappa).
  rewrite Znearest_IZR in C4.
  cbn zeta in C4.
  rewrite Z.leb_refl in C4.
  change (0 <=? Int64.max_unsigned)%Z with true in C4.
  rewrite andb_true_r in C4.
  rewrite andb_true_r in C4.
  rewrite C3F in C4.
  rewrite C4.
  reflexivity.
Qed.
                                         
Definition smallb_approx_range := 4400000000000%Z.
Lemma step1_div_longu_correct :
    forall a b,
    (1 < Int64.unsigned b <= smallb_thresh)%Z ->
    exists (q : int64),
      (step1_div_longu (Vlong a) (Vlong b)) = Vlong q /\
        (Z.abs (Int64.unsigned a - Int64.unsigned b*Int64.unsigned q) <= smallb_approx_range)%Z.
Proof.
  intros a b b_RANGE.

  pose proof (Int64.unsigned_range a) as a_RANGE.
  change Int64.modulus with 18446744073709551616%Z in a_RANGE.
  set (a' := Int64.unsigned a) in *.
  set (b' := Int64.unsigned b) in *.

  destruct (Z_le_gt_dec a' 0).
  { assert (a' = 0%Z) as ZERO by lia.
    exists Int64.zero.
    rewrite ZERO.
    rewrite Int64.unsigned_zero.
    replace (Z.abs (0 - b' * 0))%Z with 0%Z by lia.
    replace a with Int64.zero; cycle 1.
    {
      unfold a' in ZERO.
      unfold Int64.zero.
      rewrite <- ZERO.
      apply Int64.repr_unsigned.
    }
    split.
    { apply step1_div_longu_a0.
      lia.
    }
    unfold smallb_approx_range.
    lia.
  }

  unfold step1_div_longu.
  assert (1 < b')%Z as b_NOT01 by lia.   
  destruct (step1_real_div_longu_correct a b b_NOT01) as (q & C1E & C1R & C1F & C1S).
  rewrite C1E. cbn.
  unfold Float.to_longu_ne.
  pose proof (ZofB_ne_range_correct 53 1024 q 0 Int64.max_unsigned) as C2.
  rewrite C1F in C2.

  
  assert (1 <= IZR a' <= 18446744073709551615)%R as a_RANGE'.
  { split; apply IZR_le; lia. }
  assert (2 <= IZR b' <= IZR smallb_thresh)%R as b_RANGE'.
  { split; apply IZR_le; lia. }
  assert (1 <= IZR b' <= IZR smallb_thresh)%R as b_RANGE'' by lra.
  pose proof (step1_smallb_real (IZR a') (IZR b') a_RANGE' b_RANGE'') as DELTA.
  fold a' in C1R.
  fold b' in C1R.
  rewrite <- C1R in DELTA.

  assert (0 <= B2R _ _ q)%R as q_NONNEG.
  { apply Bsign_false_nonneg. assumption. }
  cbn in C2.
  rewrite Zle_bool_true in C2; cycle 1.
  { apply Znearest_IZR_le. assumption. }
  assert (B2R _ _ q <= 9223376000000000000)%R as q_SMALL.
  { replace (B2R _ _ q) with
      ((IZR a' / IZR b') + (B2R _ _ q * IZR b' - IZR a') / IZR b')%R; cycle 1.
    { field. lra. }
    unfold smallb_approx_real_range in DELTA.
    unfold smallb_thresh in b_RANGE'.
    set (y := (B2R 53 1024 q * IZR b' - IZR a')%R) in *.
    gappa.
  }
  rewrite Zle_bool_true in C2; cycle 1.
  { apply Znearest_le_IZR. lra. }
  cbn in C2.

  change Int64.max_unsigned with 18446744073709551615%Z.
  rewrite C2.
  cbn.

  econstructor. split. reflexivity.
  rewrite Int64.unsigned_repr; cycle 1.
  { split.
    - apply Znearest_IZR_le. lra.
    - apply Znearest_le_IZR.
      change Int64.max_unsigned with 18446744073709551615%Z.
      lra.
  }
  apply le_IZR.
  rewrite abs_IZR.
  unfold smallb_approx_real_range, smallb_approx_range, smallb_thresh in *.
  rewrite minus_IZR.
  rewrite mult_IZR.
  set (q_r := B2R 53 1024 q) in *.
  assert (Rabs (IZR (ZnearestE q_r) - q_r) <= / 2)%R as NEAR
      by apply Znearest_imp2.
  set (q_i := IZR (ZnearestE q_r)) in *.
  replace  (IZR a' - IZR b' * q_i)%R with
    (-(IZR b' * (q_i - q_r)) - (q_r * IZR b' - IZR a'))%R by ring.
  set (delta1 := (q_i - q_r)%R) in *.
  set (delta2 := (q_r * IZR b' - IZR a')%R) in *.
  gappa.
Qed.

(*
Lemma range_up_le :
  forall a x b,
    (IZR a <= IZR x <= (IZR b) - 1)%R ->
    (a <= x < b)%Z.
Proof.
  intros until b. intro RANGE.
  split.
  { apply le_IZR. lra. }
  assert (x <= b-1)%Z.
  { apply le_IZR. rewrite minus_IZR. lra. }
  lia.
Qed.
 *)

(*
Lemma find_quotient:
  forall (a b : Z)
         (b_POS : (b > 0)%Z)
         (invb : R)
         (eps : R)
         (invb_EPS : (Rabs (invb * IZR b - 1) <= eps)%R)
         (SMALL : (IZR (Z.abs a)*eps < IZR b / 2)%R),
    (a / b)%Z =
      let q := ZnearestE ((IZR a) * invb) in
      if (b*q >? a)%Z
      then (q-1)%Z
      else q.
Proof.
  intros.
  rewrite <- Rabs_Zabs in SMALL.
  set (q := ZnearestE (IZR a * invb)%R).
  cbn zeta.
  set (b' := IZR b) in *.
  set (a' := IZR a) in *.
  assert (1 <= b')%R as b_POS'.
  { apply IZR_le.
    lia.
  }
  
  pose proof (Znearest_imp2 (fun x : Z => negb (Z.even x)) (a' * invb)) as NEAR.
  fold q in NEAR.
  set (q' := IZR q) in *.
  assert ((Rabs a') * Rabs (invb * b' - 1) <= (Rabs a') * eps)%R as S1.
  { apply Rmult_le_compat_l.
    apply Rabs_pos.
    assumption. }
  rewrite <- Rabs_mult in S1.
  assert ((Rabs b') * Rabs (q' - a' * invb) <= (Rabs b') / 2)%R as S2.
  { apply Rmult_le_compat_l.
    apply Rabs_pos.
    assumption. }
  rewrite <- Rabs_mult in S2.
  rewrite (Rabs_right b') in S2 by lra.
  pose proof (Rabs_triang (a' * (invb * b' - 1))
                          (b' * (q' - a' * invb))) as TRIANGLE.
  replace (a' * (invb * b' - 1) + b' * (q' - a' * invb))%R with
    (b' * q' - a')%R in TRIANGLE by ring.
  assert (Rabs (b' * q' - a') < b')%R as DELTA by lra.
  
  pose proof (Zgt_cases (b * q) a)%Z as CASE.
  destruct (_ >? _)%Z.
  { unfold b', q', a' in DELTA.
    rewrite <- mult_IZR in DELTA.
    rewrite <- minus_IZR in DELTA.
    rewrite Rabs_Zabs in DELTA.
    apply lt_IZR in DELTA.
    rewrite Z.abs_eq in DELTA by lia.
    
    apply Zdiv_unique with (b := (a - (q-1)*b)%Z).
    ring.
    split; lia.
  }
  
  rewrite <- Rabs_Ropp in DELTA.
  unfold b', q', a' in DELTA.
  rewrite <- mult_IZR in DELTA.
  rewrite <- minus_IZR in DELTA.
  rewrite <- opp_IZR in DELTA.
  rewrite Rabs_Zabs in DELTA.
  apply lt_IZR in DELTA.
  rewrite Z.abs_eq in DELTA by lia.

  apply Zdiv_unique with (b := (a - q*b)%Z).
  ring.
  
  split; lia.
Qed.
 *)

Lemma find_quotient:
  forall (a b : Z)
         (b_POS : (0 < b)%Z)
         (qr : R)
         (GAP : (Rabs (IZR a / IZR b - qr) < /2)%R),
    (a / b)%Z =
      let q := ZnearestE qr in
      if (b*q >? a)%Z
      then (q-1)%Z
      else q.
Proof.
  intros.
  set (q := ZnearestE qr).
  cbn zeta.
  set (b' := IZR b) in *.
  set (a' := IZR a) in *.
  assert (1 <= b')%R as b_POS'.
  { apply IZR_le.
    lia.
  }
  
  pose proof (Znearest_imp2 (fun x : Z => negb (Z.even x)) qr) as ROUND.
  fold q in ROUND.
  set (q' := IZR q) in *.
  
  pose proof (Rabs_triang (a' / b' - qr)
                          (qr - q'))%R as TRIANGLE.
  replace ((a' / b' - qr) +  (qr - q'))%R with
    (a' / b' - q')%R in TRIANGLE by ring.
  rewrite <- Rabs_Ropp in ROUND.
  replace (- (q' - qr))%R with (qr - q')%R in ROUND by ring.
  assert (Z.abs (a - b*q) < b)%Z as DELTA.
  { apply lt_IZR.
    rewrite <- Rabs_Zabs.
    rewrite minus_IZR.
    rewrite mult_IZR.
    fold a' q' b'.
    apply Rmult_lt_reg_r with (r := (/b')%R).
    { apply Rinv_0_lt_compat. lra. }
    rewrite Rinv_r by lra.
    replace (/ b')%R with (/ Rabs(b'))%R ; cycle 1.
    { f_equal.
      apply Rabs_pos_eq. lra. }
    rewrite <- Rabs_Rinv by lra.
    rewrite <- Rabs_mult.
    replace ((a' - b' * q') * / b')%R with (a'/b' - q')%R by (field ; lra).
    lra.
  }
   
  pose proof (Zgt_cases (b * q) a)%Z as CASE.
  destruct (_ >? _)%Z.
  { apply Zdiv_unique with (b := (a - (q-1)*b)%Z).
    ring.
    split; lia.
  }

  apply Zdiv_unique with (b := (a - q*b)%Z).
  ring.
  split; lia.
Qed.

Definition step2_real_div_long a b :=
  Val.mulf (Val.maketotal (Val.floatoflong a)) (approx_inv_longu b).

Definition smalla_thresh := 34184372088832%Z.

Lemma step2_real_div_long_smalla_correct :
    forall (a b : int64)
           (a_SMALL : (Z.abs (Int64.signed a) <= smalla_thresh)%Z)
           (b_NOT0 : (0 < Int64.unsigned b)%Z),
    exists (q : float),
      (step2_real_div_long (Vlong a) (Vlong b)) = Vfloat q /\
        (Rabs ((B2R _ _ q) - (IZR (Int64.signed a)) / (IZR (Int64.unsigned b))) <= (32767/65536))%R /\
      is_finite _ _ q = true.
Proof.
  intros.
  unfold step2_real_div_long.
  destruct (approx_inv_longu_correct_rel b b_NOT0) as (f & C0E & C0F & C0R).
  rewrite C0E.
  econstructor.
  split. reflexivity.
  Local Transparent Float.of_long.
  unfold Float.mul, Float.of_long.
  set (re := (@eq_refl Datatypes.comparison Lt)) in *.
  pose proof (Int64.unsigned_range b) as b_RANGE.
  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
  set (a' := Int64.signed a) in *.
  set (b' := Int64.unsigned b) in *.
  assert (1 <= IZR b' <= 18446744073709551615)%R as b_RANGE'.
  { split; apply IZR_le; lia.
  }
  assert(Rabs (IZR a') <= IZR smalla_thresh)%R as a_RANGE'.
  { rewrite Rabs_Zabs.
    apply IZR_le.
    assumption.
  }
  assert (- 2 ^ 53 <= a' <= 2 ^ 53)%Z as SILLY.
  { unfold smalla_thresh in a_SMALL.
    apply Z.abs_le.
    lia.
  } 
  destruct (BofZ_exact 53 1024 re re (Int64.signed a) SILLY) as (C1R & C1F & C1S).
  fold a' in C1R, C1F, C1S.
  pose proof (Bmult_correct 53 1024 re re Float.binop_nan mode_NE (BofZ 53 1024 re re a') f) as C2.
  rewrite Rlt_bool_true in C2 ; cycle 1.
  { clear C2.
    apply Rabs_relax with (b := bpow radix2 53).
    { apply bpow_lt. lia. }
    cbn.
    rewrite C1R.
    unfold approx_inv_rel_thresh in C0R.
    replace (B2R 53 1024 f) with
      ((1/IZR b') * ((IZR b' * B2R 53 1024 f - 1) + 1))%R ; cycle 1.
    { field. lra. }
    unfold smalla_thresh in *.
    gappa.
  }
  rewrite C0F in C2.
  rewrite C1R in C2.
  rewrite C1F in C2.
  rewrite C1S in C2.
  cbn in C2.
  destruct C2 as (C2R & C2F & _).
  split.
  2: exact C2F.
  rewrite C2R.
  replace (IZR a' * (B2R 53 1024 f))%R with
    ((IZR a'/IZR b') * ((IZR b' * B2R 53 1024 f - 1) + 1))%R ; cycle 1.
  { field. lra. }
  set (delta1 := (IZR b' * B2R 53 1024 f - 1)%R) in *.
  set (q1 := (IZR a' / IZR b' * (delta1 + 1))%R).
  replace (rd q1) with (((rd q1) - q1) + q1)%R by ring.
  set (delta2 := ((rd q1) - q1)%R).
  unfold q1.
  replace (delta2 + IZR a' / IZR b' * (delta1 + 1) - IZR a' / IZR b')%R with
    (delta2 + (IZR a' / IZR b') * delta1)%R by ring.
  unfold delta2.
  unfold q1.
  unfold approx_inv_rel_thresh in *.
  unfold smalla_thresh in *.
  gappa.
Qed.

Definition step2_div_long' a b :=
  Val.maketotal (Val.longoffloat_ne (step2_real_div_long a b)).

Definition step2_div_long a b :=
  let q := step2_div_long' a b in
  Val.select (Val.cmpl_bool Cgt (Val.subl (Val.mull q b) a) (Vlong Int64.zero))
             (Val.addl q (Vlong Int64.mone)) q Tlong.

Lemma Znearest_lub :
  forall choice (n : Z) (x : R), (IZR n <= x)%R -> (n <= Znearest choice x)%Z.
Proof.
  intros until x. intro BND.
  pose proof (Zfloor_lub n x BND).
  pose proof (Znearest_ge_floor choice x).
  lia.
Qed.

Lemma Znearest_glb :
  forall choice (n : Z) (x : R), (x <= IZR n)%R -> (Znearest choice x <= n)%Z.
Proof.
  intros until x. intro BND.
  pose proof (Zceil_glb n x BND).
  pose proof (Znearest_le_ceil choice x).
  lia.
Qed.

Lemma function_ite :
  forall {A B : Type} (fn : A->B) (b : bool) (x y: A),
    fn (if b then x else y) = (if b then fn x else fn y).
Proof.
  intros.
  destruct b; reflexivity.
Qed.

Lemma normalize_ite :
  forall ty (b : bool) x y,
    Val.normalize (if b then x else y) ty =
      (if b then Val.normalize x ty else Val.normalize y ty).
Proof.
  intros.
  destruct b; reflexivity.
Qed.


Lemma int64_mul_signed_unsigned:
    forall x y : int64,
      Int64.mul x y = Int64.repr (Int64.signed x * Int64.unsigned y).
Proof.
  intros.
  unfold Int64.mul.
  apply Int64.eqm_samerepr.
  apply Int64.eqm_mult.
  - apply Int64.eqm_sym.
    apply Int64.eqm_signed_unsigned.
  - apply Int64.eqm_refl.
Qed.

Lemma int64_eqm_signed_repr:
  forall z : Z, Int64.eqm z (Int64.signed (Int64.repr z)).
Proof.
  intros.
  apply Int64.eqm_trans with (y := Int64.unsigned (Int64.repr z)).
  - apply Int64.eqm_unsigned_repr.
  - apply Int64.eqm_sym.
    apply Int64.eqm_signed_unsigned.
Qed.

Lemma signed_repr_sub:
  forall x y,
    Int64.repr (Int64.signed (Int64.repr x) - y) =
    Int64.repr (x - y).
Proof.
  intros.
  apply Int64.eqm_samerepr.
  apply Int64.eqm_sub.
  - apply Int64.eqm_sym.
    apply int64_eqm_signed_repr.
  - apply Int64.eqm_refl.
Qed.
                                                                        
Lemma step2_div_long_smalla_correct :
    forall a b
      (a_SMALL : (Z.abs (Int64.signed a) <= smalla_thresh)%Z)
      (b_NOT0 : (0 < Int64.unsigned b)%Z)
      (b_NOT_VERY_BIG : (Int64.unsigned b <= Int64.max_signed)%Z),
      step2_div_long (Vlong a) (Vlong b) = Vlong (Int64.repr (Int64.signed a / Int64.unsigned b))%Z.
Proof.
  intros.
  pose proof (Int64.unsigned_range b) as b_RANGE.
  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
  set (a' := (Int64.signed a)) in *.
  set (b' := (Int64.unsigned b)) in *.
  assert (Rabs (IZR a') <= IZR smalla_thresh)%R as a_RANGE'.
  { rewrite Rabs_Zabs.
    apply IZR_le.
    assumption.
  }
  assert (1 <= IZR b' <= 18446744073709551615)%R as b_RANGE'.
  { split; apply IZR_le; lia.
  } 
  destruct (step2_real_div_long_smalla_correct a b a_SMALL b_NOT0) as (q & C1R & C1E & C1F).
  fold a' b' in C1E.
  assert ((Int64.min_signed <=? ZnearestE (B2R 53 1024 q))=true)%Z as q_LOW.
  { apply Zle_imp_le_bool.
    change Int64.min_signed with (-9223372036854775808)%Z.
    apply Znearest_lub.
    set (q' :=  B2R 53 1024 q) in *.
    replace q' with (IZR a' / IZR b' + (q' - IZR a' / IZR b'))%R by ring.
    unfold smalla_thresh in a_RANGE'.
    gappa.
  }
  assert ((ZnearestE (B2R 53 1024 q) <=? Int64.max_signed)=true)%Z as q_HIGH.
  { apply Zle_imp_le_bool.
    change Int64.max_signed with (9223372036854775807)%Z.
    apply Znearest_glb.
    set (q' :=  B2R 53 1024 q) in *.
    replace q' with (IZR a' / IZR b' + (q' - IZR a' / IZR b'))%R by ring.
    unfold smalla_thresh in a_RANGE'.
    gappa.
  }
  unfold step2_div_long, step2_div_long'.
  rewrite C1R.
  cbn.
  unfold Float.to_long_ne.
  rewrite (ZofB_ne_range_correct _ _ q Int64.min_signed Int64.max_signed).
  rewrite C1F.
  rewrite q_LOW.
  rewrite q_HIGH.
  cbn.
  rewrite normalize_ite.
  cbn.
  rewrite <- (function_ite Vlong).
  f_equal.
  unfold Int64.lt.
  set (q' :=  B2R 53 1024 q) in *.
  fold a'.
  assert (Int64.signed (Int64.repr (ZnearestE q')) = ZnearestE q') as q_SIGNED.
  { apply Int64.signed_repr.
    split; lia.
  }
  rewrite Int64.add_signed.
  rewrite q_SIGNED.
  rewrite Int64.signed_mone.
  rewrite Int64.signed_zero.
  rewrite <- (function_ite Int64.repr).
  f_equal.
  replace  (ZnearestE q' + -1)%Z with (ZnearestE q' - 1)%Z by ring.
  
  set (q'' :=  (ZnearestE q')) in *.
  Check Int64.sub_signed.
  fold a'.
  rewrite int64_mul_signed_unsigned.
  rewrite q_SIGNED.
  fold b'.

  rewrite Int64.sub_signed.
  fold a'.
  rewrite signed_repr_sub.

  assert ((Rabs (IZR a' / IZR b' - q') < / 2)%R) as HALF.
  { replace (IZR a' / IZR b' - q')%R with
      (-(q' - IZR a' / IZR b'))%R by ring.
    rewrite Rabs_Ropp.
    lra.
  }
  pose proof (find_quotient a' b' b_NOT0 q' HALF) as QUOTIENT.
  fold q'' in QUOTIENT.
  cbn zeta in QUOTIENT.

  assert (b' <> 0)%Z as NONZ by lia.
  pose proof (Zmod_eq_full a' b' NONZ) as MOD.
  assert (b' > 0)%Z as b_GT0 by lia.
  pose proof (Z_mod_lt a' b' b_GT0) as MOD_LT.
  destruct (Z_lt_dec a' (b' * q'')) as [LT | GE].
  { replace (b' * q'' >? a')%Z with true in QUOTIENT by lia.
    replace q'' with (1 + (a' / b'))%Z by lia.
    replace ((1 + a' / b') * b' - a')%Z
      with (-(a' - a' / b' * b')+b')%Z by ring.
    rewrite <- MOD.
    rewrite Int64.signed_repr; cycle 1.
    { change Int64.min_signed with (-9223372036854775808)%Z in *.
      change Int64.max_signed with (9223372036854775807)%Z in *.
      lia.
    }
    rewrite zlt_true by lia.
    ring.
  }
  replace (b' * q'' >? a')%Z with false in QUOTIENT by lia.
  rewrite <- QUOTIENT.
  replace (a' / b' * b' - a')%Z with (-(a' - a' / b' * b'))%Z by ring.
  rewrite <- MOD.
  rewrite Int64.signed_repr ; cycle 1.
  { change Int64.min_signed with (-9223372036854775808)%Z in *.
    change Int64.max_signed with (9223372036854775807)%Z in *.
    lia.
  }
  rewrite zlt_false by lia.
  reflexivity.
Qed.

Definition twostep_div_longu a b :=
  let q1 := step1_div_longu a b in
  let q2 := step2_div_long (Val.subl a (Val.mull b q1)) b in
  Val.addl q1 q2.

Lemma unsigned_repr_sub :
  forall a b,
    Int64.repr (a - b) = Int64.repr (a - Int64.unsigned (Int64.repr b)).
Proof.
  intros.
  apply Int64.eqm_samerepr.
  apply Int64.eqm_sub.
  - apply Int64.eqm_refl.
  - apply Int64.eqm_unsigned_repr.
Qed.

Lemma unsigned_repr_add :
  forall a b,
    Int64.repr (a + b) = Int64.repr (a + Int64.unsigned (Int64.repr b)).
Proof.
  intros.
  apply Int64.eqm_samerepr.
  apply Int64.eqm_add.
  - apply Int64.eqm_refl.
  - apply Int64.eqm_unsigned_repr.
Qed.
    
Lemma twostep_div_longu_smallb_correct :
    forall a b
      (b_RANGE : (1 < Int64.unsigned b <= smallb_thresh)%Z),
          (twostep_div_longu (Vlong a) (Vlong b)) =
            (Val.maketotal (Val.divlu (Vlong a) (Vlong b))).
Proof.
  intros.
  unfold twostep_div_longu.
  destruct (step1_div_longu_correct a b b_RANGE) as (q1 & C1R & C1E).
  rewrite C1R.
  set (q1' := Int64.unsigned q1) in *.
  set (b' := Int64.unsigned b) in *.
  set (a' := Int64.unsigned a) in *.
  assert ( Z.abs (Int64.signed (Int64.sub a (Int64.mul b q1))) <= smalla_thresh)%Z as r1_SMALL.
  { unfold smalla_thresh, smallb_approx_range in *.
    unfold Int64.sub, Int64.mul.
    fold q1' b' a'.
    rewrite <- unsigned_repr_sub.
    rewrite Int64.signed_repr ; cycle 1.
    { change Int64.min_signed with (-9223372036854775808)%Z.
      change Int64.max_signed with (9223372036854775807)%Z.
      lia.
    }
    lia.
  }
  assert (0 < b')%Z as b_NOT0 by lia.
  assert (b' <= Int64.max_signed)%Z as b_NOTBIG.
  { change Int64.max_signed with (9223372036854775807)%Z.
    unfold smallb_thresh in b_RANGE.
    lia.
  }
  cbn.
  rewrite (step2_div_long_smalla_correct (Int64.sub a (Int64.mul b q1)) b r1_SMALL b_NOT0 b_NOTBIG).
  unfold Int64.add, Int64.sub, Int64.mul, Int64.divu.
  fold q1' b' a'.
  rewrite <- unsigned_repr_sub.
  rewrite <- unsigned_repr_add.
  rewrite Int64.signed_repr ; cycle 1.
  {
    change Int64.min_signed with (-9223372036854775808)%Z.
    change Int64.max_signed with (9223372036854775807)%Z.
    unfold smallb_approx_range in *.
    lia.
  }
  rewrite Z.add_comm.
  rewrite <- Z.div_add by lia.
  replace (a' - b' * q1' + q1' * b')%Z with a' by ring.
  rewrite Int64.eq_false ; cycle 1.
  { intro Z. unfold b' in b_NOT0. rewrite Z in b_NOT0.
    rewrite Int64.unsigned_zero in b_NOT0.
    lia.
  }
  reflexivity.
Qed.


Lemma step2_real_div_long_bigb_correct :
    forall (a b : int64)
           (b_BIG : ((Int64.unsigned b) > smallb_thresh)%Z),
    exists (q : float),
      (step2_real_div_long (Vlong a) (Vlong b)) = Vfloat q /\
        (Rabs ((B2R _ _ q) - (IZR (Int64.signed a)) / (IZR (Int64.unsigned b))) <= (32767/65536))%R /\
      is_finite _ _ q = true.
Proof.
  intros.
  unfold step2_real_div_long.
  assert (0 < Int64.unsigned b)%Z as b_NOT0 by (unfold smallb_thresh in *; lia).
  destruct (approx_inv_longu_correct_rel b b_NOT0) as (f & C0E & C0F & C0R).
  rewrite C0E.
  econstructor.
  split. reflexivity.
  Local Transparent Float.of_long.
  unfold Float.mul, Float.of_long.
  set (re := (@eq_refl Datatypes.comparison Lt)) in *.
  pose proof (Int64.unsigned_range b) as b_RANGE.
  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
  pose proof (Int64.signed_range a) as a_RANGE.
  set (a' := Int64.signed a) in *.
  set (b' := Int64.unsigned b) in *.
  assert (IZR (1 + smallb_thresh) <= IZR b' <= 18446744073709551615)%R as b_RANGE'.
  { split; apply IZR_le; lia.
  }
  assert(IZR Int64.min_signed <= IZR a' <= IZR Int64.max_signed)%R as a_RANGE'.
  { split; apply IZR_le; lia.
  }
  change Int64.min_signed with (-9223372036854775808)%Z in a_RANGE'.
  change Int64.max_signed with (9223372036854775807)%Z in a_RANGE'.
  pose proof (BofZ_correct 53 1024 re re a') as C1.
  rewrite Rlt_bool_true in C1 ; cycle 1.
  { clear C1.
    apply Rabs_relax with (b := bpow radix2 64).
    { apply bpow_lt; lia. }
    cbn.
    gappa.
  }
  cbn in C1.
  destruct C1 as (C1R & C1F & C1S).

  unfold smallb_thresh in b_RANGE'; cbn in b_RANGE'.

  pose proof (Bmult_correct 53 1024 re re Float.binop_nan mode_NE (BofZ 53 1024 re re a') f) as C2.
  rewrite Rlt_bool_true in C2 ; cycle 1.
  { clear C2.
    apply Rabs_relax with (b := bpow radix2 53).
    { apply bpow_lt. lia. }
    cbn.
    rewrite C1R.
    unfold approx_inv_rel_thresh in C0R.
    replace (B2R 53 1024 f) with
      ((1/IZR b') * ((IZR b' * B2R 53 1024 f - 1) + 1))%R ; cycle 1.
    { field.  lra. }
    gappa.
  }
  rewrite C0F in C2.
  rewrite C1R in C2.
  rewrite C1F in C2.
  rewrite C1S in C2.
  cbn in C2.
  destruct C2 as (C2R & C2F & _).
  split.
  2: exact C2F.
  rewrite C2R.
  set (f' := (B2R 53 1024 f)) in *.  
  replace (rd(rd (IZR a') * f') - IZR a' / IZR b')%R with
    ((rd(rd (IZR a') * f') - IZR a' * f') + IZR a' / IZR b' * (IZR b' * f' - 1))%R ; cycle 1.
  { field. lra. }
  unfold approx_inv_rel_thresh in *.
  gappa.
Qed.
                                                                        
Lemma step2_div_long_bigb_correct :
    forall a b
           (b_BIG : ((Int64.unsigned b) > smallb_thresh)%Z)
           (b_NOT_TOO_BIG : ((Int64.unsigned b) <= Int64.max_signed)%Z),
      step2_div_long (Vlong a) (Vlong b) = Vlong (Int64.repr (Int64.signed a / Int64.unsigned b))%Z.
Proof.
  intros.
  pose proof (Int64.unsigned_range b) as b_RANGE.
  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
  pose proof (Int64.signed_range a) as a_RANGE.
  set (a' := (Int64.signed a)) in *.
  set (b' := (Int64.unsigned b)) in *.
  assert (IZR (1 + smallb_thresh) <= IZR b' <= 18446744073709551615)%R as b_RANGE'.
  { split; apply IZR_le; lia.
  }
  assert(IZR Int64.min_signed <= IZR a' <= IZR Int64.max_signed)%R as a_RANGE'.
  { split; apply IZR_le; lia.
  }
  unfold smallb_thresh in *; cbn in b_RANGE'.
  change Int64.min_signed with (-9223372036854775808)%Z in *.
  change Int64.max_signed with (9223372036854775807)%Z in *.
  assert (0 < b')%Z as b_NOT0 by lia.

  destruct (step2_real_div_long_bigb_correct a b b_BIG) as (q & C1R & C1E & C1F).
  fold a' b' in C1E.
  assert ((Int64.min_signed <=? ZnearestE (B2R 53 1024 q))=true)%Z as q_LOW.
  { apply Zle_imp_le_bool.
    change Int64.min_signed with (-9223372036854775808)%Z.
    apply Znearest_lub.
    set (q' :=  B2R 53 1024 q) in *.
    replace q' with (IZR a' / IZR b' + (q' - IZR a' / IZR b'))%R by ring.
    gappa.
  }
  assert ((ZnearestE (B2R 53 1024 q) <=? Int64.max_signed)=true)%Z as q_HIGH.
  { apply Zle_imp_le_bool.
    change Int64.max_signed with (9223372036854775807)%Z.
    apply Znearest_glb.
    set (q' :=  B2R 53 1024 q) in *.
    replace q' with (IZR a' / IZR b' + (q' - IZR a' / IZR b'))%R by ring.
    gappa.
  }
  unfold step2_div_long, step2_div_long'.
  rewrite C1R.
  cbn.
  unfold Float.to_long_ne.
  rewrite (ZofB_ne_range_correct _ _ q Int64.min_signed Int64.max_signed).
  rewrite C1F.
  rewrite q_LOW.
  rewrite q_HIGH.
  cbn.
  rewrite normalize_ite.
  cbn.
  rewrite <- (function_ite Vlong).
  f_equal.
  unfold Int64.lt.
  set (q' :=  B2R 53 1024 q) in *.
  fold a'.
  assert (Int64.signed (Int64.repr (ZnearestE q')) = ZnearestE q') as q_SIGNED.
  { apply Int64.signed_repr.
    split; lia.
  }
  rewrite Int64.add_signed.
  rewrite q_SIGNED.
  rewrite Int64.signed_mone.
  rewrite Int64.signed_zero.
  rewrite <- (function_ite Int64.repr).
  f_equal.
  replace  (ZnearestE q' + -1)%Z with (ZnearestE q' - 1)%Z by ring.
  
  set (q'' :=  (ZnearestE q')) in *.
  Check Int64.sub_signed.
  fold a'.
  rewrite int64_mul_signed_unsigned.
  rewrite q_SIGNED.
  fold b'.

  rewrite Int64.sub_signed.
  fold a'.
  rewrite signed_repr_sub.

  assert ((Rabs (IZR a' / IZR b' - q') < / 2)%R) as HALF.
  { replace (IZR a' / IZR b' - q')%R with
      (-(q' - IZR a' / IZR b'))%R by ring.
    rewrite Rabs_Ropp.
    lra.
  }
  pose proof (find_quotient a' b' b_NOT0 q' HALF) as QUOTIENT.
  fold q'' in QUOTIENT.
  cbn zeta in QUOTIENT.

  assert (b' <> 0)%Z as NONZ by lia.
  pose proof (Zmod_eq_full a' b' NONZ) as MOD.
  assert (b' > 0)%Z as b_GT0 by lia.
  pose proof (Z_mod_lt a' b' b_GT0) as MOD_LT.
  destruct (Z_lt_dec a' (b' * q'')) as [LT | GE].
  { replace (b' * q'' >? a')%Z with true in QUOTIENT by lia.
    replace q'' with (1 + (a' / b'))%Z by lia.
    replace ((1 + a' / b') * b' - a')%Z
      with (-(a' - a' / b' * b')+b')%Z by ring.
    rewrite <- MOD.
    rewrite Int64.signed_repr; cycle 1.
    { change Int64.min_signed with (-9223372036854775808)%Z in *.
      change Int64.max_signed with (9223372036854775807)%Z in *.
      lia.
    }
    rewrite zlt_true by lia.
    ring.
  }
  replace (b' * q'' >? a')%Z with false in QUOTIENT by lia.
  rewrite <- QUOTIENT.
  replace (a' / b' * b' - a')%Z with (-(a' - a' / b' * b'))%Z by ring.
  rewrite <- MOD.
  rewrite Int64.signed_repr ; cycle 1.
  { change Int64.min_signed with (-9223372036854775808)%Z in *.
    change Int64.max_signed with (9223372036854775807)%Z in *.
    lia.
  }
  rewrite zlt_false by lia.
  reflexivity.
Qed.