open Core

let _ = assert (Sys.word_size_in_bits = 64)
let chlist s = String.to_list s
let ( <> ) l1 l2 = not (List.equal Char.( = ) l1 l2)
let str = String.of_char_list
let prec = 53
let emax = 1024
let exp_bits = 11
let saturated_exp = Int.shift_left 1 exp_bits - 1

let rec int_bits (x : int) : bool list =
  if Int.(x < 0) then raise (Invalid_argument "Number must be nonnegative")
  else if Int.(x = 0) then []
  else
    let q = x /% 2 and r = x % 2 in
    int_bits q @ [ r = 1 ]

let int_to_positive (x : int) : Extraction.Internal.BinNums.positive =
  if x <= 0 then raise (Invalid_argument "Number must be positive")
  else
    let bits = List.tl_exn (int_bits x) in
    List.fold_left
      ~f:(fun acc digit ->
        if digit then Extraction.Internal.BinNums.Coq_xI acc
        else Extraction.Internal.BinNums.Coq_xO acc)
      ~init:Extraction.Internal.BinNums.Coq_xH bits

let int_to_z (x : int) : Extraction.Internal.BinNums.coq_Z =
  if x = 0 then Z0
  else if x < 0 then Zneg (int_to_positive (-x))
  else Zpos (int_to_positive x)

let rec int63_of_positive (x : Extraction.Internal.BinNums.positive) : Int63.t =
  let two = Int63.(succ one) in
  match x with
  | Coq_xH -> Int63.of_int_exn 1
  | Coq_xO x -> Int63.(two * int63_of_positive x)
  | Coq_xI x -> Int63.((two * int63_of_positive x) + one)

let int63_of_z (x : Extraction.Internal.BinNums.coq_Z) : Int63.t =
  match x with
  | Z0 -> Int63.zero
  | Zpos x -> int63_of_positive x
  | Zneg x -> Int63.neg (int63_of_positive x)

let int63_to_positive x = Int63.to_int_exn x |> int_to_positive

(* Conversions are the same as those in Coq's FloatOps.
   See https://github.com/coq/coq/blob/master/theories/Floats/FloatOps.v *)

let normfr_mantissa f =
  let f = Float.abs f in
  if Float.(f >= 0.5) && Float.(f < 1.) then Float.to_int (Float.ldexp f prec)
  else 0

let float_to_flocq (x : float) : Extraction.Internal.Floats.float =
  match Float.classify x with
  | Zero -> Extraction.Internal.Floats.B754_zero (Float.ieee_negative x)
  | Nan ->
      Extraction.Internal.Floats.B754_nan
        (Float.ieee_negative x, Float.ieee_mantissa x |> int63_to_positive)
  | Infinite -> Extraction.Internal.Floats.B754_infinity (Float.ieee_negative x)
  | Normal | Subnormal -> (
      let prec_z = int_to_z prec and emax_z = int_to_z emax in
      let r, exp = Float.frexp x in
      let e = int_to_z (exp - prec) and r' = int_to_z (normfr_mantissa r) in
      let shr, e' =
        Extraction.Internal.Floats.(shr_fexp prec_z emax_z r' e Coq_loc_Exact)
      in
      match shr.shr_m with
      | Zpos p -> B754_finite (Float.is_negative x, p, e')
      | Zneg _ | Z0 -> assert false)

let float_from_flocq x : float =
  match x with
  | Extraction.Internal.Floats.B754_zero s ->
      Float.create_ieee_exn ~negative:s ~mantissa:Int63.zero ~exponent:0
  | Extraction.Internal.Floats.B754_infinity s ->
      if s then Float.neg_infinity else Float.infinity
  | Extraction.Internal.Floats.B754_nan (s, m) ->
      let m_int = int63_of_positive m in
      Float.create_ieee_exn ~negative:s ~mantissa:m_int ~exponent:saturated_exp
  | Extraction.Internal.Floats.B754_finite (s, m, e) ->
      let pm = Float.of_int63 (int63_of_positive m) in
      let f = Float.ldexp pm (Int63.to_int_exn (int63_of_z e)) in
      if s then Float.neg f else f

open struct
  open Base_quickcheck

  let%expect_test "float conversion" =
    Test.run_exn
      (module Float)
      ~f:(fun x -> [%test_eq: float] (float_from_flocq (float_to_flocq x)) x)
end