Program Listing for File monomial.hpp

↰ Return to documentation for file (ellip/monomial.hpp)

#ifndef _MONOMIAL_HPP
#define _MONOMIAL_HPP

#include <cassert>
#include <initializer_list>
#include <valarray>

template <typename _Tp>
class monomial
{
    using Vec = std::valarray<_Tp>;
    using _Self = monomial<_Tp>;

  public:
    explicit monomial(size_t n)
        : _b(_Tp(0))
        , _a(_Tp(0), n)
    {
    }

    monomial(const _Tp& c, const Vec& a)
        : _b(log(c))
        , _a(a)
    {
    }

    monomial(const _Tp& c, std::initializer_list<_Tp> lst)
        : _b(log(c))
        , _a(lst)
    {
    }

    monomial(size_t n, const _Tp& c)
        : _b(log(c))
        , _a(_Tp(0), n)
    {
    }

    template <typename _Up, class Map>
    monomial(const monomial<_Up>& mon, const Map& polarity);

    monomial(size_t n, const _Tp ar[])
        : _b(log(ar[0]))
        , _a(_Tp(0), n)
    {
        for (auto i = 1; i <= n; ++i)
            _a[i - 1] = ar[i];
    }

    ~monomial() { }

    _Self& operator*=(const _Self& M)
    {
        _a += M._a;
        _b += M._b;
        return *this;
    }

    _Self& operator/=(const _Self& M)
    {
        _a -= M._a;
        _b -= M._b;
        return *this;
    }

    _Self& operator*=(const _Tp& c)
    {
        _b += log(c);
        return *this;
    }

    _Self& operator/=(const _Tp& c)
    {
        _b -= log(c);
        return *this;
    }


    _Self operator*(const _Tp& c) const
    {
        return _Self(*this) *= c;
    }

    _Self operator/(const _Tp& c) const
    {
        return _Self(*this) /= c;
    }

    void sqrt()
    {
        _a /= 2.;
        _b /= 2.;
    }

    template <typename _Up>
    _Tp lse(const std::valarray<_Up>& y) const
    {
        _Tp res = _b;
        for (auto i = 0; i != _a.size(); ++i)
        {
            if (_a[i] == _Tp(0))
                continue;
            res += _a[i] * y[i];
        }
        return res;
    }

    template <typename _Up>
    Vec lse_gradient(const _Up&) const
    {
        return _a;
    }

    template <typename _Up>
    _Tp log_exp_fvalue_with_gradient(const _Up& y, Vec& g) const
    {
        assert(_a.size() == y.size());
        _Tp res = _b;
        for (auto i = 0; i != _a.size(); ++i)
        {
            if (_a[i] == _Tp(0))
                continue;
            res += _a[i] * y[i];
        }
        g = _a; // gradient
        return res;
    }

    template <typename _Up>
    Vec lse_gradient(const _Up& y, _Tp& f) const
    {
        assert(_a.size() == y.size());
        f = lse(y);
        return _a;
    }

    void set_coeff(_Tp c)
    {
        // assert(c > 0);
        _b = log(c);
    }

  public:
    _Tp _b;
    Vec _a;
};


// Non-member functions

template <typename _Tp>
inline monomial<_Tp> sqrt(monomial<_Tp> m)
{
    m.sqrt();
    return m;
}


template <typename _Tp, typename _Up>
inline monomial<_Tp> operator/(const _Up& c, const monomial<_Tp>& m)
{
    return monomial<_Tp>(m._a.size(), _Tp(c)) / m;
}

template <typename _Tp, typename _Up>
inline monomial<_Tp> operator*(const _Up& c, const monomial<_Tp>& m)
{
    return m * _Tp(c); // assume commutative
}

template <typename _Tp>
inline monomial<_Tp> operator*(monomial<_Tp> lhs, const monomial<_Tp>& rhs)
{
    lhs *= rhs;
    return lhs;
}

template <typename _Tp>
inline monomial<_Tp> operator/(monomial<_Tp> lhs, const monomial<_Tp>& rhs)
{
    lhs /= rhs;
    return lhs;
}


#endif