C++ Logo

STD-DISCUSSION

Advanced search

Subject: Re: Problems with pow(std::complex<T>, double)
From: Peter Sommerlad (C++) (peter.cpp_at_[hidden])
Date: 2021-03-03 10:01:34


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

STD-DISCUSSION list run by std-discussion-owner@lists.isocpp.org

Older Archives on Google Groups