octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Ed Meyer
Subject: Re: round-off error in std::pow(std::complex<T>, double) in C++11
Date: Thu, 24 Jan 2013 18:02:34 -0800



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

--
Ed Meyer

reply via email to

[Prev in Thread] Current Thread [Next in Thread]