C++ Logo

STD-DISCUSSION

Advanced search

Subject: Re: [isocpp-sci] Problems with pow(std::complex<T>, double)
From: Bell, Ian H. (Fed) (ian.bell_at_[hidden])
Date: 2021-03-03 10:59:03


Python must be special casing integer exponents that can be exactly represented in double precision. For instance for an exponent *just* above 2.0, we have a similar problem in Python:

h = 1e-100
z = -0.1+1j*h
p = 2.00000001
print(pow(z, p), pow(z, p).imag/h)

yielding

(0.00999999976974149+3.1415925780626665e-10j) 3.1415925780626666e+90

which is off by ~90 orders of magnitude from the correct value of -0.200000001

-----Original Message-----
From: Mark Hoemmen <mhoemmen_at_[hidden]>
Sent: Wednesday, March 3, 2021 11:58 AM
To: sci_at_[hidden]; std-discussion_at_[hidden]
Cc: Peter Sommerlad (C++) <peter.cpp_at_[hidden]>; Bell, Ian H. (Fed) <ian.bell_at_[hidden]>
Subject: Re: [isocpp-sci] [std-discussion] Problems with pow(std::complex<T>, double)

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://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgodb
> olt.org%2Fz%2FGET64M&amp;data=04%7C01%7Cian.bell%40nist.gov%7C042f0cfb
> 34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C
> 637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjo
> iV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&amp;sdata=0vOuLB3RfmPvx
> L5I3nNXINChDJ4b1w730ta%2FbPAXDBs%3D&amp;reserved=0
>
> #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://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgod
>> bolt.org%2Fz%2Fa35G55&amp;data=04%7C01%7Cian.bell%40nist.gov%7C042f0c
>> fb34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1
>> %7C637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJ
>> QIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&amp;sdata=w%2BtYYz
>> wrGQXBb2iaXhQ9C65F8W127TnsDa7iM%2F3qrz4%3D&amp;reserved=0
>>
>> #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://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgo
>>> dbolt.org%2Fz%2Fx1Ga99&amp;data=04%7C01%7Cian.bell%40nist.gov%7C042f
>>> 0cfb34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%
>>> 7C1%7C637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDA
>>> iLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&amp;sdata=PRy
>>> 5gn7oZ7SoA7PPXgeHv4Dnhd3mPnO4UjtCQqlAIok%3D&amp;reserved=0
>>>
>>> #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://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fs
>>> inews.siam.org%2FDetails-Page%2Fdifferentiation-without-a-difference
>>> &amp;data=04%7C01%7Cian.bell%40nist.gov%7C042f0cfb34a84fc3940a08d8de
>>> 657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C6375038746097311
>>> 24%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBT
>>> iI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&amp;sdata=ThlJMhFiWnISI07qXe8pxDxob
>>> ilWxYSPVODikEJ6R%2FQ%3D&amp;reserved=0)
>>>
>>>     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://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fli
>>> sts.isocpp.org%2Fmailman%2Flistinfo.cgi%2Fstd-discussion&amp;data=04
>>> %7C01%7Cian.bell%40nist.gov%7C042f0cfb34a84fc3940a08d8de657376%7C2ab
>>> 5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C637503874609731124%7CUnknown
>>> %7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiL
>>> CJXVCI6Mn0%3D%7C3000&amp;sdata=w%2BRTwAMx2OpVfcaLP453DXKdo1iTd2%2BGW
>>> G4xW28iMK4%3D&amp;reserved=0
>>>
>>>
>>>
>>
>>
>
>


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

Older Archives on Google Groups