## A portable OCaml function to print floats

By pascal on Saturday, October 29 2011, 14:47 - Permalink

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

- they should print all the information available (doing otherwise might mask an existing bug);
- they should produce results readable enough for a person to validate the result of each test once;
- 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