std::sqrt(std::complex<T>) produces incorrect results for large real components

[cjdb]

The following program produces an incorrect complex value due to floating-point numbers being unable to store a precise enough angle when converting to polar values.

template<class T>
void sqrt()
{
  auto const x = std::complex<T>{std::numeric_limits<float>::lowest(), 1.0f};
  std::println("x = {} + i * {}", x.real(), x.imag());
  auto const s = std::sqrt(x);
  std::println("sqrt(x) = {} + i * {}", s.real(), s.imag());
}

int main()
{
  std::println("==== expected ====");
  std::println("x = {} + i * {}", std::numeric_limits<float>::lowest(), 1.0f);
  std::println("sqrt(x) = 2.7105055e-20 + i * 1.8446744e19");

  std::println("==== sqrt<float>() ====");
  sqrt<float>();

  std::println("==== sqrt<double>() ====");
  sqrt<double>();
  
  std::println("==== sqrt<long double>() ====");
  sqrt<long double>();
}

==== expected ====
x = -3.4028235e+38 + i * 1
sqrt(x) = 2.7105055e-20 + i * 1.8446744e19
==== sqrt<float>() ====
x = -3.4028235e+38 + i * 1
sqrt(x) = -806332727296 + i * 1.8446743e+19
==== sqrt<double>() ====
x = -3.4028234663852886e+38 + i * 1
sqrt(x) = 1129.537270565105 + i * 18446743523953729536
==== sqrt<long double>() ====
x = -3.4028234663852886e+38 + i * 1
sqrt(x) = -0.4626957582701517 + i * 18446743523953729536

Confirmed by @lntue.

[llvmbot]

@llvm/issue-subscribers-bug

Author: Christopher Di Bella (cjdb)

The following program produces an incorrect complex value due to floating-point numbers being unable to store a precise enough angle when converting to polar values.

template<class T>
void sqrt()
{
  auto const x = std::complex<T>{std::numeric_limits<float>::lowest(), 1.0f};
  std::println("x = {} + i * {}", x.real(), x.imag());
  auto const s = std::sqrt(x);
  std::println("sqrt(x) = {} + i * {}", s.real(), s.imag());
}

int main()
{
  std::println("==== expected ====");
  std::println("x = {} + i * {}", std::numeric_limits<float>::lowest(), 1.0f);
  std::println("sqrt(x) = 2.7105055e-20 + i * 1.8446744e19");

  std::println("==== sqrt<float>() ====");
  sqrt<float>();

  std::println("==== sqrt<double>() ====");
  sqrt<double>();
  
  std::println("==== sqrt<long double>() ====");
  sqrt<long double>();
}

==== expected ====
x = -3.4028235e+38 + i * 1
sqrt(x) = 2.7105055e-20 + i * 1.8446744e19
==== sqrt<float>() ====
x = -3.4028235e+38 + i * 1
sqrt(x) = -806332727296 + i * 1.8446743e+19
==== sqrt<double>() ====
x = -3.4028234663852886e+38 + i * 1
sqrt(x) = 1129.537270565105 + i * 18446743523953729536
==== sqrt<long double>() ====
x = -3.4028234663852886e+38 + i * 1
sqrt(x) = -0.4626957582701517 + i * 18446743523953729536

Confirmed by @lntue.