continuous_sqrt.hpp

Go to the documentation of this file.

#pragma once

#include <complex>
#include <stdexcept>

#include "pi.hpp"


namespace waveblocks {
    namespace math {
        template<class T>
        class ContinuousSqrt
        {
        private:
            std::complex<T> sqrt_;

            T state_;

            bool empty_;

        public:
            ContinuousSqrt()
                : sqrt_()
                , state_()
                , empty_(true)
            { }

            ContinuousSqrt(std::complex<T> sqrt)
                : sqrt_(sqrt)
                , state_(std::arg(sqrt))
                , empty_(false)
            { }

            static T continuate(T ref, T arg)
            {
                const T PI = pi<T>();
                const T RANGE = 0.25*PI; // 0.5*pi allows all inputs

                // determine, how long one needs to rotate the reference angle counter-clock-wise
                // to hit the angle of the 1st root in the domain = [-2*pi;2*pi].
                T rot = arg - ref;

                // force rotation into domain [-pi;pi]
                while (rot >= PI) {
                    rot -= 2.0*PI;
                }
                while (rot < -PI) {
                    rot += 2.0*PI;
                }

                if (-RANGE < rot && rot < RANGE) {
                    return arg;
                }
                else if (rot > PI-RANGE || rot < -PI+RANGE) {
                    if (arg > 0.0)
                        return arg - PI;
                    else
                        return arg + PI;
                }
                else {
                    throw std::runtime_error("continuous_sqrt: too large deviation between computed square root and reference solution");
                }
            }

            std::complex<T> operator()(std::complex<T> input)
            {
                if (empty_) {
                    state_ = 0.5*std::arg(input); // choose principal solution
                    empty_ = false;
                } else {
                    state_ = continuate(state_, 0.5*std::arg(input) );
                }

                sqrt_ = std::polar(std::sqrt(std::abs(input)), state_);

                return sqrt_;
            }

            std::complex<T> operator()() const
            {
                return sqrt_;
            }
            T get_state(void) const
            {
                return state_;
            }
        };
    }
}