Verifying numerical precision with Frama-C's value analysis — part 2
Pascal Cuoq - 10th Feb 2011In this sequel to a previous post about numerical precision and the value analysis 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.
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.
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;
}
You can also download sin1.c
The tools we have are the same as last time. The first step is to run a basic value analysis with command frama-c -val sin1.c -main interp_sin. Precision should not be expected from this first run:
[value] Values for function interp_sin:
i ∈ [0..90]
r ∈ [-90. .. 90.]
result ∈ [-179. .. 181.]
But let us pay attention to another important piece of information in the log:
sin1.c:99:[kernel] warning: accessing out of bounds index. assert (0 ≤ (int )(i+1)) ∧
((int )(i+1) < 91);
The index i is computed from the argument degrees and the value analysis was able to take advantage of the precondition limiting the value of degrees to 90.0 but the function actually accesses table_sin[i+1]. 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 double table_sin[90] = ...).
Digression: for the mistake in declaring the size of
table_sinas 90 any serious compiler would have been able to spot that something was wrong because initial values were provided for 91 cells (e.g.sin_interp.c:92: warning: excess elements in array initializer). 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 use the value analysis in addition to your compiler's warnings. Frama-C's value analysis would still warn about accessing atable_sinof size 90 out of bounds though.
Since interp_sin's contract is that it is called with argument degrees at most 90.0 one simple fix is to add one element in table_sin:
double table_sin[92] =
...
0.999848 /* 89 */
1.000000 /* 90 */
1.000000 /* work around overflow at i+1 for degrees=90.0 */ } ;
This does not make the value analysis any more precise but at least the alarm disappears.
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.
/*@ 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;
}
This function also contains computations that can benefit from option -subdivide-float-var. This option described in last post automatically refines the precision of the analysis without requiring any annotations. The command-line thus grows to:
frama-c -main interp_sin sin2.c -slevel 100 -val -subdivide-float-var 50 -float-relative
It is normal for the result corresponding to the large undivided range from 4.25 to 90.0 to be imprecise:
[value] Called Frama_C_show_each_d([4.25 ++ 85.75] [4..90] [-85.75 ++ 171.75]
[-79.698667 ++ 159.769407])
The interesting part is what happens for the other intervals:
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]
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.
Analyzing this issue it turns out the reason for the irreducible imprecision is that an input interval such as [2.75 .. 3.0] for degrees means that i can be either 2 or 3 which means that r gets approximated. Interval arithmetics cannot see that r computed as degrees - i is always positive. The solution is to subdivide smarter so that all the values inside one same interval for degrees produce the same value for i:
/*@ 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 ; */
If you are using version 20110201 or later 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:
sin3.c:98:[value] Assertion got status valid.
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:
[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])
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 -float-hex. 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.
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])
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 double immediately below 4.25. Similarly and more importantly 0x1.fffffffffffffp1 is the first double immediately below 4.0 and the value analysis knows that the truncation to int of this number as that of any number in [0x1.e0p1 .. 0x1.fffffffffffffp1] (in decimal [3.75 .. 4.0]) is 3.
The next step in reproducing the method from part 1 is to use option -all-rounding-modes 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:
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])
Unfortunately option -all-rounding-modes 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 degrees < 4.0 to degrees <= 0x1.fffffffffffffp1.
There is no convenient way around this limitation. You can write degrees <= 0x1.fffffffffffffp1 in the assertion instead of degrees < 4.0 but then the assertion will not be provable with option -all-rounding-modes for exactly the same reason. With extended doubles or for that matter reals there exists numbers between 0x1.fffffffffffffp1 and 4.0.
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 the mailing list. 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.
Many thanks to Zaynah Dargaye and David Delmas for their remarks.