Complex log implementations might want to compute the magnitude as
sqrt(real(x)*real(x) + imag(x)*imag(x)). This results in loss of
accuracy if one of real(x) or imag(x) (the latter in this case) is small
relative to the other. Using std::hypot instead could help. I'm a bit
busy at the moment but I can study this in more detail if you would
like, though I'm sure there are plenty of more qualified experts here :
- ) .

mfh

On 3/3/2021 09:01, Peter Sommerlad (C++) via Sci wrote:
> 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_lists.isocpp.org
>>> <mailto:std-discussion_at_lists.isocpp.org>> 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].org
>>> <mailto:Std-Discussion_at_[hidden]>
>>> https://lists.isocpp.org/mailman/listinfo.cgi/std-discussion
>>>
>>>
>>>
>>
>>
>
>