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

[Marc Glisse <marc dot glisse at inria dot fr>]

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

> On 23 January 2013 10:21, Marc Glisse <marc.glisse@inria.fr> wrote:
> > 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.
>
> I am not very good at reading standardese, but what does "effectively
> cast" mean here?

The effect is the same as if they were cast (easiest way to satisfy this 
requirement is to cast indeed).

> > Maybe the libm function doesn't round well then?
>
> If you have floats or doubles as your type, I don't think you can
> avoid the rounding error.

If pow is a primitive function, avoiding rounding errors seems doable 
(not easy but possible). One issue is that the standard defines pow(x,y) 
as exp(y*log(x)), and I don't know if they mean it in a mathematical 
sense (infinite precision) or in an implementation sense (compute log(x), 
round, multiply with y, round, etc). I hope it is the former. The code in 
<complex> that uses log, exp, polar instead of forwarding to the libc pow 
function might need to be protected by _GLIBCXX_FAST_MATH.

> It seems to me that having the int overload, truncating its precision
> if necessary, and then performing computations with integers (e.g. by
> binary exponentiation), should not break anything else.

If we do we probably want that code to work not just for int but also 
unsigned short, etc. But then should we also re-introduce special code for 
pow(double,int)? I think forwarding to the generic libm function is not 
that bad (it often detects int arguments).

I guess I'll shut up and wait for the real libstdc++ experts to speak up 
:-)

--
Marc Glisse