:github_url: https://github.com/svenevs/exhale-companion


.. _program_listing_file_ellip_monomial.hpp:

Program Listing for File monomial.hpp
=====================================

|exhale_lsh| :ref:`Return to documentation for file <file_ellip_monomial.hpp>` (``ellip/monomial.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #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