MultiIndex.hpp

#ifndef __Fiber_MultiIndex_HPP
#define __Fiber_MultiIndex_HPP
 
#include <eagle/vecmath.hpp>
#include <eagle/Assignment.hpp>
#include <eagle/Operator.hpp>
#include <eagle/FixedArray.hpp>
#include <eagle/DomainVector.hpp>
 
#include "Index.hpp"
 
#include <array>
#include <initializer_list>
 
namespace Fiber
{
 
using   Eagle::FixedArray;
 
using   Eagle::Mult;
using   Eagle::Div;
using   Eagle::Add;
using   Eagle::Sub;
 
struct  CreateFromLinear {};
 
struct  Power2Alignment {};
 
using   Dims_t = int;
 
template <Dims_t Dims>
        class MultiIndex;
 
 
template <Dims_t Dims>
struct  MultiIndexIterator : MultiIndex<Dims>
{
        const MultiIndex<Dims> Size;
        bool                Iteratable;
 
        MultiIndexIterator(const MultiIndex<Dims>&MaxIndex)
        : Size(MaxIndex)
        {
                Iteratable = MaxIndex.size()>0;
        }
 
        MultiIndexIterator(const MultiIndex<Dims>&CurrentIndex, const MultiIndex<Dims>&MaxIndex)
        : MultiIndex<Dims>(CurrentIndex)
        , Size(MaxIndex)
        {}
 
        void    operator++()
        {
                Iteratable = this->inc(Size);
        }
 
        MultiIndex<Dims>&operator*()
        {
                return *this;
        }
 
        bool    operator!=(const MultiIndexIterator&)
        {
                return Iteratable;
        }
};
 
 
template <>
        class MultiIndex<1>
{
        index_t idx;
 
protected:
        typedef index_t Index;
        typedef Index value_type;
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Add&) noexcept
        : idx(M.idx + D.maxidx() )
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Sub&) noexcept
        : idx(M.idx - D.maxidx() )
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Mult&) noexcept
        : idx(M.idx * D.maxidx() )
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Div&) noexcept
        : idx(M.idx / D.maxidx() )
        {}
 
        constexpr MultiIndex(unsigned int bits, const ::Eagle::BinaryAnd&) noexcept
        : idx( (bits & 1) ? 1 : 0 )
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const ::Eagle::BinaryAnd&) noexcept
        : idx( M.idx > D.maxidx() ? M.idx - D.maxidx() : D.maxidx() - M.idx )
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const Power2Alignment&) noexcept
        : idx( RoundUpToNextHighestPowerOf2( M.idx) )
        {}
 
        constexpr MultiIndex(const std::array<index_t, 1>&A)
        : idx( A[0] )
        {}
 
 
public:
        constexpr MultiIndex() noexcept
        : idx(0)
        {}
 
        constexpr MultiIndex(const Index&I) noexcept
        : idx(I)
        {}
 
        constexpr MultiIndex(const MultiIndex&Midx, const Index&I) noexcept
        : idx(I)
        {}
 
        constexpr MultiIndex(const Index&lin_idx, const MultiIndex&Dimens, const CreateFromLinear&)
        : idx(lin_idx % Dimens.maxidx() )
        {}
 
        MultiIndexIterator<1> begin() const
        {
                return MultiIndexIterator<1>(*this);
        }
 
        MultiIndexIterator<1> end() const
        {
                return MultiIndexIterator<1>(*this, *this);
        }
 
 
 
static  int     log2(index_t N) noexcept
        {
                if (N&1) return  0;
                else     return -1;
        }
 
        constexpr int getSignedOrientation(const MultiIndex&B) const
        {
                if (idx<B.idx) return +1;
                if (idx>B.idx) return -1;
                return 0;
        }
 
        constexpr const Index&operator[](int) const noexcept
        {       return idx;      }
 
        constexpr Index&operator[](int) noexcept
        {       return idx;      }
 
 
        constexpr const Index&maxidx() const noexcept
        {       return idx; }
 
        constexpr Index&maxidx() noexcept
        {       return idx;             }
 
        constexpr Index size() const noexcept
        {       return idx;             }
 
 
        constexpr MultiIndex&operator+=(const MultiIndex&D) noexcept
        {
                idx += D.idx;
                return *this;
        }
 
        constexpr MultiIndex&operator-=(const MultiIndex&D) noexcept
        {
                idx -= D.idx;
                return *this;
        }
 
        constexpr MultiIndex operator+(const MultiIndex&D) noexcept
        {
                return MultiIndex(idx + D.idx);
        }
 
 
 
        friend constexpr MultiIndex clamp(const MultiIndex&M, const MultiIndex&D) noexcept
        {
        MultiIndex retval;
                retval.idx      = clamp(M.idx, D.idx);
                return retval;
        }
 
        friend constexpr MultiIndex clamp_m1(const MultiIndex&M, const MultiIndex&D) noexcept
        {
        MultiIndex retval;
                retval.idx      = clamp(M.idx, D.idx-1);
                return retval;
        }
 
 
        bool operator==(const MultiIndex&D) const noexcept
        {       return idx == D.idx;    }
 
        bool operator!=(const MultiIndex&D) const noexcept
        {
                return idx != D.maxidx();
        }
 
        bool    inc(const MultiIndex&Dimens, const MultiIndex&Increment) noexcept
        {
                idx += Increment.maxidx();
                if (idx >= Dimens.maxidx() )
                {
                        idx = Index(0);
                        return false;
                }
                return true;
        }
 
        bool    inc(const MultiIndex&Dimens) noexcept
        {
                idx++;
                if (idx >= Dimens.maxidx() )
                {
                        idx = Index(0);
                        return false;
                }
                return true;
        }
 
        Index   linear(const MultiIndex&/*Dimens*/) const
        {       return idx;                             }
 
        bool    isWithin(const MultiIndex&Range) const noexcept
        {
//              return maxidx()>=0 && maxidx() < Range.maxidx();
                return                maxidx() < Range.maxidx();
        }
 
//               MSVC                           GCC             G++                   SGI C++
#if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM) || defined(__clang__)
 
        std::ostream&print(std::ostream&os, char sep) const
        //{     return os << IntValueTrait<Index>::getValue(idx);       }
        {       return os << idx;       }
 
        friend std::ostream&operator<<(std::ostream&os, const MultiIndex&MI)
        {
                return MI.print(os, 'x');
        }
#endif
 
 
template <class Functor, Dims_t SuperDims>
static  bool ForEachRecursion(Functor&F, const MultiIndex<SuperDims>&SuperIndex,
                              const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment)
        {
                for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
                {
                MultiIndex<SuperDims+1> SupraIndex(SuperIndex, idx);
                        if (!F(SupraIndex) )
                                return false;
                }
                return true;
        }
};
 
 
template <Dims_t Dims>
        class MultiIndex : public MultiIndex<Dims-1>
{
        typedef MultiIndex<Dims-1>  base_t;
 
        index_t idx;
 
public:
        typedef index_t value_type;
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Add&) noexcept
        : base_t(M,D, Add() )
        , idx(M.idx + D.maxidx())
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Sub&) noexcept
        : base_t(M,D, Sub() )
        , idx(M.idx - D.maxidx())
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Mult&) noexcept
        : base_t(M,D, Mult() )
        , idx(M.idx * D.maxidx())
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const Div&)
        : base_t(M,D, Div() )
        , idx(M.idx / D.maxidx())
        {}
 
        constexpr MultiIndex(unsigned int bits, const ::Eagle::BinaryAnd&) noexcept
        : base_t(bits, ::Eagle::BinaryAnd() )
        , idx( bits & (1<<(Dims-1)) ? 1 : 0 )
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const MultiIndex&D, const ::Eagle::BinaryAnd&) noexcept
        : base_t(M,D, ::Eagle::BinaryAnd() )
        , idx( M.idx > D.maxidx() ? M.idx - D.maxidx() : D.maxidx() - M.idx )
        {}
 
        constexpr MultiIndex(const MultiIndex&M, const Power2Alignment&) noexcept
        : base_t(M, Power2Alignment() )
        , idx( RoundUpToNextHighestPowerOf2( M.idx) )
        {}
 
public:
        constexpr MultiIndex(const std::array<index_t, Dims>&A)
        {
                for(int i=0; i<Dims; i++)
                        (*this)[i] = A[i];
        }
 
        constexpr MultiIndex(const std::initializer_list<index_t>&A)
        {
        Dims_t  i = 0;
                for(auto I : A)
                {
                        (*this)[i] = I;
                        i++;
                        if (i>=Dims)
                                break;
                }
        }
 
        MultiIndexIterator<Dims> begin() const
        {
                return MultiIndexIterator<Dims>(*this);
        }
 
        MultiIndexIterator<Dims> end() const
        {
                return MultiIndexIterator<Dims>(*this, *this);
        }
 
 
 
public:
        friend MultiIndex distance(const MultiIndex&M0, const MultiIndex&M1)
        {
                return MultiIndex(M0,M1, ::Eagle::BinaryAnd() );
        }
 
        friend MultiIndex bitalign(const MultiIndex&M) noexcept
        {
                return MultiIndex(M, Power2Alignment() );
        }
 
static  MultiIndex BitIndex(unsigned int bits) noexcept
        {
                return MultiIndex(bits, ::Eagle::BinaryAnd() );
        }
 
static  MultiIndex Axis(unsigned int orientation) noexcept
        {
                return BitIndex( 1 << orientation );
        }
 
        index_t BitIndex() const
        {
                return linear( MultiIndex(2) );
        }
 
static  int     log2(index_t N) noexcept
        {
                if (N & (1<<(Dims-1) ) )
                    return Dims-1;
 
                return base_t::log2(N);
        }
 
        index_t Orientation() const noexcept
        {
                return log2( BitIndex() );
        }
 
        enum
        {
                dims = Dims,
 
                SIZE = Dims
        };
 
        constexpr int   getSignedOrientation(const MultiIndex&B) const
        {
                if (idx<B.idx) return +Dims;
                if (idx>B.idx) return -int(Dims);
                return subidx().getSignedOrientation(B.subidx() );
        }
 
        constexpr MultiIndex() noexcept
        : idx(0)
        {}
 
        explicit constexpr MultiIndex(const index_t&I) noexcept
        : base_t(I), idx(I)
        {}
 
        constexpr MultiIndex(const MultiIndex<Dims-1>&Midx, const index_t&I) noexcept
        : base_t(Midx), idx(I)
        {}
 
        constexpr MultiIndex(const MultiIndex<Dims-1>&m, const index_t&Slice, int SliceDirection)
        {
                for(int i=0, k=0; i<Dims; i++)
                {
                        if (i==SliceDirection)
                                (*this)[i] = Slice;
                        else
                        {
                                (*this)[i] = m[k];
                                k++;
                        }
                }
        }
 
        constexpr MultiIndex(const index_t&lin_idx, const MultiIndex&Dimens, const CreateFromLinear&)
        : base_t( lin_idx, Dimens, CreateFromLinear() )
        , idx( (lin_idx / Dimens.subidx().size()) % Dimens.maxidx() )
        {}
 
        constexpr const MultiIndex<Dims-1>&subidx() const noexcept
        {       return *this;   }
 
        constexpr MultiIndex<Dims-1>&subidx() noexcept
        {       return *this;   }
 
        constexpr const index_t&operator[](std::size_t  i) const
        {       return (i>=Dims-1) ? idx : subidx()[i];         }
 
        constexpr index_t&operator[](std::size_t i) noexcept
        {       return (i>=Dims-1) ? idx : subidx()[i];         }
 
 
        Eagle::Assignment<MultiIndex,0> operator=(const index_t&i)
        {
                return Eagle::Assignment<MultiIndex, 0>(*this, i);
        }
 
        constexpr const index_t&maxidx() const noexcept
        {       return idx; }
 
        constexpr index_t&maxidx() noexcept
        {       return idx; }
 
        index_t size() const noexcept
        {
                return idx*subidx().size();
        }
 
        constexpr bool operator!=(const MultiIndex&D) const noexcept
        {
                return idx != D.maxidx() || subidx() != D.subidx() ;
        }
 
        constexpr bool operator==(const MultiIndex&M) const noexcept
        {
                return idx == M.maxidx() && subidx() == M.subidx() ;
        }
 
        MultiIndex&operator+=(const MultiIndex&D) noexcept
        {
                idx += D.maxidx();
                subidx() += D.subidx();
                return *this;
        }
 
        MultiIndex&operator-=(const MultiIndex&D) noexcept
        {
                idx -= D.maxidx();
                subidx() -= D.subidx();
                return *this;
        }
 
/*
        friend MultiIndex clamp(const MultiIndex&M, const MultiIndex&D)
        {
        MultiIndex retval;
                retval.idx      = clamp(M.idx     , D.idx);
                retval.subidx() = clamp(M.subidx(), D.subidx());
                return retval;
        }
 
        friend MultiIndex clamp_m1(const MultiIndex&M, const MultiIndex&D)
        {
        MultiIndex retval;
                retval.idx      = clamp_m1(M.idx     , D.idx-1);
                retval.subidx() = clamp_m1(M.subidx(), D.subidx());
                return retval;
        }
*/
 
        bool    inc(const MultiIndex&Dimens) noexcept
        {
                if ( subidx().inc(Dimens) )
                        return true;
 
                idx ++;
                if ( idx >= Dimens.maxidx() )
                {
                        idx = index_t(0);
                        return false;
                }
                return true;
        }
 
        bool    inc(const MultiIndex&Dimens, const MultiIndex&Increment) noexcept
        {
                if (subidx().inc( Dimens, Increment ) )
                        return true;
 
                idx += Increment.maxidx();
                if ( idx >= Dimens.maxidx() )
                {
                        idx = index_t(0);
                        return false;
                }
                return true;
        }
 
inline  index_t linear(const MultiIndex&Dimens) const
#if     1
        ;
#else
        {
#ifndef NDEBUG
                if (!(idx >= 0 && idx < Dimens.maxidx() ))
                {
                        MemCore::breakpoint();
                        assert(idx >= 0 && idx < Dimens.maxidx() );
                }
#endif
                /*
                 Pretty sure, this formula can be optimized.
                 The size() call on the subdimension recursively multiplies
                 its elements, but this could be moved into the summation/
                 multiplication recursive call.
                */
        index_t L = Dimens.subidx().size();
                return idx * L + subidx().linear( Dimens.subidx() );
        }
#endif
 
 
//               MSVC                           GCC             G++                   SGI C++
#if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM) || defined(__clang__)
 
        std::ostream&print(std::ostream&os, char sep) const
        {
                base_t::print(os, sep);
                return os << sep << idx;
        }
 
        friend std::ostream&operator<<(std::ostream&os, const MultiIndex&MI)
        {
                return MI.print(os, 'x');
        }
#endif
 
        bool    isWithin(const MultiIndex&Range) const noexcept
        {
        // return maxidx()>=0 && maxidx() < Range.maxidx() &&
                return  maxidx() < Range.maxidx() &&
                        subidx().isWithin( Range );
        }
 
        bool    operator<(const MultiIndex&Range) const noexcept
        {
                return isWithin( Range );
        }
 
        bool    operator>(const MultiIndex&Range) const noexcept
        {
                return !isWithin( Range );
        }
 
        MultiIndex operator+(const MultiIndex&D) const noexcept
        {       return MultiIndex(*this, D, Add() );    }
 
        MultiIndex operator-(const MultiIndex&D) const noexcept
        {       return MultiIndex(*this, D, Sub() );    }
 
        MultiIndex operator*(const MultiIndex&D) const noexcept
        {       return MultiIndex(*this, D, Mult() );   }
 
        MultiIndex operator/(const MultiIndex&D) const noexcept
        {       return MultiIndex(*this, D, Div() );    }
 
 
        MultiIndex operator&(int i) const
        {
        const MultiIndex&M = *this;
                if (i==0) return M;
 
                if (i<0)
                {
                        if (M[-i-1] == 0)
                                return M;
 
                        return M - MultiIndex::Axis(-i-1);
                }
                else
                        return M + MultiIndex::Axis( i-1);
        }
 
        friend MultiIndex operator|(const MultiIndex&M, int i)
        {
                if (i==0) return M;
 
                if (i<0)
                        return M + MultiIndex::Axis(-i-1);
                else
                {
                        if (M[-i-1] == 0)
                                return M;
 
                        return M - MultiIndex::Axis( i-1);
                }
        }
 
protected:
 
template <class Functor, Dims_t SuperDims>
static  bool ForEachRecursion(Functor&F, const MultiIndex<SuperDims>&SuperIndex,
                              const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment)
        {
                for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
                {
                MultiIndex<SuperDims+1> SupraIndex(SuperIndex, idx);
                        if (!MultiIndex<Dims-1>::ForEachRecursion(F,SupraIndex, Start, End, Increment) )
                                return false;
                }
                return true;
        }
 
public:
 
template <class Functor>
static  bool ForEach(Functor&F, const MultiIndex&Start, const MultiIndex&End, const MultiIndex&Increment = MultiIndex(1) )
        {
                for(index_t idx = Start.idx; idx < End.idx; idx += Increment.idx)
                {
                MultiIndex<1> SupraIndex(idx);
                        if (!MultiIndex<Dims-1>::ForEachRecursion(F, SupraIndex, Start, End, Increment) )
                                return false;
                }
                return true;
        }
};
 
 
inline constexpr MultiIndex<1> MIndex(index_t i0)
{
        return  {i0};
//MultiIndex<1> m = i0;
//      return m;
}
 
inline constexpr MultiIndex<2> MIndex(index_t i0, index_t i1)
{
        return {i0, i1};
 
//MultiIndex<2> m; m = i0, i1;
//      return m;
}
 
inline constexpr MultiIndex<3> MIndex(index_t i0, index_t i1, index_t i2)
{
        return  {i0, i1, i2};
//MultiIndex<3> m; m = i0, i1, i2;
//      return m;
}
 
inline constexpr MultiIndex<4> MIndex(index_t i0, index_t i1, index_t i2, index_t i3)
{
        return  {i0, i1, i2, i3};
//MultiIndex<4> m; m = i0, i1, i2, i3;
//      return m;
}
 
//               MSVC                           GCC             G++                   SGI C++
#if     defined(_IOSTREAM_) || defined(_GLIBCXX_OSTREAM) || _CPP_IOSTREAM || defined(__SGI_STL_IOSTREAM) || defined(__clang__)
template <Dims_t Dims>
  inline std::ostream&operator<<(std::ostream&os, const MultiIndex<Dims>&M)
{
        return M.print( os<<'[', ',') << ']' ;
}
#endif
 
 
template <Dims_t Dims>
inline constexpr FixedArray<double, Dims> div(const MultiIndex<Dims>&Divisor, const MultiIndex<Dims>&Dividend)
{
FixedArray<double, Dims> result;
 
        for(Dims_t i=0; i<Dims; i++)
                result[i] = double(Divisor[i])/Dividend[i];
 
        return result;
}
 
 
template <class VectorType, class Domain, class scalar_t>
Eagle::DomainVector<VectorType, Domain, scalar_t>
CellSize(const Eagle::DomainVector<VectorType, Domain, scalar_t>&V,
         const MultiIndex<VectorType::SIZE>&Dividend,
         const MultiIndex<VectorType::SIZE>&Stride = MultiIndex<VectorType::SIZE>(1) )
{
Eagle::DomainVector<VectorType, Domain, scalar_t>
        result = V;
 
constexpr Dims_t DIMS = VectorType::SIZE;
 
        for(Dims_t i=0; i<DIMS; i++)
        {
        index_t divisor = Dividend[i];
 
                if (divisor > Stride[i])
                        result[i] /= divisor - Stride[i];
 
        }
 
        return result;
}
 
template <class VectorType, class Domain, class scalar_t>
Eagle::DomainVector<VectorType, Domain, scalar_t>
CellSize0(const Eagle::DomainVector<VectorType, Domain, scalar_t>&V,
          const MultiIndex<VectorType::SIZE>&Dividend,
          const MultiIndex<VectorType::SIZE>&Stride = MultiIndex<VectorType::SIZE>(1) )
{
Eagle::DomainVector<VectorType, Domain, scalar_t>
        result = V;
 
constexpr Dims_t DIMS = VectorType::SIZE;
 
        for(Dims_t i=0; i<DIMS; i++)
        {
        index_t divisor = Dividend[i];
 
                if (divisor > Stride[i])
                        result[i] /= divisor - Stride[i];
                else
                        result[i] = 0;
 
        }
 
        return result;
}
 
template <Dims_t Dims>
MultiIndex<Dims> isLargerThan(const MultiIndex<Dims>&I, const MultiIndex<Dims>&CompareValue = MultiIndex<Dims>(0),
                           const MultiIndex<Dims>&Yes = MultiIndex<Dims>(1),
                           const MultiIndex<Dims>&No  = MultiIndex<Dims>(0) )
{
MultiIndex<Dims> retval = No;
        for(Dims_t i=0; i<Dims; i++)
        {
                if (I[i]>CompareValue[i])
                        retval[i] = Yes[i];
        }
        return retval;
}
 
 
 
template <Dims_t Dims>
inline MultiIndex<Dims> operator+(index_t C, const MultiIndex<Dims>&I)
{
        return MultiIndex<Dims>(C) + I;
}
 
 
template <Dims_t Dims>
inline MultiIndex<Dims> operator+(const MultiIndex<Dims>&I, index_t C)
{
        return I + MultiIndex<Dims>(C);
}
 
template <Dims_t Dims, class Int>
inline MultiIndex<Dims> operator*(const MultiIndex<Dims>&MIn, Int Factor)
{
MultiIndex<Dims> M = MIn;
 
        for(Dims_t i=0; i<Dims; i++)
                M[i] *= Factor;
 
        return M;
}
 
template <Dims_t Dims, class Int>
inline MultiIndex<Dims>&operator*=(MultiIndex<Dims>&M, Int Factor)
{
        for(Dims_t i=0; i<Dims; i++)
                M[i] *= Factor;
 
        return M;
}
 
template <Dims_t Dims, class Int>
inline MultiIndex<Dims> operator/(const MultiIndex<Dims>&MIn, Int Factor)
{
MultiIndex<Dims> M = MIn;
 
        for(Dims_t i=0; i<Dims; i++)
                M[i] /= Factor;
 
        return M;
}
 
template <Dims_t Dims, class Int>
inline MultiIndex<Dims>&operator/=(MultiIndex<Dims>&M, Int Factor)
{
        for(Dims_t i=0; i<Dims; i++)
                M[i] /= Factor;
 
        return M;
}
 
 
template <Dims_t Dims>
inline FixedArray<double, Dims> operator*(const FixedArray<double, Dims>&a, const MultiIndex<Dims>&I)
{
FixedArray<double,Dims> m;
        for(Dims_t i=0; i<Dims; i++)
                m[i] = a[i]*I[i];
 
        return m;
}
 
template <Dims_t Dims>
inline FixedArray<double, Dims> operator*(const MultiIndex<Dims>&I, const FixedArray<double, Dims>&a)
{
FixedArray<double,Dims> m;
        for(Dims_t i=0; i<Dims; i++)
                m[i] = I[i]*a[i];
 
        return m;
}
 
template <Dims_t Dims>
inline FixedArray<double, Dims> operator/(const FixedArray<double, Dims>&a, const MultiIndex<Dims>&I)
{
FixedArray<double,Dims> m;
        for(Dims_t i=0; i<Dims; i++)
                m[i] = a[i]/I[i];
 
        return m;
}
 
template <Dims_t Dims>
inline FixedArray<double, Dims>&operator/=(FixedArray<double, Dims>&a, const MultiIndex<Dims>&I)
{
        for(Dims_t i=0; i<Dims; i++)
                a[i] /= I[i];
 
        return a;
}
 
template <Dims_t Dims>
inline FixedArray<double, Dims>&operator*=(FixedArray<double, Dims>&a, const MultiIndex<Dims>&I)
{
        for(Dims_t i=0; i<Dims; i++)
                a[i] *= I[i];
 
        return a;
}
 
 
 
template <Dims_t Dims>
inline FixedArray<int, Dims> MultiIndexToFixedArray(const MultiIndex<Dims>&Value)
{
FixedArray<int, Dims> result;
 
        for(Dims_t i=0; i<Dims; i++)
                result[i] = Value[i];
 
        return result;
}
 
template <Dims_t Dims>
inline MultiIndex<Dims>  FixedArrayToMultiIndex(const FixedArray<int, Dims>&Value)
{
MultiIndex<Dims> result;
 
        for(Dims_t i=0; i<Dims; i++)
                result[i] = Value[i];
 
        return result;
}
 
template <Dims_t Dims>
inline MultiIndex<Dims>  minMultiIndex(const MultiIndex<Dims>&A, const MultiIndex<Dims>&B)
{
MultiIndex<Dims>  result;
 
        for(Dims_t i=0; i<Dims; i++)
        {
                result[i] = A[i]<B[i] ? A[i] : B[i];
        }
        return result;
}
 
 
template <Dims_t Dims>
inline MultiIndex<Dims>  maxMultiIndex(const MultiIndex<Dims>&A, const MultiIndex<Dims>&B)
{
MultiIndex<Dims>  result;
 
        for(Dims_t i=0; i<Dims; i++)
        {
                result[i] = A[i]>B[i] ? A[i] : B[i];
        }
        return result;
}
 
 
 
#if 0
template <Dims_t Dims>
        inline FixedArray<double, Dims> Modulo(double*a, double range=1.0)
        {
        FixedArray<double, Dims> m;
                for(Dims_t i=0; i<Dims; i++)
                        m[i] = fmod(a[i] , range);
                return m;
        }
 
 
template <Dims_t Dims>
        inline FixedArray<double, Dims> Modulo(double*a, double range=1.0)
        {
        FixedArray<double, Dims> m;
                for(Dims_t i=0; i<Dims; i++)
                        m[i] = fmod(a[i] , range);
                return m;
        }
 
template <Dims_t Dims>
        inline FixedArray<double, Dims> Modulo(const FixedArray<double, Dims>&a, const FixedArray<double, Dims>&range)
        {
        MultiIndex<Dims, double> m;
                for(Dims_t i=0; i<Dims; i++)
                        m[i] = fmod(a[i] , range[i]);
                return m;
        }
#endif
 
} // namespace Fiber
 
namespace std
{
 
template <Fiber::Dims_t Dims>
inline  string to_string(const Fiber::MultiIndex<Dims>&M)
{
Fiber::Dims_t   i = 0;
 
string  retval = "[";
 
        for(; i<Dims-1; i++)
                retval += std::to_string(M[i]) + "x";
 
        retval += std::to_string(M[i]) + "]";
 
        return retval;
}
 
} // namespace std
 
 
#ifdef  ENABLE_MEMCORE_ASSERT
#include <memcore/Verbose.hpp>
#endif
 
namespace Fiber
{
 
template <Dims_t Dims>
inline  index_t MultiIndex<Dims>::linear(const MultiIndex<Dims>&Dimens) const
{
#ifndef NDEBUG
        if (!(idx >= 0 && idx < Dimens.maxidx() ))
        {
#ifdef  ENABLE_MEMCORE_ASSERT
                Verbose(0) << "Index " << idx << " beyond allowed maximum " << Dimens.maxidx() << " out of " << Dimens;
                Assert(idx >= 0 && idx < Dimens.maxidx() );
#endif
                assert(idx >= 0 && idx < Dimens.maxidx() );
        }
#endif
        /*
          Pretty sure, this formula can be optimized.
          The size() call on the subdimension recursively multiplies
          its elements, but this could be moved into the summation/
          multiplication recursive call.
        */
index_t L = Dimens.subidx().size();
return idx * L + subidx().linear( Dimens.subidx() );
}
 
 
}
 
 
#endif  // __Fiber_MultiIndex_HPP