C++ Logo

std-discussion

Advanced search

Re: Problems with pow(std::complex<T>, double)

From: Peter Sommerlad (C++) <"Peter>
Date: Wed, 3 Mar 2021 17:01:34 +0100
Taking the generic definition of complex pow function I can confirm that
the pow() implementation is carrying the same error. May be, what you
attempt is just beyond reasonable precision to expect from floating point.

I am CC-in numerics SG6 to query the numerics experts there.

https://godbolt.org/z/GET64M

#include <iostream>
#include <complex>
#include <cmath>

int main()

{
      using namespace std::literals;

      auto x = -0.1+ 1e-100i;
      auto y = 2.0+0i;
      auto z = std::exp(y*std::log(x));
      auto zpow = std::pow(x,y);

      std::cout << z << zpow << std::pow(-0.1+ 1e-100i, 2) << ":" << x *
x <<std::endl;

}


Peter.

Peter Sommerlad (C++) via Std-Discussion wrote on 03.03.21 16:48:
> using an integral 2nd argument to pow() solves the issue. May be
> python is optimizing by checking that 2.0 is actually integral.
>
>
> https://godbolt.org/z/a35G55
>
> #include <iostream>
> #include <complex>
> #include <cmath>
>
> int main()
> {
> using namespace std::literals;
> auto x = -0.1+ 1e-100i;
> std::cout << std::pow(-0.1+ 1e-100i, 2) << ":" << x * x <<std::endl;
> }
>
>
> But I am far from being a numerics or python expert.
>
> Regards
> Peter.
>
> will wray via Std-Discussion wrote on 03.03.21 16:29:
>> Confirm; something fishy here.
>> Curious if the unqualified 'pow' was a culprit, I qualified as std::pow
>> and still confirm the same bad result on gcc, clang, icc and msvc latest.
>>
>> The equivalent C program has the same output in gcc and clang
>> https://godbolt.org/z/x1Ga99
>>
>> #include "stdio.h"
>> #include "complex.h"
>> #include "math.h"
>>
>> int main() {
>> double complex x = -0.1 + I * 1e-100;
>> double complex p = cpow(x, 2.0);
>> printf("%g,%g\n",creal(p),cimag(p));
>> }
>>
>> Implies:
>> (1) cut-n-paste code of same numerical precision 'bug' across
>> implementations, or
>> (2) precision 'bug' in the specification of complex pow
>>
>> If 2. then it may be hard to correct it as a defect as it will change
>> existing results?
>>
>> Not sure if there's still a dedicated numerics SG -
>> SG19 ML has some remit for numerics, including automatic differentiation
>> SG14 also does numerics
>> or, this may also be in the remit of the new Joint C and C++ Study Group
>> (there has been recent work on complex in C)
>>
>>
>>
>>
>> On Wed, Mar 3, 2021 at 9:05 AM Bell, Ian H. (Fed) via Std-Discussion
>> <std-discussion_at_[hidden]
>> <mailto:std-discussion_at_[hidden]>> wrote:
>>
>> The recent “rediscovery” of complex step derivatives
>> (https://sinews.siam.org/Details-Page/differentiation-without-a-difference)
>>
>> has made numerical differentiation as accurate as any other
>> evaluation in double precision arithmetic. To fully make use of
>> this technique, all functions must accept either complex or double
>> arguments. In principle that is no problem for C++. In practice,
>> serious problems occur in some cases.____
>>
>> __ __
>>
>> Here first is a simple example in Python of when things go right.
>> The derivative of x^2.0 is 2.0*x, so the derivative of x^2.0 at
>> x=-0.1 should be dy/dx=-0.2. In Python, no problem to use complex
>> step derivatives to evaluate:____
>>
>> __ __
>>
>> h = 1e-100____
>>
>> z = -0.1+1j*h____
>>
>> print(pow(z, 2.0), pow(z, 2.0).imag/h, (z*z).imag/h)____
>>
>> __ __
>>
>> gives____
>>
>> __ __
>>
>> (0.010000000000000002-2e-101j) -0.2 -0.2____
>>
>> __ __
>>
>> On the contrary, the same example in C++ ____
>>
>> __ __
>>
>> #include <iostream>____
>>
>> #include <complex>____
>>
>> #include <cmath>____
>>
>> __ __
>>
>> int main()____
>>
>> {____
>>
>> std::cout << pow(std::complex<double>(-0.1, 1e-100), 2.0) <<
>> std::endl; ____
>>
>> }____
>>
>> __ __
>>
>> gives____
>>
>> __ __
>>
>> (0.01,-2.44929e-18)____
>>
>> __ __
>>
>> __ __
>>
>> I believe the problem has to do with the handling of the branch-cut
>> of the log function. In any case, this demonstrates a result that
>> is silently in error by 83 orders of magnitude! Had I multiplied the
>> complex step by itself rather than pow(z,2.0), I would have obtained
>> the correct result.____
>>
>> __ __
>>
>> I realize that I am probing an uncomfortable part of the complex
>> plane, but I wonder if this could be handled more like Python, to
>> minimize surprises for complex step derivative approaches?____
>>
>> __ __
>>
>> Ian____
>>
>> __ __
>>
>> -- Std-Discussion mailing list
>> Std-Discussion_at_[hidden]
>> <mailto:Std-Discussion_at_[hidden]>
>> https://lists.isocpp.org/mailman/listinfo.cgi/std-discussion
>>
>>
>>
>
>


-- 
Peter Sommerlad
Better Software: Consulting, Training, Reviews
Modern, Safe & Agile C++
peter.cpp_at_[hidden]
+41 79 432 23 32

Received on 2021-03-03 10:01:41