g++ Bug? Why w/ and w/o "O3" options gives different result(g++3.3.3) for this program?

[Peng Yu]

I used the following command lines to compile the progrogram below.
However, the output of the exe files is different. The w/o O3 exe
produced the corrected result.

Is this indicate that there is a bug with g++3.3.3's O3 option? Or I
shouldn't use O3 in this case?

I'm very appreciated that if you will help me.

Thanks,
Peng

g++-3.3 -c  -O3  -o main.o main.cc
g++-3.3   main.o   -o main
result:
temp = (0.0134963,0)
temp = (0.0132417,0)

g++-3.3 -c   -o main.o main.cc
g++-3.3   main.o   -o main
result:
temp = (0.0140056,0)
temp = (0.0140056,0)


The program:
/////////////////////////////////////////////
#include <cmath>
#include <iostream>
#include <iterator>
#include <complex>
#include <algorithm>
#include <iomanip>

#define NA       0.5
#define s        0.5 //for coherent light this is 0.0
#define INTEGRAL_PIECES 25

#define N        18//34 //should be greater than (1+s)*2 / Delta_f

#define L_x      5//10 //in the unit of \lambda\mathit{NA}
#define L_y      5//10 //in the unit of \lambda\mathit{NA}
#define Delta_f  (1./L_x)
#define Delta_g  (1./L_y)

#define Delta_z  0.


using namespace std;

std::complex<double> scrK(double f, double g, double delta_z){
  if((f * f + g * g) < 1.)
    return exp(std::complex<double>(0, M_PI * delta_z * NA * (f * f +
g * g)));
  else if((f * f + g * g) == 1.)
    return exp(std::complex<double>(0, M_PI * delta_z * NA * (f * f +
g * g))) / 2.;
  else
    return 0;
}

double J_minus_0(double f, double g){
  if((f * f + g * g) < s * s )
    return 1 / (M_PI * s * s);
  else if((f * f + g * g) == s * s )
    return 1 / (2 * M_PI * s * s);
  else
    return 0;
}

std::complex<double> compute_scrT(double f_prime, double g_prime,
double f_double_prime, double g_double_prime, double delta_z){
  double f_range_left[3];
  double f_range_right[3];

  double g_range_down[3];
  double g_range_up[3];

  f_range_left[0] = - f_prime        - 1;
  f_range_left[1] = - f_double_prime - 1;
  f_range_left[2] =                  - s;

  f_range_right[0] = - f_prime        + 1;
  f_range_right[1] = - f_double_prime + 1;
  f_range_right[2] =                  + s;

  g_range_down[0] = - g_prime        - 1;
  g_range_down[1] = - g_double_prime - 1;
  g_range_down[2] =                  - s;

  g_range_up[0] = - g_prime        + 1;
  g_range_up[1] = - g_double_prime + 1;
  g_range_up[2] =                  + s;

  std::sort(f_range_left, f_range_left + 3);
  std::sort(f_range_right, f_range_right + 3);
  std::sort(g_range_down, g_range_down + 3);
  std::sort(g_range_up, g_range_up + 3);

  if((f_range_left[2] > f_range_right[0]) || (g_range_down[2] >
g_range_up[0]))
    return 0;
  else{
    double delta_freq = double(s) / double(INTEGRAL_PIECES);
    int f_index_left  = int(ceil(f_range_left[2] / delta_freq));
    int f_index_right = int(floor(f_range_right[2] / delta_freq));
    int g_index_down  = int(ceil(g_range_down[2] / delta_freq));
    int g_index_up    = int(floor(g_range_up[2] / delta_freq));


    double f, g;
    std::complex<double> result = 0;
    for(int i = f_index_left; i <= f_index_right; i ++){
      for(int j = g_index_down; j <= g_index_up; j ++){
        f = i * delta_freq;
        g = j * delta_freq;
        result += J_minus_0(f ,g) * scrK(f + f_prime, g + g_prime,
delta_z) * conj(scrK(f + f_double_prime, g + g_double_prime,
delta_z));
      }
    }
    return result * delta_freq * delta_freq;
  }
}

int main(int argc, char *argv[])
{
        int i1, i2, i3, i4;
        i1 = 6, i2 = 7, i3 = 8, i4 = 16;
        std::complex<double> temp = compute_scrT((i1 - N/2)*Delta_f,
(i2 - N/2)*Delta_g, (i3 - N/2)*Delta_f, (i4 - N/2)*Delta_g, 0);
        std::cout << "temp = " << temp << std::endl;
        i1 = 6, i2 = 7, i3 = 8, i4 = 16;
        temp = compute_scrT((i3 - N/2)*Delta_f, (i4 - N/2)*Delta_g,
(i1 - N/2)*Delta_f, (i2 - N/2)*Delta_g, 0);
        std::cout << "temp = " << temp << std::endl;
}