Re: wrong result in complex pow

From: Gero Peterhoff (g.peterhoff_at_[hidden])
Date: 2021-04-18 18:08:36

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