dichotomy

dicho.ml

(* x + a = b

   xl + au >= bl *)

let au = 4.0
let bl = float_of_string "0x0.0000000000001p0"

let neg_zero = Int64.float_of_bits 0x8000_0000_0000_0000L
let pos_zero = Int64.float_of_bits 0x0000_0000_0000_0000L

let fsucc f = Int64.(float_of_bits @@ succ @@ bits_of_float f)
let fpred f = Int64.(float_of_bits @@ pred @@ bits_of_float f)

let rm = 0x0000_0000_0000_0001L
let lm = 0x4000_0000_0000_0000L

let on_bit f pos =
  let mask = Int64.(shift_right lm pos) in
  Int64.(float_of_bits @@ logor (bits_of_float f) mask)

let off_bit f pos =
  let mask = Int64.(shift_right rm pos) in
  Int64.(float_of_bits @@ logxor (bits_of_float f) mask)

let ppf f = Printf.printf "%Lx" (Int64.bits_of_float f)

let off_current_neg f a b = f +. a <= b
let off_current_pos f a b = f +. a >= b

(* no closure allocation *)
let rec dichotomy cmp f a target i =
  if i >= 63 then f else
    let nf = on_bit f i in
    if cmp nf a target then
      dichotomy cmp f a target (i + 1)
    else
      dichotomy cmp nf a target (i + 1)

let find_neg a target = dichotomy off_current_neg neg_zero a target 0
let find_pos a target = dichotomy off_current_pos pos_zero a target 0

let ppd d = Printf.printf "%LX\n" d