Different array access methods cause different floating point calculation results

Different array access methods cause different floating point calculation results

[Florian Hackenberger]

Hi!

Maybe this is considered normal, but I experience odd behaviour, when 
calculating X*X + Y*Y using direct array access compared to array access via 
functions returning references to array elements. Here is what I'm trying to 
do (reduced to the bare minimum to show the bug):

const double& getState(unsigned int index) {
                return y[index];
}

double y[2];
int main(int, char**) {
        y[0] = 1.0000002232024669535093153172056190669536590576171875;
        y[1] = 
6.24999999999999945247790289482026082623633556067943572998046875e-07;

        double temp1 = getState(0)*getState(0) + getState(1)*getState(1);
        double temp2 = y[0]*y[0] + y[1]*y[1];

        if(temp2 != temp1) {
                return -1;
        }
	return 0;
}

When compiling with:
g++ -g -O0 -o test test.cpp
and running it, I obtain the following results:
temp1 = 1.0000004464053742214701969714951701462268829345703125
temp2 = 1.000000446405374443514801896526478230953216552734375

Of course I know about floating point rounding errors, but when calculating an 
expression twice using the very same input values, I would expect the result 
to be the same.

The results vary, depending on the optimisation I choose for compiling:
-O0:
1.0000004464053742214701969714951701462268829345703125
1.000000446405374443514801896526478230953216552734375
-O1:
1.000000446405374443514801896526478230953216552734375
1.0000004464053742214701969714951701462268829345703125
-O2:
1.000000446405374443514801896526478230953216552734375
1.0000004464053742214701969714951701462268829345703125
-O3:
1.0000004464053742214701969714951701462268829345703125
1.0000004464053742214701969714951701462268829345703125

When using -O0 the compiler produces the following assembler code:
    0x08048618 <main+0>:    push   %ebp
    0x08048619 <main+1>:    mov    %esp,%ebp
    0x0804861b <main+3>:    sub    $0x18,%esp
    0x0804861e <main+6>:    and    $0xfffffff0,%esp
    0x08048621 <main+9>:    sub    $0x10,%esp
    0x08048624 <main+12>:   movl   $0x3bea5b53,0x8049a00
    0x0804862e <main+22>:   movl   $0x3ff00000,0x8049a04
    0x08048638 <main+32>:   movl   $0x88e368f0,0x8049a08
    0x08048642 <main+42>:   movl   $0x3ea4f8b5,0x8049a0c
    0x0804864c <main+52>:   movl   $0x64,0x8049978
    0x08048656 <main+62>:   movl   $0x77d4be65,0x4(%esp)
    0x0804865e <main+70>:   movl   $0x3ff00000,0x8(%esp)
    0x08048666 <main+78>:   movl   $0x8049970,(%esp)
    0x0804866d <main+85>:   call   0x80484b8 <_ZNSolsEd@plt>
    0x08048672 <main+90>:   mov    %eax,(%esp)
    0x08048675 <main+93>:   call   0x80484d8 
<_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_@plt>
    0x0804867a <main+98>:   movl   $0x77d4be65,0x4(%esp)
    0x08048682 <main+106>:  movl   $0x3ff00000,0x8(%esp)
    0x0804868a <main+114>:  movl   $0x8049970,(%esp)
    0x08048691 <main+121>:  call   0x80484b8 <_ZNSolsEd@plt>
    0x08048696 <main+126>:  mov    %eax,(%esp)
    0x08048699 <main+129>:  call   0x80484d8 
<_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_@plt>
    0x0804869e <main+134>:  movl   $0x77d4be65,0x4(%esp)
    0x080486a6 <main+142>:  movl   $0x3ff00000,0x8(%esp)
    0x080486ae <main+150>:  movl   $0x8049970,(%esp)
    0x080486b5 <main+157>:  call   0x80484b8 <_ZNSolsEd@plt>
    0x080486ba <main+162>:  mov    %eax,(%esp)
    0x080486bd <main+165>:  call   0x80484d8 
<_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_@plt>
    0x080486c2 <main+170>:  xor    %eax,%eax
    0x080486c4 <main+172>:  leave  
    0x080486c5 <main+173>:  ret    

I would consider this behaviour a compiler bug, but maybe someone can explain 
it, before I file a bug report. BTW I have tried g++ (GCC) 4.0.3 (Ubuntu 
4.0.3-1ubuntu5) and g++ (GCC) 4.0.2 20051125 (Red Hat 4.0.2-8), the results 
are the same.

Regards,
	Florian

-- 
Florian Hackenberger
student @
University of Technology
Graz, Austria
florian@hackenberger.at
www.hackenberger.at

Re: Different array access methods cause different floating point calculation results

[Brian Budge]

Without having time to analyze your program in depth, my guess would
be that in one case temporaries are stored in memory, and the other
they are sticking around in registers.  When they are stored to
memory, some precision can be lost because, for example, the x87 fpu
has 80 bits of precision, but when it is stored to memory, it is
stored as a 64 bit double.

If you're running on an x86 platform with sse, you might try the -msee
-mfpmath=sse options when compiling, which theoretically should use
the 64 bit sse registers instead of the fpu stack.

  Brian

On 9/7/06, Florian Hackenberger <f.hackenberger@chello.at> wrote:
> Hi!
>
> Maybe this is considered normal, but I experience odd behaviour, when
> calculating X*X + Y*Y using direct array access compared to array access via
> functions returning references to array elements. Here is what I'm trying to
> do (reduced to the bare minimum to show the bug):
>
> const double& getState(unsigned int index) {
>                 return y[index];
> }
>
> double y[2];
> int main(int, char**) {
>         y[0] = 1.0000002232024669535093153172056190669536590576171875;
>         y[1] =
> 6.24999999999999945247790289482026082623633556067943572998046875e-07;
>
>         double temp1 = getState(0)*getState(0) + getState(1)*getState(1);
>         double temp2 = y[0]*y[0] + y[1]*y[1];
>
>         if(temp2 != temp1) {
>                 return -1;
>         }
>         return 0;
> }
>
> When compiling with:
> g++ -g -O0 -o test test.cpp
> and running it, I obtain the following results:
> temp1 = 1.0000004464053742214701969714951701462268829345703125
> temp2 = 1.000000446405374443514801896526478230953216552734375
>
> Of course I know about floating point rounding errors, but when calculating an
> expression twice using the very same input values, I would expect the result
> to be the same.
>
> The results vary, depending on the optimisation I choose for compiling:
> -O0:
> 1.0000004464053742214701969714951701462268829345703125
> 1.000000446405374443514801896526478230953216552734375
> -O1:
> 1.000000446405374443514801896526478230953216552734375
> 1.0000004464053742214701969714951701462268829345703125
> -O2:
> 1.000000446405374443514801896526478230953216552734375
> 1.0000004464053742214701969714951701462268829345703125
> -O3:
> 1.0000004464053742214701969714951701462268829345703125
> 1.0000004464053742214701969714951701462268829345703125
>
> When using -O0 the compiler produces the following assembler code:
>     0x08048618 <main+0>:    push   %ebp
>     0x08048619 <main+1>:    mov    %esp,%ebp
>     0x0804861b <main+3>:    sub    $0x18,%esp
>     0x0804861e <main+6>:    and    $0xfffffff0,%esp
>     0x08048621 <main+9>:    sub    $0x10,%esp
>     0x08048624 <main+12>:   movl   $0x3bea5b53,0x8049a00
>     0x0804862e <main+22>:   movl   $0x3ff00000,0x8049a04
>     0x08048638 <main+32>:   movl   $0x88e368f0,0x8049a08
>     0x08048642 <main+42>:   movl   $0x3ea4f8b5,0x8049a0c
>     0x0804864c <main+52>:   movl   $0x64,0x8049978
>     0x08048656 <main+62>:   movl   $0x77d4be65,0x4(%esp)
>     0x0804865e <main+70>:   movl   $0x3ff00000,0x8(%esp)
>     0x08048666 <main+78>:   movl   $0x8049970,(%esp)
>     0x0804866d <main+85>:   call   0x80484b8 <_ZNSolsEd@plt>
>     0x08048672 <main+90>:   mov    %eax,(%esp)
>     0x08048675 <main+93>:   call   0x80484d8
> <_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_@plt>
>     0x0804867a <main+98>:   movl   $0x77d4be65,0x4(%esp)
>     0x08048682 <main+106>:  movl   $0x3ff00000,0x8(%esp)
>     0x0804868a <main+114>:  movl   $0x8049970,(%esp)
>     0x08048691 <main+121>:  call   0x80484b8 <_ZNSolsEd@plt>
>     0x08048696 <main+126>:  mov    %eax,(%esp)
>     0x08048699 <main+129>:  call   0x80484d8
> <_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_@plt>
>     0x0804869e <main+134>:  movl   $0x77d4be65,0x4(%esp)
>     0x080486a6 <main+142>:  movl   $0x3ff00000,0x8(%esp)
>     0x080486ae <main+150>:  movl   $0x8049970,(%esp)
>     0x080486b5 <main+157>:  call   0x80484b8 <_ZNSolsEd@plt>
>     0x080486ba <main+162>:  mov    %eax,(%esp)
>     0x080486bd <main+165>:  call   0x80484d8
> <_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_@plt>
>     0x080486c2 <main+170>:  xor    %eax,%eax
>     0x080486c4 <main+172>:  leave
>     0x080486c5 <main+173>:  ret
>
> I would consider this behaviour a compiler bug, but maybe someone can explain
> it, before I file a bug report. BTW I have tried g++ (GCC) 4.0.3 (Ubuntu
> 4.0.3-1ubuntu5) and g++ (GCC) 4.0.2 20051125 (Red Hat 4.0.2-8), the results
> are the same.
>
> Regards,
>         Florian
>
> --
> Florian Hackenberger
> student @
> University of Technology
> Graz, Austria
> florian@hackenberger.at
> www.hackenberger.at
>

Re: Different array access methods cause different floating point calculation results

[Andrew Haley]

Florian Hackenberger writes:
 > Hi!
 > 
 > Maybe this is considered normal,

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323

Andrew.

Re: Different array access methods cause different floating point calculation results

[Florian Hackenberger]

On Thursday 07 September 2006 17:28, Andrew Haley wrote:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
I see, please forgive my ignorance and thank you all for the explanations / 
pointers to explanations!

Regards,
	Florian
-- 
Florian Hackenberger
student @
University of Technology
Graz, Austria
florian@hackenberger.at
www.hackenberger.at

Re: Different array access methods cause different floating point calculation results

[Florian Hackenberger]

On Thursday 07 September 2006 17:28, Andrew Haley wrote:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
Maybe somebody could add it to the GCC FAQ?

Regards,
	Florian
-- 
Florian Hackenberger
student @
University of Technology
Graz, Austria
florian@hackenberger.at
www.hackenberger.at

Re: Different array access methods cause different floating point calculation results

[Andrew Haley]

Florian Hackenberger writes:
 > On Thursday 07 September 2006 17:28, Andrew Haley wrote:
 > > http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
 > Maybe somebody could add it to the GCC FAQ?

I can do that.  Do you really think it would help?

Andrew.

Re: Different array access methods cause different floating point calculation results

[Florian Hackenberger]

On Friday 08 September 2006 11:19, Andrew Haley wrote:
>  > Maybe somebody could add it to the GCC FAQ?
> I can do that.  Do you really think it would help?
Yes indeed, at least for me it would have helped. I read the FAQ before 
posting here. Although I'm not sure if all the other users out there are 
reading FAQ's, it would at least catch a few cases.

Regards,
	Florian

-- 
Florian Hackenberger
student @
University of Technology
Graz, Austria
florian@hackenberger.at
www.hackenberger.at