9

I have two similar programs one in C++ and another in D.

The compilation is on on Windows7 64bit, to 64bit binaries.

C++ version, VS 2013:

#include <iostream>
#include <string>

int main(int argc, char* argv[])
{
    float eps = 1.0f;
    float f = 0.0f;
    while (f + eps != f)
        f += 1.0f;

    std::cout << "eps = " + std::to_string(eps) + ", max_f = " + std::to_string(f) << std::endl;
    return 0;
}

D version, DMD v2.066.1:

import std.stdio;
import std.conv;

int main(string[] argv)
{
    float eps = 1.0f;
    float f = 0.0f;
    while (f + eps != f)
        f += 1.0f;

    writeln("eps = " ~ to!string(eps) ~ ", max_f = " ~ to!string(f));
    return 0;
}

C++ version works as expected and finds that f + e == f when f = 16777216.

But D version hungs forever. When I put breakpoint I see that in D version f also 16777216 (after running for some time) and Watch window (I use VisualD) shows that (f + e != f) is 'false' so the loop should be terminate but it's not the case during runtime.

I think assembly could give the answer but I'm not very good with it.

I'm new to D, so it should be the case that I misused the language/compiler (compiled with DMD just as 'dmd test.d' without additional options and also from VS with VisualD with default options). Any ideas what could be wrong with D version of the program? Thanks!

Disassembly:

C++:

000000013F7D1410  mov         rax,rsp  
000000013F7D1413  push        rbp  
000000013F7D1414  lea         rbp,[rax-5Fh]  
000000013F7D1418  sub         rsp,0E0h  
000000013F7D141F  mov         qword ptr [rbp+17h],0FFFFFFFFFFFFFFFEh  
000000013F7D1427  mov         qword ptr [rax+8],rbx  
000000013F7D142B  movaps      xmmword ptr [rax-18h],xmm6  
000000013F7D142F  xorps       xmm1,xmm1  
    float eps = 1.0f;
    float f = 0.0f;
000000013F7D1432  movss       xmm6,dword ptr [__real@3f800000 (013F7D67E8h)]  
000000013F7D143A  nop         word ptr [rax+rax]  
        f += 1.0f;
000000013F7D1440  addss       xmm1,xmm6  
    while (f + eps != f)
000000013F7D1444  movaps      xmm0,xmm1  
000000013F7D1447  addss       xmm0,xmm6  
000000013F7D144B  ucomiss     xmm0,xmm1  
000000013F7D144E  jp          main+30h (013F7D1440h)  
000000013F7D1450  jne         main+30h (013F7D1440h)  

D:

000000013F761002  mov         ebp,esp  
000000013F761004  sub         rsp,50h  
{
    float eps = 1.0f;
000000013F761008  xor         eax,eax  
000000013F76100A  mov         dword ptr [rbp-50h],eax  
000000013F76100D  movss       xmm0,dword ptr [rbp-50h]  
000000013F761012  movss       dword ptr [f],xmm0  
    float f = 0.0f;
    while (f + eps != f)
        f += 1.0f;
000000013F761017  movss       xmm1,dword ptr [__NULL_IMPORT_DESCRIPTOR+1138h (013F7C3040h)]  
000000013F76101F  movss       xmm2,dword ptr [f]  
000000013F761024  addss       xmm2,xmm1  
000000013F761028  movss       dword ptr [f],xmm2  
000000013F76102D  fld         dword ptr [f]  
000000013F761030  fadd        dword ptr [__NULL_IMPORT_DESCRIPTOR+1138h (013F7C3040h)]  
000000013F761036  fld         dword ptr [f]  
000000013F761039  fucomip     st,st(1)  
000000013F76103B  fstp        st(0)  
000000013F76103D  jne         D main+17h (013F761017h)  
000000013F76103F  jp          D main+17h (013F761017h)  

Summary

Accept harold's answer that program behavior is due to the mixed FPU and SSE usage.

Here's a summary what happens in D assembly snippet. In fact the loop will run forever.

SSE behaves strictly according to IEEE-754 when f reaches 16777216.0 and we add 1.0 to this value (f += 1.0f) we still obtain 16777216.0 in xmm2 register, then we store it to memory.

(f + eps != f) expression is computed on the FPU. Since FPU registers have enough precision (f+eps) results in 16777217.0. If we stored this result back to memory into float variable then we would get expected value 16777216.0 (since 16777217.0 is not represented as float). And (f + eps != f) would be 'false' and loop would terminate. But we do not store any numbers back to memory and perform comparison on the FPU (since we have both operands). It means that we compare one number that is computed strictly according to IEEE-754 (f) and another that is computed with 80bit accuracy (f+eps). 16777216.0 != 16777217.0 and the loop runs forever.

I'm not an expert in this area but for me it looks like that doing floating point with SSE instructions is more robust as was demonstrated in C++ version of the program.

Update

I had a discussion on the D forum http://forum.dlang.org/thread/[email protected]

It turned out that program behaves correctly - it's according to the language specification that intermediate calculations can be performed with higher accuracy.

The robust implementation for any D compiler is:

import std.stdio;
int main()
{
    const float eps = 1.0f;
    const float step = 1.0;

    float f = 0.0f;
    float fPlusEps = f + eps;
    while (f != fPlusEps)
    {
        f += step;
        fPlusEps = f + eps;
    }
    writeln("eps = ", eps, ", max_f = ", f);
    return 0;
}
5
  • This might have something to do with optimization and floating point math compiler options. Try playing around with them? Also try to write the loop in different ways (just to make sure any optimization doesn't break it). Commented Feb 7, 2015 at 10:14
  • Yes, maybe dmd has some different floating point defaults than VC but I performed very simple operations and expected that this will work the same as in C++. Commented Feb 7, 2015 at 10:31
  • Show the disassembly. In particular, is one of them really using float sized floats while the other is using 80bit extended precision? It would take a long time (not forever, but too long) to break out of that loop with 80bit extended precision.
    – user555045
    Commented Feb 7, 2015 at 10:52
  • 1
    Just few more notes. I'm testing here agains IEEE-754 specification. It should work correctly for any implementation. 80 bit conversion it's only impementation details of x86. When we got results back in memory it should be according to specification. I don't think we should think much about such low level details. And yes, I agree we should never compare results FP computations with == or != but in this case I'm testing floating points behavior according to IEEE-754, so I believe it's valid for this case. Commented Feb 7, 2015 at 11:00
  • 2
    Unfortunately compilers often don't work that way, and they'll go ahead and use 80bit precision without the extra steps needed to chop the result back its proper precision. As you say, it should be according to the specification when you get it back in memory, but you probably don't get it back in memory, it probably stays in a register as an optimization. Anyway this is easy to check: see the disassembly.
    – user555045
    Commented Feb 7, 2015 at 11:14

2 Answers 2

6

Mixed FPU and SSE code, that's .. really strange. I see absolutely no reason to implement it this way.

But they have, and the result is that f + eps != f is evaluated with 80bit extended precision, while
f += 1.0f is evaluated using 32bit floats.

That means the loop can never end, since f will stop going up before the value that makes
f + eps != f false (which, in 80bit precision, is huge) is reached.

3
  • Thansk for checking the assembly!. I'll check if I can do something with DMD compiler. Another option is to compile with GDC. Commented Feb 7, 2015 at 11:27
  • I figured out that the loop runs forever - it stucks with the same values for each iteration. I updated the question with more detailed explanation of the assembly snippet. Commented Feb 7, 2015 at 14:32
  • 4
    Compiling with GDC or LDC yields the expected result. As explained in the answer, the cause is the 80-bit precision. The way DMD does it's codegen is rather special, the backend is originally from Digital Mars C, and it was designed to compile code fast, rather than with the maximum final performance.
    – Orvid King
    Commented Feb 7, 2015 at 20:23
1

Trying to break a loop with != or == with floating point values is looking for troubles.

The different behavior is mot likely due to the float to double to 80-bits internal floating point conversion compiler may adopt when passing values to the FPU.

When extending the mantissa, in particular- some compilers or optimizer can decide to let the less significant bit "random" instead of zeroed. So 1.0f, when given to the FPU may become 1.000000000000000000000012134432 that -according to a float- precision, is still 1.0, but wen 1.000000000000000000000012134432 and 1.000000000000000000000089544455 (the two tail are random) are compared by the FPU, look different.

You should verify how C++ and D compiler treat the floating point extension/reduction and eventually configure the appropriate switches: if the two compilers are not from the same manufacturer, thay had probably made different choices for their respective defaults.

18
  • 2
    @harold: To make the code run faster: the FPU use internally 80 bit representaion. If you provide less, and don't reset the tail you need less communication bandwidth between CPU and FPU. It doesn't have much sense in understanding if if has more or less trouble: if you do a sequence of FADD and FMUL, you will in any case, sooner or later, abandon the "exact result arena". It is using == (or !=) the real "trouble". If I have to design a new language nowadays, I'll deny == for floating points. (And write an article titled "== considered harmful") Commented Feb 7, 2015 at 10:55
  • 1
    Floating point works correctly for integer values until some maximum value (and I try to find it in the above code). Up to that limit addition, multiplication, subrraction and exact comparison should work correctly and produce exactly integer values (i.e. zero after floating point). That's how integer arithmetic works in JavaScript without any issues which uses doubles for all numerical calculations. Commented Feb 7, 2015 at 11:06
  • 2
    @EmilioGaravaglia you can't not reset the tail. Instructions that load floating point numbers do that automatically.
    – user555045
    Commented Feb 7, 2015 at 11:23
  • 2
    A compiler that puts garbage in the low bits of a floating-point number is a broken compiler. No compiler has ever done that.
    – tmyklebu
    Commented Feb 7, 2015 at 19:47
  • 3
    @EmilioGaravaglia: Go read the IEEE standard and tell me where it says arithmetic on floating-point numbers is permitted to produce garbage results. Then tell me where the C or C++ standards say a strict implementation of the IEEE standard is required. Then tell me which compilers actually implement the C or C++ standards, especially as they pertain to floating-point math. Inserting garbage into the significand of a floating-point number is a heinous sin that will break a lot of carefully-written code, which is why compilers don't and have never done that.
    – tmyklebu
    Commented Feb 7, 2015 at 22:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.