libalgebra-demo/libalgebra/tensor_basis.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)

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




//  tensor_basis.h


// Include once wrapper
#ifndef DJC_COROPA_LIBALGEBRA_TENSORBASISH_SEEN
#define DJC_COROPA_LIBALGEBRA_TENSORBASISH_SEEN


#include "_tensor_basis.h"
#include "constpower.h"
#include "basis_traits.h"


template<typename SCA, DEG n_letters, DEG max_degree>
class tensor_basis
{
public:
        typedef _tensor_basis<n_letters, max_degree> KEY;
        const KEY empty_key;
        typedef std::map<KEY, SCA> MAP;
private:
        DEG _size;
public:

        tensor_basis(void)
                : _size ((ConstPower < n_letters, max_degree + 1 > ::ans - 1) / (n_letters - 1)) {}
public:
        inline KEY keyofletter(LET letter) const
        {
                return KEY(letter);
        }
        inline bool letter(const KEY& k) const
        {
                return k.size() == 1;
        }
        inline LET getletter(const KEY& k) const
        {
                return k.FirstLetter();
        }
        inline KEY lparent(const KEY& k) const
        {
                return k.lparent();
        }
        inline KEY rparent(const KEY& k) const
        {
                return k.rparent();
        }
        inline DEG degree(const KEY& k) const
        {
                return k.size();
        }
        inline DEG size(void) const
        {
                return _size;
        }
        inline KEY begin(void) const
        {
                return empty_key;
        }
        inline KEY end(void) const
        {
                KEY result; // empty key.
                result.push_back(0); // invalid key.    
                return result;
        }
        inline KEY nextkey(const KEY& k) const
        {
                KEY::size_type i;
                for (i = k.size() - 1; i >= 0; --i)
                        if (k[i] < n_letters)
                        {
                                KEY result(k);
                                result[i] += 1;
                                return result;
                        }
                return end();
        }
        std::string key2string(const KEY& k) const
        {
                std::ostringstream oss;

                if (k.size() > 0)
                {
                        oss << k.FirstLetter();
                        KEY kk = k.rparent();
                        for (unsigned i = 1; i < k.size(); ++i)
                        {
                                oss << "," << kk.FirstLetter();
                                kk = kk.rparent();
                        }
                }
                return oss.str();
        }
};


template<typename SCA, typename RAT, DEG n_letters, DEG max_degree>
class free_tensor_basis : public tensor_basis<SCA, n_letters, max_degree>,
                                                  public basis_traits<With_Degree, n_letters, max_degree>
{
public:
        typedef tensor_basis<SCA, n_letters, max_degree> TBASIS;
        typedef RAT RATIONAL;
        typedef free_tensor<SCA, RAT, n_letters, max_degree> TENSOR;
public:
        free_tensor_basis(void) {}
public:

        static inline TENSOR prod(const KEY& k1, const KEY& k2)
        {
                static SCA one(+1);
                TENSOR result;
                if ((max_degree == 0) || (k1.size() + k2.size() <= max_degree))
                {
                        KEY concat = k1 * k2;
                        result[concat] = one;
                }
                return result;
        }
        inline friend
                std::ostream& operator<<(std::ostream& os, const std::pair<free_tensor_basis*, KEY>& t)
        {
                return os << (t.first)->key2string(t.second);
        }
};


template<typename SCA, typename RAT, DEG n_letters, DEG max_degree>
class shuffle_tensor_basis : public tensor_basis<SCA, n_letters, max_degree>,
                                                         public basis_traits<With_Degree, n_letters, max_degree>
{
public:
        typedef tensor_basis<SCA, n_letters, max_degree> TBASIS;
        typedef typename TBASIS::KEY KEY;
        typedef typename TBASIS::MAP MAP;
        typedef RAT RATIONAL;
        typedef shuffle_tensor<SCA, RAT, n_letters, max_degree> TENSOR;
public:
        shuffle_tensor_basis(void) {}
public:

        inline const TENSOR& prod(const KEY& k1, const KEY& k2) const
        {               
                static boost::recursive_mutex table_access;
                // get exclusive recursive access for the thread 
                boost::lock_guard<boost::recursive_mutex> lock(table_access); 

                typedef std::map<std::pair<KEY, KEY>, TENSOR> TABLE_T;
                static TABLE_T table;
                typename TABLE_T::iterator it;
                std::pair<KEY, KEY> p(std::min(k1, k2), std::max(k1, k2));
                it = table.find(p);
                if (it == table.end())
                        return table[p] = _prod(k1, k2);
                else
                        return it->second;
        }
        inline friend
                std::ostream& operator<<(std::ostream& os, const std::pair<shuffle_tensor_basis*, KEY>& t)
        {
                return os << (t.first)->key2string(t.second);
        }
private:
        
        inline void tensor_add_mul(const TENSOR& a, const TENSOR& b, TENSOR& ans) const
        {
                sparse_vector<typename free_tensor_basis<SCA, RAT, n_letters, max_degree>,MAP>
                ::const_iterator i, j;
                for (i = a.begin(); i != a.end(); ++i)
                        for (j = b.begin(); j != b.end(); ++j)
                                ans.add_scal_prod(free_tensor_prod(i->first, j->first),
                                        i->second * j->second);
        }

        static inline TENSOR free_tensor_prod(const KEY& k1, const KEY& k2)
        {
                static SCA one(+1);
                TENSOR result;
                if ((max_degree == 0) || (k1.size() + k2.size() <= max_degree))
                {
                        KEY concat = k1 * k2;
                        result[concat] = one;
                }
                return result;
        }

        inline TENSOR _prod(const KEY& k1, const KEY& k2) const
        {
                TENSOR result;
                unsigned i, j;
                static SCA one(+1);

                if ((max_degree == 0) || (k1.size() + k2.size() <= max_degree))

                {
                        if (k1.size() == 0)
                        {
                                result[k2] = one;
                                return result;
                        }
                        if (k2.size() == 0)
                        {
                                result[k1] = one;
                                return result;
                        }
                        // k1.size() >= 1 and k2.size() >= 1
                        tensor_add_mul((TENSOR)k1.lparent(), prod(k1.rparent(), k2), result);
                        tensor_add_mul((TENSOR)k2.lparent(), prod(k1, k2.rparent()), result);
                }
                return result;
        }
};
//#endif


// Include once wrapper
#endif // DJC_COROPA_LIBALGEBRA_TENSORBASISH_SEEN

//EOF.