Different array access methods cause different floating point calculation results

[Florian Hackenberger <f dot hackenberger at chello dot at>]

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