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

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




//  lie_basis.h


// Include once wrapper
#ifndef DJC_COROPA_LIBALGEBRA_LIEBASISH_SEEN
#define DJC_COROPA_LIBALGEBRA_LIEBASISH_SEEN


class hall_basis
{
public:
        typedef DEG KEY;
        typedef std::pair<KEY, KEY> PARENT;
protected:
        std::vector<PARENT> hall_set;
        std::map<PARENT, KEY> reverse_map;
        std::vector<DEG> degrees;
        std::string letters;
        std::map<LET, KEY> ltk;
        DEG curr_degree;
public:
        hall_basis(DEG n_letters)
                : curr_degree(0)
        {
                // We put something at position 0 to make indexing more logical
                degrees.push_back(0);
                PARENT p(0,0);
                hall_set.push_back(p);
        
                for (LET c = 1; c <= n_letters; ++c)
                        letters += (char) c;
        
                // We add the letters to the basis, starting from position 1.
                std::string::size_type i;
                for (i = 1; i <= letters.size(); ++i)
                {
                        PARENT parents(0,i);
                        hall_set.push_back(parents); // at [i]
                        degrees.push_back(1); // at [i]
                        ltk[letters[i - 1]] = (LET) i;
                }
                curr_degree = 1;
                // To construct the rest of the basis now, call growup(max_degree) here.
        }

        inline void growup(DEG desired_degree)
        {
                for (DEG d = curr_degree + 1; d <= desired_degree; ++d)
                {
                        KEY bound = (KEY)hall_set.size();
                        for (KEY i = 1; i <= bound; ++i)
                                for (KEY j = i + 1; j <= bound; ++j)
                                        if ((degrees[i] + degrees[j] == d) && (hall_set[j].first <= i))
                                        {
                                                PARENT parents(i, j);
                                                hall_set.push_back(parents);  // at position max_key.
                                                degrees.push_back(d);         // at position max_key.
                                                reverse_map[parents] = (KEY) hall_set.size() - 1;
                                        }
                        ++curr_degree;
                }
        }
        inline DEG degree(const KEY& k) const
        {
                return degrees[k];
        }
        inline KEY keyofletter(LET letter) const
        {
                return ltk.find(letter)->second;
        }
        inline KEY lparent(const KEY& k) const
        {
                return hall_set[k].first;
        }
        inline KEY rparent(const KEY& k) const
        {
                return hall_set[k].second;
        }
        inline bool letter(const KEY& k) const
        {
                return ((k > 0) && (k <= letters.size()));
        }
        inline LET getletter(const KEY& k) const
        {
                return letters[k - 1];
        }
        inline KEY begin(void) const
        {
                return 1;
        }
        inline KEY end(void) const
        {
                return 0;
        }
        inline KEY nextkey(const KEY& k) const
        {
                if (k < (hall_set.size() - 1))
                        return (k + 1);
                else
                        return 0;
        }
        inline DEG keypos(const KEY& k) const
        {
                return k;
        }
        inline DEG size(void) const
        {
                return ( (KEY) hall_set.size() - 1);
        }
        inline friend std::ostream& operator<<(std::ostream& os, hall_basis& b)
        {       
                for (KEY k = b.begin(); k != b.end(); k = b.nextkey(k))
                        os << b.key2string(k) << ' ';
                return os;
        }

        
        inline const std::string& key2string(const KEY& k) const
        {
                static boost::recursive_mutex table_access;
                // get exclusive recursive access for the thread 
                boost::lock_guard<boost::recursive_mutex> lock(table_access); 

                static std::vector<std::string> table;

                if (k > table.size())
                {
                        for (KEY i = (KEY) table.size() + 1; i <= k; ++i)
                                table.push_back(_key2string(i));
                }
                return table[k - 1];
        }
private:
        std::string _key2string(const KEY& k) const
        {
                std::ostringstream oss;
                if (k > 0)
                {
                        if (letter(k))
                                oss << getletter(k);
                        else
                        {
                                oss << '[';
                                oss << key2string(lparent(k));
                                oss << ',';
                                oss << key2string(rparent(k));
                                oss << ']';
                        }
                }
                return oss.str();
        }
};


template<typename SCA, typename RAT, DEG n_letters, DEG max_degree>
class lie_basis : public hall_basis,
                                  public basis_traits<With_Degree, n_letters, max_degree>
{
public:
        typedef hall_basis::KEY KEY;
        typedef std::map<KEY, SCA> MAP;
        typedef lie<SCA, RAT, n_letters, max_degree> LIE;
        typedef RAT RATIONAL;
public:
        lie_basis(void)
                : hall_basis(n_letters) {}

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

                static std::map<PARENT, LIE> table;
                static typename std::map<PARENT, LIE>::iterator it;
                PARENT p(k1, k2);
                it = table.find(p);
                if (it == table.end())
                        return table[p] = _prod(k1, k2);
                else
                        return it->second;
        }

        LIE replace(const KEY& k, const std::vector<LET>& s, const std::vector<LIE*>& v, std::map<KEY, LIE>& table)
        {
                typename std::map<KEY, LIE>::iterator it;
                it = table.find(k);
                if (it != table.end())
                        return it->second;
                else
                {
                        if (letter(k))
                        {
                                typename std::vector<LET>::size_type i;
                                for (i = 0; i < s.size(); ++i)
                                        if (s[i] == getletter(k))
                                                return table[k] = *(v[i]);
                                return (table[k] = (LIE)k);
                        }
                        else
                                return (table[k]
                                                = replace(lparent(k), s, v, table)
                                                * replace(rparent(k), s, v, table));
                }
        }
private:
        LIE _prod(const KEY& k1, const KEY& k2)
        {       
                static LIE empty;
                // [A,A] = 0.
                if (k1 == k2)
                        return empty;
                // if index(A) > index(B) we use [A,B] = -[B,A] 
                if (k1 > k2)
                        return -prod(k2, k1);
                //
                DEG target_degree = degrees[k1] + degrees[k2];
                if ((max_degree > 0) && (target_degree > max_degree))
                        return empty; // degree truncation
                // We grow up the basis up to the desired degree.
                growup(target_degree);
                // We look up for the desired product in our basis.
                PARENT parents(k1, k2);
                typename std::map<PARENT, KEY>::const_iterator it;
                it = reverse_map.find(parents);
                if (it != reverse_map.end())
                {
                        // [k1,k2] exists in the basis.
                        LIE result(it->second);
                        return result;
                }
                else
                        // [k1,k2] does not exists in the basis.
                {
                        // Since k1 <= k2, k2 is not a letter because if it was a letter, 
                        // then also k1, which is impossible since [k1,k2] is not in the basis.
                        // We use Jacobi: [k1,k2] = [k1,[k3,k4]]] = [[k1,k3],k4]-[[k1,k4],k3] 
                        KEY k3(lparent (k2));
                        KEY k4(rparent (k2));
                        LIE result(prod(k1, k3) * (LIE)k4);
                        return result.sub_mul(prod(k1, k4), (LIE)k3);
                }
        }
        inline friend std::ostream& operator<<(std::ostream& os, const std::pair<lie_basis*, KEY>& t)
        {
                return os << t.first->key2string(t.second);
        }
};


// Include once wrapper
#endif // DJC_COROPA_LIBALGEBRA_LIEBASISH_SEEN

//EOF.