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

>

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