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 pow(const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE& x,
const complex& a)
{
constexpr BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE
zero_v{0};
const complex
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 sqrt(const complex& 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(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(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 (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 sqrt(const complex& x)
{
// branchfree
return std::pow(x, boost::math::constants::half());
// 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())},
t{(std::abs(x.imag ()) * boost::math::constants::half()) / s};
return (x.real() >= zero_v) ?
complex{s, std::copysign (t, x.imag())}:
complex{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())};
return complex{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())},
t{(std::abs(x.imag()) * boost::math::constants::half()) / s};
return (x.real() > zero_v) ?
complex{s, std::copysign(t, x.imag())}:
complex{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 cbrt(const complex& 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
