C++ Logo

std-discussion

Advanced search

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

From: will wray <wjwray_at_[hidden]>
Date: Wed, 3 Mar 2021 10:29:52 -0500
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]> 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]
> https://lists.isocpp.org/mailman/listinfo.cgi/std-discussion
>

Received on 2021-03-03 09:30:12