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 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 *)

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

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.