From: Christopher Kormanyos (e_float_at_[hidden])
Date: 2021-04-19 15:28:39
> i see that the pow-bug> (https://github.com/boostorg/math/issues/506)> is also included in 1.76
Hi Gero,
Thank you for following up on this.
You are correct.
The state of that issue is *open*, which means,unfortunately, that I have not yet fixed it.
The intent is to gather outstanding relatedpoints, bug-lets and the like and fix theseen-masse in a more dedicated driveforward.
The time for me to drive forward on thisis rapidly approaching. The deal is that1.76 was full of many larger organizationsland code-technical efforts in both Mathas well as Multiprecision, followed bycontinued cleanup to this very day.
Forward mothion on this issue is plannedforthcoming.
kind regards, Chris
On Sunday, April 18, 2021, 8:08:39 PM GMT+2, Gero Peterhoff <g.peterhoff_at_[hidden]> wrote:
Hello Chris,
i see that the pow-bug (https://github.com/boostorg/math/issues/506) is also included in 1.76. I think regardless of special values ​​(zero, nan, inf) this should be fixed:
inline complex <BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE> pow(const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE& x,
                                const complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>& a)
{
   constexpr BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE
      zero_v{0};
   const complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>
      log_i_zero{std::log(std::abs(x)), std::atan2(zero_v, x)};
   return std::exp(a * log_i_zero);
}
There is also an error in your sqrt:
inline complex <BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE> sqrt(const complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>& x)
{
   using std::fabs;
   using std::sqrt;
   if (x.real() > 0)
   {
      const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE s = sqrt((fabs(x.real()) + std::abs (x)) / 2);
      return complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>(s, x.imag()
/ (s * 2));
   }
   else if (x.real() < 0)
   {
      const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE s = sqrt((fabs(x.real()) + std::abs(x)) / 2);
      const bool imag_is_neg = (x.imag() < 0);
      return complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>(fabs(x.imag()) / (s * 2), (imag_is_neg ? -s : s));
   }
   else
   {
       const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE sqrt_xi_half =
sqrt(x.imag() / 2);
      return complex <BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE> (sqrt_xi_half, sqrt_xi_half);
   }
}
As you can see this gives wrong results for real==0 and imag<0 -> std::sqrt from negative number. There are several ways to remedy the error:
inline complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE> sqrt(const complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>& x)
{
   // branchfree
   return std::pow(x, boost::math::constants::half<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>());
   // constexpr BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE
   //    zero_v{0};
   /*
   // 2 branches
   const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE
      s{std::sqrt((std::abs(x.real()) + std::abs(x)) * boost::math::constants::half<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>())},
      t{(std::abs(x.imag ()) * boost::math::constants::half<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>()) / s};
   return (x.real() >= zero_v)   ?
      complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>{s, std::copysign (t, x.imag())}:
      complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>{t, std::copysign (s, x.imag())};
   */
   /*
   // 3 branches
   if (x.real() == zero_v)
   {
      const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE
         s{std::sqrt(std::abs(x.imag()) * boost::math::constants::half<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>())};
      return complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>{s, std::copysign(s, x.imag())};
   }
   else
   {
      const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE
         s{std::sqrt((std::abs(x.real()) + std::abs(x)) * boost::math::constants::half<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>())},
         t{(std::abs(x.imag()) * boost::math::constants::half<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>()) / s};
      return (x.real() > zero_v)   ?
         complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>{s, std::copysign(t, x.imag())}:
         complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>{t, std::copysign(s, x.imag())};
   }
   */
}
I think pow(x, 0.5) is accurate enough - it is also the fastest. According to the same scheme, cbrt can also be implemented, which is missing so far:
inline complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE> cbrt(const complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>& x)
{
   constexpr BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE
      three_v{3};
   return std::exp(std::log(x) / three_v);
}
Something else. I am currently working on an extended math lib, in which I provide many functions that were previously missing. But i don't know how to implement these functions (like atan2): acot2, atanh2 and acoth2. Do you know the algorithms or has information about them?
thx
Gero