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)

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



#pragma once
#include "implimentation_types.h"
#include "constlog2.h"

template <unsigned No_Letters, unsigned DEPTH>
class
_tensor_basis
{

private:

        _tensor_basis(const double base)
                                : _word(base)
        {
        }

        double _word;

        static const unsigned uBitsInLetter = ConstLog2 < No_Letters - 1 > ::ans + 1;
        static const unsigned uMaxSizeAlphabet = (1 << uBitsInLetter);
        static const unsigned uMaxWordLength = 52 / uBitsInLetter;


public:

        typedef alg::LET LET;
        

        _tensor_basis(void)
                : _word(1.)
        {
                assert(DEPTH <= uMaxWordLength);
        }

        ~_tensor_basis(void)
        {
        }

        inline _tensor_basis& push_back(const _tensor_basis& rhs)
        {
                int iExponent;
                frexp(rhs._word, &iExponent);
                double dPowerOfTwo = ldexp(.5, iExponent);
                _word = (_word * dPowerOfTwo + rhs._word) - dPowerOfTwo;
                return *this;
        }

        inline _tensor_basis operator* (const _tensor_basis& rhs) const
        {
                int iExponent;
                frexp(rhs._word, &iExponent);
                double dPowerOfTwo = ldexp(.5, iExponent);
                return (_word * dPowerOfTwo + rhs._word) - dPowerOfTwo;
        }


        inline bool operator < (const _tensor_basis & rhs) const
        {
                assert(size() <= DEPTH);
                return _word < rhs._word;
        }

        _tensor_basis(const LET uLetter)
                : _word (uMaxSizeAlphabet + (uLetter - 1)%uMaxSizeAlphabet)
        {
                assert(uMaxSizeAlphabet >= uLetter && uLetter > 0);
        }

        inline unsigned size () const
        {
                int iExponent;
                frexp(_word, &iExponent);
                assert((iExponent - 1) % uBitsInLetter == 0);
                return (iExponent - 1) / uBitsInLetter;

        }

        inline LET FirstLetter() const
        {
                static const double dShiftPlus1(uMaxSizeAlphabet * 2);
                static const double dShift(uMaxSizeAlphabet);
                static const double dMinusShift = 1 / dShift;

                assert(size() > 0);
                int iExponent;
                double dMantissa = frexp(_word, &iExponent);
                double ans;
                modf(dMantissa * dShiftPlus1, &ans);
                return LET(ans - dShift) + 1;
        }

        inline _tensor_basis lparent() const
        {
                static const double dShiftPlus1(uMaxSizeAlphabet * 2);
                static const double dShift(uMaxSizeAlphabet);
                static const double dMinusShift = 1 / dShift;

                assert(size() > 0);
                int iExponent;
                double dMantissa = frexp(_word, &iExponent);
                double ans;
                modf(dMantissa * dShiftPlus1, &ans);
                return ans;
        }
        
        inline _tensor_basis rparent() const
        {
                static const double dShiftPlus1(uMaxSizeAlphabet * 2);
                static const double dShift(uMaxSizeAlphabet);
                static const double dMinusShift = 1 / dShift;

                assert(size() > 0);
                int iExponent;
                double dMantissa = frexp(_word, &iExponent);
                double ans;
                double dPowerOfTwo = ldexp(.5, iExponent - uBitsInLetter);
                return (modf(dMantissa * dShiftPlus1, &ans) + 1.) * dPowerOfTwo;
        }

        static _tensor_basis end()
        {
                return std::numeric_limits<double>::infinity();
        }

        friend std::ostream & operator << (std::ostream & os, const typename _tensor_basis < No_Letters, DEPTH > & word)
        {
                int iExponent;
                unsigned count = word.size();
                double dNormalised = frexp(word._word, &iExponent) * 2. - 1.;
                os << "(";
                while (count > 0)
                {
                        double letter;
                        dNormalised = modf(dNormalised * uMaxSizeAlphabet, &letter);
                        os << letter + 1.;
                        --count;
                        if (count != 0)
                                os << ",";
                }
                return os << ")";
        }
};