libalgebra-demo/libalgebra/polynomials.h

/* *************************************************************

Copyright 2010 Terry Lyons, Stephen Buckley, Djalil Chafai, 
Greg Gyurkó and Arend Janssen. 

Distributed under the terms of the GNU General Public License, 
Version 3. (See accompanying file License.txt)

************************************************************* */




//  polynomials.h


// Include once wrapper
#ifndef DJC_COROPA_LIBALGEBRA_POLYNOMIALSH_SEEN
#define DJC_COROPA_LIBALGEBRA_POLYNOMIALSH_SEEN


template<typename SCA, typename RAT>
class poly :
public algebra<poly_basis<SCA, RAT> >
{
public:
        typedef poly_basis<SCA, RAT> BASIS;
        typedef typename BASIS::KEY KEY;
        typedef sparse_vector<BASIS> VECT;
        typedef algebra<BASIS> ALG;
        typedef typename ALG::iterator iterator;
        typedef typename ALG::const_iterator const_iterator;
public:
  poly(void) {}
  poly(const poly& p) : ALG(p) {}
  poly(const ALG& a) : ALG(a) {}        
  poly(const VECT& v) : ALG(v) {}       
  explicit poly(const SCA& s) : ALG(poly::basis.empty_key, s) {}
  explicit poly(const KEY& k) : ALG(k) {}
  explicit poly(LET letter, const SCA& s)
        : ALG(VECT::basis.keyofletter(letter), s) {}
public:
  inline __DECLARE_BINARY_OPERATOR(poly,*,*=,SCA)
  inline __DECLARE_BINARY_OPERATOR(poly,/,/=,RAT)
  inline __DECLARE_BINARY_OPERATOR(poly,*,*=,poly)
  inline __DECLARE_BINARY_OPERATOR(poly,+,+=,poly)
  inline __DECLARE_BINARY_OPERATOR(poly,-,-=,poly)
  inline __DECLARE_UNARY_OPERATOR(poly,-,-,ALG)
  inline SCA eval(const std::map<LET, SCA>& values) const
        {
                SCA result(VECT::zero);
                for (const_iterator i = VECT::begin(); i != VECT::end(); ++i)
                        result += VECT::basis.eval_key(i->first, values) * i->second;
                return result;
        }

public:
  inline static poly prediff(const KEY& k1, const typename LET& k2)
        {
                typename KEY::iterator it;
                typename
                KEY k(k1);
                it = k.find(k2);
                poly result; //zero
                if (it != k.end())
                {
                        if (it->second == 1)
                        {
                                k.erase(it);
                                result = poly(k);
                        }
                        else
                        {
                                SCA coeff = (it->second)--;
                                poly temp1(k);
                                poly temp2(coeff);
                                result = temp2 * temp1;
                        }
                }
                return result;
        }
public:
                
  inline static poly diff(const poly& p1, const typename LET& k2)
        {
                poly result;
                const_iterator it;
                for (it = p1.begin(); it != p1.end(); ++it)
                {
                        result += poly(it->second) * prediff(it->first, k2);
                }
                return result;
        }

  inline friend poly exp(const poly& arg, DEG max_degree = 3)
        {
                static KEY kunit;
                poly result(kunit);
                for (DEG i = max_degree; i >= 1; --i)
                {
                        result.mul_scal_div(arg, (RAT)i);
                        result += (poly)kunit;
                }
                return result;
        }

  inline friend poly log(const poly& arg, DEG max_degree = 3)
        {
                static KEY kunit;
                poly tunit(kunit);
                poly x(arg);
                iterator it = x.find(kunit);
                if (it != x.end())
                        x.erase(it);
                poly result;
                for (DEG i = max_degree; i >= 1; --i)
                {
                        if (i % 2 == 0)
                                result.sub_scal_div(tunit, (RAT)i);
                        else
                                result.add_scal_div(tunit, (RAT)i);
                        result *= x;
                }
                return result;
        }
};

// Include once wrapper
#endif // DJC_COROPA_LIBALGEBRA_POLYNOMIALSH_SEEN

//EOF.