Why does this code fail when an int is used but not when a double is used (possible compiler bug)?

Question

Given the following very simple program in GCC 7.3.1:

//exponentTest.cc
#include <complex> 
#include <iostream>        

int main(int argc, char **argv)
{
  std::complex<double> iCpx = std::complex<double>(0,1); 
  double baseline[3] = {-0.1, 0, 0};
  
  double theta = atan2(baseline[1], baseline[0]); 
  std::cout << "Theta: " << theta << std::endl;

  for(size_t m =1; m < 10; m++)
    {
      std::complex<double> thetExp = exp(-1 * m * theta * iCpx); 
      std::cout << "m(" << m << "): " << thetExp << std::endl;
    }
  return 0; 
}

When theta is PI, this code gives an incorrect result:

[stix@localhost ~]$ gcc exponentTest.cc -o expTest -lm -lstdc++
[stix@localhost ~]$ ./expTest
Theta: 3.14159
m(1): (-0.963907,0.26624)
m(2): (-0.963907,0.26624)
m(3): (-0.963907,0.26624)
m(4): (-0.963907,0.26624)
m(5): (-0.963907,0.26624)
m(6): (-0.963907,0.26624)
m(7): (-0.963907,0.26624)
m(8): (-0.963907,0.26624)
m(9): (-0.963907,0.26624)

The correct answer should, however, be alternating +-1.

After a bit of wailing and gnashing of teeth, I tracked it down to the use of the size_t in the for() loop. I replaced it with an unsigned int, and the unit test worked, but my wider codebase where this is used didn't. Frustrated, I ultimately explicitly cast the inputs to std::exp to be double:

//Improved exponentTest.cc
#include <complex> 
#include <iostream> 



int main(int argc, char **argv)
{
  std::complex<double> iCpx = std::complex<double>(0,1); 
  double baseline[3] = {-0.1, 0, 0};
  
  double theta = atan2(baseline[1], baseline[0]); 
  std::cout << "Theta: " << theta << std::endl;

  for(size_t m = 1; m < 10; m++)
    {
      std::complex<double> thetExp = exp(-1 * (double)m * theta * iCpx); 
      std::cout << "m(" << m << "): " << thetExp << std::endl;
    }
  return 0; 
}

This version consistently gives the correct result:

[stix@localhost ~]$ gcc exponentTest.cc -o expTest -lm -lstdc++
[stix@localhost ~]$ ./expTest
Theta: 3.14159
m(1): (-1,-1.22465e-16)
m(2): (1,2.44929e-16)
m(3): (-1,-3.67394e-16)
m(4): (1,4.89859e-16)
m(5): (-1,-6.12323e-16)
m(6): (1,7.34788e-16)
m(7): (-1,-8.57253e-16)
m(8): (1,9.79717e-16)
m(9): (-1,-1.10218e-15)

So problem solved. However, I don't know why the problem occurred in the first place, and it somewhat smells of a compiler bug. The worst part is that the problem is inconsistent without the explicit cast of m in the exp() function; sometimes it provides correct results, other times it doesn't.

My understanding is that the C++ standard requires that a compiler automatically cast all ints to doubles in something like the line:

double = int * double * int * double;

But the compiler is clearly not doing this.

Is this a compiler bug or is there some gotcha I'm not thinking about in mixing doubles and ints?

Edit: Per one of the comments the exponent line should throw a warning since an unsigned type is being multiplied by -1, however, that is not the case:

[stix@localhost ~]$ gcc exponentTest.cc -o expTest -lm -lstdc++ -Wall -Wextra -pedantic-errors
exponentTest.cc: In function ‘int main(int, char**)’:
exponentTest.cc:6:14: warning: unused parameter ‘argc’ [-Wunused-parameter]
 int main(int argc, char **argv)
              ^~~~
exponentTest.cc:6:27: warning: unused parameter ‘argv’ [-Wunused-parameter]
 int main(int argc, char **argv)
                           ^~~~

Relevant software versions (I'm using devtoolset-7 for GCC 7):

[stix@localhost ~]$ gcc --version
gcc (GCC) 7.3.1 20180303 (Red Hat 7.3.1-5)
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

[stix@localhost ~]$ cat /etc/redhat-release
CentOS Linux release 7.9.2009 (Core)
[stix@localhost ~]$ rpm -qa | grep devtoolset
devtoolset-7-gcc-plugin-devel-7.3.1-5.16.el7.x86_64
devtoolset-7-binutils-2.28-11.el7.x86_64
devtoolset-7-gcc-7.3.1-5.16.el7.x86_64
devtoolset-7-gcc-gfortran-7.3.1-5.16.el7.x86_64
devtoolset-7-runtime-7.1-4.el7.x86_64
devtoolset-7-libquadmath-devel-7.3.1-5.16.el7.x86_64
devtoolset-7-gcc-gdb-plugin-7.3.1-5.16.el7.x86_64
devtoolset-7-libstdc++-devel-7.3.1-5.16.el7.x86_64
devtoolset-7-gcc-c++-7.3.1-5.16.el7.x86_64

Answer 1

The multiplication operator (*) has left-to-right associativity. Thus, in your argument expression:

exp(-1 * m * theta * iCpx)

The -1 * m is performed first and is done using the size_t type. This will result in a very large number (due to -1 being interpreted as an unsigned constant) which, when converted/cast to a double will almost certainly lose precision.

Adding some diagnostics to your loop reveals the problem in all its glory:

for (size_t m = 1; m < 10; m++) {
    std::cout << (-1 * m) << " " << (-1 * (double)m) << std::endl; // Diagnostic!
    std::complex<double> thetExp = exp(-1 * m * theta * iCpx);
    std::cout << "m(" << m << "): " << thetExp << std::endl;
}

Output:

Theta: 3.14159
18446744073709551615 -1
m(1): (-0.963907,0.26624)
18446744073709551614 -2
m(2): (-0.963907,0.26624)
18446744073709551613 -3
m(3): (-0.963907,0.26624)
18446744073709551612 -4
m(4): (-0.963907,0.26624)
18446744073709551611 -5
m(5): (-0.963907,0.26624)
18446744073709551610 -6
m(6): (-0.963907,0.26624)
18446744073709551609 -7
m(7): (-0.963907,0.26624)
18446744073709551608 -8
m(8): (-0.963907,0.26624)
18446744073709551607 -9
m(9): (-0.963907,0.26624)

Explicitly casting m to a double forces that initial multiplication to be (correctly) performed in double-precision arithmetic. (Note: casting m to a (signed) int will also work.)

---

Note that clang-cl warns about this:

warning : implicit conversion changes signedness: 'int' to 'unsigned long long' [-Wsign-conversion]

Answer 2

The simplest fix is to change 1 to 1.0. That forces the computations to be done with doubles:

std::complex<double> thetExp = exp(-1.0 * m * theta * iCpx);

Why is that? Let's take a look at the failing expression:

-1 * m * theta * iCpx

The types here are:

int * size_t * double * std::complex<double>

C++ doesn't look at all the types and pick the "highest" to promote everything to. Instead, it looks at the binary operations one by one, left-to-right, as if there were parentheses grouping them like:

(((int * size_t) * double) * std::complex<double>)

You get in trouble because int * size_t is performed first. Integer promotion rules apply and the int is converted to size_t since on your platform size_t is larger. This means a signed 32-bit integer is being converted to a 64-bit unsigned integer.

Try this and you can see the problem:

std::cout << (size_t) -1 << "\n";

It prints:

18446744073709551615

When you changed -1 * m to -1 * (double) m that got rid of the signed/unsigned problem. -1 was promoted to a double which is simply -1.0.