Analyzing single-precision floating-point constants

The previous post on this blog points out how subtle just floating-point constants can be. In a previous previous post exactly one year ago I was already alluding to these difficulties in the context of Frama-C's front-end.

Programming a code analyzer that can answer question 5

How should the program below be parsed to have a chance to be analyzed correctly?

main(){
  float f1 = 1.01161128282547f;
  float f2 = 1.01161128282547;
  return f1 == f2;
}

In order to avoid the very double-rounding that could make variable f1 appear equal to f2 the method I delineated one year ago was to call strtof() to parse the single-precision constant. An analyzer that carelessly parses the initializer of f1 with a call to double-precision function strtod() would wrongly conclude that both variables are double-rounded and therefore equal. In OCaml it is easy to slip because by default strtod() is the only interfaced conversion function.

I have now tried the "interface and call strtof()" method and in this post I report on the lesson I learned. It fits in one sentence.

Calling strtof() from OCaml in order to avoid double rounding issues only works if the function strtof() is correct on the host system in the first place.

The diagram and the program explaining the problem are almost identical from the previous post's question 5:

#include <stdlib.h>
#include <stdio.h>
main(){
  float f1 = strtof("1.01161128282547"  NULL);
  float f2 = strtod("1.01161128282547"  NULL);
  printf("f1:%.9g\
f2:%.9g"  f1  f2);
}

       ~1.01161122                    ~1.01161128                     ~1.01161134
            +------------------------------+------------------------------+
           f2                                ^                            f1
                                       original number

This program when executed should display:

$ ./a.out
f1:1.01161134
f2:1.01161122

It is expected for the program to print 1.01161122 for f2 as the string in the program is converted to a double near 1.01161128 and then to a float (the float near 1.01161122 is chosen). For variable f1 calling function strtof() to convert the string directly to a float should result in the float near 1.01161134 being computed. Indeed this is what the program does on Mac OS X. However:

$ ./a.out
f1:1.01161122
f2:1.01161122
$ uname -a
Linux 2.6.34.6-xxxx-grs-ipv6-64 #3 SMP Fri Sep 17 16:06:38 UTC 2010 x86_64 GNU/Linux

This is a bug in Glibc.

Note that the standard only requires that conversion functions strtod() strtof() and conversions done by the C compiler when parsing must produce a result within one ULP of the exact value in an implementation-defined manner. However I do not think that "trying to return the nearest float but failing in some difficult cases and returning the second nearest" qualifies as "implementation-defined". There's an implementation all right but it's not much of a definition.

An analyzer works on the source code. It is natural during the elaboration of the Abstract Syntax Tree to convert strings from the program's source code into floating-point numbers and the natural way to do this is to carefully call the host's function strtof() of single-precision constants and strtod() on double-precision constants. Unfortunately if the analyzer is running on GNU/Linux all this care is for nothing: the AST will contain wrong floating-point values (wrong in the sense of being different from the values the compiler will use for the same constants) and therefore the analyzer will for instance wrongly conclude that the program in question 5 in last post returns 1.

Mitigation technique

A simple technique to mitigate this issue is to analyze the program from question 5 less precisely taking into account rounding issues so as not to say anything false. This is what the value analysis in Nitrogen does:

$ ~/frama-c-Nitrogen-20111001/bin/toplevel.opt -val q5.c
[...]
          f1 ∈ [1.01161122322 .. 1.01161134243]
          f2 ∈ [1.01161122322 .. 1.01161134243]
          __retres ∈ {0; 1}

The sets computed for f1 and f2 are both over-approximated but correct. As a result of the approximation the analyzer is unsure whether the program returns 0 or 1. This is a bit unsatisfactory because the analysis is imprecise on a simple and deterministic program. But at least it is sound.

The solution: write your own strtof()

According to this post on the Exploring Binary blog it is easy to write your own correctly rounded decimal to binary function as long as you do not worry too much about performance. In the context of static analysis parsing time is small with respect to the time taken by subsequent treatments and it is reassuring to know that the values in the AST are the correct values regardless of bugs in the host's libraries. Following this principle I wrote Frama-C's own Caml conversion function based on the outline there. Preliminary tests indicate that this function is adequate for this usage. It will in all likelihood be used for both single and double precision constants in the next Frama-C release.

This means that the value analysis in the Oxygen release will be able to give precise and correct answers to all questions in the floating-point quiz from last post. No "maybe it is and maybe it ain't" answers as Nitrogen gives on this quiz and no incorrect answers like very old Frama-C versions gave. Here is the development version answering question 5:

$ bin/toplevel.opt -val q5.c
[...]
warning: Floating-point constant 1.01161128282547f is not represented exactly.
  Will use 0x1.02f8f60000000p0. See documentation for option -warn-decimal-float
[...]
          f1 ∈ 1.01161134243
          f2 ∈ 1.01161122322
          __retres ∈ {0}

Small approximations soon snowball into large ones so I think it is important to be able to give the most precise value to the simplest expressions of the C language constants.

Should we worry about being too precise ?

The problem encountered with strtof() also highlights the relative fuzziness of the standard. If Frama-C is too precise it may omit the actual behavior of the program as compiled by the compiler. This would be unfortunate. Recent versions of GCC go to great lengths to offer correct rounding to nearest (using the MPFR library). Another compiler may itself use strtof() for parsing. For covering both the possibilities that the program is compiled with GCC or this hypothetical other compiler Nitrogen's imprecise handling of single-precision constants was better.

I would say that a Frama-C plug-in may wish to offer both modes through a command-line option. The value analysis probably will eventually. Regardless it is also important that Frama-C always provides consistent results regardless of the host platform so it was necessary to implement a correctly rounded string -> float function even if we then add and subtract an ULP here and there in order to capture all possible compilations.

Conclusion

In the introduction I self-centeredly linked to two earlier posts from this blog. The post from one year ago showed one difficulty as an example of why single-precision floating-point was not supported precisely in Frama-C in general. The post from last week showed several contrived examples that can only be handled with a precise understanding of floating-point —including but not limited to the issue described one year ago— in the form of a quiz. What I find interesting is that there were unforeseen difficulties in the implementation of even the foreseen difficulty of getting the right value for a single-precision constant. And to reiterate the current development version of the value analysis correctly answers all five questions in the quiz.

This post was improved by suggestions from Fabrice Derepas and Florent Kirchner.