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

[Gabriel Dos Reis]

On Thu, Jan 24, 2013 at 8:02 PM, Ed Meyer <address@hidden> wrote:
>
>
> On Thu, Jan 24, 2013 at 7:24 AM, Gabriel Dos Reis
> <address@hidden> wrote:
>>
>> On Wed, Jan 23, 2013 at 3:55 PM, Marc Glisse <address@hidden> wrote:
>> > On Wed, 23 Jan 2013, Jordi Gutiérrez Hermoso wrote:
>> >
>> >> On 23 January 2013 10:21, Marc Glisse <address@hidden> 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
>>
>> I don't know what else will break if the old function is restored.
>> I am in favour of it.
>> In the interest of full disclosure, I've always considered the
>> effort to mimic C99's <complex.h> and <tgmath.h> to be misguided.
>>
>> -- Gaby
>
>
> On the other hand I think it is a bad idea to write code that depends on
> a particular implementation or rounding behaviour

How *concretely* to do you propose to avoid that?  Are you suggesting
something to the programmers or to the implementers?

-- Gaby