A convert-to-floating-point function at the frontier of formal verification

The function in this post to the musl mailing list is a good candidate for formal verification. It should just barely be possible to verify the absence of undefined behaviors with Frama-C's value analysis and it should just barely be possible to verify it functionally with one of the Jessie or Wp plug-ins.

Context

Musl is a libc project. Key quote: “musl is lightweight fast simple free and strives to be correct in the sense of standards-conformance and safety”.

The function highlighted in this post converts the decimal representation of a number to a long double in a correctly rounded manner. The long double type is not very portable and therefore I do not expect any Frama-C plug-in to support it but after replacing all instances of long double with double and choosing the corresponding LDBL_MANT_DIG and LDBL_MAX_EXP the function becomes a plain decimal-to-double conversion function that it should be possible to verify formally with Frama-C.

Motivation

The formal verification I am encouraging the valiant reader to attempt would be useful too! This is a brand-new implementation. I was looking for one such decimal-to-binary function at about the time I was writing this post. At the time there were about four categories of Open Source implementations:

1. the very wrong version available in Ruby and Tcl to compile on platforms that lack strtod(). This one can return a result off by 15 ULPs or so;
2. the slightly wrong and somewhat complicated version (Glibc);
3. the correctly rounded but non-portable and terribly complicated version (David M. Gay's implementation);
4. and the naïve version (the pseudo-code from Rick Regan that I ended up basing my own function on).

The new function I am suggesting to formally verify nicely fills the gap between 3. and 4. above. It is quite clean and portable but from a quick code review I expect it to be faster than say my own implementation. It is written to be correctly rounded according to the FPU's rounding mode. It seems to use base-10⁹ multi-precision integers. I am still trying to make sense of parts of it.