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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
---
layout: post
author: Pascal Cuoq
date: 2011-02-10 09:11 +0200
categories: floating-point value
format: xhtml
title: "Verifying numerical precision with Frama-C's value analysis — part 2"
summary:
---
{% raw %}
<p>In this sequel to <a href="/index.php?post/2011/01/22/Verifying-numerical-precision-with-Frama-C-s-value-analysis" hreflang="en">a previous post about numerical precision and the value analysis</a> we tackle another extremely common implementation technique the linear interpolation table. In an attempt to make things less boring the approximated function this time is a sine and takes as its argument a double representing an angle in degrees. The interpolation table contains values for each integral angle between 0 and 90.</p>
<p>As before we are focusing on verification so the code already exists and probably works fine too. Our goal here is to double-check that it works. We make up the specification as we go along which is the main difference with real life where the specification would have been defined before the function was written.</p>
<pre>double table_sin[91] =
{ 0.000000 /* 0 */
0.017452 /* 1 */
0.034899 /* 2 */
0.052336 /* 3 */
...
0.998630 /* 87 */
0.999391 /* 88 */
0.999848 /* 89 */
1.000000 /* 90 */ } ;
/*@ requires 0.0 <= degrees <= 90.0 ; */
double interp_sin(double degrees)
{
int i = (int) degrees;
double r = degrees - i;
double result = table_sin[i] * (1-r) + table_sin[i+1] * r;
return result;
}
</pre>
<p>You can also <a href="/assets/img/blog/imported-posts/sin1.c">download sin1.c</a></p>
<p>The tools we have are the same as last time. The first step is to run a basic value analysis with command <code>frama-c -val sin1.c -main interp_sin</code>. Precision should not be expected from this first run:</p>
<pre>[value] Values for function interp_sin:
i ∈ [0..90]
r ∈ [-90. .. 90.]
result ∈ [-179. .. 181.]
</pre>
<p>But let us pay attention to another important piece of information in the log:</p>
<pre>sin1.c:99:[kernel] warning: accessing out of bounds index. assert (0 ≤ (int )(i+1)) ∧
((int )(i+1) < 91);
</pre>
<p>The index <code>i</code> is computed from the argument <code>degrees</code> and the value analysis was able to take advantage of the precondition limiting the value of <code>degrees</code> to 90.0 but the function actually accesses <code>table_sin[i+1]</code>. So we have found a bug which was not our initial goal. Since this was an honest mistake I made when preparing this post I thought I would leave it in. Do not complain for I did spare you the previous mistake I made (I initially wrote <code>double table_sin[90] = ...</code>).</p>
<blockquote><p>Digression: for the mistake in declaring the size of <code>table_sin</code> as 90 any serious compiler would have been able to spot that something was wrong because initial values were provided for 91 cells (e.g. <code>sin_interp.c:92: warning: excess elements in array initializer</code>). Frama-C does not make any effort to duplicate work that has already been done many times over in compilers so you should not necessarily expect it to warn about local inconsistencies such as this one. If you want the best of both worlds <strong>use the value analysis in addition to your compiler's warnings</strong>. Frama-C's value analysis would still warn about accessing a <code>table_sin</code> of size 90 out of bounds though.</p>
</blockquote>
<p>Since <code>interp_sin</code>'s contract is that it is called with argument <code>degrees</code> at most 90.0 one simple fix is to add one element in <code>table_sin</code>:</p>
<pre>double table_sin[92] =
...
0.999848 /* 89 */
1.000000 /* 90 */
1.000000 /* work around overflow at i+1 for degrees=90.0 */ } ;
</pre>
<p>This does not make the value analysis any more precise but at least the alarm disappears.</p>
<p>The next step is to gain precision by sub-dividing input intervals. We subdivide the first few degrees in quarter-degree intervals and for clarity leave the largest part of the input range ([4.25 .. 90.0]) undivided.</p>
<pre>/*@ requires 0.0 <= degrees <= 90.0 ; */
double interp_sin(double degrees)
{
/*@ assert
degrees <= 0.25 ||
0.25 <= degrees <= 0.5 ||
0.5 <= degrees <= 0.75 ||
...
3.75 <= degrees <= 4.0 ||
4.0 <= degrees <= 4.25 ||
4.25 <= degrees ; */
int i = (int) degrees;
double r = degrees - i;
double result = table_sin[i] * (1-r) + table_sin[i+1] * r;
Frama_C_show_each_d(degrees i r result);
return result;
}
</pre>
<p><a href="/assets/img/blog/imported-posts/sin2.c">Download sin2.c</a></p>
<p>This function also contains computations that can benefit from option <code>-subdivide-float-var</code>. This option <a href="/index.php?post/2011/01/22/Verifying-numerical-precision-with-Frama-C-s-value-analysis" hreflang="en">described in last post</a> automatically refines the precision of the analysis without requiring any annotations. The command-line thus grows to:</p>
<pre>frama-c -main interp_sin sin2.c -slevel 100 -val -subdivide-float-var 50 -float-relative
</pre>
<p>It is normal for the result corresponding to the large undivided range from 4.25 to 90.0 to be imprecise:</p>
<pre>[value] Called Frama_C_show_each_d([4.25 ++ 85.75] [4..90] [-85.75 ++ 171.75]
[-79.698667 ++ 159.769407])
</pre>
<p>The interesting part is what happens for the other intervals:</p>
<pre>degrees i r result
[4. ++ 0.25] {4; } [0. ++ 0.25] [0.06975478125 ++ 0.00435161819458]
[3.75 ++ 0.25] {3; 4; } [-0.25 ++ 1.25] [0.043630998695 ++ 0.0435276801381]
[3.5 ++ 0.25] {3; } [0.5 ++ 0.25] [0.0610458515625 ++ 0.004355198349]
[3.25 ++ 0.25] {3; } [0.25 ++ 0.25] [0.0566908515625 ++ 0.004355198349]
[3. ++ 0.25] {3; } [0. ++ 0.25] [0.0523358515625 ++ 0.004355198349]
[2.75 ++ 0.25] {2; 3; } [-0.25 ++ 1.25] [0.0261847499594 ++ 0.0435715837503]
[2.5 ++ 0.25] {2; } [0.5 ++ 0.25] [0.0436174958397 ++ 0.0043592552011]
[2.25 ++ 0.25] {2; } [0.25 ++ 0.25] [0.0392582458397 ++ 0.0043592552011]
[2. ++ 0.25] {2; } [0. ++ 0.25] [0.0348989958397 ++ 0.0043592552011]
[1.75 ++ 0.25] {1; 2; } [-0.25 ++ 1.25] [0.008731 ++ 0.0436050104082]
[1.5 ++ 0.25] {1; } [0.5 ++ 0.25] [0.026175499998 ++ 0.00436175000263]
[1.25 ++ 0.25] {1; } [0.25 ++ 0.25] [0.021813749998 ++ 0.00436175000263]
[1. ++ 0.25] {1; } [0. ++ 0.25] [0.017451999998 ++ 0.00436175000263]
[0.75 ++ 0.25] {0; 1; } [-0.25 ++ 1.25] [-0.00872475 ++ 0.0436237500051]
[0.5 ++ 0.25] {0; } [0.5 ++ 0.25] [0.008726 ++ 0.004363]
[0.25 ++ 0.25] {0; } [0.25 ++ 0.25] [0.004363 ++ 0.004363]
[-0. ++ 0.25] {0; } [-0. ++ 0.25] [0. ++ 0.004363]
</pre>
<p>Every few input intervals the result interval is ten times wider than usual. This could become a problem. Indeed if you try further to reduce the width of of one of the problematic input intervals you will notice that there always remains one output interval that refuses to get slimmer.</p>
<p>Analyzing this issue it turns out the reason for the irreducible imprecision is that an input interval such as [2.75 .. 3.0] for <code>degrees</code> means that <code>i</code> can be either 2 or 3 which means that <code>r</code> gets approximated. Interval arithmetics cannot see that <code>r</code> computed as <code>degrees - i</code> is always positive. The solution is to <strong>subdivide smarter</strong> so that all the values inside one same interval for <code>degrees</code> produce the same value for <code>i</code>:</p>
<pre> /*@ assert
degrees < 0.25 ||
0.25 <= degrees < 0.5 ||
0.5 <= degrees < 0.75 ||
...
3.75 <= degrees < 4.0 ||
4.0 <= degrees < 4.25 ||
4.25 <= degrees ; */
</pre>
<p>If you are using <strong>version 20110201 or later</strong> the value analysis is able to guarantee that just as before the assertion is always true and thus that using it does not incur any additional proof obligation:</p>
<pre>sin3.c:98:[value] Assertion got status valid.
</pre>
<p>With the new assertion all the result intervals are uniformly narrow with a width less than 0.005 for all the few intervals we are currently looking at:</p>
<pre>[value] Called Frama_C_show_each_d([4. ++ 0.25] {4; } [0. ++ 0.25]
[0.06975478125 ++ 0.00435161819458])
[value] Called Frama_C_show_each_d([3.75 ++ 0.25] {3; } [0.75 ++ 0.25]
[0.0654008515625 ++ 0.004355198349])
[value] Called Frama_C_show_each_d([3.5 ++ 0.25] {3; } [0.5 ++ 0.25]
[0.0610458515625 ++ 0.004355198349])
...
[value] Called Frama_C_show_each_d([0.25 ++ 0.25] {0; } [0.25 ++ 0.25]
[0.004363 ++ 0.004363])
[value] Called Frama_C_show_each_d([-0. ++ 0.25] {0; } [-0. ++ 0.25]
[0. ++ 0.004363])
</pre>
<p>The input sub-intervals appear the same as before in the log but they are subtly different. The best way to see the difference is to request floating-point intervals to be displayed exactly to the bit with option <code>-float-hex</code>. It takes some time getting used to hexadecimal floating-point numbers; if you wish do not worry about it and skip to the part about Emacs Lisp at the end of the post.</p>
<pre>frama-c -main interp_sin sin3.c -slevel 100 -val -subdivide-float-var 50 -float-hex
...
[value] Called Frama_C_show_each_d([0x1.0000000000000p2 .. 0x1.0ffffffffffffp2]
{4; }
[0x0.0000000000000p-1022 .. 0x1.fffffffffffe0p-3]
[0x1.1db73083558a7p-4 .. 0x1.2f8a31209edbep-4])
[value] Called Frama_C_show_each_d([0x1.e000000000000p1 .. 0x1.fffffffffffffp1]
{3; }
[0x1.8000000000000p-1 .. 0x1.ffffffffffffcp-1]
[0x1.0be1c36976bc2p-4 .. 0x1.1db8851116a8bp-4])
[value] Called Frama_C_show_each_d([0x1.c000000000000p1 .. 0x1.dffffffffffffp1]
{3; }
[0x1.0000000000000p-1 .. 0x1.7fffffffffffcp-1]
[0x1.f4166e008e9b4p-5 .. 0x1.0be1f8a7e73a3p-4])
[value] Called Frama_C_show_each_d([0x1.a000000000000p1 .. 0x1.bffffffffffffp1]
{3; }
[0x1.0000000000000p-2 .. 0x1.ffffffffffff8p-2]
[0x1.d069552e2fbe3p-5 .. 0x1.f416d87d6f976p-5])
</pre>
<p>The first of these input intervals is when displayed exactly [0x1.0p2 .. 0x1.0ffffffffffffp2]. With the assertion used in sin2.c the corresponding interval was [0x1.0p2 .. 0x1.10p2] or in decimal [4. .. 4.25]. The number 0x1.0ffffffffffffp2 is the first <code>double</code> immediately below 4.25. Similarly and more importantly 0x1.fffffffffffffp1 is the first <code>double</code> immediately below 4.0 and the value analysis knows that the truncation to <code>int</code> of this number as that of any number in [0x1.e0p1 .. 0x1.fffffffffffffp1] (in decimal [3.75 .. 4.0]) is 3.</p>
<p>The next step in reproducing the method from <a href="/index.php?post/2011/01/22/Verifying-numerical-precision-with-Frama-C-s-value-analysis" hreflang="en">part 1</a> is to use option <code>-all-rounding-modes</code> to capture all possible FPU rounding modes extra precision introduced by the compiler for temporary results and incidentally to get result intervals that capture what the program would have computed if it had used reals instead of doubles:</p>
<pre>frama-c -main interp_sin sin3.c -slevel 100 -val -subdivide-float-var 50 -all-rounding-modes -float-relative
...
[value] Called Frama_C_show_each_d([3.75 ++ 0.25] {3; 4; } [-0.25 ++ 1.25]
[0.043630998695 ++ 0.0435276801381])
</pre>
<p>Unfortunately option <code>-all-rounding-modes</code> brings back the imprecision that we had just found a way around. Indeed by instructing the analyzer to take into account the possibility that the compiler might be using 80-bit extended doubles for representing variables typed as doubles in the program we are preventing it from reasoning its way from <code>degrees < 4.0</code> to <code>degrees <= 0x1.fffffffffffffp1</code>.</p>
<p>There is no convenient way around this limitation. You can write <code>degrees <= 0x1.fffffffffffffp1</code> in the assertion instead of <code>degrees < 4.0</code> but then the assertion will not be provable with option <code>-all-rounding-modes</code> for exactly the same reason. With extended doubles or for that matter reals there exists numbers between 0x1.fffffffffffffp1 and 4.0.</p>
<p>There remains the possibility to compare each input output interval pair to a reference function known to be increasing on [0.0 .. 90.0]. Since this post is getting long I will leave this as an optional exercise as well as actually finishing to sub-divide the interval [4.25 .. 90.0]. Spam forces us to limit the possibility to comment on this blog but should these exercises lead to any remarks they will always be welcome on <a href="http://lists.gforge.inria.fr/cgi-bin/mailman/listinfo/frama-c-discuss" hreflang="en">the mailing list</a>. One most welcome remark from someone familiar with Emacs Lisp would be a ready-made macro to generate subdivision assertions very quickly in the Emacs text editor.</p>
<p>Many thanks to Zaynah Dargaye and David Delmas for their remarks.</p>
<p>In this sequel to <a href="/index.php?post/2011/01/22/Verifying-numerical-precision-with-Frama-C-s-value-analysis" hreflang="en">a previous post about numerical precision and the value analysis</a> we tackle another extremely common implementation technique the linear interpolation table. In an attempt to make things less boring the approximated function this time is a sine and takes as its argument a double representing an angle in degrees. The interpolation table contains values for each integral angle between 0 and 90.</p>
<p>As before we are focusing on verification so the code already exists and probably works fine too. Our goal here is to double-check that it works. We make up the specification as we go along which is the main difference with real life where the specification would have been defined before the function was written.</p>
<pre>double table_sin[91] =
{ 0.000000 /* 0 */
0.017452 /* 1 */
0.034899 /* 2 */
0.052336 /* 3 */
...
0.998630 /* 87 */
0.999391 /* 88 */
0.999848 /* 89 */
1.000000 /* 90 */ } ;
/*@ requires 0.0 <= degrees <= 90.0 ; */
double interp_sin(double degrees)
{
int i = (int) degrees;
double r = degrees - i;
double result = table_sin[i] * (1-r) + table_sin[i+1] * r;
return result;
}
</pre>
<p>You can also <a href="/assets/img/blog/imported-posts/sin1.c">download sin1.c</a></p>
<p>The tools we have are the same as last time. The first step is to run a basic value analysis with command <code>frama-c -val sin1.c -main interp_sin</code>. Precision should not be expected from this first run:</p>
<pre>[value] Values for function interp_sin:
i ∈ [0..90]
r ∈ [-90. .. 90.]
result ∈ [-179. .. 181.]
</pre>
<p>But let us pay attention to another important piece of information in the log:</p>
<pre>sin1.c:99:[kernel] warning: accessing out of bounds index. assert (0 ≤ (int )(i+1)) ∧
((int )(i+1) < 91);
</pre>
<p>The index <code>i</code> is computed from the argument <code>degrees</code> and the value analysis was able to take advantage of the precondition limiting the value of <code>degrees</code> to 90.0 but the function actually accesses <code>table_sin[i+1]</code>. So we have found a bug which was not our initial goal. Since this was an honest mistake I made when preparing this post I thought I would leave it in. Do not complain for I did spare you the previous mistake I made (I initially wrote <code>double table_sin[90] = ...</code>).</p>
<blockquote><p>Digression: for the mistake in declaring the size of <code>table_sin</code> as 90 any serious compiler would have been able to spot that something was wrong because initial values were provided for 91 cells (e.g. <code>sin_interp.c:92: warning: excess elements in array initializer</code>). Frama-C does not make any effort to duplicate work that has already been done many times over in compilers so you should not necessarily expect it to warn about local inconsistencies such as this one. If you want the best of both worlds <strong>use the value analysis in addition to your compiler's warnings</strong>. Frama-C's value analysis would still warn about accessing a <code>table_sin</code> of size 90 out of bounds though.</p>
</blockquote>
<p>Since <code>interp_sin</code>'s contract is that it is called with argument <code>degrees</code> at most 90.0 one simple fix is to add one element in <code>table_sin</code>:</p>
<pre>double table_sin[92] =
...
0.999848 /* 89 */
1.000000 /* 90 */
1.000000 /* work around overflow at i+1 for degrees=90.0 */ } ;
</pre>
<p>This does not make the value analysis any more precise but at least the alarm disappears.</p>
<p>The next step is to gain precision by sub-dividing input intervals. We subdivide the first few degrees in quarter-degree intervals and for clarity leave the largest part of the input range ([4.25 .. 90.0]) undivided.</p>
<pre>/*@ requires 0.0 <= degrees <= 90.0 ; */
double interp_sin(double degrees)
{
/*@ assert
degrees <= 0.25 ||
0.25 <= degrees <= 0.5 ||
0.5 <= degrees <= 0.75 ||
...
3.75 <= degrees <= 4.0 ||
4.0 <= degrees <= 4.25 ||
4.25 <= degrees ; */
int i = (int) degrees;
double r = degrees - i;
double result = table_sin[i] * (1-r) + table_sin[i+1] * r;
Frama_C_show_each_d(degrees i r result);
return result;
}
</pre>
<p><a href="/assets/img/blog/imported-posts/sin2.c">Download sin2.c</a></p>
<p>This function also contains computations that can benefit from option <code>-subdivide-float-var</code>. This option <a href="/index.php?post/2011/01/22/Verifying-numerical-precision-with-Frama-C-s-value-analysis" hreflang="en">described in last post</a> automatically refines the precision of the analysis without requiring any annotations. The command-line thus grows to:</p>
<pre>frama-c -main interp_sin sin2.c -slevel 100 -val -subdivide-float-var 50 -float-relative
</pre>
<p>It is normal for the result corresponding to the large undivided range from 4.25 to 90.0 to be imprecise:</p>
<pre>[value] Called Frama_C_show_each_d([4.25 ++ 85.75] [4..90] [-85.75 ++ 171.75]
[-79.698667 ++ 159.769407])
</pre>
<p>The interesting part is what happens for the other intervals:</p>
<pre>degrees i r result
[4. ++ 0.25] {4; } [0. ++ 0.25] [0.06975478125 ++ 0.00435161819458]
[3.75 ++ 0.25] {3; 4; } [-0.25 ++ 1.25] [0.043630998695 ++ 0.0435276801381]
[3.5 ++ 0.25] {3; } [0.5 ++ 0.25] [0.0610458515625 ++ 0.004355198349]
[3.25 ++ 0.25] {3; } [0.25 ++ 0.25] [0.0566908515625 ++ 0.004355198349]
[3. ++ 0.25] {3; } [0. ++ 0.25] [0.0523358515625 ++ 0.004355198349]
[2.75 ++ 0.25] {2; 3; } [-0.25 ++ 1.25] [0.0261847499594 ++ 0.0435715837503]
[2.5 ++ 0.25] {2; } [0.5 ++ 0.25] [0.0436174958397 ++ 0.0043592552011]
[2.25 ++ 0.25] {2; } [0.25 ++ 0.25] [0.0392582458397 ++ 0.0043592552011]
[2. ++ 0.25] {2; } [0. ++ 0.25] [0.0348989958397 ++ 0.0043592552011]
[1.75 ++ 0.25] {1; 2; } [-0.25 ++ 1.25] [0.008731 ++ 0.0436050104082]
[1.5 ++ 0.25] {1; } [0.5 ++ 0.25] [0.026175499998 ++ 0.00436175000263]
[1.25 ++ 0.25] {1; } [0.25 ++ 0.25] [0.021813749998 ++ 0.00436175000263]
[1. ++ 0.25] {1; } [0. ++ 0.25] [0.017451999998 ++ 0.00436175000263]
[0.75 ++ 0.25] {0; 1; } [-0.25 ++ 1.25] [-0.00872475 ++ 0.0436237500051]
[0.5 ++ 0.25] {0; } [0.5 ++ 0.25] [0.008726 ++ 0.004363]
[0.25 ++ 0.25] {0; } [0.25 ++ 0.25] [0.004363 ++ 0.004363]
[-0. ++ 0.25] {0; } [-0. ++ 0.25] [0. ++ 0.004363]
</pre>
<p>Every few input intervals the result interval is ten times wider than usual. This could become a problem. Indeed if you try further to reduce the width of of one of the problematic input intervals you will notice that there always remains one output interval that refuses to get slimmer.</p>
<p>Analyzing this issue it turns out the reason for the irreducible imprecision is that an input interval such as [2.75 .. 3.0] for <code>degrees</code> means that <code>i</code> can be either 2 or 3 which means that <code>r</code> gets approximated. Interval arithmetics cannot see that <code>r</code> computed as <code>degrees - i</code> is always positive. The solution is to <strong>subdivide smarter</strong> so that all the values inside one same interval for <code>degrees</code> produce the same value for <code>i</code>:</p>
<pre> /*@ assert
degrees < 0.25 ||
0.25 <= degrees < 0.5 ||
0.5 <= degrees < 0.75 ||
...
3.75 <= degrees < 4.0 ||
4.0 <= degrees < 4.25 ||
4.25 <= degrees ; */
</pre>
<p>If you are using <strong>version 20110201 or later</strong> the value analysis is able to guarantee that just as before the assertion is always true and thus that using it does not incur any additional proof obligation:</p>
<pre>sin3.c:98:[value] Assertion got status valid.
</pre>
<p>With the new assertion all the result intervals are uniformly narrow with a width less than 0.005 for all the few intervals we are currently looking at:</p>
<pre>[value] Called Frama_C_show_each_d([4. ++ 0.25] {4; } [0. ++ 0.25]
[0.06975478125 ++ 0.00435161819458])
[value] Called Frama_C_show_each_d([3.75 ++ 0.25] {3; } [0.75 ++ 0.25]
[0.0654008515625 ++ 0.004355198349])
[value] Called Frama_C_show_each_d([3.5 ++ 0.25] {3; } [0.5 ++ 0.25]
[0.0610458515625 ++ 0.004355198349])
...
[value] Called Frama_C_show_each_d([0.25 ++ 0.25] {0; } [0.25 ++ 0.25]
[0.004363 ++ 0.004363])
[value] Called Frama_C_show_each_d([-0. ++ 0.25] {0; } [-0. ++ 0.25]
[0. ++ 0.004363])
</pre>
<p>The input sub-intervals appear the same as before in the log but they are subtly different. The best way to see the difference is to request floating-point intervals to be displayed exactly to the bit with option <code>-float-hex</code>. It takes some time getting used to hexadecimal floating-point numbers; if you wish do not worry about it and skip to the part about Emacs Lisp at the end of the post.</p>
<pre>frama-c -main interp_sin sin3.c -slevel 100 -val -subdivide-float-var 50 -float-hex
...
[value] Called Frama_C_show_each_d([0x1.0000000000000p2 .. 0x1.0ffffffffffffp2]
{4; }
[0x0.0000000000000p-1022 .. 0x1.fffffffffffe0p-3]
[0x1.1db73083558a7p-4 .. 0x1.2f8a31209edbep-4])
[value] Called Frama_C_show_each_d([0x1.e000000000000p1 .. 0x1.fffffffffffffp1]
{3; }
[0x1.8000000000000p-1 .. 0x1.ffffffffffffcp-1]
[0x1.0be1c36976bc2p-4 .. 0x1.1db8851116a8bp-4])
[value] Called Frama_C_show_each_d([0x1.c000000000000p1 .. 0x1.dffffffffffffp1]
{3; }
[0x1.0000000000000p-1 .. 0x1.7fffffffffffcp-1]
[0x1.f4166e008e9b4p-5 .. 0x1.0be1f8a7e73a3p-4])
[value] Called Frama_C_show_each_d([0x1.a000000000000p1 .. 0x1.bffffffffffffp1]
{3; }
[0x1.0000000000000p-2 .. 0x1.ffffffffffff8p-2]
[0x1.d069552e2fbe3p-5 .. 0x1.f416d87d6f976p-5])
</pre>
<p>The first of these input intervals is when displayed exactly [0x1.0p2 .. 0x1.0ffffffffffffp2]. With the assertion used in sin2.c the corresponding interval was [0x1.0p2 .. 0x1.10p2] or in decimal [4. .. 4.25]. The number 0x1.0ffffffffffffp2 is the first <code>double</code> immediately below 4.25. Similarly and more importantly 0x1.fffffffffffffp1 is the first <code>double</code> immediately below 4.0 and the value analysis knows that the truncation to <code>int</code> of this number as that of any number in [0x1.e0p1 .. 0x1.fffffffffffffp1] (in decimal [3.75 .. 4.0]) is 3.</p>
<p>The next step in reproducing the method from <a href="/index.php?post/2011/01/22/Verifying-numerical-precision-with-Frama-C-s-value-analysis" hreflang="en">part 1</a> is to use option <code>-all-rounding-modes</code> to capture all possible FPU rounding modes extra precision introduced by the compiler for temporary results and incidentally to get result intervals that capture what the program would have computed if it had used reals instead of doubles:</p>
<pre>frama-c -main interp_sin sin3.c -slevel 100 -val -subdivide-float-var 50 -all-rounding-modes -float-relative
...
[value] Called Frama_C_show_each_d([3.75 ++ 0.25] {3; 4; } [-0.25 ++ 1.25]
[0.043630998695 ++ 0.0435276801381])
</pre>
<p>Unfortunately option <code>-all-rounding-modes</code> brings back the imprecision that we had just found a way around. Indeed by instructing the analyzer to take into account the possibility that the compiler might be using 80-bit extended doubles for representing variables typed as doubles in the program we are preventing it from reasoning its way from <code>degrees < 4.0</code> to <code>degrees <= 0x1.fffffffffffffp1</code>.</p>
<p>There is no convenient way around this limitation. You can write <code>degrees <= 0x1.fffffffffffffp1</code> in the assertion instead of <code>degrees < 4.0</code> but then the assertion will not be provable with option <code>-all-rounding-modes</code> for exactly the same reason. With extended doubles or for that matter reals there exists numbers between 0x1.fffffffffffffp1 and 4.0.</p>
<p>There remains the possibility to compare each input output interval pair to a reference function known to be increasing on [0.0 .. 90.0]. Since this post is getting long I will leave this as an optional exercise as well as actually finishing to sub-divide the interval [4.25 .. 90.0]. Spam forces us to limit the possibility to comment on this blog but should these exercises lead to any remarks they will always be welcome on <a href="http://lists.gforge.inria.fr/cgi-bin/mailman/listinfo/frama-c-discuss" hreflang="en">the mailing list</a>. One most welcome remark from someone familiar with Emacs Lisp would be a ready-made macro to generate subdivision assertions very quickly in the Emacs text editor.</p>
<p>Many thanks to Zaynah Dargaye and David Delmas for their remarks.</p>
{% endraw %}