Re: [lmi] MinGW-w64 anomaly?

[Vadim Zeitlin]

On Sun, 18 Dec 2016 14:28:46 +0000 Greg Chicares <address@hidden> wrote:

GC> Vadim--Could you please take a look at this? It's an extremely simplified
GC> test case, but it's important for the way lmi performs rounding.
GC> 
GC> Apply the patch below [0] and run:
GC>   make unit_tests unit_test_targets=sandbox_test.exe
GC> 
GC> Cross-compiling for msw on debian with MinGW-w64 gcc-4.9.1, I see:

[Let me start with a digression not affecting lmi directly, so please feel
 free to skip it completely]

 Interestingly enough, this discrepancy doesn't appear if -std=c++11 option
is not used: I started testing this outside of lmi initially and saw the
same results in all 4 cases. I had to add -std=c++11 to see the same thing
as in lmi unit test. This was surprising, so I decided to check why was it
the case and it turned out that <cmath> has a std::pow() overload for
"(long double, int)" using __builtin_powil() in pre-C++11 mode, but uses
__builtin_powl() for all types of arguments for long double in C++11. And
according to http://en.cppreference.com/w/cpp/numeric/math/pow this is
correct, however I think it might be one of a few cases in which this
reference is wrong. It's quite understandable in this particular case
because it seems that the standard went through several iterations and, at
some point, these functions were indeed removed (see e.g. the proposal at
http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2011/n3286.html#550).
But the answer at http://stackoverflow.com/a/5627278/15275, written by
Howard Hinnant himself, says that this overload finally hasn't been
removed, so MinGW-w64 shouldn't have this "#if __cplusplus < 201103L" check
in its <cmath> -- even though it does in the version we're using and still
seems to be the case in the very latest version of libstdc++, see
https://github.com/gcc-mirror/gcc/blob/6514c52372bb36e555843b0419d71cf55eda0409/libstdc++-v3/include/c_global/cmath#L395
So I opened https://gcc.gnu.org/bugzilla/show_bug.cgi?id=78851 to see if it
could/should be fixed.

[End of the digression]


GC> Running sandbox_test:
GC> f2():
GC>    0.0000000001000000000000000000018377930496 1/10^10
GC>    0.0000000001000000000000000000397031165002  10^-10

 Just for the reference, here are the binary (i.e. exact) representations
of these numbers in 128 bit IEEE format:

 1.b7cdfd9d7bdba b7e0000000005b1p-34
 1.b7cdfd9d7bdba b8a0000000002b5p-34

where I added a space after the first 52 bits which is all we have in the
in-memory double representation. This at least confirms that both values
are rounded to the same 1.b7cdfd9d7bdbbp-34 when using 64 bit doubles.

GC> 
GC> f3():
GC>    0.0000000001000000000000000000018377930496 1/10^10
GC>    0.0000000001000000000000000000018377930496  10^-10
GC> 
GC> What concerns me is the second of those four numbers: it's less accurate
GC> than the other three. If an 80-byte register is being spilled to 64-byte
GC> storage, I'd expect the result to be less accurate than this. I just
GC> can't see any way to explain it.

 FWIW I could confirm, using my recently found knowledge of
-fexcess-precision option, that this is not related to spilling x87
registers to memory. If I convert your example to the sensibly equivalent C
version:
---------------------------------- >8 --------------------------------------
#define __USE_MINGW_ANSI_STDIO 1
#include <math.h>
#include <stdio.h>

void f2(long double r)
{
    int n = -10;
    printf
        ("f2():\n"
         "%.40Lf 1/10^10\n"
         "%.40Lf  10^-10\n"
         "\n",
         1.0L / powl(r, -n),
                powl(r,  n)
        );
}

void f3()
{
    long double r = 10.0L;
    int n = -10;
    printf
        ("f3():\n"
         "%.40Lf 1/10^10\n"
         "%.40Lf  10^-10\n"
         "\n",
         1.0L / powl(r, -n),
                powl(r,  n)
        );
}

int main()
{
    f2(10.0L);
    f3();
    return 0;
}
---------------------------------- >8 --------------------------------------
and compile it with -fexcess-precision=standard (and still with -O2, of
course), it doesn't change anything and I still get the same results.

 But in the C version it's very simple to see what goes on however: f3()
basically just prints the constant value, while f2() calls _powl()
function. So what you see is basically that the results of compile-time
computation are not exactly the same as the results of run-time
computation. And I can even explain this further: I'm pretty sure that
compile-time results are computed using the same __builtin_powil() that is
used in C++98 code and it's not really surprising that its much simpler
algorithm using just multiplications and inversion gives different results
from the general-purpose powl(). Notice that __builtin_powil() is, of
course, fully inlined and you get efficiently computed x^-10 by
successively computing, using the single multiplication instruction, x,
x^2, x^3, x^5, x^10 and then inverting it using just one other instruction,
so it's clearly not something powl() can do -- it's already nice that it
gives the same result as compile-time calculation for x^10 itself.

 BTW, if you add -fwhole-program (or -flto, but I can't examine its results
as easily as the assembler code) option to gcc command line, then f2()
would be reduced to the same trivial code (and both functions will be
inlined directly into main() too).


 In the C++ version it's a bit less obvious to understand what goes on, in
part because there is so much more code (you really appreciate the benefits
of printf() compared to iostream when examining disassembly of the code
using them: the C function contains ~10 instructions, the C++ one contains
~100 of them), but it's the same thing, i.e. f3() just outputs a literal
while f2() computes it using _powl().


GC> I can work around this by taking the reciprocal of the product of N tens,
GC> which seems always to be as accurate as the best case with std::pow(),
GC> but I'd rather avoid that.

 I'm not sure what exactly is worth working around here. Both results are
correct for the precision of standard double type and you would get the
same output if you stored them in a "volatile double". But if you really
want to get best performance and identical output in both cases, you could
use __builtin_powil() manually in f2().

 I hope this answers your question, please let me know if anything is still
unclear.

 Regards,
VZ