[Bug libstdc++/90638] Wrong results of pow(complex<T> , T1) function when the T1 type differs from T and from int

[Igor.Smirnov at cern dot ch]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=90638

--- Comment #7 from Igor Smirnov <Igor.Smirnov at cern dot ch> ---
There can be still problems with big integer powers of complex numbers in the
latest versions.
I just don't see a single correct answer.
As above the expected correct responses are marked by "should be" and the
received by "actually".
Notations of version started from "C++" are from the above mentioned site:

---------------------------------------------------------------


#include <iostream>
#include <complex>
#include <limits.h>

int main(void)
{
  double phi = M_PI / LONG_MAX;
  std::cout<<phi<<'\n';
  // Should be 3.40612e-19
  std::complex< double > var0 = std::polar(1.0, phi);
  std::cout<<var0 << '\n';
  // Should be (1,3.40612e-19)

  long lp1 = LONG_MAX;
  unsigned long lp2 = (unsigned long)lp1+1;
  std::cout<<lp1<<' '<<lp2<<'\n';
  // Should be 2147483647 2147483648

  std::cout<<' '<<pow(var0,  lp1)
           <<' '<<pow(var0,  lp2)
           <<'\n'<<'\n';
  // Should be (-1, zero or very small number) (-1,-3.40612e-19)
  // Actually (1,-3.40612e-19) (1,0) by 4.8.5
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++ 10.0.0.1
  // (1,-3.40612e-19) (1,0) by C++03 10.0.0.1
  // (1,-3.40612e-19) (1,0) by C++03(GNU) 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++11 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++11(GNU) 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++ 14 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++ 14(GNU) 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++ 17 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++17(GNU) 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++2a 10.0.0.1
  // (-1,1.22465e-16) (-1,1.22465e-16) by C++2a(GNU) 10.0.0.1

  long double phid = M_PI / LONG_MAX;
  std::cout<<phid<<'\n';
  // Should be 3.40612e-19
  std::complex< long double > vard0 = std::polar((long double)1.0, phid);
  std::cout<<vard0 << '\n';
  // Should be (1,3.40612e-19)


  std::cout<<' '<<pow(vard0,  lp1)
           <<' '<<pow(vard0,  lp2)
           <<'\n'<<'\n';
  // Should be (-1, zero or very small number) (-1,-3.40612e-19)
  // Actually (1,-3.40612e-19) (1,0) at 4.8.5
  // (-1,1.22898e-16) (-1,1.22465e-16) by C++ 10.0.0.1
  // (1,-3.40612e-19) (1,0) by C++03 10.0.0.1
  // (1,-3.40612e-19) (1,0) C++03(GNU) 10.0.0.1
  // etc.
  // (-1,1.22898e-16) (-1,1.22465e-16) by C++ 17 10.0.0.1
  // (-1,1.22898e-16) (-1,1.22465e-16) by C++17(GNU) 10.0.0.1
  // (-1,1.22898e-16) (-1,1.22465e-16) by C++2a(GNU) 10.0.0.1


  double phi1 = M_PI / ULONG_MAX;
  std::cout<<phi1<<'\n';
  // Should be 1.70306e-19
  std::complex< double > var1 = std::polar(1.0, phi1);
  std::cout<<var1 <<'\n';
  // Should be (1,1.70306e-19)

  unsigned long lp3 = ULONG_MAX;
  std::cout<<lp3<<'\n';
  // Should be 18446744073709551615

  std::cout<<' '<<pow(var1,  lp3)
           <<'\n'<<'\n';
  // Should be (-1, zero or very small number)
  // Actually (1,-1.70306e-19) by 4.8.5
  // (-1,1.22465e-16) by C++ 10.0.0.1
  // (1,-1.70306e-19) by C++03 10.0.0.1
  // etc.
  // (-1,1.22465e-16) by C++ 17 10.0.0.1
}

-------------------------------------------------------------
Example of the output by 4.8.5:

3.40612e-19
(1,3.40612e-19)
9223372036854775807 9223372036854775808
 (1,-3.40612e-19) (1,0)

3.40612e-19
(1,3.40612e-19)
 (1,-3.40612e-19) (1,0)

1.70306e-19
(1,1.70306e-19)
18446744073709551615
 (1,-1.70306e-19)

Other versions are in comments.