Search
Very strange problem with comparing floating point numbers
Sep 30, 2012
Ivan Agafonov
Sep 30, 2012
Andrej Mitrovic
Sep 30, 2012
Ali Çehreli
Sep 30, 2012
jerro
Oct 01, 2012
Ivan Agafonov
Oct 01, 2012
jerro
Sep 30, 2012
Tommi
Sep 30, 2012
Ali Çehreli
Sep 30, 2012
Andrej Mitrovic
Sep 30, 2012
monarch_dodra
Sep 30, 2012
so
Sep 30, 2012
monarch_dodra
Sep 30, 2012
jerro
Sep 30, 2012
monarch_dodra
Sep 30, 2012
jerro
Oct 01, 2012
Maxim Fomin
Oct 01, 2012
monarch_dodra
Oct 01, 2012
Maxim Fomin
Oct 01, 2012
monarch_dodra
Oct 02, 2012
Maxim Fomin
```// Tell me about this sutation, may be it is a bug?

import std.math;
import std.stdio;

struct Vector(int size)
{
union
{
float[size] array = 0;
struct
{
static if (size == 2) float x, y;
static if (size == 3) float x, y, z;
static if (size == 4) float x, y, z, w;
}
}

@property float lengthsqr()
{
static if (size == 2) return x*x + y*y;
static if (size >= 3) return x*x + y*y + z*z;
}

@property float length() { return sqrt(lengthsqr()); }
@property float length2() { float tmp = sqrt(lengthsqr()); return tmp; }
}

void main()
{
auto a = Vector!4([1, 2, 3, 1]);
auto a3 = Vector!3([1, 2, 3]);
assert (a.lengthsqr == 14);

auto alen = a.length; auto a3len = a3.length;

// all of this prints the same number: 0x1.deea2p+1
writefln("%a, %a", alen, a3len);
writefln("%a, %a", a.length, a3.length);
writefln("%a, %a", a.length2, a3.length2);

// passes
assert (alen == a3len);
assert (a.length2 == a3.length2);
assert (cast(real)a.length == cast(real)a3.length);

// all of this fails!!!
assert (a.length == a.length); // This is really shocking
assert (a.length == a3.length);
}
```
```On 9/30/12, Ivan Agafonov <armadil@yandex.ru> wrote:

Reduced:

import std.stdio;
import std.math;

@property float getFloat()
{
return sqrt(1.1);
}

void main()
{
writeln(getFloat == getFloat);  // fail
}

Dissasembly:
__Dmain:; Function begin, communal
enter   12, 0                                   ; 0000 _ C8, 000C, 00
call    _D4test8getFloatFNdZf                   ; 0004 _ E8,
00000000(rel)
fstp    dword [ebp-0CH]                         ; 0009 _ D9. 5D, F4
call    _D4test8getFloatFNdZf                   ; 000C _ E8,
00000000(rel)
fld     dword [ebp-0CH]                         ; 0011 _ D9. 45, F4
fxch    st1                                     ; 0014 _ D9. C9
fucompp                                         ; 0016 _ DA. E9
fnstsw  ax                                      ; 0018 _ DF. E0
sahf                                            ; 001A _ 9E
mov     eax, 1                                  ; 001B _ B8, 00000001
jpe     ?_033                                   ; 0020 _ 7A, 02
jz      ?_034                                   ; 0022 _ 74, 02
?_033:  xor     eax, eax                                ; 0024 _ 31. C0
?_034:  call    _D3std5stdio14__T7writelnTbZ7writelnFbZv; 0026 _ E8,
00000000(rel)
xor     eax, eax                                ; 002B _ 31. C0
leave                                           ; 002D _ C9
ret                                             ; 002E _ C3
; __Dmain End of function

Now compare to this which doesn't fail:

void main()
{
float fx1 = getFloat;
float fx2 = getFloat;
writeln(fx1 == fx2);  // pass
}

Dissasembly:
__Dmain:; Function begin, communal
enter   8, 0                                    ; 0000 _ C8, 0008, 00
call    _D4test8getFloatFNdZf                   ; 0004 _ E8,
00000000(rel)
fstp    dword [ebp-8H]                          ; 0009 _ D9. 5D, F8
call    _D4test8getFloatFNdZf                   ; 000C _ E8,
00000000(rel)
fstp    dword [ebp-4H]                          ; 0011 _ D9. 5D, FC
fld     dword [ebp-8H]                          ; 0014 _ D9. 45, F8
fld     dword [ebp-4H]                          ; 0017 _ D9. 45, FC
fucompp                                         ; 001A _ DA. E9
fnstsw  ax                                      ; 001C _ DF. E0
sahf                                            ; 001E _ 9E
mov     eax, 1                                  ; 001F _ B8, 00000001
jpe     ?_033                                   ; 0024 _ 7A, 02
jz      ?_034                                   ; 0026 _ 74, 02
?_033:  xor     eax, eax                                ; 0028 _ 31. C0
?_034:  call    _D3std5stdio14__T7writelnTbZ7writelnFbZv; 002A _ E8,
00000000(rel)
xor     eax, eax                                ; 002F _ 31. C0
leave                                           ; 0031 _ C9
ret                                             ; 0032 _ C3
; __Dmain End of function

In the first sample the ASM uses fstp for LHS and fld for RHS (http://en.wikipedia.org/wiki/X86_instruction_listings).

The second one uses fstp twice, then fld twice. I don't know, maybe this could be a bug. My ASM is weak but I see some difference here..
```
```On 09/29/2012 06:48 PM, Andrej Mitrovic wrote:

> The second one uses fstp twice, then fld twice. I don't know, maybe
> this could be a bug. My ASM is weak but I see some difference here..

I saw similar differences but my asm is even weaker. :)

I am writing to confirm that I can reproduce the problem with -m32 on a 64-bit system. The code works as expected without -m32.

Ali

```
```> The second one uses fstp twice, then fld twice. I don't know, maybe
> this could be a bug.

You're right the lack of one fst/fld in the first case is a bug. x87 floating point registers are 80 bit. This:

fstp    dword [ebp-0CH]

Converts the value in ST0 to single precision float and stores it to memory (and pops ST0). When it is later loaded with fld, it is not the same as before storing since some precision is lost (because the D code compares floats and not reals, this is the correct behavior). In the first example, this storing and loading only happens for the first function call. For the second call the value is returned in ST0 and stays in x87 registers until it is compared with fucompp so it is not truncated as the result of the first function call was. That's why the compared values are not equal.

```
```On Sunday, 30 September 2012 at 01:48:04 UTC, Andrej Mitrovic wrote:
>
> Dissasembly:
> __Dmain:; Function begin, communal
>         enter   12, 0                                   ; 0000 _ C8, 000C, 00
>         call    _D4test8getFloatFNdZf                   ; 0004 _ E8,
> 00000000(rel)
> ...

Can I tell DMD to produce the assembly, or what did you do to get that?
```
```On 09/30/2012 04:06 AM, Tommi wrote:
> On Sunday, 30 September 2012 at 01:48:04 UTC, Andrej Mitrovic wrote:
>>
>> Dissasembly:
>> __Dmain:; Function begin, communal
>> enter 12, 0 ; 0000 _ C8, 000C, 00
>> call _D4test8getFloatFNdZf ; 0004 _ E8,
>> 00000000(rel)
>> ...
>
> Can I tell DMD to produce the assembly, or what did you do to get that?

For a foo.d, after compiling the program and generating foo.o, the two options on Linux that I know of:

1) obj2asm that comes with dmd:

\$ obj2asm  foo.o > foo.asm

2) objdump that comes with at least my Linux distribution:

\$ objdump -d foo.o > foo.asm

Ali
```
```On Sunday, 30 September 2012 at 01:29:24 UTC, Ivan Agafonov wrote:
>
> [SNIP]
> 	// all of this fails!!!
> 	assert (a.length == a.length); // This is really shocking
> 	assert (a.length == a3.length);
> [SNIP]

This is just a fact of life regarding how floating point types work. Here is a well documented explanation. It pertains to C++, but applies.

http://www.parashift.com/c++-faq/floating-point-arith2.html

As a rule of thumb, NEVER use opEqual with floating point types aniways. You need to use some sort of comparison with leway for error, such as std.math.approxEqual.

```
```On Sunday, 30 September 2012 at 17:07:19 UTC, monarch_dodra wrote:
> On Sunday, 30 September 2012 at 01:29:24 UTC, Ivan Agafonov wrote:
>>
>> [SNIP]
>> 	// all of this fails!!!
>> 	assert (a.length == a.length); // This is really shocking
>> 	assert (a.length == a3.length);
>> [SNIP]
>
> This is just a fact of life regarding how floating point types work. Here is a well documented explanation. It pertains to C++, but applies.
>
> http://www.parashift.com/c++-faq/floating-point-arith2.html
>
> As a rule of thumb, NEVER use opEqual with floating point types aniways. You need to use some sort of comparison with leway for error, such as std.math.approxEqual.

Floating point types are trouble enough without these
optimization failures.
There are many unsolved problems, things like approxEqual are far
from answering them. Whatever the justifications they come up
with, "a.len == a.len" failure is IMO unacceptable, an opEqual
like this must not fail.

A suggestion: do what i do and have this in your config files.

alias real evil;
```
```On Sunday, 30 September 2012 at 18:31:17 UTC, so wrote:
> On Sunday, 30 September 2012 at 17:07:19 UTC, monarch_dodra wrote:
>> On Sunday, 30 September 2012 at 01:29:24 UTC, Ivan Agafonov wrote:
>>>
>>> [SNIP]
>>> 	// all of this fails!!!
>>> 	assert (a.length == a.length); // This is really shocking
>>> 	assert (a.length == a3.length);
>>> [SNIP]
>>
>> This is just a fact of life regarding how floating point types work. Here is a well documented explanation. It pertains to C++, but applies.
>>
>> http://www.parashift.com/c++-faq/floating-point-arith2.html
>>
>> As a rule of thumb, NEVER use opEqual with floating point types aniways. You need to use some sort of comparison with leway for error, such as std.math.approxEqual.
>
> Floating point types are trouble enough without these
> optimization failures.
> There are many unsolved problems, things like approxEqual are far
> from answering them. Whatever the justifications they come up
> with, "a.len == a.len" failure is IMO unacceptable, an opEqual
> like this must not fail.
>
> A suggestion: do what i do and have this in your config files.
>
> alias real evil;

I don't really agree with that. floating point operations are just inexact, regardless of optimizations. That's how they work, period.

Either you can work with inexact results, and you use them, or you can't, and don't. Banks don't use floating point types for exactly this reason. You have to know what you are getting into before you begin.

The real troubles really only start when you start using floating point type, but you expect exact results.

```
```On 9/30/12, Tommi <tommitissari@hotmail.com> wrote:
> On Sunday, 30 September 2012 at 01:48:04 UTC, Andrej Mitrovic wrote:
>>
>> Dissasembly:
>> __Dmain:; Function begin, communal
>>         enter   12, 0                                   ; 0000
>> _ C8, 000C, 00
>>         call    _D4test8getFloatFNdZf                   ; 0004
>> _ E8,
>> 00000000(rel)
>> ...
>
> Can I tell DMD to produce the assembly, or what did you do to get that?
>

I use objconv (can run on win32 only methinks) on an .obj file on
win32 via a batch file:
objconv -fnasm %~nx1 %~n1_disasm.asm  && %~n1_disasm.asm
http://www.agner.org/optimize/objconv.zip
```
« First   ‹ Prev
1 2