round-off error in std::pow(std::complex<T>, double) in C++11

[Marc Glisse]

On Wed, 23 Jan 2013, Jordi Gutiérrez Hermoso wrote:

> Greetings, fellow GNU hackers.
>
> Consider the following program,
>
>    #include <iostream>
>    #include <iomanip>
>    #include <cmath>
>    #include <complex>
>
>    int main(){
>      using namespace std;
>      cout << setprecision(16) << pow(complex<double>(0,1), 2) << endl;
>    }
>
> Compile it with and without -std=c++0x:
>
>    jordi@Iris:~$ g++ foo.c++ -o foo && ./foo
>    (-1,0)
>
>    jordi@Iris:~$ g++ foo.c++ -o foo -std=c++0x && ./foo
>    (-1,1.224646799147353e-16)
>
> The round-off error introduced by C++11 is partly due to the following
> discussion:
>
>    http://stackoverflow.com/questions/5627030/why-was-stdpowdouble-int-removed-from-c11
>
> For some reason, the pow(complex<T>, int) overload seems to be
> missing in C++11, but it doesn't seem like the standard says it should
> be. Can it be restored in GNU libstdc++ for C++11?

Here is what [cmplx.over] says:

"Function template pow shall have additional overloads sufficient to 
ensure, for a call with at least one argument of type complex<T>:
1. If either argument has type complex<long double> or type long double, 
then both arguments are effectively cast to complex<long double>.
2. Otherwise, if either argument has type complex<double>, double, or an 
integer type, then both arguments are effectively cast to complex<double>.
3. Otherwise, if either argument has type complex<float> or float, then 
both arguments are effectively cast to complex<float>."

So it looks like we are forced to convert 2 to a complex<double> and then 
call pow on that. Maybe the libm function doesn't round well then?

It should be possible to have a pow(complex<T>,int) overload (as-if), the 
question is what else is needed so it doesn't break stuff? (maybe nothing)

The change in the standard comes from the fact that the listed overloads 
were ambiguous for many calls.

-- 
Marc Glisse