A portable OCaml function to print floats

For printing the results of automated Frama-C tests, we need printing functions with the following properties:

1. they should print all the information available (doing otherwise might mask an existing bug);
2. they should produce results readable enough for a person to validate the result of each test once;
3. and their output should not differ depending on the execution platform.

For tests that manipulates floating-point data, we need a floating-point printing function with the above properties. If you start with the OCaml function Format.printf and print enough digits to make condition 1 true, then a number of issues arise with condition 3. For a while, we used hexadecimal in all tests outputs. For instance, the value for d in the program d = 3.14; was printed as d ∈ 0x1.91eb851eb851fp1. In this notation, the digits after 1. are in base 16, and p1 means \multiplied by 2^1" by similarity with e1 meaning "multiplied by 10^1". Users were not fond of this notation. I don't know why. It is accepted for input by C99 compilers and can be printed with the %a format code in printf(). It is easy to see at a glance that the d is between 3 and 4 closer to 3 because 0x1.8p1 is 3 and 0x2.0p1 is 4.

In order to improve condition 2 I wrote my own floating-point pretty-printing function. The goal here is not to have something very sophisticated but to make sure we do not miss a bug in a test output for the silly reason of having had ten fingers for a long time. For the same program d = 3.14; the new decimal pretty-printing function produces d ∈ 3.1400000000000001 which is what could be expected considering that the number 3.14 cannot be represented exactly as a double. More important is that if the program is modified to read d = 0x1.91eb851eb8520p1; changing only the least significant bit in the floating-point representation of d the output becomes d ∈ 3.1400000000000005 showing that there is a difference.

I have recently improved my improvement so I thought I might as well provide the function here (but you may have seen a previous version through other channels). The function I wrote only uses 64-bit ints and double-precision floats so that it would be straightforward to translate to most other programming languages.

It works by first isolating the sign exponent and mantissa of the floating-point number to print. The mantissa ranges over 0 .. 0xFFFFFFFFFFFFF and we would like a sequence of decimal digits that ranges over 0 .. 9999999999999998 instead. For this purpose the mantissa is bitwise or-ed into the representation of the double-precision float 1.0 giving a result that ranges over 1.0 .. 2.0 (excluded). The basic idea is then to subtract 1.0 to multiply by 1e16 to convert to an integer and to print with width 16 and with leading zeroes (the %016Ld format).

The function is at the bottom of this post. The recent improvement is this test:

      let d  re = 
	if d >= 1.5 
	then d -. 1.5  5000000000000000L 
	else d -. 1.0  0L 
      in 
      ... 

Without these lines the function was a bit unsatisfactory when the fractional part was above 0.5. The printed decimal number was a strictly increasing function of the actual floating-point number which was enough to ensure condition 1 but some numbers that were exactly representable in both the floating-point double-precision format and as decimal numbers were printed with a decimal representation other than that representing exactly the number. It was bugging me although I doubt it ever bugged someone else. This is fixed by the above code which regains a vital bit of mantissa when the fractional part of the number to print is greater than 0.5. Even with this fix the function is not correctly rounded but it doesn't claim to be. It is only the closest we have bothered to come to satisfying conditions 1-3 simultaneously.

let double_norm = Int64.shift_left 1L 52 
let double_mask = Int64.pred double_norm 
let pretty_normal ~use_hex fmt f = 
  let i = Int64.bits_of_float f in 
  let s = 0L <> (Int64.logand Int64.min_int i) in 
  let i = Int64.logand Int64.max_int i in 
  let exp = Int64.to_int (Int64.shift_right_logical i 52) in 
  let man = Int64.logand i double_mask in 
  let s = if s then "-" else "" in 
  let firstdigit  exp = 
    if exp <> 0 
    then 1  (exp - 1023) 
    else 0  -1022 
  in 
  if not use_hex 
  then begin 
      let firstdigit  man  exp = 
	if 0 <= exp && exp <= 12 
	then begin 
	    Int64.to_int 
	      (Int64.shift_right_logical 
		  (Int64.logor man double_norm) 
		  (52 - exp))  
            Int64.logand (Int64.shift_left man exp) double_mask  
            0 
	  end 
	else firstdigit  man  exp 
      in 
      let d = 
	Int64.float_of_bits 
	  (Int64.logor 0x3ff0000000000000L man) 
      in 
      let d  re = 
	if d >= 1.5 
	then d -. 1.5  5000000000000000L 
	else d -. 1.0  0L 
      in 
      let d = d *. 1e16 in 
      let decdigits = Int64.add re (Int64.of_float d) in 
      if exp = 0 
      then 
	Format.fprintf fmt "%s%d.%016Ld" 
	  s 
	  firstdigit 
	  decdigits 
      else 
	Format.fprintf fmt "%s%d.%016Ld*2^%d" 
	  s 
	  firstdigit 
	  decdigits 
	  exp 
    end 
  else 
    Format.fprintf fmt "%s0x%d.%013Lxp%d" 
      s 
      firstdigit 
      man 
      exp