Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
---
layout: post
author: Pascal Cuoq
date: 2011-10-29 12:47 +0200
categories: OCaml floating-point
format: xhtml
title: "A portable OCaml function to print floats"
summary:
---
{% raw %}
<p>For printing the results of automated Frama-C tests, we need printing functions with the following properties:</p>
<ol>
<li>they should print all the information available (doing otherwise might mask an existing bug);</li>
<li>they should produce results readable enough for a person to validate the result of each test once;</li>
<li>and their output should not differ depending on the execution platform.</li>
</ol>
<p>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 <code>Format.printf</code> 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 <code>d</code> in the program <code>d = 3.14;</code> was printed as <code>d ∈ 0x1.91eb851eb851fp1</code>. In this notation, the digits after <code>1.</code> are in base 16, and <code>p1</code> means \multiplied by <code>2^1</code>" by similarity with <code>e1</code> meaning "multiplied by <code>10^1</code>". 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 <code>%a</code> format code in <code>printf()</code>. It is easy to see at a glance that the <code>d</code> is between <code>3</code> and <code>4</code> closer to <code>3</code> because <code>0x1.8p1</code> is <code>3</code> and <code>0x2.0p1</code> is <code>4</code>.</p>
<p>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 <code>d = 3.14;</code> the new decimal pretty-printing function produces <code>d ∈ 3.1400000000000001</code> which is what could be expected considering that the number <code>3.14</code> cannot be represented exactly as a <code>double</code>. More important is that if the program is modified to read <code>d = 0x1.91eb851eb8520p1;</code> changing only the least significant bit in the floating-point representation of <code>d</code> the output becomes <code>d ∈ 3.1400000000000005</code> showing that there is a difference.</p>
<p>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.</p>
<p>It works by first isolating the sign exponent and mantissa of the floating-point number to print. The mantissa ranges over <code>0 .. 0xFFFFFFFFFFFFF</code> and we would like a sequence of decimal digits that ranges over <code>0 .. 9999999999999998</code> instead. For this purpose the mantissa is bitwise or-ed into the representation of the double-precision float <code>1.0</code> giving a result that ranges over <code>1.0 .. 2.0</code> (excluded). The basic idea is then to subtract <code>1.0</code> to multiply by <code>1e16</code> to convert to an integer and to print with width 16 and with leading zeroes (the <code>%016Ld</code> format).</p>
<p>The function is at the bottom of this post. The recent improvement is this test:</p>
<pre> let d re =
if d >= 1.5
then d -. 1.5 5000000000000000L
else d -. 1.0 0L
in
...
</pre>
<p>Without these lines the function was a bit unsatisfactory when the fractional part was above <code>0.5</code>. 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.</p>
<pre>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
</pre>" <p>For printing the results of automated Frama-C tests, we need printing functions with the following properties:</p>
<ol>
<li>they should print all the information available (doing otherwise might mask an existing bug);</li>
<li>they should produce results readable enough for a person to validate the result of each test once;</li>
<li>and their output should not differ depending on the execution platform.</li>
</ol>
<p>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 <code>Format.printf</code> 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 <code>d</code> in the program <code>d = 3.14;</code> was printed as <code>d ∈ 0x1.91eb851eb851fp1</code>. In this notation, the digits after <code>1.</code> are in base 16, and <code>p1</code> means \multiplied by <code>2^1</code>" by similarity with <code>e1</code> meaning "multiplied by <code>10^1</code>". 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 <code>%a</code> format code in <code>printf()</code>. It is easy to see at a glance that the <code>d</code> is between <code>3</code> and <code>4</code> closer to <code>3</code> because <code>0x1.8p1</code> is <code>3</code> and <code>0x2.0p1</code> is <code>4</code>.</p>
<p>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 <code>d = 3.14;</code> the new decimal pretty-printing function produces <code>d ∈ 3.1400000000000001</code> which is what could be expected considering that the number <code>3.14</code> cannot be represented exactly as a <code>double</code>. More important is that if the program is modified to read <code>d = 0x1.91eb851eb8520p1;</code> changing only the least significant bit in the floating-point representation of <code>d</code> the output becomes <code>d ∈ 3.1400000000000005</code> showing that there is a difference.</p>
<p>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.</p>
<p>It works by first isolating the sign exponent and mantissa of the floating-point number to print. The mantissa ranges over <code>0 .. 0xFFFFFFFFFFFFF</code> and we would like a sequence of decimal digits that ranges over <code>0 .. 9999999999999998</code> instead. For this purpose the mantissa is bitwise or-ed into the representation of the double-precision float <code>1.0</code> giving a result that ranges over <code>1.0 .. 2.0</code> (excluded). The basic idea is then to subtract <code>1.0</code> to multiply by <code>1e16</code> to convert to an integer and to print with width 16 and with leading zeroes (the <code>%016Ld</code> format).</p>
<p>The function is at the bottom of this post. The recent improvement is this test:</p>
<pre> let d re =
if d >= 1.5
then d -. 1.5 5000000000000000L
else d -. 1.0 0L
in
...
</pre>
<p>Without these lines the function was a bit unsatisfactory when the fractional part was above <code>0.5</code>. 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.</p>
<pre>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
</pre>"
{% endraw %}