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

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




//  algebra.h


// Include once wrapper
#ifndef DJC_COROPA_LIBALGEBRA_ALGEBRAH_SEEN
#define DJC_COROPA_LIBALGEBRA_ALGEBRAH_SEEN


template<class BASIS>
class algebra : public sparse_vector<BASIS>
{
public:
        typedef sparse_vector<BASIS> VECT;
        typedef typename VECT::iterator iterator;
        typedef typename VECT::const_iterator const_iterator;
        typedef typename VECT::KEY KEY;
        typedef typename VECT::SCALAR SCALAR;   
        typedef typename VECT::RATIONAL RATIONAL;       

        static const DEG MAX_DEGREE = BASIS::MAX_DEGREE;

        template <class Transform, unsigned DEPTH1>
inline void triangularbufferedmultiplyandcombine(const algebra& rhs, algebra& result, Transform fn) const
{
        // create buffer to avoid unnecessary calls to MAP inside loop
        std::vector<std::pair<KEY, SCALAR> > buffer(rhs.begin(), rhs.end());
        const_iterator i;

        // segment buffer according to the degree of the elements
        std::vector<std::vector<std::pair<KEY, SCALAR> >:: const_iterator>
                iterators(DEPTH1 + 1,buffer.end());
        unsigned deg = 0;
        for (std::vector<std::pair<KEY, SCALAR> >:: const_iterator j0 = buffer.begin(); j0 != buffer.end(); j0++)
        {
                DEG d = basis.degree(j0->first);
                assert(d >= deg && d <= DEPTH1); // order assumed to respect degree
                while (deg < d)
                        iterators[deg++] = j0;
                // deg == d
        }

        std::vector<std::pair<KEY, SCALAR> >:: const_iterator j, jEnd;
        for (i = VECT::begin(); i != VECT::end(); ++i)
        {
                const KEY& k = i->first;
                unsigned rhdegree = DEPTH1 - basis.degree(k);
                std::vector<std::pair<KEY, SCALAR> >:: const_iterator&
                        jEnd = iterators[rhdegree];
                for (j = buffer.begin(); j != jEnd; ++j)
                        result.add_scal_prod(VECT::basis.prod(i->first, j->first),
                                fn(i->second * j->second));
        }

}

        template <class Transform>
inline void squarebufferedmultiplyandcombine(const algebra& rhs, algebra& result, Transform fn) const
{       
        // create buffer to avoid unnecessary calls to MAP inside loop
        std::vector<std::pair<KEY, SCALAR> > buffer(rhs.begin(), rhs.end());
        const_iterator i;

        // DEPTH1 == 0
        std::vector<std::pair<KEY, SCALAR> >:: const_iterator j;
        for (i = VECT::begin(); i != VECT::end(); ++i)
        {
                for (j = buffer.begin(); j != buffer.end(); ++j)
                        result.add_scal_prod(VECT::basis.prod(i->first, j->first),
                                fn(i->second * j->second));
        }

}

        // Transformations

        struct scalar_passthrough
        {
                SCALAR operator()(const SCALAR& arg)
                {
                        return arg;
                }
        };
        struct scalar_minus
        {
                SCALAR operator()(const SCALAR& arg)
                {
                        return -arg;
                }
        };
        struct scalar_post_mult
        {
        private:
                SCALAR mFactor;
        public:
                scalar_post_mult(const SCALAR Factor = one)
                        : mFactor(Factor) {}
                SCALAR operator()(const SCALAR arg)
                {
                        return arg * mFactor;
                }
        };
        struct rational_post_div
        {
        private:
                SCALAR mFactor;
        public:
                rational_post_div(const RATIONAL Factor = one)
                        : mFactor(one / Factor) {}
                SCALAR operator()(const SCALAR arg)
                {
                        return arg * mFactor;
                }
        };
        template <unsigned DEPTH1>
inline void bufferedmultiplyandadd(const algebra& rhs, algebra& result) const
{
        scalar_passthrough fn;
        triangularbufferedmultiplyandcombine<scalar_passthrough, DEPTH1>(rhs, result, fn);
}
        template <>
        inline void bufferedmultiplyandadd<0>(const algebra& rhs, algebra& result) const
        {
                scalar_passthrough fn;
                squarebufferedmultiplyandcombine(rhs, result, fn);
        }
        template <unsigned DEPTH1>
inline void bufferedmultiplyandsub(const algebra& rhs, algebra& result) const
{
        scalar_minus fn;
        triangularbufferedmultiplyandcombine<scalar_minus, DEPTH1>(rhs, result, fn);
}

        template <>
        inline void bufferedmultiplyandsub<0>(const algebra& rhs, algebra& result) const
        {
                scalar_minus fn;
                squarebufferedmultiplyandcombine(rhs, result, fn);
        }

        struct wrapscalar
        {
                const SCALAR& hidden;
                wrapscalar(const SCALAR& s)
                        : hidden(s) {}
        };
        struct wraprational
        {
                const RATIONAL& hidden;
                wraprational(const RATIONAL& s)
                        : hidden(s) {}
        };
        template <unsigned DEPTH1>
inline void bufferedmultiplyandsmult(const wrapscalar& ss, const algebra& rhs, algebra& result) const
{
        scalar_post_mult fn(ss.hidden);
        triangularbufferedmultiplyandcombine<scalar_post_mult, DEPTH1>(rhs, result, fn);
}

        template <>
        inline void bufferedmultiplyandsmult<0>(const wrapscalar& ss, const algebra& rhs, algebra& result) const
        {
                scalar_post_mult fn(ss.hidden);
                squarebufferedmultiplyandcombine(rhs, result, fn);
        }

        template <unsigned DEPTH1>
inline void bufferedmultiplyandsdiv(const algebra& rhs, const wraprational& ss, algebra& result) const
{
        rational_post_div fn(ss.hidden);
        triangularbufferedmultiplyandcombine<rational_post_div, DEPTH1>(rhs, result, fn);
}

        template <>
        inline void bufferedmultiplyandsdiv<0>(const algebra& rhs, const wraprational& ss, algebra& result) const
        {
                rational_post_div fn(ss.hidden);
                squarebufferedmultiplyandcombine(rhs, result, fn);
        }

public:

        algebra(void) {}
        algebra(const algebra& a)
                : VECT(a) {}
        algebra(const VECT& v)
                : VECT(v) {}
        explicit algebra(const KEY& k, const SCALAR& s = VECT::one)
                : VECT(k, s) {}
public: 
        inline algebra& operator*=(const SCALAR& s)
        {
                VECT::operator *= (s);
                return *this;
        }
        inline algebra& operator/=(const RATIONAL& s)
        {
                VECT::operator /= (s);
                return *this;
        }
        inline __DECLARE_BINARY_OPERATOR(algebra,*,*=,SCALAR)
                inline __DECLARE_BINARY_OPERATOR(algebra,/,/=,SCALAR)
                inline __DECLARE_BINARY_OPERATOR(algebra,+,+=,algebra)
                inline __DECLARE_BINARY_OPERATOR(algebra,-,-=,algebra)
                inline __DECLARE_UNARY_OPERATOR(algebra,-,-,VECT);
        inline algebra& operator*=(const algebra& rhs)
        {
                algebra result;
                bufferedmultiplyandadd<MAX_DEGREE>(rhs, result);
                swap(result);
                return *this;
        }
        inline __DECLARE_BINARY_OPERATOR(algebra,*,*=,algebra);
        inline algebra& add_mul(const algebra& a, const algebra& b)
        {
                a.bufferedmultiplyandadd<MAX_DEGREE>(b, *this);
                return *this;
        }
        inline algebra& sub_mul(const algebra& a, const algebra& b)
        {
                a.bufferedmultiplyandsub<MAX_DEGREE>(b, *this);
                return *this;
        }
        inline algebra& mul_scal_prod(const algebra& rhs, const SCALAR& s)
        {
                algebra result;
                bufferedmultiplyandsmult(rhs, wrapscalar(s), result);
                swap(result);
                return *this;
        }
        inline algebra& mul_scal_div(const algebra& rhs, const RATIONAL& s)
        {
                algebra result;
                bufferedmultiplyandsdiv<MAX_DEGREE>(rhs, wraprational(s), result);
                swap(result);
                return *this;
        }
        inline friend algebra commutator(const algebra& a, const algebra& b)
        { // Returns a * b - b * a
                algebra result;
                a.bufferedmultiplyandadd<MAX_DEGREE>(b, result);
                b.bufferedmultiplyandsub<MAX_DEGREE>(a, result);
                return result;
        }
        inline algebra truncate(const DEG min, const DEG max) const
        {
                algebra result;
                const_iterator i;
                for (i = VECT::begin(); i != VECT::end(); ++i)
                        if ((VECT::basis.degree(i->first) >= min) && (VECT::basis.degree(i->first) <= max))
                                result[i->first] = i->second;
                return result;
        }
        inline DEG degree(void) const
        {
                DEG result(0);
                const_iterator i;
                for (i = VECT::begin(); i != VECT::end(); ++i)
                        result = std::max(result, VECT::basis.degree(i->first));
                return result;
        }
};

// Include once wrapper
#endif // DJC_COROPA_LIBALGEBRA_ALGEBRAH_SEEN

//EOF.