Date: Wed, 3 Mar 2021 16:59:03 +0000
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&data=04%7C01%7Cian.bell%40nist.gov%7C042f0cfb
> 34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C
> 637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjo
> iV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=0vOuLB3RfmPvx
> L5I3nNXINChDJ4b1w730ta%2FbPAXDBs%3D&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&data=04%7C01%7Cian.bell%40nist.gov%7C042f0c
>> fb34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1
>> %7C637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJ
>> QIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=w%2BtYYz
>> wrGQXBb2iaXhQ9C65F8W127TnsDa7iM%2F3qrz4%3D&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&data=04%7C01%7Cian.bell%40nist.gov%7C042f
>>> 0cfb34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%
>>> 7C1%7C637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDA
>>> iLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=PRy
>>> 5gn7oZ7SoA7PPXgeHv4Dnhd3mPnO4UjtCQqlAIok%3D&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
>>> &data=04%7C01%7Cian.bell%40nist.gov%7C042f0cfb34a84fc3940a08d8de
>>> 657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C6375038746097311
>>> 24%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBT
>>> iI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=ThlJMhFiWnISI07qXe8pxDxob
>>> ilWxYSPVODikEJ6R%2FQ%3D&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&data=04
>>> %7C01%7Cian.bell%40nist.gov%7C042f0cfb34a84fc3940a08d8de657376%7C2ab
>>> 5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C637503874609731124%7CUnknown
>>> %7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiL
>>> CJXVCI6Mn0%3D%7C3000&sdata=w%2BRTwAMx2OpVfcaLP453DXKdo1iTd2%2BGW
>>> G4xW28iMK4%3D&reserved=0
>>>
>>>
>>>
>>
>>
>
>
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&data=04%7C01%7Cian.bell%40nist.gov%7C042f0cfb
> 34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C
> 637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjo
> iV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=0vOuLB3RfmPvx
> L5I3nNXINChDJ4b1w730ta%2FbPAXDBs%3D&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&data=04%7C01%7Cian.bell%40nist.gov%7C042f0c
>> fb34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1
>> %7C637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJ
>> QIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=w%2BtYYz
>> wrGQXBb2iaXhQ9C65F8W127TnsDa7iM%2F3qrz4%3D&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&data=04%7C01%7Cian.bell%40nist.gov%7C042f
>>> 0cfb34a84fc3940a08d8de657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%
>>> 7C1%7C637503874609731124%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDA
>>> iLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=PRy
>>> 5gn7oZ7SoA7PPXgeHv4Dnhd3mPnO4UjtCQqlAIok%3D&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
>>> &data=04%7C01%7Cian.bell%40nist.gov%7C042f0cfb34a84fc3940a08d8de
>>> 657376%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C6375038746097311
>>> 24%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBT
>>> iI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=ThlJMhFiWnISI07qXe8pxDxob
>>> ilWxYSPVODikEJ6R%2FQ%3D&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&data=04
>>> %7C01%7Cian.bell%40nist.gov%7C042f0cfb34a84fc3940a08d8de657376%7C2ab
>>> 5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C637503874609731124%7CUnknown
>>> %7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiL
>>> CJXVCI6Mn0%3D%7C3000&sdata=w%2BRTwAMx2OpVfcaLP453DXKdo1iTd2%2BGW
>>> G4xW28iMK4%3D&reserved=0
>>>
>>>
>>>
>>
>>
>
>
Received on 2021-03-03 10:59:08