Re: wrong result in complex pow

From: Christopher Kormanyos (e_float_at_[hidden])
Date: 2021-04-19 15:30:51

 > acot2, atanh2 and acoth2. Do you know> the algorithms or has information about them?
I might have some notes from previousprojects. Please allow me to get backto you on this particular query.
Thanks again Gero!
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