## 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\nf2:%.9g\n", 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.

1. On Monday, January 2 2012, 19:26 by Rick Regan

Pascal,

"Recent versions of GCC go to great lengths to offer correct rounding to nearest (using the MPFR library)."

On the last version of gcc I checked (4.4.3) it still got some conversions wrong, despite its use of MPFR. Bug report http://gcc.gnu.org/bugzilla/show_bu... addresses this -- do you know if it has been fixed?

Rick

2. On Monday, January 2 2012, 22:28 by pascal

Hello Rick.

I am glad you found this blog. I enjoy yours, both for the simple algorithm that, translated to code, portably solved my own problem and for the explanations of David Gay's dtoa().

I have to admit that I was bluffing when I wrote that GCC had correct rounding because it used MPFR. I later looked at the Changelog entry, and the MPFR dependency was initially introduced to solve another problem: "The GCC middle-end has been integrated with the MPFR library. This allows GCC to evaluate and replace at compile-time calls to built-in math functions having constant arguments with their mathematically equivalent results". See also comment 4, Joseph S. Myers dated 2007-02 in the bug report you sent in.

With both GCC versions I have, 4.4.3-4ubuntu5 and SVN snapshot 172652, hosted on amd64 and therefore with no reason to use 80-bit extended doubles, all examples in that thread (the initial one and your two constants) are handled correctly. Indeed, with both these compilers, your entire examples.c is handled correctly. On what host platform do you observe the rounding bugs?

3. On Tuesday, January 3 2012, 03:07 by Rick Regan

Thanks for checking out the examples.

I ran on 32-bit/x87 FPU. 64-bit/SSE seems not to have the problem. I had wanted to look into why this happens for these examples but hadn't gotten to it (maybe I will now).

I've been reading through some of your floating-point related posts (looking forward to more).

4. On Tuesday, January 3 2012, 18:49 by Rick Regan

I reran my examples on my 32 bit machine, but updated with Ubuntu version 10.10, gcc version 4.4.5, and eglibc 2.12.1. I get the same incorrect results I got previously. I don't know why 32-bit and 64-bit should differ, since I didn't think real.c (gcc's conversion routine) used IEEE floating-point. It will have to remain a mystery for now.

(BTW, I did not author the referenced gcc bug report).

5. On Tuesday, January 3 2012, 19:37 by Rick Regan

On another topic: I don't see how the double rounding error could occur in your example, 1.01161128282547. It converts to 1.0000001011111000111101010000000000000000000000000000001... In both cases that should round to 1.0000001011111000111101, and print as 1.01161122. But being the curious type I tried this under Visual C++ and Linux gcc: Visual C++ prints both as 1.01161122; gcc (hard coded constants) prints 1.01161134 and 1.01161122; glibc (strtof() and strtod()) prints both as 1.01161122 (I am running in the same environment as per my comment above.) I can't explain gcc's answer.

6. On Tuesday, January 3 2012, 19:53 by Rick Regan

I'm mixed up. This IS a double rounding error. gcc is getting it right, and the others are wrong. Sorry for the confusion. (I bet you're sorry I found your blog now :).)

7. On Wednesday, January 4 2012, 13:59 by pascal

> gcc is getting it right, and the others are wrong. Sorry for the confusion. (I bet you're sorry I found your blog now :).)

And Frama-C's value analysis is getting it right, too :) Well, in the development version.

Actually, I am very happy that you found VC++ had issues too, because this is the sort of thing I would just assume compilers get right.

I looked at what real.c looks like in GCC's SVN. Indeed, it does not rely on the host's floating-point. My previously unsubstantiated claim about GCC "going to great lengths..." holds, although the MPFR part doesn't. They tried to make it correct. I still do not understand how there could be / have been bugs there, but the code in real.c is complex enough to have bugs in 32-bit that do not appear in 64-bit.