Double dichotomy to determine set of x such that x + a == b, given value of a and b

doubledich.ml

(* double dichotomy to determine floating point range for x
   in ``x + a == b``, given value of ``a`` and ``b`` *)
type t =
  | Range of float * float
  | Single of float
  | Empty

let to_string = function
  | Range (l, u) -> Printf.sprintf "{%.16e ... %.16e}" l u
  | Single f -> Printf.sprintf "{%.16e}" f
  | Empty -> "empty range"

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 m1 = 0x4000_0000_0000_0000L
let m2 = 0xBFFF_FFFF_FFFF_FFFFL

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

let off_bit f pos =
  let mask = Int64.(shift_right m2 pos) in
  Int64.(float_of_bits @@ logand (bits_of_float f) mask)

let dichotomy restore toggle init a b =
  let rec approach f i =
    if i >= 63 then f else
      let tf = toggle f i in
      if restore tf a b then approach f (i + 1) else approach tf (i + 1)
  in approach init 0

let upper_neg = dichotomy (fun f a b -> f +. a <= b) on_bit neg_zero

let lower_neg = dichotomy (fun f a b -> f +. a < b) on_bit neg_zero

let upper_pos = dichotomy (fun f a b -> f +. a > b) on_bit pos_zero

let lower_pos = dichotomy (fun f a b -> f +. a >= b) on_bit pos_zero

let pos_critical_point = 9.9792015476736e+291
  (* first value that can be added with a pos normalish
     to overflow to infinity *)
let neg_critical_point = -9.9792015476736e+291

let dump a b =
  Printf.printf "upper_neg : %.16e\n" (upper_neg a b);
  Printf.printf "lower_neg : %.16e\n" (lower_neg a b);
  Printf.printf "upper_pos : %.16e\n" (upper_pos a b);
  Printf.printf "lower_pos : %.16e\n" (lower_pos a b)

let range a b =
  try begin
  let l, u =
    if a = b then
      lower_neg a b, upper_pos a b else
    if a < b then begin
      if a +. max_float < b then failwith "empty range" else
      if b = infinity then
        fsucc (lower_pos a infinity), max_float
      else
        fsucc (lower_pos a b), upper_pos a b
    end
    else begin
      if a -. max_float > b then failwith "empty range" else
      if b = neg_infinity then
        (-.max_float), fsucc (upper_neg a neg_infinity)
      else
        lower_neg a b, fsucc (upper_neg a b)
    end in
  if l = u && l <> 0. then Single l else
  if l = u then Range (l, u) else
  if l > u then begin
    if l +. a = b then Single l else
    if u +. a = b then Single u else
      Empty
  end else
    Range (l, u)
  end
  with _ -> Empty

let max3 = max_float /. 3.

let check_invariant a b =
  match range a b with
  | Range (l, u) ->
    let b1 = l +. a = b in
    let b2 = u +. a = b in
    let b3 = (if l > 0. then (fpred l) else (fsucc l)) +. a <> b in
    let b4 =
      if u <> max_float then
        (if u > 0. then (fsucc u) else (fpred u)) +. a <> b
      else
        true in
    if not (b1 && b2 && b3 && b4) then failwith "Err1" else ()
  | Single s ->
    if (not ((fsucc s) +. a <> b) && ((fpred s) +. a <> b)) then
      failwith "Err2" else ()
  | Empty -> ()

let test a b = range a b |> to_string

let normal_test () =
  check_invariant 1.0 1.0;
  check_invariant 1.0 1.1;
  check_invariant 1.0 1.2;
  check_invariant 1.0 1.4;
  check_invariant 1.0 1.5;
  check_invariant 0.0 1.3;
  check_invariant 0.0 1.4;
  check_invariant 0.0 1.5;
  check_invariant 0.0 1.6;
  check_invariant 0.0 1.7;
  check_invariant max_float infinity;
  check_invariant (max_float /. 3.) infinity;
  check_invariant (max_float /. 2.) infinity;
  check_invariant (max_float /. 2.) max_float;
  check_invariant (max_float /. 2.) (max_float /. 3.)

let () = Random.self_init ()

let random_pos () =
  match Random.int 20 with
  | 0 -> min_float
  | 1 -> max_float
  | 2 -> Random.float min_float
  | 3 | 4 -> Random.float max_float
  | 5 -> Random.float 2e-308
  | 6 -> 2e-308
  | 7 | 8 | 9 -> Random.float 100.
  | 11 | 12 -> 0.0
  | 13 -> max_float /. 2.
  | 14 -> max_float /. 3.
  | _ -> Random.float 1_000_00.

let random_pair () =
  let a = random_pos () in
  let b = if Random.int 10 < 2 then infinity else random_pos () in
  if Random.int 5 < 1 then a, a else a, b

let random_test () =
  for i = 0 to 1_000_000_000 do
    if i mod 1_000_00 = 0 then
      print_endline @@ Printf.sprintf "Passed %d tests\n" i;
    let f1, f2 = random_pair () in
    try
      check_invariant f1 f2
    with _ ->
      Printf.printf "f1: %.16e\nf2: %.16e\n" f1 f2;
      dump f1 f2; assert false
  done

let () = normal_test ()
let () = random_test ()