digitalmars.D.learn - Very strange problem with comparing floating point numbers
- Ivan Agafonov (47/47) Sep 29 2012 // Tell me about this sutation, may be it is a bug?
- Andrej Mitrovic (70/71) Sep 29 2012 Reduced:
- =?UTF-8?B?QWxpIMOHZWhyZWxp?= (5/7) Sep 29 2012 I saw similar differences but my asm is even weaker. :)
- jerro (13/16) Sep 29 2012 You're right the lack of one fst/fld in the first case is a bug.
- Ivan Agafonov (2/18) Sep 30 2012 Can you or anyone report this bug? I don't know how to do this.
- jerro (11/32) Sep 30 2012 I'm not really sure that this is a bug anymore. Apparently c++
- Tommi (4/12) Sep 30 2012 Can I tell DMD to produce the assembly, or what did you do to get
- =?UTF-8?B?QWxpIMOHZWhyZWxp?= (8/17) Sep 30 2012 For a foo.d, after compiling the program and generating foo.o, the two
- Andrej Mitrovic (5/18) Sep 30 2012 I use objconv (can run on win32 only methinks) on an .obj file on
- monarch_dodra (8/14) Sep 30 2012 This is just a fact of life regarding how floating point types
- so (9/25) Sep 30 2012 Floating point types are trouble enough without these
- monarch_dodra (10/40) Sep 30 2012 I don't really agree with that. floating point operations are
- jerro (7/10) Sep 30 2012 It is true that they are inexact, but inexact and
- monarch_dodra (5/11) Sep 30 2012 Technically (AFAIK), IEEE754 does need require reproducibility,
- jerro (14/16) Sep 30 2012 It actually requires more than that:
- Maxim Fomin (13/16) Oct 01 2012 It is possible to compare exactly floating point types by binary
- monarch_dodra (6/23) Oct 01 2012 I think that what you are comparing here is the functions (the
- Maxim Fomin (59/88) Oct 01 2012 http://dpaste.dzfl.pl/1f94c0b1
- monarch_dodra (12/66) Oct 01 2012 Hum, yes, I guess I was wrong about the comparison of functions.
- Maxim Fomin (5/15) Oct 01 2012 It looks like dmd uses x87 comparison instructions which are
// 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);
}
Sep 29 2012
On 9/30/12, Ivan Agafonov <armadil yandex.ru> wrote:// Tell me about this sutation, may be it is a bug?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..
Sep 29 2012
=?UTF-8?B?QWxpIMOHZWhyZWxp?= <acehreli yahoo.com> 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
Sep 29 2012
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.
Sep 29 2012
On Sunday, 30 September 2012 at 06:20:56 UTC, jerro wrote:Can you or anyone report this bug? I don't know how to do this.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.
Sep 30 2012
On Monday, 1 October 2012 at 04:10:25 UTC, Ivan Agafonov wrote:On Sunday, 30 September 2012 at 06:20:56 UTC, jerro wrote:I'm not really sure that this is a bug anymore. Apparently c++ does this like this too (http://www.parashift.com/c++-faq/floating-point-arith2.html). I do think it would be more useful if the result of floating point comparison would be defined in a cases like this, but for floating point operation in general it makes sense to keep intermediate results in registers. So I don't know whether introducing a special case for this is worth it. If you care about this, you should probably make a thread on digitalmars.D, where more people will see it and comment on it.Can you or anyone report this bug? I don't know how to do this.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.
Sep 30 2012
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?
Sep 30 2012
On 09/30/2012 04:06 AM, Tommi wrote:On Sunday, 30 September 2012 at 01:48:04 UTC, Andrej Mitrovic wrote: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 AliDissasembly: __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?
Sep 30 2012
On 9/30/12, Tommi <tommitissari hotmail.com> wrote:On Sunday, 30 September 2012 at 01:48:04 UTC, Andrej Mitrovic wrote: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.zipDissasembly: __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?
Sep 30 2012
On Sunday, 30 September 2012 at 01:29:24 UTC, Ivan Agafonov wrote:// Tell me about this sutation, may be it is a bug? [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.
Sep 30 2012
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: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;// Tell me about this sutation, may be it is a bug? [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.
Sep 30 2012
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: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 Sunday, 30 September 2012 at 01:29:24 UTC, Ivan Agafonov wrote: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;// Tell me about this sutation, may be it is a bug? [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.
Sep 30 2012
I don't really agree with that. floating point operations are just inexact, regardless of optimizations. That's how they work, period.It is true that they are inexact, but inexact and non-deterministic are not the same thing. Floating point operations are deterministic. Doing the same operation twice on the same machine with the same rounding mode and parameters will give you the same result. This should also be true when you do those operations using D, and using == on the two results should return true.
Sep 30 2012
On Sunday, 30 September 2012 at 20:47:41 UTC, jerro wrote:Technically (AFAIK), IEEE754 does need require reproducibility, ergo determinism. You can open an ER requesting the ability to specify the FP behavior want (like rounding behavior), the way vstudio does it.I don't really agree with that. floating point operations are just inexact, regardless of optimizations. That's how they work, period.It is true that they are inexact, but inexact and non-deterministic are not the same thing. Floating point operations are deterministic.
Sep 30 2012
Technically (AFAIK), IEEE754 does need require reproducibility, ergo determinism.It actually requires more than that: "Algebraic operations covered by IEEE 754, namely + , - , * , / , square root and Binary <-> Decimal Conversion with rare exceptions, must be Correctly Rounded to the precision of the operationâs destination unless the programmer has specified a rounding other than the default. If it does not Overflow, a correctly rounded operationâs error cannot exceed half the gap between adjacent floating-point numbers astride the operationâs ideal ( unrounded ) result. Half-way cases are rounded to Nearest Even, which means that the neighbor with last digit 0 is chosen." I don't know if implementation conform exactly to IEEE 754, but I doubt there is any commonly used implementation that isn't deterministic.
Sep 30 2012
On Sunday, 30 September 2012 at 17:07:19 UTC, monarch_dodra wrote: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.It is possible to compare exactly floating point types by binary comparison, if it provides some benefits. import std.stdio; import std.math; property float getFloat() { return sqrt(1.1); } void main() { writeln(getFloat is getFloat); // doesn't fail }
Oct 01 2012
On Monday, 1 October 2012 at 11:36:43 UTC, Maxim Fomin wrote:On Sunday, 30 September 2012 at 17:07:19 UTC, monarch_dodra wrote:I think that what you are comparing here is the functions (the address), and not the results of the call. Try writeln(getFloat() is getFloat()); //*May* fail Also, "is" works like opEqual on built in types, AFAIK, it doesn't use any "binary" magic or anything like that.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.It is possible to compare exactly floating point types by binary comparison, if it provides some benefits. import std.stdio; import std.math; property float getFloat() { return sqrt(1.1); } void main() { writeln(getFloat is getFloat); // doesn't fail }
Oct 01 2012
2012/10/1 monarch_dodra <monarchdodra gmail.com>:On Monday, 1 October 2012 at 11:36:43 UTC, Maxim Fomin wrote:http://dpaste.dzfl.pl/1f94c0b1 It works with -m32 too. _Dmain: 0x0806d0e4 <+0>: push %ebp 0x0806d0e5 <+1>: mov %esp,%ebp 0x0806d0e7 <+3>: sub $0x10,%esp 0x0806d0ea <+6>: push %esi 0x0806d0eb <+7>: push %edi 0x0806d0ec <+8>: call 0x806d0d4 <_D4test8getFloatFNdZf> 0x0806d0f1 <+13>: fstps -0x10(%ebp) 0x0806d0f4 <+16>: lea -0x10(%ebp),%esi 0x0806d0f7 <+19>: call 0x806d0d4 <_D4test8getFloatFNdZf> 0x0806d0fc <+24>: fstps -0xc(%ebp) 0x0806d0ff <+27>: lea -0xc(%ebp),%edi 0x0806d102 <+30>: mov $0x4,%ecx 0x0806d107 <+35>: xor %eax,%eax 0x0806d109 <+37>: repz cmpsb %es:(%edi),%ds:(%esi) 0x0806d10b <+39>: je 0x806d112 <_Dmain+46> 0x0806d10d <+41>: sbb %eax,%eax 0x0806d10f <+43>: sbb $0xffffffff,%eax 0x0806d112 <+46>: neg %eax 0x0806d114 <+48>: sbb %eax,%eax 0x0806d116 <+50>: inc %eax 0x0806d117 <+51>: call 0x806d164 <_D3std5stdio14__T7writelnTbZ7writelnFbZv> 0x0806d11c <+56>: call 0x806d0d4 <_D4test8getFloatFNdZf> 0x0806d121 <+61>: fstps -0x8(%ebp) 0x0806d124 <+64>: lea -0x8(%ebp),%esi 0x0806d127 <+67>: call 0x806d0d4 <_D4test8getFloatFNdZf> 0x0806d12c <+72>: fstps -0x4(%ebp) 0x0806d12f <+75>: lea -0x4(%ebp),%edi 0x0806d132 <+78>: mov $0x4,%ecx 0x0806d137 <+83>: xor %eax,%eax 0x0806d139 <+85>: repz cmpsb %es:(%edi),%ds:(%esi) 0x0806d13b <+87>: je 0x806d142 <_Dmain+94> 0x0806d13d <+89>: sbb %eax,%eax 0x0806d13f <+91>: sbb $0xffffffff,%eax 0x0806d142 <+94>: neg %eax 0x0806d144 <+96>: sbb %eax,%eax 0x0806d146 <+98>: inc %eax 0x0806d147 <+99>: call 0x806d164 <_D3std5stdio14__T7writelnTbZ7writelnFbZv> 0x0806d14c <+104>: call 0x806d0d4 <_D4test8getFloatFNdZf> 0x0806d151 <+109>: sub $0x4,%esp 0x0806d154 <+112>: fstps (%esp) 0x0806d157 <+115>: call 0x806d588 <_D3std5stdio14__T7writelnTfZ7writelnFfZv> 0x0806d15c <+120>: xor %eax,%eax 0x0806d15e <+122>: pop %edi 0x0806d15f <+123>: pop %esi 0x0806d160 <+124>: leave 0x0806d161 <+125>: retOn Sunday, 30 September 2012 at 17:07:19 UTC, monarch_dodra wrote:I think that what you are comparing here is the functions (the address), and not the results of the call. Try writeln(getFloat() is getFloat()); //*May* failAs 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.It is possible to compare exactly floating point types by binary comparison, if it provides some benefits. import std.stdio; import std.math; property float getFloat() { return sqrt(1.1); } void main() { writeln(getFloat is getFloat); // doesn't fail }Also, "is" works like opEqual on built in types, AFAIK, it doesn't use any "binary" magic or anything like that.I don't understand what you are trying to say. Is operator at runtime compares two objects without calling opEquals functions (if applied on user-defined types). For built-in and derived types it is similar to == operator. Although, I am suprised that TDPL and spec doesn't mention it (focused only on CT usage), there is a paragraph (http://ddili.org/ders/d.en/null_is.html) from Turkish D book which clearly shows such usage - so, I think this a valid D feature. Object comparison at low-level (repz cmpsb) means binary comparison.
Oct 01 2012
On Monday, 1 October 2012 at 13:08:07 UTC, Maxim Fomin wrote:2012/10/1 monarch_dodra <monarchdodra gmail.com>:Hum, yes, I guess I was wrong about the comparison of functions. Sorry!On Monday, 1 October 2012 at 11:36:43 UTC, Maxim Fomin wrote:http://dpaste.dzfl.pl/1f94c0b1 [SNIP]On Sunday, 30 September 2012 at 17:07:19 UTC, monarch_dodra wrote:I think that what you are comparing here is the functions (the address), and not the results of the call. Try writeln(getFloat() is getFloat()); //*May* failAs 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.It is possible to compare exactly floating point types by binary comparison, if it provides some benefits. import std.stdio; import std.math; property float getFloat() { return sqrt(1.1); } void main() { writeln(getFloat is getFloat); // doesn't fail }What I was saying is that for built in types such a floats, "is" is (should be) no different from "==". But you catch something interesting: the fact that it provides different results is (IMO), a bug. Looking at it, I'd say the bug is probably that "==" is overly sensitive to extended precision. I've filed a BR: http://d.puremagic.com/issues/show_bug.cgi?id=8745 Please feel free to add anything to it. We'll see if Walter will react to it for a more definite answer.Also, "is" works like opEqual on built in types, AFAIK, it doesn't use any "binary" magic or anything like that.I don't understand what you are trying to say. Is operator at runtime compares two objects without calling opEquals functions (if applied on user-defined types). For built-in and derived types it is similar to == operator. Although, I am suprised that TDPL and spec doesn't mention it (focused only on CT usage), there is a paragraph (http://ddili.org/ders/d.en/null_is.html) from Turkish D book which clearly shows such usage - so, I think this a valid D feature. Object comparison at low-level (repz cmpsb) means binary comparison.
Oct 01 2012
On Monday, 1 October 2012 at 21:23:31 UTC, monarch_dodra wrote:What I was saying is that for built in types such a floats, "is" is (should be) no different from "==". But you catch something interesting: the fact that it provides different results is (IMO), a bug. Looking at it, I'd say the bug is probably that "==" is overly sensitive to extended precision. I've filed a BR: http://d.puremagic.com/issues/show_bug.cgi?id=8745 Please feel free to add anything to it. We'll see if Walter will react to it for a more definite answer.It looks like dmd uses x87 comparison instructions which are inexact comparing to is. So, similarity of is and == operators on built-in and user-defined types may be subject to float/double/real exception.
Oct 01 2012