sqrt.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_CORE_COMPLEX_FUNCTIONS_SQRT_H
#define pRC_CORE_COMPLEX_FUNCTIONS_SQRT_H
 
#include <prc/core/complex/complex.hpp>
 
namespace pRC
{
    template<IsComplex T>
    static inline constexpr auto sqrt(T const &a)
    {
        using R = typename T::Type;
 
        if(a.real() == zero())
        {
            auto const t = sqrt(abs(a.imag()) / identity<R>(2));
            if(a.imag() < zero())
            {
                return Complex(t, -t);
            }
            else
            {
                return Complex(t, t);
            }
        }
        else
        {
            auto const t = sqrt(identity<R>(2) * (abs(a) + abs(a.real())));
            auto const t2 = t / identity<R>(2);
            auto const at = a.imag() / t;
 
            if(a.real() > zero())
            {
                return Complex(t2, at);
            }
            else
            {
                auto const real = abs(at);
                if(a.imag() < zero())
                {
                    return Complex(real, -t2);
                }
                else
                {
                    return Complex(real, t2);
                }
            }
        }
    }
}
#endif // pRC_CORE_COMPLEX_FUNCTIONS_SQRT_H