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
2
Answers
The simplest fix is to change
1
to1.0
. That forces the computations to be done withdouble
s:Why is that? Let’s take a look at the failing expression:
The types here are:
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:
You get in trouble because
int * size_t
is performed first. Integer promotion rules apply and theint
is converted tosize_t
since on your platformsize_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:
It prints:
When you changed
-1 * m
to-1 * (double) m
that got rid of the signed/unsigned problem.-1
was promoted to adouble
which is simply-1.0
.The multiplication operator (
*
) has left-to-right associativity. Thus, in your argument expression:The
-1 * m
is performed first and is done using thesize_t
type. This will result in a very large number (due to-1
being interpreted as an unsigned constant) which, when converted/cast to adouble
will almost certainly lose precision.Adding some diagnostics to your loop reveals the problem in all its glory:
Output:
Explicitly casting
m
to adouble
forces that initial multiplication to be (correctly) performed in double-precision arithmetic. (Note: castingm
to a (signed)int
will also work.)Note that clang-cl warns about this: