0

While trying to debug some unit tests, I found I'm getting a difference when trying to compare the results of std::pow using std::complex<double>. In the first case, I'm evaluating std::pow directly. In the second case, I'm accessing these from inside a std::vector. When I compare them, I get a difference on the order of machine epsilon when -O3 is used. The minimum case example looks like

perlmutter:~/bug> cat test.cpp
#include <complex>
#include <iostream>
#include <vector>

int main(int argc, const char * argv[]) {
    std::vector<std::complex<double>> base(1, 4.0);
    std::vector<std::complex<double>> exp(1, -0.23);

    std::cout << std::pow(std::complex(4.0), std::complex(-0.23)) - 
                 std::pow(base.at(0),        exp.at(0)) << std::endl;
    
    return 0;
}
perlmutter:~/bug> /opt/cray/pe/gcc/11.2.0/bin/c++ --version
c++ (GCC) 11.2.0 20210728 (Cray Inc.)
Copyright (C) 2021 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.

perlmutter:~/bug> /opt/cray/pe/gcc/11.2.0/bin/c++ -O3 -o test test.cpp
perlmutter:~/bug> ./test
(1.11022e-16,-0)

This is being run on an AMD EPYC 7713 64-Core Processor. I don't see this error when using float, double, or std::complex<float>.

Retesting with gcc 12.1.0 shows the same error.

(base) perlmutter:~/bug> gcc --version
gcc (GCC) 12.1.0 20220506 (HPE)
Copyright (C) 2022 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.

perlmutter:~/bug> g++ -O3 -o test test.cpp
perlmutter:~/bug> ./test
(1.11022e-16,-0)
user1139069
  • 1,505
  • 2
  • 15
  • 27
  • Unable to reproduce with gcc 12. I get `(0,0)`. – Sam Varshavchik Jan 22 '23 at 00:01
  • @SamVarshavchik Interesting because I can. What CPU did you run on? – user1139069 Jan 22 '23 at 00:07
  • AMD Ryzen 2950X, same results on an Intel Xeon, as well. – Sam Varshavchik Jan 22 '23 at 00:08
  • 4
    Looks like they differ by exactly 1 ULP: https://godbolt.org/z/q8176Tq95 `std::pow(std::complex(4.0), std::complex(-0.23))` is constant folded but `std::pow(base.at(0), exp.at(0))` calls the library `cpow` at runtime and returns a slightly different result. See this similar question: https://stackoverflow.com/q/31418209. `-frounding-math` causes gcc to call `cpow` at runtime even with the constants (since rounding might be different than during compile time) which causes them to be equal. – Artyer Jan 22 '23 at 00:37
  • 4
    It probably doesn't make much sense to unit test on the precise value. It is not specified how accurate `std::pow` is and even if the compiler and standard library specify it to be accurate up to rounding to representable values, then there are still usually two possible outcomes. – user17732522 Jan 22 '23 at 00:44
  • But e.g. glibc on x86-64 lists up to 2 ULP error for `cpow` in their documentation at https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html. I would guess that the compiler can do it better at compile-time during constant-folding. – user17732522 Jan 22 '23 at 00:49
  • I can replicate this on Intel i7-5820K with gcc 9.4 but only on `-O3` – Something Something Jan 22 '23 at 00:58
  • @user1139069 "When I compare them, I get a difference on the order of machine epsilon" --> Please post the 2 values with sufficient (17+ significant) decimal places or in some hex notation. – chux - Reinstate Monica Jan 22 '23 at 04:48

0 Answers0