should ldexp round correctly

Question

I'm a bit surprised with MSVC ldexp behavior (it happens in Visual Studio 2013, but also with all older versions at least down to 2003...).

For example:

#include <math.h>
#include <stdio.h>
int main()
{
    double g=ldexp(2.75,-1074);
    double e=ldexp(3.0,-1074);

    printf("g=%g e=%g \n",g,e);
    return 0;
}

prints

g=9.88131e-324 e=1.4822e-323

The first one g is strangely rounded...
It is 2.75 * fmin_denormalized, so i definitely expect the second result e.
If I evaluate 2.75*ldexp(1.0,-1074) I correctly get same value as e.

Are my expectations too high, or does Microsoft fail to comply with some standard?

Answer 1

While the question does not explicitly state this, I assume that the output expected by the asker is:

g=1.4822e-323 e=1.4822e-323

This is what we would expect from a C/C++ compiler that promises strict adherence to IEEE-754. The question is tagged both C and C++, I will address C99 here as that is the standard I have in hand.

In Annex F, which describes IEC 60559 floating-point arithmetic (where IEC 60559 is basically another name for IEEE-754) the C99 standard specifies:

An implementation that defines __STDC_IEC_559__ shall conform to the specifications in this annex. [...] The scalbn and scalbln functions in <math.h> provide the scalb function recommended in the Appendix to IEC 60559.

Further down in that annex, section F.9.3.6 specifies:

On a binary system, ldexp(x, exp) is equivalent to scalbn(x, exp).

The appendix referenced by the C99 standard is the appendix of the 1985 version of IEEE-754, where we find the scalb function defined as follows:

Scalb(y, N) returns y × 2^N for integral values N without computing 2^N.

scalb is defined as a multiplication with a power of two, and multiplications must be rounded correctly based on the current rounding mode according to the standard. Therefore, with a conforming C99 compiler ldexp() must return a correctly rounded result if the compiler defines __STDC_IEC_559__. In the absence of a library call setting the rounding mode, the default rounding mode "round to nearest or even" is in effect.

I do not have access to MSVC 2013, so I do not know whether it defines that symbol or not. This could even depend on a compiler flag setting, such as /fp:strict.

After tracking down my copy of the C++11 standard, I cannot find any reference to __STDC_IEC_559__ or any language about IEEE-754 bindings. According to the answer to this question this is because that support is included by referring to the C99 standard.

Answer 2

This happens because during the ldexp calculation the 2.75 gets truncated to 2, which happens because at that small of a denormalized number the bits that represent the '.75' part get shifted off the end of the representable number and disappear. Whether this is a bug or designed behavior can be debated.

When calculating 2.75*ldexp(1.0,-1074) normal rounding happens, and the 2.75 becomes 3.

EDIT: ldexp should round correctly, and this is a bug.