polynomial.hh

#ifndef VTESTBED_CORE_POLYNOMIAL_HH
#define VTESTBED_CORE_POLYNOMIAL_HH

#include <algorithm>
#include <complex>
#include <initializer_list>
#include <type_traits>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/core/math.hh>

namespace vtb {

    namespace core {

        template <class T>
        struct is_complex: public std::false_type {};

        template <class T>
        struct is_complex<std::complex<T>>: public std::true_type {};


        template <class T>
        class Polynomial {

        public:
            typedef blitz::Array<T,1> array_type;
            typedef T value_type;
            typedef value_type* pointer;
            typedef const value_type* const_pointer;
            typedef typename std::conditional<
                is_complex<T>::value,
                Polynomial<T>,
                Polynomial<std::complex<T>>
            >::type complex_polynomial;
            typedef typename std::conditional<
                is_complex<T>::value,
                blitz::Array<T,1>,
                blitz::Array<std::complex<T>,1>
            >::type complex_array;

        private:
            array_type _coefficients;

        public:

            Polynomial() = default;

            inline explicit
            Polynomial(int size, const T& value = T{0}):
            _coefficients(size) {
                _coefficients = value;
            }

            inline explicit
            Polynomial(array_type coefs):
            _coefficients{coefs} {}

            inline
            Polynomial(std::initializer_list<T> coefs):
            _coefficients(coefs.size()) {
                std::copy(coefs.begin(), coefs.end(), this->begin());
            }

            inline
            Polynomial(const Polynomial& rhs):
            _coefficients{rhs._coefficients.copy()} {}

            inline Polynomial&
            operator=(const Polynomial& rhs) {
                this->_coefficients.reference(rhs._coefficients.copy());
                return *this;
            }

            inline
            Polynomial(Polynomial&& rhs) {
                this->_coefficients.reference(rhs._coefficients);
            }

            inline Polynomial&
            operator=(Polynomial&& rhs) {
                this->swap(rhs);
                return *this;
            }

            ~Polynomial() = default;

            inline value_type
            evaluate(value_type x) const {
                value_type result{0};
                const_pointer last = this->begin()-1;
                const_pointer first = this->end()-1;
                while (first != last) {
                    result = fma(result, x, *first);
                    --first;
                }
                return result;
            }

            inline value_type
            operator()(value_type x) const {
                return this->evaluate(x);
            }

            inline const value_type&
            operator[](int i) const noexcept {
                return this->_coefficients(i);
            }

            inline value_type&
            operator[](int i) noexcept {
                return this->_coefficients(i);
            }

            inline const_pointer
            begin() const noexcept {
                return this->_coefficients.data();
            }

            inline pointer
            begin() noexcept {
                return this->_coefficients.data();
            }

            inline const_pointer
            end() const noexcept {
                return this->_coefficients.data() + this->size();
            }

            inline pointer
            end() noexcept {
                return this->_coefficients.data() + this->size();
            }

            inline int
            size() const noexcept {
                return this->_coefficients.extent(0);
            }

            inline int
            order() const noexcept {
                return std::max(0, this->size() - 1);
            }

            inline void
            swap(Polynomial& rhs) {
                blitz::swap(this->_coefficients, rhs._coefficients);
            }

            inline const array_type&
            coefficients() const noexcept {
                return this->_coefficients;
            }

            inline array_type&
            coefficients() noexcept {
                return this->_coefficients;
            }

            inline const value_type&
            coefficients(int i) const noexcept {
                return this->_coefficients(i);
            }

            inline value_type&
            coefficients(int i) noexcept {
                return this->_coefficients(i);
            }

            inline complex_polynomial
            reciprocal() const {
                complex_polynomial result(this->size());
                const int n = this->size();
                for (int i=0; i<n; ++i) {
                    result[i] = std::conj((*this)[n-1-i]);
                }
                return result;
            }

            friend inline Polynomial
            operator*(const Polynomial& lhs, const T& rhs) {
                return Polynomial(lhs._coefficients * rhs);
            }

            friend inline Polynomial
            operator*(const T& lhs, const Polynomial& rhs) {
                return Polynomial(lhs * rhs._coefficients);
            }

            friend inline Polynomial
            operator/(const Polynomial& lhs, const T& rhs) {
                return Polynomial(lhs._coefficients / rhs);
            }

            friend inline Polynomial
            operator+(const Polynomial& lhs, const Polynomial& rhs) {
                using blitz::Range;
                Polynomial result;
                if (lhs.size() < rhs.size()) {
                    result = rhs;
                    Range r{0,lhs.size()-1};
                    result._coefficients(r) += lhs._coefficients;
                } else {
                    result = lhs;
                    Range r{0,rhs.size()-1};
                    result._coefficients(r) += rhs._coefficients;
                }
                return result;
            }

            friend inline Polynomial
            operator-(const Polynomial& lhs, const Polynomial& rhs) {
                using blitz::Range;
                Polynomial result;
                if (lhs.size() < rhs.size()) {
                    result = rhs;
                    Range r{0,lhs.size()-1};
                    result._coefficients(r) -= lhs._coefficients;
                } else {
                    result = lhs;
                    Range r{0,rhs.size()-1};
                    result._coefficients(r) -= rhs._coefficients;
                }
                return result;
            }

        };

        template <class T>
        inline void
        swap(Polynomial<T>& lhs, Polynomial<T>& rhs) {
            lhs.swap(rhs);
        }

        template <class T>
        inline auto
        schur_transform(const Polynomial<T>& p)
        -> typename std::enable_if<is_complex<T>::value, Polynomial<T>>::type {
            const int n = p.size();
            Polynomial<T> result(n-1);
            const T conj_a_0 = std::conj(p[0]);
            const T a_n = p[n-1];
            for (int i=0; i<n-1; ++i) {
                result[i] = conj_a_0*p[i] - a_n*std::conj(p[n-1-i]);
            }
            return result;
        }

        template <class T>
        inline auto
        schur_transform(const Polynomial<T>& p)
        -> typename std::enable_if<!is_complex<T>::value, Polynomial<T>>::type {
            const int n = p.size();
            Polynomial<T> result(n-1);
            const T a_0 = p[0];
            const T a_n = p[n-1];
            for (int i=0; i<n-1; ++i) {
                result[i] = a_0*p[i] - a_n*p[n-1-i];
            }
            return result;
        }

        template <class T>
        inline int
        num_roots_inside_unit_disk(Polynomial<T> p) {
            using blitz::all;
            int n = p.size();
            typedef decltype(std::real(p[0])) Real;
            blitz::Array<Real,1> mask(n);
            int negative_degree = 0;
            bool negative_multiple_times = false;
            for (int i=0; i<n-1; ++i) {
                p = schur_transform<T>(p);
                Real delta = std::real(p[0]);
                mask(i) = delta;
                #if defined(VTB_DEBUG_SCHUR_COHN)
                std::clog << "delta=" << delta << std::endl;
                #endif
                if (delta < 0) {
                    if (negative_degree == 0) {
                        negative_degree = p.size();
                    } else {
                        negative_multiple_times = true;
                    }
                }
            }
            if (negative_degree == 0) {
                return 0;
            }
            if (!negative_multiple_times) {
                return negative_degree;
            }
            return -1;
        }

    }

}

#endif // vim:filetype=cpp