A portable OCaml function to print floats

Pascal Cuoq - 29th Oct 2011

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 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))
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
Pascal Cuoq
29th Oct 2011