C++ Logo

std-discussion

Advanced search

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

From: Peter Sommerlad (C++) <"Peter>
Date: Wed, 3 Mar 2021 16:48:25 +0100
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 09:48:32