diff --git a/src/anonymous_range.cc b/src/anonymous_range.cc
deleted file mode 100644
index 9363cdc..0000000
--- a/src/anonymous_range.cc
+++ /dev/null
@@ -1,13 +0,0 @@
-// -*- C++ -*-
-
-#include "anonymous_range.h"
-
-namespace MultiArrayTools
-{
-
-
-
-
-
-
-}
diff --git a/src/anonymous_range.h b/src/anonymous_range.h
index efe21c2..19bcba0 100644
--- a/src/anonymous_range.h
+++ b/src/anonymous_range.h
@@ -18,6 +18,10 @@ namespace MultiArrayTools
 
 }
 
-#include "anonymous_range.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+// ...
 
 #endif
diff --git a/src/block.cc b/src/block.cc
deleted file mode 100644
index a351545..0000000
--- a/src/block.cc
+++ /dev/null
@@ -1,216 +0,0 @@
-// -*- C++ -*-
-
-#include "block.h"
-
-namespace MultiArrayHelper
-{
-
-    /*********************
-     *   BlockBinaryOp   *
-     *********************/
-
-    template <typename T, class OpFunc, class BlockClass1, class BlockClass2>
-    BlockResult<T>
-    BlockBinaryOp<T,OpFunc,BlockClass1,BlockClass2>::operator()(const BlockClass1& arg1,
-								const BlockClass2& arg2)
-    {
-	static OpFunc f;
-	BlockResult<T> res(arg1.size());
-	assert(arg1.size() == arg2.size());
-	for(size_t i = 0; i != arg1.size(); ++i){
-	    res[i] = f(arg1[i], arg2[i]);
-	}
-	return res;
-    }
-
-    template <typename T, class OpFunc, class BlockClass>
-    BlockBinaryOpSelf<T,OpFunc,BlockClass>::BlockBinaryOpSelf(BlockResult<T>& res) : mRes(res) {}
-
-    template <typename T, class OpFunc, class BlockClass>
-    void BlockBinaryOpSelf<T,OpFunc,BlockClass>::operator()(const BlockClass& arg)
-    {
-	static OpFunc f;
-	if(mRes.size() == 0) { mRes.assign(arg.size(), static_cast<T>(0)); }
-	assert(mRes.size() == arg.size());
-	for(size_t i = 0; i != arg.size(); ++i){
-	    mRes[i] = f(mRes[i], arg[i]);
-	}
-    }
-
-    
-    /*****************
-     *   BlockBase   *
-     *****************/
-    
-    template <typename T>
-    BlockBase<T>::BlockBase(size_t size) : mSize(size) {}
-
-    template <typename T>
-    size_t BlockBase<T>::size() const
-    {
-	return mSize;
-    }
-
-    /************************
-     *   MutableBlockBase   *
-     ************************/
-    
-    template <typename T>
-    MutableBlockBase<T>::MutableBlockBase(size_t size) : BlockBase<T>(size) {}
-
-    
-    /*************
-     *   Block   *
-     *************/
-
-    template <typename T>
-    Block<T>::Block(const std::vector<T>& data,
-		    size_t begPos, size_t size, size_t stepSize) :
-	BlockBase<T>(size),
-	mData(&data),
-	mBegPtr(data.data() + begPos),
-	mStepSize(stepSize) {}
-
-    template <typename T>
-    BlockType Block<T>::type() const
-    {
-	return mStepSize == 0 ? BlockType::VALUE :
-	    ( mStepSize == 1 ? BlockType::BLOCK : BlockType::SPLIT );
-    }
-    
-    template <typename T>
-    const T& Block<T>::operator[](size_t i) const
-    {
-	
-	return *(mBegPtr + i * mStepSize);
-    }
-
-    template <typename T>
-    Block<T>& Block<T>::set(size_t npos)
-    {
-	mBegPtr = &(*mData)[npos];
-	return *this;
-    }
-
-    template <typename T>
-    size_t Block<T>::stepSize() const
-    {
-	return 1;
-    }
-    
-    /**************
-     *   MBlock   *
-     **************/
-
-    template <typename T>
-    MBlock<T>::MBlock(std::vector<T>& data,
-		      size_t begPos, size_t size, size_t stepSize) :
-	MutableBlockBase<T>(size),
-	mData(&data),
-	mBegPtr(data.data() + begPos),
-	mStepSize(stepSize) {}
-
-    template <typename T>
-    template <class BlockClass>
-    MBlock<T>& MBlock<T>::operator=(const BlockClass& in)
-    {
-	for(size_t i = 0; i != BlockBase<T>::mSize; ++i){
-	    (*this)[i] = in[i];
-	}
-	return *this;
-    }
-   
-    template <typename T>
-    BlockType MBlock<T>::type() const
-    {
-	return mStepSize == 0 ? BlockType::VALUE :
-	    ( mStepSize == 1 ? BlockType::BLOCK : BlockType::SPLIT );
-    }
-    
-    template <typename T>
-    const T& MBlock<T>::operator[](size_t i) const
-    {
-	
-	return *(mBegPtr + i * mStepSize);
-    }
-
-    template <typename T>
-    T& MBlock<T>::operator[](size_t i)
-    {
-	
-	return *(mBegPtr + i * mStepSize);
-    }
-    
-    template <typename T>
-    MBlock<T>& MBlock<T>::set(size_t npos)
-    {
-	mBegPtr = &(*mData)[npos];
-	return *this;
-    }   
-
-    template <typename T>
-    size_t MBlock<T>::stepSize() const
-    {
-	return 1;
-    }
-    
-    /*******************
-     *   BlockResult   *
-     *******************/
-
-    template <typename T>
-    BlockResult<T>::BlockResult(size_t size) :
-	MutableBlockBase<T>(size),
-	mRes(size) {}
-
-    template <typename T>
-    template <class BlockClass>
-    BlockResult<T>& BlockResult<T>::operator=(const BlockClass& in)
-    {
-	//CHECK;
-	for(size_t i = 0; i != BlockBase<T>::mSize; ++i){
-	    (*this)[i] = in[i];
-	}
-	return *this;
-    }
-    
-    template <typename T>
-    BlockType BlockResult<T>::type() const
-    {
-	return BlockType::RESULT;
-    }
-    
-    template <typename T>
-    const T& BlockResult<T>::operator[](size_t i) const
-    {
-	return mRes[i];
-    }
-
-    template <typename T>
-    T& BlockResult<T>::operator[](size_t i)
-    {
-	
-	return mRes[i];
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::set(size_t npos)
-    {
-	return *this;
-    }
-
-    template <typename T>
-    size_t BlockResult<T>::stepSize() const
-    {
-	return 1;
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::assign(size_t size, const T& val)
-    {
-	BB::mSize = size;
-	mRes.assign(BB::mSize, val);
-	return *this;
-    }
-    
-} // end namespace MultiArrayHelper
diff --git a/src/block.h b/src/block.h
index e007350..cf19632 100644
--- a/src/block.h
+++ b/src/block.h
@@ -153,6 +153,221 @@ namespace MultiArrayHelper
     
 } // end namespace MultiArrayHelper
 
-#include "block.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayHelper
+{
+
+    /*********************
+     *   BlockBinaryOp   *
+     *********************/
+
+    template <typename T, class OpFunc, class BlockClass1, class BlockClass2>
+    BlockResult<T>
+    BlockBinaryOp<T,OpFunc,BlockClass1,BlockClass2>::operator()(const BlockClass1& arg1,
+								const BlockClass2& arg2)
+    {
+	static OpFunc f;
+	BlockResult<T> res(arg1.size());
+	assert(arg1.size() == arg2.size());
+	for(size_t i = 0; i != arg1.size(); ++i){
+	    res[i] = f(arg1[i], arg2[i]);
+	}
+	return res;
+    }
+
+    template <typename T, class OpFunc, class BlockClass>
+    BlockBinaryOpSelf<T,OpFunc,BlockClass>::BlockBinaryOpSelf(BlockResult<T>& res) : mRes(res) {}
+
+    template <typename T, class OpFunc, class BlockClass>
+    void BlockBinaryOpSelf<T,OpFunc,BlockClass>::operator()(const BlockClass& arg)
+    {
+	static OpFunc f;
+	if(mRes.size() == 0) { mRes.assign(arg.size(), static_cast<T>(0)); }
+	assert(mRes.size() == arg.size());
+	for(size_t i = 0; i != arg.size(); ++i){
+	    mRes[i] = f(mRes[i], arg[i]);
+	}
+    }
+
+    
+    /*****************
+     *   BlockBase   *
+     *****************/
+    
+    template <typename T>
+    BlockBase<T>::BlockBase(size_t size) : mSize(size) {}
+
+    template <typename T>
+    size_t BlockBase<T>::size() const
+    {
+	return mSize;
+    }
+
+    /************************
+     *   MutableBlockBase   *
+     ************************/
+    
+    template <typename T>
+    MutableBlockBase<T>::MutableBlockBase(size_t size) : BlockBase<T>(size) {}
+
+    
+    /*************
+     *   Block   *
+     *************/
+
+    template <typename T>
+    Block<T>::Block(const std::vector<T>& data,
+		    size_t begPos, size_t size, size_t stepSize) :
+	BlockBase<T>(size),
+	mData(&data),
+	mBegPtr(data.data() + begPos),
+	mStepSize(stepSize) {}
+
+    template <typename T>
+    BlockType Block<T>::type() const
+    {
+	return mStepSize == 0 ? BlockType::VALUE :
+	    ( mStepSize == 1 ? BlockType::BLOCK : BlockType::SPLIT );
+    }
+    
+    template <typename T>
+    const T& Block<T>::operator[](size_t i) const
+    {
+	
+	return *(mBegPtr + i * mStepSize);
+    }
+
+    template <typename T>
+    Block<T>& Block<T>::set(size_t npos)
+    {
+	mBegPtr = &(*mData)[npos];
+	return *this;
+    }
+
+    template <typename T>
+    size_t Block<T>::stepSize() const
+    {
+	return 1;
+    }
+    
+    /**************
+     *   MBlock   *
+     **************/
+
+    template <typename T>
+    MBlock<T>::MBlock(std::vector<T>& data,
+		      size_t begPos, size_t size, size_t stepSize) :
+	MutableBlockBase<T>(size),
+	mData(&data),
+	mBegPtr(data.data() + begPos),
+	mStepSize(stepSize) {}
+
+    template <typename T>
+    template <class BlockClass>
+    MBlock<T>& MBlock<T>::operator=(const BlockClass& in)
+    {
+	for(size_t i = 0; i != BlockBase<T>::mSize; ++i){
+	    (*this)[i] = in[i];
+	}
+	return *this;
+    }
+   
+    template <typename T>
+    BlockType MBlock<T>::type() const
+    {
+	return mStepSize == 0 ? BlockType::VALUE :
+	    ( mStepSize == 1 ? BlockType::BLOCK : BlockType::SPLIT );
+    }
+    
+    template <typename T>
+    const T& MBlock<T>::operator[](size_t i) const
+    {
+	
+	return *(mBegPtr + i * mStepSize);
+    }
+
+    template <typename T>
+    T& MBlock<T>::operator[](size_t i)
+    {
+	
+	return *(mBegPtr + i * mStepSize);
+    }
+    
+    template <typename T>
+    MBlock<T>& MBlock<T>::set(size_t npos)
+    {
+	mBegPtr = &(*mData)[npos];
+	return *this;
+    }   
+
+    template <typename T>
+    size_t MBlock<T>::stepSize() const
+    {
+	return 1;
+    }
+    
+    /*******************
+     *   BlockResult   *
+     *******************/
+
+    template <typename T>
+    BlockResult<T>::BlockResult(size_t size) :
+	MutableBlockBase<T>(size),
+	mRes(size) {}
+
+    template <typename T>
+    template <class BlockClass>
+    BlockResult<T>& BlockResult<T>::operator=(const BlockClass& in)
+    {
+	//CHECK;
+	for(size_t i = 0; i != BlockBase<T>::mSize; ++i){
+	    (*this)[i] = in[i];
+	}
+	return *this;
+    }
+    
+    template <typename T>
+    BlockType BlockResult<T>::type() const
+    {
+	return BlockType::RESULT;
+    }
+    
+    template <typename T>
+    const T& BlockResult<T>::operator[](size_t i) const
+    {
+	return mRes[i];
+    }
+
+    template <typename T>
+    T& BlockResult<T>::operator[](size_t i)
+    {
+	
+	return mRes[i];
+    }
+
+    template <typename T>
+    BlockResult<T>& BlockResult<T>::set(size_t npos)
+    {
+	return *this;
+    }
+
+    template <typename T>
+    size_t BlockResult<T>::stepSize() const
+    {
+	return 1;
+    }
+
+    template <typename T>
+    BlockResult<T>& BlockResult<T>::assign(size_t size, const T& val)
+    {
+	BB::mSize = size;
+	mRes.assign(BB::mSize, val);
+	return *this;
+    }
+    
+} // end namespace MultiArrayHelper
 
 #endif
diff --git a/src/container_range.cc b/src/container_range.cc
deleted file mode 100644
index c0fbdd6..0000000
--- a/src/container_range.cc
+++ /dev/null
@@ -1,296 +0,0 @@
-// -*- C++ -*-
-
-#include "container_range.h"
-#include "pack_num.h"
-
-namespace MultiArrayTools
-{
-    namespace
-    {
-	using namespace MultiArrayHelper;
-    }
-
-    
-    /**********************
-     *   ContainerIndex   *
-     **********************/
-    /*
-    template <class... Indices>
-    ContainerIndex<Indices...>::ContainerIndex(const ContainerIndex& in) :
-	IndexInterface<std::tuple<typename Indices::MetaType...> >(in)
-    {
-	PackNum<sizeof...(Indices)-1>::copy(mIPack, in);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-    }
-
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator=(const ContainerIndex& in)
-    {
-	IndexI::operator=(in);
-	PackNum<sizeof...(Indices)-1>::copy(mIPack, in);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	return *this;
-    }
-    */
-    template <class... Indices>
-    template <class MRange>
-    ContainerIndex<Indices...>::ContainerIndex(const std::shared_ptr<MRange>& range) :
-	IndexInterface<std::tuple<typename Indices::MetaType...> >(range, 0)
-    {
-	PackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
-	PackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack);
-    }
-
-    template <class... Indices>
-    IndexType ContainerIndex<Indices...>::type() const
-    {
-	return IndexType::CONT;
-    }
-    
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator++()
-    {
-	if(mExternControl){
-	    IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	}
-	PackNum<sizeof...(Indices)-1>::pp( mIPack );
-	++IB::mPos;
-	return *this;
-    }
-
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator--()
-    {
-	if(mExternControl){
-	    IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	}
-	PackNum<sizeof...(Indices)-1>::mm( mIPack );
-	--IB::mPos;
-	return *this;
-
-    }
-
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator=(size_t pos)
-    {
-	IB::mPos = pos;
-	PackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
-	return *this;
-    }
-
-    template <class... Indices>
-    int ContainerIndex<Indices...>::pp(std::shared_ptr<IndexBase>& idxPtr)
-    {
-	int tmp = PackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtr);
-	IB::mPos += tmp;
-	return tmp;
-    }
-
-    template <class... Indices>
-    int ContainerIndex<Indices...>::mm(std::shared_ptr<IndexBase>& idxPtr)
-    {
-	int tmp = PackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtr);
-	IB::mPos -= tmp;
-	return tmp;
-    }
-    
-    template <class... Indices>
-    typename ContainerIndex<Indices...>::MetaType ContainerIndex<Indices...>::meta() const
-    {
-	MetaType metaTuple;
-	PackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
-	return metaTuple;
-    }
-
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::at(const MetaType& metaPos)
-    {
-	PackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	return *this;
-    }
-
-    template <class... Indices>
-    size_t ContainerIndex<Indices...>::dim() const
-    {
-	return sizeof...(Indices);
-    }
-
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::sync()
-    {
-	if(mExternControl){
-	    IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	    //VCHECK(id());
-	    //VCHECK(sizeof...(Indices));
-	    //assert(IB::mPos < IB::max());
-	}
-	return *this;
-    }
-    
-    template <class... Indices>
-    template <size_t N>
-    auto ContainerIndex<Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
-    {
-	return *std::get<N>( mIPack );
-    }
-
-    template <class... Indices>
-    template <size_t N>
-    auto ContainerIndex<Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
-    {
-	return std::get<N>( mIPack );
-    }
-
-    template <class... Indices>
-    std::shared_ptr<IndexBase> ContainerIndex<Indices...>::getPtr(size_t n) const
-    {
-	if(n >= sizeof...(Indices)){
-	    assert(0);
-	    // throw !!
-	}
-	ContainerIndex<Indices...> const* t = this;
-	return PackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
-    }
-
-    
-    template <class... Indices>
-    size_t ContainerIndex<Indices...>::getStepSize(size_t n) const
-    {
-	if(n >= sizeof...(Indices)){
-	    assert(0);
-	    // throw !!
-	}
-	return mBlockSizes[n+1];
-    }
-	
-    template <class... Indices>
-    bool ContainerIndex<Indices...>::first() const
-    {
-	return IB::pos() == 0;
-    }
-
-    template <class... Indices>
-    bool ContainerIndex<Indices...>::last() const
-    {
-	return IB::pos() == IB::mRangePtr->size() - 1;
-    }
-    
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator()(const std::shared_ptr<Indices>&... inds)
-    {
-	PackNum<sizeof...(Indices)-1>::swapIndices(mIPack, inds...);
-	mExternControl = true;
-	return sync();
-    }    
-
-    template <class... Indices>
-    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator()()
-    {
-	return sync();
-    }
-
-    template <class... Indices>
-    std::shared_ptr<typename ContainerIndex<Indices...>::RangeType> ContainerIndex<Indices...>::range() const
-    {
-	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
-    }
-    
-    /*****************************
-     *   ContainerRangeFactory   *
-     *****************************/
-
-    template <class... Ranges>
-    ContainerRangeFactory<Ranges...>::ContainerRangeFactory(const std::shared_ptr<Ranges>&... rs)
-    {
-	mProd = std::shared_ptr<ContainerRange<Ranges...> >( new ContainerRange<Ranges...>( rs... ) );
-    }
-    
-    template <class... Ranges>
-    ContainerRangeFactory<Ranges...>::
-    ContainerRangeFactory(const typename ContainerRange<Ranges...>::SpaceType& space)
-    {
-	mProd = std::shared_ptr<ContainerRange<Ranges...> >( new ContainerRange<Ranges...>( space ) );
-    }
-
-    template <class... Ranges>
-    std::shared_ptr<RangeBase> ContainerRangeFactory<Ranges...>::create()
-    {
-	setSelf();
-	return mProd;
-    }
-    
-    /**********************
-     *   ContainerRange   *
-     **********************/
-
-    template <class... Ranges>
-    ContainerRange<Ranges...>::ContainerRange(const std::shared_ptr<Ranges>&... rs) :
-	mSpace( std::make_tuple( rs... ) ) {}
-
-    template <class... Ranges>
-    ContainerRange<Ranges...>::ContainerRange(const SpaceType& space) : mSpace( space ) {}
-
-    template <class... Ranges>
-    size_t ContainerRange<Ranges...>::dim() const
-    {
-	return sizeof...(Ranges);
-    }
-
-    template <class... Ranges>
-    size_t ContainerRange<Ranges...>::size() const
-    {
-	return PackNum<sizeof...(Ranges)-1>::getSize(mSpace);
-    }
-
-    template <class... Ranges>
-    template <size_t N>
-    auto ContainerRange<Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
-    {
-	return *std::get<N>( mSpace );
-    }
-
-    template <class... Ranges>
-    template <size_t N>
-    auto ContainerRange<Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
-    {
-	return std::get<N>( mSpace );
-    }
-
-    template <class... Ranges>
-    const typename ContainerRange<Ranges...>::SpaceType& ContainerRange<Ranges...>::space() const
-    {
-	return mSpace;
-    }
-        
-    template <class... Ranges>
-    typename ContainerRange<Ranges...>::IndexType ContainerRange<Ranges...>::begin() const
-    {
-	ContainerIndex<typename Ranges::IndexType...>
-	    i( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
-	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-	i = 0;
-	return i;
-    }
-
-    template <class... Ranges>
-    typename ContainerRange<Ranges...>::IndexType ContainerRange<Ranges...>::end() const
-    {
-	ContainerIndex<typename Ranges::IndexType...>
-	    i( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
-	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-	i = size();
-	return i;
-    }
-    
-    template <class... Ranges>
-    std::shared_ptr<IndexBase> ContainerRange<Ranges...>::index() const
-    {
-	return std::make_shared<ContainerIndex<typename Ranges::IndexType...> >
-	    ( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
-	      ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-    }
-    
-} // end namespace MultiArrayTools
diff --git a/src/container_range.h b/src/container_range.h
index 4482567..902ce65 100644
--- a/src/container_range.h
+++ b/src/container_range.h
@@ -11,6 +11,8 @@
 #include "range_base.h"
 #include "index_base.h"
 
+#include "pack_num.h"
+
 namespace MultiArrayTools
 {
     
@@ -137,6 +139,283 @@ namespace MultiArrayTools
    
 } // end namespace MultiArrayTools
 
-#include "container_range.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+    
+    /**********************
+     *   ContainerIndex   *
+     **********************/
+
+    template <class... Indices>
+    template <class MRange>
+    ContainerIndex<Indices...>::ContainerIndex(const std::shared_ptr<MRange>& range) :
+	IndexInterface<std::tuple<typename Indices::MetaType...> >(range, 0)
+    {
+	PackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
+	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
+	PackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack);
+    }
+
+    template <class... Indices>
+    IndexType ContainerIndex<Indices...>::type() const
+    {
+	return IndexType::CONT;
+    }
+    
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator++()
+    {
+	if(mExternControl){
+	    IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	}
+	PackNum<sizeof...(Indices)-1>::pp( mIPack );
+	++IB::mPos;
+	return *this;
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator--()
+    {
+	if(mExternControl){
+	    IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	}
+	PackNum<sizeof...(Indices)-1>::mm( mIPack );
+	--IB::mPos;
+	return *this;
+
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	PackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
+	return *this;
+    }
+
+    template <class... Indices>
+    int ContainerIndex<Indices...>::pp(std::shared_ptr<IndexBase>& idxPtr)
+    {
+	int tmp = PackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtr);
+	IB::mPos += tmp;
+	return tmp;
+    }
+
+    template <class... Indices>
+    int ContainerIndex<Indices...>::mm(std::shared_ptr<IndexBase>& idxPtr)
+    {
+	int tmp = PackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtr);
+	IB::mPos -= tmp;
+	return tmp;
+    }
+    
+    template <class... Indices>
+    typename ContainerIndex<Indices...>::MetaType ContainerIndex<Indices...>::meta() const
+    {
+	MetaType metaTuple;
+	PackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
+	return metaTuple;
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::at(const MetaType& metaPos)
+    {
+	PackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
+	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+
+    template <class... Indices>
+    size_t ContainerIndex<Indices...>::dim() const
+    {
+	return sizeof...(Indices);
+    }
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::sync()
+    {
+	if(mExternControl){
+	    IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	    //VCHECK(id());
+	    //VCHECK(sizeof...(Indices));
+	    //assert(IB::mPos < IB::max());
+	}
+	return *this;
+    }
+    
+    template <class... Indices>
+    template <size_t N>
+    auto ContainerIndex<Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
+    {
+	return *std::get<N>( mIPack );
+    }
+
+    template <class... Indices>
+    template <size_t N>
+    auto ContainerIndex<Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>( mIPack );
+    }
+
+    template <class... Indices>
+    std::shared_ptr<IndexBase> ContainerIndex<Indices...>::getPtr(size_t n) const
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	ContainerIndex<Indices...> const* t = this;
+	return PackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
+    }
+
+    
+    template <class... Indices>
+    size_t ContainerIndex<Indices...>::getStepSize(size_t n) const
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	return mBlockSizes[n+1];
+    }
+	
+    template <class... Indices>
+    bool ContainerIndex<Indices...>::first() const
+    {
+	return IB::pos() == 0;
+    }
+
+    template <class... Indices>
+    bool ContainerIndex<Indices...>::last() const
+    {
+	return IB::pos() == IB::mRangePtr->size() - 1;
+    }
+    
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator()(const std::shared_ptr<Indices>&... inds)
+    {
+	PackNum<sizeof...(Indices)-1>::swapIndices(mIPack, inds...);
+	mExternControl = true;
+	return sync();
+    }    
+
+    template <class... Indices>
+    ContainerIndex<Indices...>& ContainerIndex<Indices...>::operator()()
+    {
+	return sync();
+    }
+
+    template <class... Indices>
+    std::shared_ptr<typename ContainerIndex<Indices...>::RangeType> ContainerIndex<Indices...>::range() const
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
+    }
+    
+    /*****************************
+     *   ContainerRangeFactory   *
+     *****************************/
+
+    template <class... Ranges>
+    ContainerRangeFactory<Ranges...>::ContainerRangeFactory(const std::shared_ptr<Ranges>&... rs)
+    {
+	mProd = std::shared_ptr<ContainerRange<Ranges...> >( new ContainerRange<Ranges...>( rs... ) );
+    }
+    
+    template <class... Ranges>
+    ContainerRangeFactory<Ranges...>::
+    ContainerRangeFactory(const typename ContainerRange<Ranges...>::SpaceType& space)
+    {
+	mProd = std::shared_ptr<ContainerRange<Ranges...> >( new ContainerRange<Ranges...>( space ) );
+    }
+
+    template <class... Ranges>
+    std::shared_ptr<RangeBase> ContainerRangeFactory<Ranges...>::create()
+    {
+	setSelf();
+	return mProd;
+    }
+    
+    /**********************
+     *   ContainerRange   *
+     **********************/
+
+    template <class... Ranges>
+    ContainerRange<Ranges...>::ContainerRange(const std::shared_ptr<Ranges>&... rs) :
+	mSpace( std::make_tuple( rs... ) ) {}
+
+    template <class... Ranges>
+    ContainerRange<Ranges...>::ContainerRange(const SpaceType& space) : mSpace( space ) {}
+
+    template <class... Ranges>
+    size_t ContainerRange<Ranges...>::dim() const
+    {
+	return sizeof...(Ranges);
+    }
+
+    template <class... Ranges>
+    size_t ContainerRange<Ranges...>::size() const
+    {
+	return PackNum<sizeof...(Ranges)-1>::getSize(mSpace);
+    }
+
+    template <class... Ranges>
+    template <size_t N>
+    auto ContainerRange<Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
+    {
+	return *std::get<N>( mSpace );
+    }
+
+    template <class... Ranges>
+    template <size_t N>
+    auto ContainerRange<Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
+    {
+	return std::get<N>( mSpace );
+    }
+
+    template <class... Ranges>
+    const typename ContainerRange<Ranges...>::SpaceType& ContainerRange<Ranges...>::space() const
+    {
+	return mSpace;
+    }
+        
+    template <class... Ranges>
+    typename ContainerRange<Ranges...>::IndexType ContainerRange<Ranges...>::begin() const
+    {
+	ContainerIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+
+    template <class... Ranges>
+    typename ContainerRange<Ranges...>::IndexType ContainerRange<Ranges...>::end() const
+    {
+	ContainerIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = size();
+	return i;
+    }
+    
+    template <class... Ranges>
+    std::shared_ptr<IndexBase> ContainerRange<Ranges...>::index() const
+    {
+	return std::make_shared<ContainerIndex<typename Ranges::IndexType...> >
+	    ( std::dynamic_pointer_cast<ContainerRange<Ranges...> >
+	      ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+    }
+    
+} // end namespace MultiArrayTools
+
 
 #endif
diff --git a/src/helper_tools.h b/src/helper_tools.h
index db85569..c0f41ee 100644
--- a/src/helper_tools.h
+++ b/src/helper_tools.h
@@ -13,9 +13,9 @@ namespace MultiArrayTools
     
 }
 
-// =====================
-// ===    C O D E    ===
-// =====================
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
 
 namespace MultiArrayTools
 {
diff --git a/src/index_base.cc b/src/index_base.cc
deleted file mode 100644
index 6d31d79..0000000
--- a/src/index_base.cc
+++ /dev/null
@@ -1,67 +0,0 @@
-// -*- C++ -*-
-#include "index_base.h"
-
-namespace MultiArrayTools
-{
-    /*****************
-     *   IndexBase   *
-     *****************/
-
-    IndexBase::IndexBase(const std::shared_ptr<RangeBase>& range,
-			 size_t pos) : mRangePtr(range),
-				       mPos(pos)
-    {
-	mId = indexId();
-    }
-    
-    bool IndexBase::operator==(const IndexBase& in) const
-    {
-	return in.mPos == mPos and in.mRangePtr.get() == mRangePtr.get();
-    }
-    
-    bool IndexBase::operator!=(const IndexBase& in) const
-    {
-	return in.mPos != mPos or in.mRangePtr.get() != mRangePtr.get();
-    }
-    
-    size_t IndexBase::pos() const
-    {
-	return mPos;
-    }
-
-    size_t IndexBase::max() const
-    {
-	return mRangePtr->size();
-    }
-    /*
-    bool IndexBase::locked() const
-    {
-	return mLocked;
-    }
-    
-    IndexBase& IndexBase::lock(std::shared_ptr<IndexBase>& idx)
-    {
-	mLocked = (idx.get() == this);
-	return *this;
-    }
-    */
-
-    std::shared_ptr<RangeBase> IndexBase::rangePtr() const
-    {
-	return mRangePtr;
-    }
-    
-    IndexBase::operator size_t() const
-    {
-	return pos();
-    }
-    
-    /**********************
-     *   IndexInterface   *	     
-     **********************/
-
-    template <typename MetaType>
-    IndexInterface<MetaType>::IndexInterface(const std::shared_ptr<RangeBase>& rangePtr, size_t pos) :
-	IndexBase(rangePtr, pos) {}
-    
-}
diff --git a/src/index_base.h b/src/index_base.h
index a45f904..0cce418 100644
--- a/src/index_base.h
+++ b/src/index_base.h
@@ -90,6 +90,73 @@ namespace MultiArrayTools
     
 }
 
-#include "index_base.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    /*****************
+     *   IndexBase   *
+     *****************/
+
+    IndexBase::IndexBase(const std::shared_ptr<RangeBase>& range,
+			 size_t pos) : mRangePtr(range),
+				       mPos(pos)
+    {
+	mId = indexId();
+    }
+    
+    bool IndexBase::operator==(const IndexBase& in) const
+    {
+	return in.mPos == mPos and in.mRangePtr.get() == mRangePtr.get();
+    }
+    
+    bool IndexBase::operator!=(const IndexBase& in) const
+    {
+	return in.mPos != mPos or in.mRangePtr.get() != mRangePtr.get();
+    }
+    
+    size_t IndexBase::pos() const
+    {
+	return mPos;
+    }
+
+    size_t IndexBase::max() const
+    {
+	return mRangePtr->size();
+    }
+    /*
+    bool IndexBase::locked() const
+    {
+	return mLocked;
+    }
+    
+    IndexBase& IndexBase::lock(std::shared_ptr<IndexBase>& idx)
+    {
+	mLocked = (idx.get() == this);
+	return *this;
+    }
+    */
+
+    std::shared_ptr<RangeBase> IndexBase::rangePtr() const
+    {
+	return mRangePtr;
+    }
+    
+    IndexBase::operator size_t() const
+    {
+	return pos();
+    }
+    
+    /**********************
+     *   IndexInterface   *	     
+     **********************/
+
+    template <typename MetaType>
+    IndexInterface<MetaType>::IndexInterface(const std::shared_ptr<RangeBase>& rangePtr, size_t pos) :
+	IndexBase(rangePtr, pos) {}
+    
+}
 
 #endif
diff --git a/src/manipulator.cc b/src/manipulator.cc
deleted file mode 100644
index 43f8a60..0000000
--- a/src/manipulator.cc
+++ /dev/null
@@ -1,56 +0,0 @@
-// -*- C++ -*-
-
-#include "manipulator.h"
-
-namespace MultiArrayTools
-{
-
-    template <typename T>
-    void ManipulatorBase<T>::setup(std::vector<T>& targetRef, size_t begin, size_t end)
-    {
-	mBegin = begin;
-	mEnd = end;
-	mTargetPtr = &targetRef;
-    }
-
-    template <typename T>
-    void ManipulatorBase<T>::reset()
-    {
-	mLastPos = 0;
-    }
-
-    template <typename T>
-    BinReader<T>::BinReader(const std::string& fname, size_t readBegin):
-	mFs(fname, std::fstream::in | std::fstream::binary),
-	mReadBegin(readBegin) {}
-
-    template <typename T>
-    BinReader<T>::~BinReader()
-    {
-	mFs.close();
-    }
-    
-    template <typename T>
-    void BinReader<T>::execute()
-    {
-	const size_t readSize = MB::mEnd - MB::mBegin;
-	const size_t blockSize = sizeof(T) * readSize;
-	char* buffer = new char[blockSize];
-
-	mFs.seekg(mReadBegin + MB::mLastPos, mFs.beg);
-	mFs.read(buffer, blockSize);
-	MB::mLastPos += blockSize;
-	
-	T* content = reinterpret_cast<T*>( buffer );
-
-	if(MB::mTargetPtr->size() < MB::mEnd){
-	    assert(0);
-	    // throw
-	}
-	
-	std::copy(&content[0], &content[readSize], MB::mTargetPtr->begin() + MB::mBegin);
-
-	delete[] buffer;	
-    }
-
-}
diff --git a/src/manipulator.h b/src/manipulator.h
deleted file mode 100644
index f884d0d..0000000
--- a/src/manipulator.h
+++ /dev/null
@@ -1,56 +0,0 @@
-// -*- C++ -*-
-
-#ifndef __manipulator_h__
-#define __manipulator_h__
-
-#include <cstdlib>
-#include <vector>
-#include <string>
-#include <fstream>
-
-#include "base_def.h"
-
-namespace MultiArrayTools
-{
-
-    // some kind of input stream or sth similar...
-    template <typename T>
-    class ManipulatorBase
-    {
-    public:
-	virtual ~ManipulatorBase() {}
-	
-	virtual void setup(std::vector<T>& targetRef, size_t begin, size_t end);
-	virtual void execute(size_t pos) = 0;
-	virtual void reset();
-	
-    protected:
-
-	size_t mLastPos = 0;
-	size_t mBegin = 0;
-	size_t mEnd = 0;
-	std::vector<T>* mTargetPtr = nullptr;
-    };
-
-    template <typename T>
-    class BinReader : public ManipulatorBase<T>
-    {
-    public:
-	typedef ManipulatorBase<T> MB;
-
-	BinReader(const std::string& fname, size_t readBegin = 0);
-	virtual ~BinReader();
-	
-	
-	virtual void execute() override;
-
-    private:
-	size_t mReadBegin;
-	std::fstream mFs;
-    };
-
-}
-
-#include "manipulator.cc"
-
-#endif
diff --git a/src/multi_array.cc b/src/multi_array.cc
deleted file mode 100644
index 58df5ec..0000000
--- a/src/multi_array.cc
+++ /dev/null
@@ -1,595 +0,0 @@
-// -*- C++ -*-
-
-#include "multi_array.h"
-
-namespace MultiArrayTools
-{
-    /**************************************
-     *  MultiArrayBase::const_iterator    *	     
-     **************************************/
-
-    template <typename T, class... SRanges>
-    MultiArrayBase<T,SRanges...>::const_iterator::const_iterator(const MultiArrayBase<T,SRanges...>& ma):
-	mMAPtr(&ma), mPos(0) { }
-
-    template <typename T, class... SRanges>
-    MultiArrayBase<T,SRanges...>::const_iterator::const_iterator(const MultiArrayBase<T,SRanges...>& ma,
-							     const typename CRange::IndexType& index):
-	mMAPtr(&ma), mPos(index.pos()) { }
-    
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::const_iterator::operator==(const const_iterator& it) const
-    {
-	return mMAPtr == it.mMAPtr and mPos == it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::const_iterator::operator!=(const const_iterator& it) const
-    {
-	return mMAPtr != it.mMAPtr or mPos != it.mPos;
-    }
-    
-    template <typename T, class... SRanges>
-    const T& MultiArrayBase<T,SRanges...>::const_iterator::operator*() const
-    {
-	return mMAPtr->data()[mPos];
-    }
-
-    template <typename T, class... SRanges>
-    T const* MultiArrayBase<T,SRanges...>::const_iterator::operator->() const
-    {
-	return &mMAPtr->data()[mPos];
-    }
-    
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator++()
-    {
-	++mPos;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator++(int)
-    {
-	const_iterator tmp(*this);
-	++mPos;
-	return tmp;
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator--()
-    {
-	--mPos;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator--(int)
-    {
-	const_iterator tmp(*this);
-	--mPos;
-	return tmp;
-    }
-    
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator+=(int diff)
-    {
-	mPos += diff;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator-=(int diff)
-    {
-	mPos -= diff;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator+(int num) const
-    {
-	const_iterator tmp(*this);
-	tmp += num;
-	return tmp;
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator-(int num) const
-    {
-	const_iterator tmp(*this);
-	tmp -= num;
-    }
-    
-    template <typename T, class... SRanges>
-    int MultiArrayBase<T,SRanges...>::const_iterator::operator-(const const_iterator& it) const
-    {
-	return mPos - it.mPos;
-    }
-    
-    template <typename T, class... SRanges>
-    const T& MultiArrayBase<T,SRanges...>::const_iterator::operator[](int num) const
-    {
-	return *(operator+(num));
-    }
-    
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::const_iterator::operator<(const const_iterator& it) const
-    {
-	return mPos < it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::const_iterator::operator>(const const_iterator& it) const
-    {
-	return mPos > it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::const_iterator::operator<=(const const_iterator& it) const
-    {
-	return mPos <= it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::const_iterator::operator>=(const const_iterator& it) const
-    {
-	return mPos >= it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::IndexType
-    MultiArrayBase<T,SRanges...>::const_iterator::index() const
-    {
-	auto i = mMAPtr->beginIndex();
-	i = mPos;
-	return i;
-    }
-    
-    /**********************
-     *  MultiArrayBase    *	     
-     **********************/
-
-    template <typename T, class... SRanges>
-    MultiArrayBase<T,SRanges...>::MultiArrayBase(const std::shared_ptr<SRanges>&... ranges)
-    {
-	ContainerRangeFactory<SRanges...> crf(ranges...);
-	mRange = std::dynamic_pointer_cast<ContainerRange<SRanges...> >( crf.create() );
-    }
-    
-    template <typename T, class... SRanges>
-    size_t MultiArrayBase<T,SRanges...>::size() const
-    {
-	return mRange->size();
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::begin() const
-    {
-	return const_iterator(*this, beginIndex());
-    }
-    
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::end() const
-    {
-	return const_iterator(*this, endIndex());
-    }
-    
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::IndexType
-    MultiArrayBase<T,SRanges...>::beginIndex() const
-    {
-	return mRange->begin();
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::IndexType
-    MultiArrayBase<T,SRanges...>::endIndex() const
-    {
-	return mRange->end();
-    }
-
-    template <typename T, class... SRanges>
-    const std::shared_ptr<typename MultiArrayBase<T,SRanges...>::CRange>&
-    MultiArrayBase<T,SRanges...>::range() const
-    {
-	return mRange;
-    }
-
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::isConst() const
-    {
-	return true;
-    }
-    
-    template <typename T, class... SRanges>
-    ConstOperationRoot<T,SRanges...>
-    MultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
-    {
-	return ConstOperationRoot<T,SRanges...>(*this, inds...);
-    }
-    
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::isInit() const
-    {
-	return mInit;
-    }
-
-    /****************************************
-     *  MutableMultiArrayBase::iterator     *	     
-     ****************************************/
-    
-    template <typename T, class... SRanges>
-    MutableMultiArrayBase<T,SRanges...>::iterator::iterator(MutableMultiArrayBase<T,SRanges...>& ma):
-	mMAPtr(&ma), mPos(0)
-    { }
-    
-    template <typename T, class... SRanges>
-    MutableMultiArrayBase<T,SRanges...>::iterator::iterator(MutableMultiArrayBase<T,SRanges...>& ma,
-							const typename CRange::IndexType& index):
-	mMAPtr(&ma), mPos(index.pos())
-    { }
-    
-    template <typename T, class... SRanges>
-    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator==(const iterator& it) const
-    {
-	return mMAPtr == it.mMAPtr and mPos == it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator!=(const iterator& it) const
-    {
-	return mMAPtr != it.mMAPtr or mPos != it.mPos;
-    }
-    
-    template <typename T, class... SRanges>
-    const T& MutableMultiArrayBase<T,SRanges...>::iterator::operator*() const
-    {
-	return mMAPtr->data()[mPos];
-    }
-
-    template <typename T, class... SRanges>
-    T const* MutableMultiArrayBase<T,SRanges...>::iterator::operator->() const
-    {
-	return &mMAPtr->data()[mPos];
-    }
-
-    template <typename T, class... SRanges>
-    T& MutableMultiArrayBase<T,SRanges...>::iterator::operator*()
-    {
-	return mMAPtr->data()[mPos];
-    }
-
-    template <typename T, class... SRanges>
-    T* MutableMultiArrayBase<T,SRanges...>::iterator::operator->()
-    {
-	return &mMAPtr->data()[mPos];
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator++()
-    {
-	++mPos;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator++(int)
-    {
-	iterator tmp(*this);
-	++mPos;
-	return tmp;
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator--()
-    {
-	--mPos;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator--(int)
-    {
-	iterator tmp(*this);
-	--mPos;
-	return tmp;
-    }
-    
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator+=(int diff)
-    {
-	mPos += diff;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator-=(int diff)
-    {
-	mPos -= diff;
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator+(int num) const
-    {
-	iterator tmp(*this);
-	tmp += num;
-	return tmp;
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator-(int num) const
-    {
-	iterator tmp(*this);
-	tmp -= num;
-    }
-    
-    template <typename T, class... SRanges>
-    int MutableMultiArrayBase<T,SRanges...>::iterator::operator-(const iterator& it) const
-    {
-	return mPos - it.mPos;
-    }
-    
-    template <typename T, class... SRanges>
-    const T& MutableMultiArrayBase<T,SRanges...>::iterator::operator[](int num) const
-    {
-	return *(operator+(num));
-    }
-
-    template <typename T, class... SRanges>
-    T& MutableMultiArrayBase<T,SRanges...>::iterator::operator[](int num)
-    {
-	return *(operator+(num));
-    }
-    
-    template <typename T, class... SRanges>
-    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator<(const iterator& it) const
-    {
-	return mPos < it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator>(const iterator& it) const
-    {
-	return mPos > it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator<=(const iterator& it) const
-    {
-	return mPos <= it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator>=(const iterator& it) const
-    {
-	return mPos >= it.mPos;
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::IndexType
-    MutableMultiArrayBase<T,SRanges...>::iterator::index() const
-    {
-	auto i = mMAPtr->beginIndex();
-	i = mPos;
-	return i;
-    }    
-    
-    /******************************
-     *  MutableMultiArrayBase     *	     
-     ******************************/
-
-    template <typename T, class... SRanges>
-    MutableMultiArrayBase<T,SRanges...>::MutableMultiArrayBase(const std::shared_ptr<SRanges>&... ranges) :
-	MultiArrayBase<T,SRanges...>(ranges...) {}
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::begin()
-    {
-	return iterator(*this, MAB::beginIndex());
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::end()
-    {
-	return iterator(*this, MAB::endIndex());
-    }
-    
-    template <typename T, class... SRanges>
-    bool MutableMultiArrayBase<T,SRanges...>::isConst() const
-    {
-	return false;
-    }
-
-    template <typename T, class... SRanges>
-    OperationRoot<T,SRanges...>
-    MutableMultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds)
-    {
-	return OperationRoot<T,SRanges...>(*this, inds...);
-    }
-
-    template <typename T, class... SRanges>
-    ConstOperationRoot<T,SRanges...>
-    MutableMultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
-    {
-	return ConstOperationRoot<T,SRanges...>(*this, inds...);
-    }
-
-    
-    /*******************
-     *  MultiArray     *	     
-     *******************/
-
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges) :
-	MutableMultiArrayBase<T,SRanges...>(ranges...),
-	mCont(MAB::mRange->size())
-    {
-	MAB::mInit = true;
-    }
-
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec) :
-	MutableMultiArrayBase<T,SRanges...>(ranges...),
-	mCont(vec)
-    {
-	MAB::mInit = true;
-	if(mCont.size() > MAB::mRange->size()){
-	    mCont.erase(mCont.begin() + MAB::mRange->size(), mCont.end());
-	}
-    }
-    
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec) :
-	MutableMultiArrayBase<T,SRanges...>(ranges...),
-	mCont(vec)
-    {
-	MAB::mInit = true;
-	if(mCont.size() > MAB::mRange->size()){
-	    mCont.erase(mCont.begin() + MAB::mRange->size(), mCont.end());
-	}
-    }
-    /*
-    template <typename T, class... SRanges>
-    template <class Range2, class Range3>
-    MultiArray<T,SRanges...>::MultiArray(const MultiArray<MultiArray<T,Range2>,Range3> in) :
-	MutableMultiArrayBase<T,SRanges...>(merge(in.range(), in[ in.beginIndex() ].range()))
-	// assert that Range2 has always same extension
-    {
-	MAB::mInit = true;
-	mCont.clear();
-	for(auto i = in.beginIndex(); i != in.endIndex(); ++i){
-	    mCont.insert(mCont.end(), in[i].mCont.begin(), in[i].mCont.end());
-	}
-	assert(mCont.size() == MAB::mRange->size());
-    }
-    */
-    /*
-    template <typename T, class... SRanges>
-    template <class Range2, class Range3>
-    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator=(const MultiArray<MultiArray<T,Range2>,Range3> in)
-    {
-	MAB::mRange.reset(new Range(merge(in.range(), in[ in.beginIndex() ].range())));
-	// assert that Range2 has always same extension
-	mCont.clear();
-	for(auto i = in.beginIndex(); i != in.endIndex(); ++i){
-	    mCont.insert(mCont.end(), in[i].mCont.begin(), in[i].mCont.end());
-	}
-	assert(mCont.size() == MAB::mRange->size());
-	return *this;
-	} */   
-
-    template <typename T, class... SRanges>
-    T& MultiArray<T,SRanges...>::operator[](const typename CRange::IndexType& i)
-    {
-	return mCont[ i.pos() ];
-    }
-
-    template <typename T, class... SRanges>
-    const T& MultiArray<T,SRanges...>::operator[](const typename CRange::IndexType& i) const
-    {
-	return mCont[ i.pos() ];
-    }
-
-    template <typename T, class... SRanges>
-    T& MultiArray<T,SRanges...>::at(const typename CRange::IndexType::MetaType& meta)
-    {
-	return mCont[ MAB::beginIndex().at(meta).pos() ];
-    }
-
-    template <typename T, class... SRanges>
-    const T& MultiArray<T,SRanges...>::at(const typename CRange::IndexType::MetaType& meta) const
-    {
-	return mCont[ MAB::beginIndex().at(meta).pos() ];
-    }
-	
-    template <typename T, class... SRanges>
-    bool MultiArray<T,SRanges...>::isConst() const
-    {
-	return false;
-    }
-    
-    template <typename T, class... SRanges>
-    bool MultiArray<T,SRanges...>::isSlice() const
-    {
-	return false;
-    }
-
-    template <typename T, class... SRanges>
-    template <class... SRanges2>
-    MultiArray<T,SRanges2...> MultiArray<T,SRanges...>::format(const std::shared_ptr<SRanges2>&... nrs)
-    {
-	return MultiArray<T,SRanges2...>( nrs... , std::move(mCont) );
-    }
-    
-    template <typename T, class... SRanges>
-    const T* MultiArray<T,SRanges...>::data() const
-    {
-	return mCont.data();
-    }
-
-    template <typename T, class... SRanges>
-    T* MultiArray<T,SRanges...>::data()
-    {
-	return mCont.data();
-    }
-
-    template <typename T, class... SRanges>
-    const std::vector<T>& MultiArray<T,SRanges...>::datav() const
-    {
-	return mCont;
-    }
-
-    template <typename T, class... SRanges>
-    std::vector<T>& MultiArray<T,SRanges...>::datav()
-    {
-	return mCont;
-    }
-
-    
-    /*
-    template <typename T, class... SRanges>
-    void MultiArray<T,SRanges...>::manipulate(ManipulatorBase<T>& mb,
-					 const typename Range::IndexType& manBegin,
-					 const typename Range::IndexType& manEnd)
-    {
-	mb.setup(mCont, manBegin.pos(), manEnd.pos());
-	mb.execute();
-    }
-    */
-
-
-    /****************************
-     *   FunctionalMultiArray   *	     
-     ****************************/
-
-    /*
-    template <typename T, class Range, class Function>
-    FunctionalMultiArray<T,Range,Function>::FunctionalMultiArray(const Range& range) :
-	MultiArrayBase<T,SRanges...>(range), mFunc() {}
-    */
-    template <typename T, class Function, class... SRanges>
-    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges,
-								  const Function& func) :
-	MultiArrayBase<T,SRanges...>(ranges...), mFunc(func) {}
-
-    template <typename T, class Function, class... SRanges>
-    const T& FunctionalMultiArray<T,Function,SRanges...>::operator[](const typename CRange::IndexType& i) const
-    {
-	mVal = mFunc(i);
-	return mVal;
-    }
-
-    template <typename T, class Function, class... SRanges>
-    bool FunctionalMultiArray<T,Function,SRanges...>::isConst() const
-    {
-	return true;
-    }
-
-    template <typename T, class Function, class... SRanges>
-    bool FunctionalMultiArray<T,Function,SRanges...>::isSlice() const
-    {
-	return false;
-    }
-    
-}
diff --git a/src/multi_array.h b/src/multi_array.h
index 4e37e70..90bedcd 100644
--- a/src/multi_array.h
+++ b/src/multi_array.h
@@ -264,6 +264,600 @@ namespace MultiArrayTools
     
 }
 
-#include "multi_array.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    /**************************************
+     *  MultiArrayBase::const_iterator    *	     
+     **************************************/
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::const_iterator::const_iterator(const MultiArrayBase<T,SRanges...>& ma):
+	mMAPtr(&ma), mPos(0) { }
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::const_iterator::const_iterator(const MultiArrayBase<T,SRanges...>& ma,
+							     const typename CRange::IndexType& index):
+	mMAPtr(&ma), mPos(index.pos()) { }
+    
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator==(const const_iterator& it) const
+    {
+	return mMAPtr == it.mMAPtr and mPos == it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator!=(const const_iterator& it) const
+    {
+	return mMAPtr != it.mMAPtr or mPos != it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MultiArrayBase<T,SRanges...>::const_iterator::operator*() const
+    {
+	return mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T const* MultiArrayBase<T,SRanges...>::const_iterator::operator->() const
+    {
+	return &mMAPtr->data()[mPos];
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator++()
+    {
+	++mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator++(int)
+    {
+	const_iterator tmp(*this);
+	++mPos;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator--()
+    {
+	--mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator--(int)
+    {
+	const_iterator tmp(*this);
+	--mPos;
+	return tmp;
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator+=(int diff)
+    {
+	mPos += diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator& MultiArrayBase<T,SRanges...>::const_iterator::operator-=(int diff)
+    {
+	mPos -= diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator+(int num) const
+    {
+	const_iterator tmp(*this);
+	tmp += num;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::const_iterator::operator-(int num) const
+    {
+	const_iterator tmp(*this);
+	tmp -= num;
+    }
+    
+    template <typename T, class... SRanges>
+    int MultiArrayBase<T,SRanges...>::const_iterator::operator-(const const_iterator& it) const
+    {
+	return mPos - it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MultiArrayBase<T,SRanges...>::const_iterator::operator[](int num) const
+    {
+	return *(operator+(num));
+    }
+    
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator<(const const_iterator& it) const
+    {
+	return mPos < it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator>(const const_iterator& it) const
+    {
+	return mPos > it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator<=(const const_iterator& it) const
+    {
+	return mPos <= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::const_iterator::operator>=(const const_iterator& it) const
+    {
+	return mPos >= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::const_iterator::index() const
+    {
+	auto i = mMAPtr->beginIndex();
+	i = mPos;
+	return i;
+    }
+    
+    /**********************
+     *  MultiArrayBase    *	     
+     **********************/
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::MultiArrayBase(const std::shared_ptr<SRanges>&... ranges)
+    {
+	ContainerRangeFactory<SRanges...> crf(ranges...);
+	mRange = std::dynamic_pointer_cast<ContainerRange<SRanges...> >( crf.create() );
+    }
+    
+    template <typename T, class... SRanges>
+    size_t MultiArrayBase<T,SRanges...>::size() const
+    {
+	return mRange->size();
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::begin() const
+    {
+	return const_iterator(*this, beginIndex());
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::const_iterator MultiArrayBase<T,SRanges...>::end() const
+    {
+	return const_iterator(*this, endIndex());
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::beginIndex() const
+    {
+	return mRange->begin();
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::endIndex() const
+    {
+	return mRange->end();
+    }
+
+    template <typename T, class... SRanges>
+    const std::shared_ptr<typename MultiArrayBase<T,SRanges...>::CRange>&
+    MultiArrayBase<T,SRanges...>::range() const
+    {
+	return mRange;
+    }
+
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::isConst() const
+    {
+	return true;
+    }
+    
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	return ConstOperationRoot<T,SRanges...>(*this, inds...);
+    }
+    
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::isInit() const
+    {
+	return mInit;
+    }
+
+    /****************************************
+     *  MutableMultiArrayBase::iterator     *	     
+     ****************************************/
+    
+    template <typename T, class... SRanges>
+    MutableMultiArrayBase<T,SRanges...>::iterator::iterator(MutableMultiArrayBase<T,SRanges...>& ma):
+	mMAPtr(&ma), mPos(0)
+    { }
+    
+    template <typename T, class... SRanges>
+    MutableMultiArrayBase<T,SRanges...>::iterator::iterator(MutableMultiArrayBase<T,SRanges...>& ma,
+							const typename CRange::IndexType& index):
+	mMAPtr(&ma), mPos(index.pos())
+    { }
+    
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator==(const iterator& it) const
+    {
+	return mMAPtr == it.mMAPtr and mPos == it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator!=(const iterator& it) const
+    {
+	return mMAPtr != it.mMAPtr or mPos != it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MutableMultiArrayBase<T,SRanges...>::iterator::operator*() const
+    {
+	return mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T const* MutableMultiArrayBase<T,SRanges...>::iterator::operator->() const
+    {
+	return &mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T& MutableMultiArrayBase<T,SRanges...>::iterator::operator*()
+    {
+	return mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    T* MutableMultiArrayBase<T,SRanges...>::iterator::operator->()
+    {
+	return &mMAPtr->data()[mPos];
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator++()
+    {
+	++mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator++(int)
+    {
+	iterator tmp(*this);
+	++mPos;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator--()
+    {
+	--mPos;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator--(int)
+    {
+	iterator tmp(*this);
+	--mPos;
+	return tmp;
+    }
+    
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator+=(int diff)
+    {
+	mPos += diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator& MutableMultiArrayBase<T,SRanges...>::iterator::operator-=(int diff)
+    {
+	mPos -= diff;
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator+(int num) const
+    {
+	iterator tmp(*this);
+	tmp += num;
+	return tmp;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::iterator::operator-(int num) const
+    {
+	iterator tmp(*this);
+	tmp -= num;
+    }
+    
+    template <typename T, class... SRanges>
+    int MutableMultiArrayBase<T,SRanges...>::iterator::operator-(const iterator& it) const
+    {
+	return mPos - it.mPos;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& MutableMultiArrayBase<T,SRanges...>::iterator::operator[](int num) const
+    {
+	return *(operator+(num));
+    }
+
+    template <typename T, class... SRanges>
+    T& MutableMultiArrayBase<T,SRanges...>::iterator::operator[](int num)
+    {
+	return *(operator+(num));
+    }
+    
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator<(const iterator& it) const
+    {
+	return mPos < it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator>(const iterator& it) const
+    {
+	return mPos > it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator<=(const iterator& it) const
+    {
+	return mPos <= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::iterator::operator>=(const iterator& it) const
+    {
+	return mPos >= it.mPos;
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::IndexType
+    MutableMultiArrayBase<T,SRanges...>::iterator::index() const
+    {
+	auto i = mMAPtr->beginIndex();
+	i = mPos;
+	return i;
+    }    
+    
+    /******************************
+     *  MutableMultiArrayBase     *	     
+     ******************************/
+
+    template <typename T, class... SRanges>
+    MutableMultiArrayBase<T,SRanges...>::MutableMultiArrayBase(const std::shared_ptr<SRanges>&... ranges) :
+	MultiArrayBase<T,SRanges...>(ranges...) {}
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::begin()
+    {
+	return iterator(*this, MAB::beginIndex());
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::iterator MutableMultiArrayBase<T,SRanges...>::end()
+    {
+	return iterator(*this, MAB::endIndex());
+    }
+    
+    template <typename T, class... SRanges>
+    bool MutableMultiArrayBase<T,SRanges...>::isConst() const
+    {
+	return false;
+    }
+
+    template <typename T, class... SRanges>
+    OperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds)
+    {
+	return OperationRoot<T,SRanges...>(*this, inds...);
+    }
+
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	return ConstOperationRoot<T,SRanges...>(*this, inds...);
+    }
+
+    
+    /*******************
+     *  MultiArray     *	     
+     *******************/
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges) :
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mCont(MAB::mRange->size())
+    {
+	MAB::mInit = true;
+    }
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec) :
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mCont(vec)
+    {
+	MAB::mInit = true;
+	if(mCont.size() > MAB::mRange->size()){
+	    mCont.erase(mCont.begin() + MAB::mRange->size(), mCont.end());
+	}
+    }
+    
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec) :
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mCont(vec)
+    {
+	MAB::mInit = true;
+	if(mCont.size() > MAB::mRange->size()){
+	    mCont.erase(mCont.begin() + MAB::mRange->size(), mCont.end());
+	}
+    }
+    /*
+    template <typename T, class... SRanges>
+    template <class Range2, class Range3>
+    MultiArray<T,SRanges...>::MultiArray(const MultiArray<MultiArray<T,Range2>,Range3> in) :
+	MutableMultiArrayBase<T,SRanges...>(merge(in.range(), in[ in.beginIndex() ].range()))
+	// assert that Range2 has always same extension
+    {
+	MAB::mInit = true;
+	mCont.clear();
+	for(auto i = in.beginIndex(); i != in.endIndex(); ++i){
+	    mCont.insert(mCont.end(), in[i].mCont.begin(), in[i].mCont.end());
+	}
+	assert(mCont.size() == MAB::mRange->size());
+    }
+    */
+    /*
+    template <typename T, class... SRanges>
+    template <class Range2, class Range3>
+    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator=(const MultiArray<MultiArray<T,Range2>,Range3> in)
+    {
+	MAB::mRange.reset(new Range(merge(in.range(), in[ in.beginIndex() ].range())));
+	// assert that Range2 has always same extension
+	mCont.clear();
+	for(auto i = in.beginIndex(); i != in.endIndex(); ++i){
+	    mCont.insert(mCont.end(), in[i].mCont.begin(), in[i].mCont.end());
+	}
+	assert(mCont.size() == MAB::mRange->size());
+	return *this;
+	} */   
+
+    template <typename T, class... SRanges>
+    T& MultiArray<T,SRanges...>::operator[](const typename CRange::IndexType& i)
+    {
+	return mCont[ i.pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    const T& MultiArray<T,SRanges...>::operator[](const typename CRange::IndexType& i) const
+    {
+	return mCont[ i.pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    T& MultiArray<T,SRanges...>::at(const typename CRange::IndexType::MetaType& meta)
+    {
+	return mCont[ MAB::beginIndex().at(meta).pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    const T& MultiArray<T,SRanges...>::at(const typename CRange::IndexType::MetaType& meta) const
+    {
+	return mCont[ MAB::beginIndex().at(meta).pos() ];
+    }
+	
+    template <typename T, class... SRanges>
+    bool MultiArray<T,SRanges...>::isConst() const
+    {
+	return false;
+    }
+    
+    template <typename T, class... SRanges>
+    bool MultiArray<T,SRanges...>::isSlice() const
+    {
+	return false;
+    }
+
+    template <typename T, class... SRanges>
+    template <class... SRanges2>
+    MultiArray<T,SRanges2...> MultiArray<T,SRanges...>::format(const std::shared_ptr<SRanges2>&... nrs)
+    {
+	return MultiArray<T,SRanges2...>( nrs... , std::move(mCont) );
+    }
+    
+    template <typename T, class... SRanges>
+    const T* MultiArray<T,SRanges...>::data() const
+    {
+	return mCont.data();
+    }
+
+    template <typename T, class... SRanges>
+    T* MultiArray<T,SRanges...>::data()
+    {
+	return mCont.data();
+    }
+
+    template <typename T, class... SRanges>
+    const std::vector<T>& MultiArray<T,SRanges...>::datav() const
+    {
+	return mCont;
+    }
+
+    template <typename T, class... SRanges>
+    std::vector<T>& MultiArray<T,SRanges...>::datav()
+    {
+	return mCont;
+    }
+
+    
+    /*
+    template <typename T, class... SRanges>
+    void MultiArray<T,SRanges...>::manipulate(ManipulatorBase<T>& mb,
+					 const typename Range::IndexType& manBegin,
+					 const typename Range::IndexType& manEnd)
+    {
+	mb.setup(mCont, manBegin.pos(), manEnd.pos());
+	mb.execute();
+    }
+    */
+
+
+    /****************************
+     *   FunctionalMultiArray   *	     
+     ****************************/
+
+    /*
+    template <typename T, class Range, class Function>
+    FunctionalMultiArray<T,Range,Function>::FunctionalMultiArray(const Range& range) :
+	MultiArrayBase<T,SRanges...>(range), mFunc() {}
+    */
+    template <typename T, class Function, class... SRanges>
+    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges,
+								  const Function& func) :
+	MultiArrayBase<T,SRanges...>(ranges...), mFunc(func) {}
+
+    template <typename T, class Function, class... SRanges>
+    const T& FunctionalMultiArray<T,Function,SRanges...>::operator[](const typename CRange::IndexType& i) const
+    {
+	mVal = mFunc(i);
+	return mVal;
+    }
+
+    template <typename T, class Function, class... SRanges>
+    bool FunctionalMultiArray<T,Function,SRanges...>::isConst() const
+    {
+	return true;
+    }
+
+    template <typename T, class Function, class... SRanges>
+    bool FunctionalMultiArray<T,Function,SRanges...>::isSlice() const
+    {
+	return false;
+    }
+    
+}
 
 #endif
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
deleted file mode 100644
index bad3d9f..0000000
--- a/src/multi_array_operation.cc
+++ /dev/null
@@ -1,377 +0,0 @@
-// -*- C++ -*-
-
-#include "multi_array_operation.h"
-
-namespace MultiArrayTools
-{
-    namespace
-    {
-	using namespace MultiArrayHelper;
-    }
-
-    void seekIndexInst(std::shared_ptr<IndexBase> i, std::vector<std::shared_ptr<IndexBase> >& ivec)
-    {
-	for(size_t inum = 0; inum != i->rangePtr()->dim(); ++inum){
-	    auto ii = i->getPtr(inum);
-	    if(ii->type() == IndexType::MULTI or
-	       ii->type() == IndexType::CONT){
-		seekIndexInst(ii, ivec);
-	    }
-	    ivec.push_back(ii);
-	}
-    }
-
-    BTSS getBlockType(std::shared_ptr<IndexBase> i,
-		      std::shared_ptr<IndexBase> j,
-		      bool first, size_t higherStepSize)
-    {
-	// returning BlockType and step size is redundant (change in the future)
-	// stepSize == 0 => VALUE
-	// stepSize == 1 => BLOCK
-	// stepSize > 1 => SPLIT :)
-	BTSS out(BlockType::VALUE, 0);
-	size_t lastNum = i->rangePtr()->dim();
-	for(size_t inum = 0; inum != lastNum; ++inum){
-	    auto ii = i->getPtr(inum);
-	    if(ii == j){
-		
-		if(inum == lastNum - 1 and first){
-		    out = BTSS(BlockType::BLOCK, 1);
-		}
-		else {
-		    first = false;
-		    out = BTSS(BlockType::SPLIT, i->getStepSize(inum) * higherStepSize + out.second);
-		}
-		continue;
-	    }
-	    
-	    if(ii->type() == IndexType::MULTI or
-	       ii->type() == IndexType::CONT){
-
-		BTSS tmp = getBlockType(ii, j, inum == lastNum - 1, i->getStepSize(inum) * higherStepSize);
-		if(tmp.first != BlockType::VALUE){
-		    out = tmp;
-		}
-	    }
-	}
-	return out;
-    }
-
-    template <typename T>
-    std::shared_ptr<Block<T> > makeBlock(const std::vector<T>& vec, size_t stepSize, size_t blockSize)
-    {
-	return std::make_shared<Block<T> >(vec, 0, blockSize, stepSize);
-    }
-
-    template <typename T>
-    std::shared_ptr<MBlock<T> > makeBlock(std::vector<T>& vec, size_t stepSize, size_t blockSize)
-    {
-	return std::make_shared<MBlock<T> >(vec, 0, blockSize, stepSize);
-    }
-
-    size_t getBTNum(const std::vector<BTSS>& mp, BlockType bt)
-    {
-	size_t out = 0;
-	for(auto& xx: mp){
-	    if(xx.first == bt){
-		++out;
-	    }
-	}
-	return out;
-    }
-
-    void minimizeAppearanceOfType(std::map<std::shared_ptr<IndexBase>, std::vector<BTSS> >& mp,
-				  BlockType bt)
-    {
-	size_t minNum = getBTNum( mp.begin()->second, bt );
-	for(auto& mm: mp){
-	    size_t tmp = getBTNum( mm.second, bt );
-	    if(tmp < minNum){
-		minNum = tmp;
-	    }
-	}
-
-	for(auto mit = mp.begin(); mit != mp.end(); ){
-	    size_t tmp = getBTNum( mit->second, bt );
-	    if(tmp > minNum){
-		mit = mp.erase(mit);
-	    }
-	    else {
-		++mit;
-	    }
-	}
-
-    }
-    
-    template <class OpClass>
-    std::shared_ptr<IndexBase> seekBlockIndex(std::shared_ptr<IndexBase> ownIdx,
-					      const OpClass& second)
-    {
-	std::vector<std::shared_ptr<IndexBase> > ivec;
-	seekIndexInst(ownIdx, ivec);
-	std::map<std::shared_ptr<IndexBase>, std::vector<BTSS> > mp;
-	
-	for(auto& xx: ivec){
-	    mp[xx] = second.block(xx);
-	}
-	// seek minimal number of VALUEs => guarantees absence of conflicting blocks
-	minimizeAppearanceOfType(mp, BlockType::VALUE);
-	// seek mininmal number of SPLITs => maximal vectorization possible
-	minimizeAppearanceOfType(mp, BlockType::SPLIT);
-	
-	return mp.begin()->first;
-    }
-
-    /***************************
-     *   OperationTemplate     *
-     ***************************/
-
-    template <typename T, class OperationClass>
-    OperationTemplate<T,OperationClass>::OperationTemplate(OperationClass* oc) : mOc(oc) {}
-    
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationTemplate<T,OperationClass>::operator+(const Second& in) const
-	-> Operation<T,std::plus<T>,OperationClass,Second>
-    {
-    	return Operation<T,std::plus<T>,OperationClass,Second>(*mOc, in);
-    }
-
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationTemplate<T,OperationClass>::operator-(const Second& in) const
-	-> Operation<T,std::minus<T>,OperationClass,Second>
-    {
-    	return Operation<T,std::minus<T>,OperationClass,Second>(*mOc, in);
-    }
-    
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationTemplate<T,OperationClass>::operator*(const Second& in) const
-	-> Operation<T,std::multiplies<T>,OperationClass,Second>
-    {
-    	return Operation<T,std::multiplies<T>,OperationClass,Second>(*mOc, in);
-    }
-
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationTemplate<T,OperationClass>::operator/(const Second& in) const
-	-> Operation<T,std::divides<T>,OperationClass,Second>
-    {
-    	return Operation<T,std::divides<T>,OperationClass,Second>(*mOc, in);
-    }
-
-    template <typename T, class OperationClass>
-    template <class IndexType>
-    auto OperationTemplate<T,OperationClass>::c(std::shared_ptr<IndexType>& ind) const
-	-> Contraction<T,OperationClass,IndexType>
-    {
-	return Contraction<T,OperationClass,IndexType>(*mOc, ind);
-    }
-
-    
-    /*************************
-     *   OperationMaster     *
-     *************************/
-
-    template <typename T, class OpClass, class... Ranges>
-    OperationMaster<T,OpClass,Ranges...>::
-    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
-		    std::shared_ptr<typename CRange::IndexType>& index) :
-	mSecond(second), mArrayRef(ma), mIndex()
-    {
-	MultiRangeFactory<Ranges...> mrf( index->range() );
-	std::shared_ptr<MultiRange<Ranges...> > mr =
-	    std::dynamic_pointer_cast<MultiRange<Ranges...> >( mrf.create() );
-	mIndex = std::make_shared<IndexType>( mr->begin() );
-	(*mIndex) = *index;
-
-	auto blockIndex = seekBlockIndex( mIndex, second);
-
-	block(blockIndex);
-	second.block(blockIndex);
-	
-	for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndex) ){
-	    get() = mSecond.get();
-	}
-    }
-    
-    template <typename T, class OpClass, class... Ranges>
-    MBlock<T>& OperationMaster<T,OpClass,Ranges...>::get()
-    {
-	block();
-	return *mBlockPtr;
-    }
-    
-    template <typename T, class OpClass, class... Ranges>
-    const Block<T>& OperationMaster<T,OpClass,Ranges...>::get() const
-    {
-	block();
-	return *mBlockPtr;
-    }
-
-    template <typename T, class OpClass, class... Ranges>
-    std::vector<BTSS> OperationMaster<T,OpClass,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
-    {
-	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
-	mBlockPtr = makeBlock(mArrayRef.datav(), btv[0].second, blockIndex->max());
-	return btv;
-    }
-
-    template <typename T, class OpClass, class... Ranges>
-    const OperationMaster<T,OpClass,Ranges...>& OperationMaster<T,OpClass,Ranges...>::block() const
-    {
-	mBlockPtr->set( mIndex->pos() );
-	return *this;
-    }
-    
-    /****************************
-     *   ConstOperationRoot     *
-     ****************************/
-
-    template <typename T, class... Ranges>
-    ConstOperationRoot<T,Ranges...>::
-    ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
-		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	OperationTemplate<T,ConstOperationRoot<T,Ranges...> >(this),
-	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
-    {
-	(*mIndex)(indices...);
-    }
-
-    template <typename T, class... Ranges>
-    const Block<T>& ConstOperationRoot<T,Ranges...>::get() const
-    {
-	block();
-	return *mBlockPtr;
-    }
-
-    template <typename T, class... Ranges>
-    std::vector<BTSS> ConstOperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
-    {
-	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
-	mBlockPtr = makeBlock(mArrayRef.datav(), btv[0].second, blockIndex->max());
-	return btv;
-    }
-
-    template <typename T, class... Ranges>
-    const ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::block() const
-    {
-	mBlockPtr->set( (*mIndex)().pos() );
-	return *this;
-    }
-    
-    /***********************
-     *   OperationRoot     *
-     ***********************/
-    
-    template <typename T, class... Ranges>
-    OperationRoot<T,Ranges...>::
-    OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
-		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	OperationTemplate<T,OperationRoot<T,Ranges...> >(this),
-	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
-    {
-	(*mIndex)(indices...);
-    }
-
-    template <typename T, class... Ranges>
-    template <class OpClass>
-    OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
-    {
-	return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex);
-    }
-
-    template <typename T, class... Ranges>
-    const MBlock<T>& OperationRoot<T,Ranges...>::get() const
-    {
-	block();
-	return *mBlockPtr;
-    }
-    
-    template <typename T, class... Ranges>
-    MBlock<T>& OperationRoot<T,Ranges...>::get()
-    {
-	block();
-	return *mBlockPtr;
-    }
-
-    template <typename T, class... Ranges>
-    std::vector<BTSS> OperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
-    {
-	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
-	mBlockPtr = makeBlock(mArrayRef.datav(), btv[0].second, blockIndex->max());
-	return btv;
-    }
-
-    template <typename T, class... Ranges>
-    const OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::block() const
-    {
-	mBlockPtr->set( (*mIndex)().pos() );
-	return *this;
-    }
-    
-    /*******************
-     *   Operation     *
-     *******************/
-    
-    template <typename T, class OpFunction, class... Ops>
-    Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
-	OperationTemplate<T,Operation<T,OpFunction,Ops...> >(this),
-	mOps(ops...) {}
-
-    template <typename T, class OpFunction, class... Ops>
-    const BlockResult<T>& Operation<T,OpFunction,Ops...>::get() const
-    {
-	mRes = std::move( PackNum<sizeof...(Ops)-1>::template unpackArgs<T,OpFunction>(mOps) );
-	return mRes;
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    std::vector<BTSS> Operation<T,OpFunction,Ops...>::block(const std::shared_ptr<IndexBase> blockIndex) const
-    {
-	std::vector<BTSS> btv;
-	PackNum<sizeof...(Ops)-1>::makeBlockTypeVec(btv, mOps, blockIndex);
-	return btv;
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    const Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::block() const
-    {
-	//mBlockPtr->set( mIndex->pos() );
-	return *this;
-    }
-
-    /*********************
-     *   Contraction     *
-     *********************/
-
-    template <typename T, class Op, class IndexType>
-    Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
-	OperationTemplate<T,Contraction<T,Op,IndexType> >(this),
-	mOp(op),
-	mInd(ind) {}
-
-    template <typename T, class Op, class IndexType>
-    const BlockResult<T>& Contraction<T,Op,IndexType>::get() const
-    {
-	BlockBinaryOpSelf<T,std::plus<T>,decltype(mOp.get())> f(mRes);
-	for(*mInd = 0; mInd->pos() != mInd->max(); ++(*mInd)){
-	    f(mOp.get());
-	}
-	return mRes;
-    }
-    
-    template <typename T, class Op, class IndexType>
-    std::vector<BTSS> Contraction<T,Op,IndexType>::block(const std::shared_ptr<IndexBase> blockIndex) const
-    {
-	return mOp.block(blockIndex);
-    }
-
-    template <typename T, class Op, class IndexType>
-    const Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::block() const
-    {
-	return *this;
-    }
-    
-}
diff --git a/src/multi_array_operation.h b/src/multi_array_operation.h
index 6e17fa0..90537b7 100644
--- a/src/multi_array_operation.h
+++ b/src/multi_array_operation.h
@@ -251,6 +251,382 @@ namespace MultiArrayTools
     
 }
 
-#include "multi_array_operation.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+
+    void seekIndexInst(std::shared_ptr<IndexBase> i, std::vector<std::shared_ptr<IndexBase> >& ivec)
+    {
+	for(size_t inum = 0; inum != i->rangePtr()->dim(); ++inum){
+	    auto ii = i->getPtr(inum);
+	    if(ii->type() == IndexType::MULTI or
+	       ii->type() == IndexType::CONT){
+		seekIndexInst(ii, ivec);
+	    }
+	    ivec.push_back(ii);
+	}
+    }
+
+    BTSS getBlockType(std::shared_ptr<IndexBase> i,
+		      std::shared_ptr<IndexBase> j,
+		      bool first, size_t higherStepSize)
+    {
+	// returning BlockType and step size is redundant (change in the future)
+	// stepSize == 0 => VALUE
+	// stepSize == 1 => BLOCK
+	// stepSize > 1 => SPLIT :)
+	BTSS out(BlockType::VALUE, 0);
+	size_t lastNum = i->rangePtr()->dim();
+	for(size_t inum = 0; inum != lastNum; ++inum){
+	    auto ii = i->getPtr(inum);
+	    if(ii == j){
+		
+		if(inum == lastNum - 1 and first){
+		    out = BTSS(BlockType::BLOCK, 1);
+		}
+		else {
+		    first = false;
+		    out = BTSS(BlockType::SPLIT, i->getStepSize(inum) * higherStepSize + out.second);
+		}
+		continue;
+	    }
+	    
+	    if(ii->type() == IndexType::MULTI or
+	       ii->type() == IndexType::CONT){
+
+		BTSS tmp = getBlockType(ii, j, inum == lastNum - 1, i->getStepSize(inum) * higherStepSize);
+		if(tmp.first != BlockType::VALUE){
+		    out = tmp;
+		}
+	    }
+	}
+	return out;
+    }
+
+    template <typename T>
+    std::shared_ptr<Block<T> > makeBlock(const std::vector<T>& vec, size_t stepSize, size_t blockSize)
+    {
+	return std::make_shared<Block<T> >(vec, 0, blockSize, stepSize);
+    }
+
+    template <typename T>
+    std::shared_ptr<MBlock<T> > makeBlock(std::vector<T>& vec, size_t stepSize, size_t blockSize)
+    {
+	return std::make_shared<MBlock<T> >(vec, 0, blockSize, stepSize);
+    }
+
+    size_t getBTNum(const std::vector<BTSS>& mp, BlockType bt)
+    {
+	size_t out = 0;
+	for(auto& xx: mp){
+	    if(xx.first == bt){
+		++out;
+	    }
+	}
+	return out;
+    }
+
+    void minimizeAppearanceOfType(std::map<std::shared_ptr<IndexBase>, std::vector<BTSS> >& mp,
+				  BlockType bt)
+    {
+	size_t minNum = getBTNum( mp.begin()->second, bt );
+	for(auto& mm: mp){
+	    size_t tmp = getBTNum( mm.second, bt );
+	    if(tmp < minNum){
+		minNum = tmp;
+	    }
+	}
+
+	for(auto mit = mp.begin(); mit != mp.end(); ){
+	    size_t tmp = getBTNum( mit->second, bt );
+	    if(tmp > minNum){
+		mit = mp.erase(mit);
+	    }
+	    else {
+		++mit;
+	    }
+	}
+
+    }
+    
+    template <class OpClass>
+    std::shared_ptr<IndexBase> seekBlockIndex(std::shared_ptr<IndexBase> ownIdx,
+					      const OpClass& second)
+    {
+	std::vector<std::shared_ptr<IndexBase> > ivec;
+	seekIndexInst(ownIdx, ivec);
+	std::map<std::shared_ptr<IndexBase>, std::vector<BTSS> > mp;
+	
+	for(auto& xx: ivec){
+	    mp[xx] = second.block(xx);
+	}
+	// seek minimal number of VALUEs => guarantees absence of conflicting blocks
+	minimizeAppearanceOfType(mp, BlockType::VALUE);
+	// seek mininmal number of SPLITs => maximal vectorization possible
+	minimizeAppearanceOfType(mp, BlockType::SPLIT);
+	
+	return mp.begin()->first;
+    }
+
+    /***************************
+     *   OperationTemplate     *
+     ***************************/
+
+    template <typename T, class OperationClass>
+    OperationTemplate<T,OperationClass>::OperationTemplate(OperationClass* oc) : mOc(oc) {}
+    
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator+(const Second& in) const
+	-> Operation<T,std::plus<T>,OperationClass,Second>
+    {
+    	return Operation<T,std::plus<T>,OperationClass,Second>(*mOc, in);
+    }
+
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator-(const Second& in) const
+	-> Operation<T,std::minus<T>,OperationClass,Second>
+    {
+    	return Operation<T,std::minus<T>,OperationClass,Second>(*mOc, in);
+    }
+    
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator*(const Second& in) const
+	-> Operation<T,std::multiplies<T>,OperationClass,Second>
+    {
+    	return Operation<T,std::multiplies<T>,OperationClass,Second>(*mOc, in);
+    }
+
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationTemplate<T,OperationClass>::operator/(const Second& in) const
+	-> Operation<T,std::divides<T>,OperationClass,Second>
+    {
+    	return Operation<T,std::divides<T>,OperationClass,Second>(*mOc, in);
+    }
+
+    template <typename T, class OperationClass>
+    template <class IndexType>
+    auto OperationTemplate<T,OperationClass>::c(std::shared_ptr<IndexType>& ind) const
+	-> Contraction<T,OperationClass,IndexType>
+    {
+	return Contraction<T,OperationClass,IndexType>(*mOc, ind);
+    }
+
+    
+    /*************************
+     *   OperationMaster     *
+     *************************/
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::
+    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+		    std::shared_ptr<typename CRange::IndexType>& index) :
+	mSecond(second), mArrayRef(ma), mIndex()
+    {
+	MultiRangeFactory<Ranges...> mrf( index->range() );
+	std::shared_ptr<MultiRange<Ranges...> > mr =
+	    std::dynamic_pointer_cast<MultiRange<Ranges...> >( mrf.create() );
+	mIndex = std::make_shared<IndexType>( mr->begin() );
+	(*mIndex) = *index;
+
+	auto blockIndex = seekBlockIndex( mIndex, second);
+
+	block(blockIndex);
+	second.block(blockIndex);
+	
+	for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndex) ){
+	    get() = mSecond.get();
+	}
+    }
+    
+    template <typename T, class OpClass, class... Ranges>
+    MBlock<T>& OperationMaster<T,OpClass,Ranges...>::get()
+    {
+	block();
+	return *mBlockPtr;
+    }
+    
+    template <typename T, class OpClass, class... Ranges>
+    const Block<T>& OperationMaster<T,OpClass,Ranges...>::get() const
+    {
+	block();
+	return *mBlockPtr;
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    std::vector<BTSS> OperationMaster<T,OpClass,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
+    {
+	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
+	mBlockPtr = makeBlock(mArrayRef.datav(), btv[0].second, blockIndex->max());
+	return btv;
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    const OperationMaster<T,OpClass,Ranges...>& OperationMaster<T,OpClass,Ranges...>::block() const
+    {
+	mBlockPtr->set( mIndex->pos() );
+	return *this;
+    }
+    
+    /****************************
+     *   ConstOperationRoot     *
+     ****************************/
+
+    template <typename T, class... Ranges>
+    ConstOperationRoot<T,Ranges...>::
+    ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
+		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	OperationTemplate<T,ConstOperationRoot<T,Ranges...> >(this),
+	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
+    {
+	(*mIndex)(indices...);
+    }
+
+    template <typename T, class... Ranges>
+    const Block<T>& ConstOperationRoot<T,Ranges...>::get() const
+    {
+	block();
+	return *mBlockPtr;
+    }
+
+    template <typename T, class... Ranges>
+    std::vector<BTSS> ConstOperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
+    {
+	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
+	mBlockPtr = makeBlock(mArrayRef.datav(), btv[0].second, blockIndex->max());
+	return btv;
+    }
+
+    template <typename T, class... Ranges>
+    const ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::block() const
+    {
+	mBlockPtr->set( (*mIndex)().pos() );
+	return *this;
+    }
+    
+    /***********************
+     *   OperationRoot     *
+     ***********************/
+    
+    template <typename T, class... Ranges>
+    OperationRoot<T,Ranges...>::
+    OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	OperationTemplate<T,OperationRoot<T,Ranges...> >(this),
+	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
+    {
+	(*mIndex)(indices...);
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
+    {
+	return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex);
+    }
+
+    template <typename T, class... Ranges>
+    const MBlock<T>& OperationRoot<T,Ranges...>::get() const
+    {
+	block();
+	return *mBlockPtr;
+    }
+    
+    template <typename T, class... Ranges>
+    MBlock<T>& OperationRoot<T,Ranges...>::get()
+    {
+	block();
+	return *mBlockPtr;
+    }
+
+    template <typename T, class... Ranges>
+    std::vector<BTSS> OperationRoot<T,Ranges...>::block(const std::shared_ptr<IndexBase> blockIndex) const
+    {
+	std::vector<BTSS> btv(1, getBlockType(mIndex, blockIndex, true) );
+	mBlockPtr = makeBlock(mArrayRef.datav(), btv[0].second, blockIndex->max());
+	return btv;
+    }
+
+    template <typename T, class... Ranges>
+    const OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::block() const
+    {
+	mBlockPtr->set( (*mIndex)().pos() );
+	return *this;
+    }
+    
+    /*******************
+     *   Operation     *
+     *******************/
+    
+    template <typename T, class OpFunction, class... Ops>
+    Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
+	OperationTemplate<T,Operation<T,OpFunction,Ops...> >(this),
+	mOps(ops...) {}
+
+    template <typename T, class OpFunction, class... Ops>
+    const BlockResult<T>& Operation<T,OpFunction,Ops...>::get() const
+    {
+	mRes = std::move( PackNum<sizeof...(Ops)-1>::template unpackArgs<T,OpFunction>(mOps) );
+	return mRes;
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    std::vector<BTSS> Operation<T,OpFunction,Ops...>::block(const std::shared_ptr<IndexBase> blockIndex) const
+    {
+	std::vector<BTSS> btv;
+	PackNum<sizeof...(Ops)-1>::makeBlockTypeVec(btv, mOps, blockIndex);
+	return btv;
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    const Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::block() const
+    {
+	//mBlockPtr->set( mIndex->pos() );
+	return *this;
+    }
+
+    /*********************
+     *   Contraction     *
+     *********************/
+
+    template <typename T, class Op, class IndexType>
+    Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
+	OperationTemplate<T,Contraction<T,Op,IndexType> >(this),
+	mOp(op),
+	mInd(ind) {}
+
+    template <typename T, class Op, class IndexType>
+    const BlockResult<T>& Contraction<T,Op,IndexType>::get() const
+    {
+	BlockBinaryOpSelf<T,std::plus<T>,decltype(mOp.get())> f(mRes);
+	for(*mInd = 0; mInd->pos() != mInd->max(); ++(*mInd)){
+	    f(mOp.get());
+	}
+	return mRes;
+    }
+    
+    template <typename T, class Op, class IndexType>
+    std::vector<BTSS> Contraction<T,Op,IndexType>::block(const std::shared_ptr<IndexBase> blockIndex) const
+    {
+	return mOp.block(blockIndex);
+    }
+
+    template <typename T, class Op, class IndexType>
+    const Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::block() const
+    {
+	return *this;
+    }
+    
+}
 
 #endif
diff --git a/src/multi_range.cc b/src/multi_range.cc
deleted file mode 100644
index 049ab48..0000000
--- a/src/multi_range.cc
+++ /dev/null
@@ -1,320 +0,0 @@
-
-#include "multi_range.h"
-#include "pack_num.h"
-
-namespace MultiArrayTools
-{
-
-    namespace
-    {
-	using namespace MultiArrayHelper;
-    }
-
-    /******************
-     *   MultiIndex   *
-     ******************/
-
-    /*
-    template <class... Indices>
-    MultiIndex<Indices...>::MultiIndex(const MultiIndex<Indices...>& in) :
-	IndexInterface<std::tuple<typename Indices::MetaType...> >(in)
-    {
-	PackNum<sizeof...(Indices)-1>::copy(mIPack, in);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-    }
-	
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(const MultiIndex<Indices...>& in)
-    {
-	IndexI::operator=(in);
-	PackNum<sizeof...(Indices)-1>::copy(mIPack, in);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	return *this;
-    }
-    */
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(ContainerIndex<Indices...>& ci)
-    {	
-	PackNum<sizeof...(Indices)-1>::copyInst(mIPack, ci);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	return *this;
-    }
-    
-    template <class... Indices>
-    template <class MRange>
-    MultiIndex<Indices...>::MultiIndex(const std::shared_ptr<MRange>& range) :
-	IndexInterface<std::tuple<typename Indices::MetaType...> >(range, 0)
-    {
-	PackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
-	PackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack); // has one more element!
-    }
-
-    template <class... Indices>
-    IndexType MultiIndex<Indices...>::type() const
-    {
-	return IndexType::MULTI;
-    }
-    
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::operator++()
-    {
-	PackNum<sizeof...(Indices)-1>::pp( mIPack );
-	++IB::mPos;
-	return *this;
-    }
-
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::operator--()
-    {
-	PackNum<sizeof...(Indices)-1>::mm( mIPack );
-	--IB::mPos;
-	return *this;
-    }
-    
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(size_t pos)
-    {
-	IB::mPos = pos;
-	PackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
-	return *this;
-    }
-
-    template <class... Indices>
-    int MultiIndex<Indices...>::pp(std::shared_ptr<IndexBase>& idxPtr)
-    {
-	int tmp = PackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtr);
-	IB::mPos += tmp;
-	return tmp;
-    }
-
-    template <class... Indices>
-    int MultiIndex<Indices...>::mm(std::shared_ptr<IndexBase>& idxPtr)
-    {
-	int tmp = PackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtr);
-	IB::mPos -= tmp;
-	return tmp;
-    }
-
-    template <class... Indices>
-    template <size_t DIR>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::up()
-    {
-	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
-	IB::mPos += PackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
-	PackNum<DIR>::pp( mIPack );
-	return *this;
-    }
-
-    template <class... Indices>
-    template <size_t DIR>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::down()
-    {
-	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
-	IB::mPos -= PackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
-	PackNum<DIR>::mm( mIPack );
-	return *this;
-    }
-    
-    template <class... Indices>
-    size_t MultiIndex<Indices...>::dim() const
-    {
-	return sizeof...(Indices);
-    }
-    
-    template <class... Indices>
-    template <size_t N>
-    auto MultiIndex<Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
-    {
-	return *std::get<N>(mIPack);
-    }
-
-    template <class... Indices>
-    template <size_t N>
-    auto MultiIndex<Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
-    {
-	return std::get<N>(mIPack);
-    }
-    
-    template <class... Indices>
-    const IndexBase& MultiIndex<Indices...>::get(size_t n) const
-    {
-	if(n >= sizeof...(Indices)){
-	    assert(0);
-	    // throw !!
-	}
-	MultiIndex<Indices...> const* t = this;
-	return PackNum<sizeof...(Indices)-1>::getIndex(*t, n);
-    }
-
-    template <class... Indices>
-    std::shared_ptr<IndexBase> MultiIndex<Indices...>::getPtr(size_t n) const
-    {
-	if(n >= sizeof...(Indices)){
-	    assert(0);
-	    // throw !!
-	}
-	MultiIndex<Indices...> const* t = this;
-	return PackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
-    }
-
-    template <class... Indices>
-    size_t MultiIndex<Indices...>::getStepSize(size_t n) const
-    {
-	if(n >= sizeof...(Indices)){
-	    assert(0);
-	    // throw !!
-	}
-	return mBlockSizes[n+1];
-    }
-    
-    template <class... Indices>
-    typename MultiIndex<Indices...>::MetaType MultiIndex<Indices...>::meta() const
-    {
-        MetaType metaTuple;
-	PackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
-	return metaTuple;
-    }
-
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::at(const MultiIndex<Indices...>::MetaType& metaPos)
-    {
-	PackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
-	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	return *this;
-    }
-
-    template <class... Indices>
-    bool MultiIndex<Indices...>::first() const
-    {
-	return IB::mPos == 0;
-    }
-
-    template <class... Indices>
-    bool MultiIndex<Indices...>::last() const
-    {
-	return IB::mPos == IB::mRangePtr->size() - 1;
-    }
-
-    template <class... Indices>
-    std::shared_ptr<typename MultiIndex<Indices...>::RangeType> MultiIndex<Indices...>::range() const
-    {
-	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
-    }
-    /*
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::lock(std::shared_ptr<IndexBase>& idx)
-    {
-	IB::mLocked = (idx.get() == this);
-	PackNum<sizeof...(Indices)-1>::lock(mIPack, idx);
-	return *this;
-	}*/
-    
-    template <class... Indices>
-    MultiIndex<Indices...>& MultiIndex<Indices...>::operator()(std::shared_ptr<Indices>&... indices)
-    {
-	PackNum<sizeof...(Indices)-1>::swapIndices(mIPack, indices...);
-	PackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, IB::mPos);
-	return *this;
-    }
-    
-    /*************************
-     *   MultiRangeFactory   *
-     *************************/
-
-    template <class... Ranges>
-    MultiRangeFactory<Ranges...>::MultiRangeFactory(const std::shared_ptr<Ranges>&... rs)
-    {
-	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( rs... ) );
-    }
-    
-    template <class... Ranges>
-    MultiRangeFactory<Ranges...>::MultiRangeFactory(const typename MultiRange<Ranges...>::SpaceType& st)
-    {
-	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( st ) );
-    }
-
-    template <class... Ranges>
-    MultiRangeFactory<Ranges...>::MultiRangeFactory(const std::shared_ptr<ContainerRange<Ranges...> >& cr)
-    {
-	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( cr->space() ) );
-    }
-    
-    template <class... Ranges>
-    std::shared_ptr<RangeBase> MultiRangeFactory<Ranges...>::create()
-    {
-	setSelf();
-	return mProd;
-    }
-    
-    /******************
-     *   MultiRange   *
-     ******************/
-
-    template <class... Ranges>
-    MultiRange<Ranges...>::MultiRange(const std::shared_ptr<Ranges>&... rs) : mSpace(std::make_tuple(rs...)) {}
-
-    template <class... Ranges>
-    MultiRange<Ranges...>::MultiRange(const SpaceType& space) : mSpace( space ) {}
-
-    template <class... Ranges>
-    template <size_t N>
-    auto MultiRange<Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
-    {
-	return *std::get<N>(mSpace);
-    }
-
-    template <class... Ranges>
-    template <size_t N>
-    auto MultiRange<Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
-    {
-	return std::get<N>(mSpace);
-    }
-
-    template <class... Ranges>
-    size_t MultiRange<Ranges...>::dim() const
-    {
-	return sdim;
-    }
-    
-    template <class... Ranges>
-    size_t MultiRange<Ranges...>::size() const
-    {
-	return PackNum<sizeof...(Ranges)-1>::getSize(mSpace);
-    }
-
-    template <class... Ranges>
-    const typename MultiRange<Ranges...>::SpaceType& MultiRange<Ranges...>::space() const
-    {
-	return mSpace;
-    }
-    
-    template <class... Ranges>
-    typename MultiRange<Ranges...>::IndexType MultiRange<Ranges...>::begin() const
-    {
-	MultiIndex<typename Ranges::IndexType...>
-	    i( std::dynamic_pointer_cast<MultiRange<Ranges...> >
-	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-	i = 0;
-	return i;
-    }
-
-    template <class... Ranges>
-    typename MultiRange<Ranges...>::IndexType MultiRange<Ranges...>::end() const
-    {
-	MultiIndex<typename Ranges::IndexType...>
-	    i( std::dynamic_pointer_cast<MultiRange<Ranges...> >
-	       ( std::shared_ptr<RangeBase>( RB::mThis )) );
-	i = size();
-	return i;
-    }
-
-    template <class... Ranges>
-    std::shared_ptr<IndexBase> MultiRange<Ranges...>::index() const
-    {
-	return std::make_shared<MultiIndex<typename Ranges::IndexType...> >
-	    ( std::dynamic_pointer_cast<MultiRange<Ranges...> >
-	      ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-    }
-}
diff --git a/src/multi_range.h b/src/multi_range.h
index bfd1285..ecfe670 100644
--- a/src/multi_range.h
+++ b/src/multi_range.h
@@ -11,6 +11,8 @@
 #include "range_base.h"
 #include "index_base.h"
 
+#include "pack_num.h"
+
 namespace MultiArrayTools
 {
     
@@ -152,6 +154,325 @@ namespace MultiArrayTools
     
 }
 
-#include "multi_range.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+
+    /******************
+     *   MultiIndex   *
+     ******************/
+
+    /*
+    template <class... Indices>
+    MultiIndex<Indices...>::MultiIndex(const MultiIndex<Indices...>& in) :
+	IndexInterface<std::tuple<typename Indices::MetaType...> >(in)
+    {
+	PackNum<sizeof...(Indices)-1>::copy(mIPack, in);
+	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+    }
+	
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(const MultiIndex<Indices...>& in)
+    {
+	IndexI::operator=(in);
+	PackNum<sizeof...(Indices)-1>::copy(mIPack, in);
+	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+    */
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(ContainerIndex<Indices...>& ci)
+    {	
+	PackNum<sizeof...(Indices)-1>::copyInst(mIPack, ci);
+	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+    
+    template <class... Indices>
+    template <class MRange>
+    MultiIndex<Indices...>::MultiIndex(const std::shared_ptr<MRange>& range) :
+	IndexInterface<std::tuple<typename Indices::MetaType...> >(range, 0)
+    {
+	PackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
+	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
+	PackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack); // has one more element!
+    }
+
+    template <class... Indices>
+    IndexType MultiIndex<Indices...>::type() const
+    {
+	return IndexType::MULTI;
+    }
+    
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator++()
+    {
+	PackNum<sizeof...(Indices)-1>::pp( mIPack );
+	++IB::mPos;
+	return *this;
+    }
+
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator--()
+    {
+	PackNum<sizeof...(Indices)-1>::mm( mIPack );
+	--IB::mPos;
+	return *this;
+    }
+    
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	PackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
+	return *this;
+    }
+
+    template <class... Indices>
+    int MultiIndex<Indices...>::pp(std::shared_ptr<IndexBase>& idxPtr)
+    {
+	int tmp = PackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtr);
+	IB::mPos += tmp;
+	return tmp;
+    }
+
+    template <class... Indices>
+    int MultiIndex<Indices...>::mm(std::shared_ptr<IndexBase>& idxPtr)
+    {
+	int tmp = PackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtr);
+	IB::mPos -= tmp;
+	return tmp;
+    }
+
+    template <class... Indices>
+    template <size_t DIR>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::up()
+    {
+	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
+	IB::mPos += PackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
+	PackNum<DIR>::pp( mIPack );
+	return *this;
+    }
+
+    template <class... Indices>
+    template <size_t DIR>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::down()
+    {
+	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
+	IB::mPos -= PackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
+	PackNum<DIR>::mm( mIPack );
+	return *this;
+    }
+    
+    template <class... Indices>
+    size_t MultiIndex<Indices...>::dim() const
+    {
+	return sizeof...(Indices);
+    }
+    
+    template <class... Indices>
+    template <size_t N>
+    auto MultiIndex<Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
+    {
+	return *std::get<N>(mIPack);
+    }
+
+    template <class... Indices>
+    template <size_t N>
+    auto MultiIndex<Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>(mIPack);
+    }
+    
+    template <class... Indices>
+    const IndexBase& MultiIndex<Indices...>::get(size_t n) const
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	MultiIndex<Indices...> const* t = this;
+	return PackNum<sizeof...(Indices)-1>::getIndex(*t, n);
+    }
+
+    template <class... Indices>
+    std::shared_ptr<IndexBase> MultiIndex<Indices...>::getPtr(size_t n) const
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	MultiIndex<Indices...> const* t = this;
+	return PackNum<sizeof...(Indices)-1>::getIndexPtr(*t, n);
+    }
+
+    template <class... Indices>
+    size_t MultiIndex<Indices...>::getStepSize(size_t n) const
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	return mBlockSizes[n+1];
+    }
+    
+    template <class... Indices>
+    typename MultiIndex<Indices...>::MetaType MultiIndex<Indices...>::meta() const
+    {
+        MetaType metaTuple;
+	PackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
+	return metaTuple;
+    }
+
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::at(const MultiIndex<Indices...>::MetaType& metaPos)
+    {
+	PackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
+	IB::mPos = PackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+
+    template <class... Indices>
+    bool MultiIndex<Indices...>::first() const
+    {
+	return IB::mPos == 0;
+    }
+
+    template <class... Indices>
+    bool MultiIndex<Indices...>::last() const
+    {
+	return IB::mPos == IB::mRangePtr->size() - 1;
+    }
+
+    template <class... Indices>
+    std::shared_ptr<typename MultiIndex<Indices...>::RangeType> MultiIndex<Indices...>::range() const
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
+    }
+    /*
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::lock(std::shared_ptr<IndexBase>& idx)
+    {
+	IB::mLocked = (idx.get() == this);
+	PackNum<sizeof...(Indices)-1>::lock(mIPack, idx);
+	return *this;
+	}*/
+    
+    template <class... Indices>
+    MultiIndex<Indices...>& MultiIndex<Indices...>::operator()(std::shared_ptr<Indices>&... indices)
+    {
+	PackNum<sizeof...(Indices)-1>::swapIndices(mIPack, indices...);
+	PackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, IB::mPos);
+	return *this;
+    }
+    
+    /*************************
+     *   MultiRangeFactory   *
+     *************************/
+
+    template <class... Ranges>
+    MultiRangeFactory<Ranges...>::MultiRangeFactory(const std::shared_ptr<Ranges>&... rs)
+    {
+	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( rs... ) );
+    }
+    
+    template <class... Ranges>
+    MultiRangeFactory<Ranges...>::MultiRangeFactory(const typename MultiRange<Ranges...>::SpaceType& st)
+    {
+	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( st ) );
+    }
+
+    template <class... Ranges>
+    MultiRangeFactory<Ranges...>::MultiRangeFactory(const std::shared_ptr<ContainerRange<Ranges...> >& cr)
+    {
+	mProd = std::shared_ptr< MultiRange<Ranges...> >( new MultiRange<Ranges...>( cr->space() ) );
+    }
+    
+    template <class... Ranges>
+    std::shared_ptr<RangeBase> MultiRangeFactory<Ranges...>::create()
+    {
+	setSelf();
+	return mProd;
+    }
+    
+    /******************
+     *   MultiRange   *
+     ******************/
+
+    template <class... Ranges>
+    MultiRange<Ranges...>::MultiRange(const std::shared_ptr<Ranges>&... rs) : mSpace(std::make_tuple(rs...)) {}
+
+    template <class... Ranges>
+    MultiRange<Ranges...>::MultiRange(const SpaceType& space) : mSpace( space ) {}
+
+    template <class... Ranges>
+    template <size_t N>
+    auto MultiRange<Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
+    {
+	return *std::get<N>(mSpace);
+    }
+
+    template <class... Ranges>
+    template <size_t N>
+    auto MultiRange<Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
+    {
+	return std::get<N>(mSpace);
+    }
+
+    template <class... Ranges>
+    size_t MultiRange<Ranges...>::dim() const
+    {
+	return sdim;
+    }
+    
+    template <class... Ranges>
+    size_t MultiRange<Ranges...>::size() const
+    {
+	return PackNum<sizeof...(Ranges)-1>::getSize(mSpace);
+    }
+
+    template <class... Ranges>
+    const typename MultiRange<Ranges...>::SpaceType& MultiRange<Ranges...>::space() const
+    {
+	return mSpace;
+    }
+    
+    template <class... Ranges>
+    typename MultiRange<Ranges...>::IndexType MultiRange<Ranges...>::begin() const
+    {
+	MultiIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<MultiRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+
+    template <class... Ranges>
+    typename MultiRange<Ranges...>::IndexType MultiRange<Ranges...>::end() const
+    {
+	MultiIndex<typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<MultiRange<Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis )) );
+	i = size();
+	return i;
+    }
+
+    template <class... Ranges>
+    std::shared_ptr<IndexBase> MultiRange<Ranges...>::index() const
+    {
+	return std::make_shared<MultiIndex<typename Ranges::IndexType...> >
+	    ( std::dynamic_pointer_cast<MultiRange<Ranges...> >
+	      ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+    }
+}
 
 #endif
diff --git a/src/range_base.cc b/src/range_base.cc
deleted file mode 100644
index 133078c..0000000
--- a/src/range_base.cc
+++ /dev/null
@@ -1,30 +0,0 @@
-
-#include "range_base.h"
-
-namespace MultiArrayTools
-{
-
-    /*************************
-     *   RangeFactoryBase    *
-     *************************/
-    
-    void RangeFactoryBase::setSelf()
-    {
-	mProd->mThis = mProd;
-    }
-    
-    /******************
-     *   RangeBase    *
-     ******************/
-
-    bool RangeBase::operator==(const RangeBase& in) const
-    {
-	return this == &in;
-    }
-    
-    bool RangeBase::operator!=(const RangeBase& in) const
-    {
-	return this != &in;
-    }
-
-}
diff --git a/src/range_base.h b/src/range_base.h
index e06f000..c977c99 100644
--- a/src/range_base.h
+++ b/src/range_base.h
@@ -81,6 +81,36 @@ namespace MultiArrayTools
     
 }
 
-#include "range_base.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+
+    /*************************
+     *   RangeFactoryBase    *
+     *************************/
+    
+    void RangeFactoryBase::setSelf()
+    {
+	mProd->mThis = mProd;
+    }
+    
+    /******************
+     *   RangeBase    *
+     ******************/
+
+    bool RangeBase::operator==(const RangeBase& in) const
+    {
+	return this == &in;
+    }
+    
+    bool RangeBase::operator!=(const RangeBase& in) const
+    {
+	return this != &in;
+    }
+
+}
 
 #endif
diff --git a/src/range_factory.cc b/src/range_factory.cc
deleted file mode 100644
index 4ab5e40..0000000
--- a/src/range_factory.cc
+++ /dev/null
@@ -1,8 +0,0 @@
-// -*- C++ -*-
-
-#include "range_factory.h"
-
-namespace MultiArrayTools
-{
-    
-} //namespace MultiArrayTools
diff --git a/src/range_factory.h b/src/range_factory.h
index daca600..99a4714 100644
--- a/src/range_factory.h
+++ b/src/range_factory.h
@@ -21,6 +21,6 @@ namespace MultiArrayTools
     
 } //namespace MultiArrayTools
 
-#include "range_factory.cc"
+// === NO TEMPLATE CODE HERE ===
 
 #endif
diff --git a/src/single_range.cc b/src/single_range.cc
deleted file mode 100644
index 7e0e98d..0000000
--- a/src/single_range.cc
+++ /dev/null
@@ -1,379 +0,0 @@
-
-#include "single_range.h"
-
-namespace MultiArrayTools
-{
-    /******************
-     *  SingleIndex   *	     
-     ******************/
-
-    template <typename U, RangeType TYPE>
-    SingleIndex<U,TYPE>::SingleIndex(const std::shared_ptr<SingleRange<U,TYPE> >& range) :
-	IndexInterface<U>(range, 0) {}
-
-    template <typename U, RangeType TYPE>
-    IndexType SingleIndex<U,TYPE>::type() const
-    {
-	return IndexType::SINGLE;
-    }
-	
-    template <typename U, RangeType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator=(size_t pos)
-    {
-	IB::mPos = pos;
-	return *this;
-    }
-    
-    template <typename U, RangeType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator++()
-    {
-	++IB::mPos;
-	return *this;
-    }
-
-    template <typename U, RangeType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator--()
-    {
-	--IB::mPos;
-	return *this;
-    }
-
-    template <typename U, RangeType TYPE>
-    int SingleIndex<U,TYPE>::pp(std::shared_ptr<IndexBase>& idxPtr)
-    {
-	++(*this);
-	return 1;
-    }
-
-    template <typename U, RangeType TYPE>
-    int SingleIndex<U,TYPE>::mm(std::shared_ptr<IndexBase>& idxPtr)
-    {
-	--(*this);
-	return 1;
-    }
-    
-    template <typename U, RangeType TYPE>
-    U SingleIndex<U,TYPE>::meta() const
-    {
-	return std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->get( IB::pos() );
-    }
-
-    template <typename U, RangeType TYPE>
-    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::at(const U& metaPos)
-    {
-	operator=( std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->getMeta( metaPos ) );
-	return *this;
-    }
-
-    template <typename U, RangeType TYPE>
-    size_t SingleIndex<U,TYPE>::dim() const
-    {
-	return 1;
-    }
-
-    template <typename U, RangeType TYPE>
-    bool SingleIndex<U,TYPE>::last() const
-    {
-	return IB::mPos == IB::mRangePtr->size() - 1;
-    }
-
-    template <typename U, RangeType TYPE>
-    bool SingleIndex<U,TYPE>::first() const
-    {
-	return IB::mPos == 0;
-    }
-
-    template <typename U, RangeType TYPE>
-    std::shared_ptr<IndexBase> SingleIndex<U,TYPE>::getPtr(size_t n) const
-    {
-	return std::shared_ptr<IndexBase>();
-    }
-
-    template <typename U, RangeType TYPE>    
-    size_t SingleIndex<U,TYPE>::getStepSize(size_t n) const
-    {
-	return 1;
-    }
-    
-    /********************
-     *   SingleRange    *
-     ********************/
-
-    template <typename U, RangeType TYPE>
-    SingleRangeFactory<U,TYPE>::SingleRangeFactory(const std::vector<U>& space)
-    {
-	mProd = std::shared_ptr<oType>( new SingleRange<U,TYPE>( space ) );
-    }
-
-    template <typename U, RangeType TYPE>
-    std::shared_ptr<RangeBase> SingleRangeFactory<U,TYPE>::create()
-    {
-	setSelf();
-	return mProd;
-    }
-    
-    /********************
-     *   SingleRange    *
-     ********************/
-    
-    template <typename U, RangeType TYPE>
-    SingleRange<U,TYPE>::SingleRange(const std::vector<U>& space) : RangeInterface<SingleIndex<U,TYPE> >(),
-	mSpace(space) {}
-    
-    template <typename U, RangeType TYPE>
-    const U& SingleRange<U,TYPE>::get(size_t pos) const
-    {
-	return mSpace[pos];
-    }
-
-    template <typename U, RangeType TYPE>
-    size_t SingleRange<U,TYPE>::getMeta(const U& metaPos) const
-    {
-	size_t cnt = 0;
-	for(auto& x: mSpace){
-	    if(x == metaPos){
-		return cnt;
-	    }
-	    ++cnt;
-	}
-	return cnt;
-    }
-
-    template <typename U, RangeType TYPE>
-    size_t SingleRange<U,TYPE>::size() const
-    {
-	return mSpace.size();
-    }
-
-    template <typename U, RangeType TYPE>
-    size_t SingleRange<U,TYPE>::dim() const
-    {
-	return 1;
-    }
-
-    template <typename U, RangeType TYPE>
-    typename SingleRange<U,TYPE>::IndexType SingleRange<U,TYPE>::begin() const
-    {
-	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
-			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-	i = 0;
-	return i;
-    }
-    
-    template <typename U, RangeType TYPE>
-    typename SingleRange<U,TYPE>::IndexType SingleRange<U,TYPE>::end() const
-    {
-	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
-			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-	i = size();
-	return i;
-    }
-
-    // put this in the interface class !!!
-    template <typename U, RangeType TYPE>
-    std::shared_ptr<IndexBase> SingleRange<U,TYPE>::index() const
-    {
-	return std::make_shared<SingleIndex<U,TYPE> >
-	    ( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
-	      ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-    }
-
-    /*
-    // specializations (not updated!!!)
-
-    SingleRange<int,RangeType::SPACE>::SingleRange(size_t ext) :
-	RangeBase<SingleIndex<int,RangeType::SPACE> >(),
-	mExt(ext) {}
-    
-    int SingleRange<int,RangeType::SPACE>::get(size_t pos) const
-    {
-	return (pos <= mExt / 2) ? static_cast<int>( pos ) : static_cast<int>( pos ) - static_cast<int>( mExt );
-    }
-
-    size_t SingleRange<int,RangeType::SPACE>::getMeta(int metaPos) const
-    {
-	return (metaPos < 0) ? metaPos + mExt : metaPos;
-    }
-
-    
-    size_t SingleRange<int,RangeType::SPACE>::size() const
-    {
-	return mExt;
-    }
-
-    
-    MultiRangeType SingleRange<int,RangeType::SPACE>::type() const
-    {
-	return MultiRangeType(RangeType::SPACE);
-    }
-    
-    
-    SingleIndex<int,RangeType::SPACE> SingleRange<int,RangeType::SPACE>::begin() const
-    {
-	return SingleIndex<int,RangeType::SPACE>(this, 0);
-    }
-
-    
-    SingleIndex<int,RangeType::SPACE> SingleRange<int,RangeType::SPACE>::end() const
-    {
-	return SingleIndex<int,RangeType::SPACE>(this, size());
-    }
-
-    // 
-    
-    SingleRange<size_t,RangeType::DISTANCE>::SingleRange(size_t ext) :
-	RangeBase<SingleIndex<size_t,RangeType::DISTANCE> >(),
-	mExt(ext) {}
-    
-    size_t SingleRange<size_t,RangeType::DISTANCE>::get(size_t pos) const
-    {
-	return pos;
-    }
-
-    size_t SingleRange<size_t,RangeType::DISTANCE>::getMeta(size_t metaPos) const
-    {
-	return metaPos;
-    }
-    
-    size_t SingleRange<size_t,RangeType::DISTANCE>::size() const
-    {
-	return mExt;
-    }
-    
-    MultiRangeType SingleRange<size_t,RangeType::DISTANCE>::type() const
-    {
-	return MultiRangeType(RangeType::DISTANCE);
-    }
-   
-    SingleIndex<size_t,RangeType::DISTANCE> SingleRange<size_t,RangeType::DISTANCE>::begin() const
-    {
-	return SingleIndex<size_t,RangeType::DISTANCE>(this, 0);
-    }
-    
-    SingleIndex<size_t,RangeType::DISTANCE> SingleRange<size_t,RangeType::DISTANCE>::end() const
-    {
-	return SingleIndex<size_t,RangeType::DISTANCE>(this, size());
-    }
-
-    // 
-    
-    SingleRange<size_t,RangeType::ENSEMBLE>::SingleRange(size_t num) :
-	RangeBase<SingleIndex<size_t,RangeType::ENSEMBLE> >(),
-	mNum(num) {}
-    
-    
-    size_t SingleRange<size_t,RangeType::ENSEMBLE>::get(size_t pos) const
-    {
-	return pos;
-    }
-
-    
-    size_t SingleRange<size_t,RangeType::ENSEMBLE>::getMeta(size_t metaPos) const
-    {
-	return metaPos;
-    }
-
-    
-    size_t SingleRange<size_t,RangeType::ENSEMBLE>::size() const
-    {
-	return mNum;
-    }
-
-    
-    MultiRangeType SingleRange<size_t,RangeType::ENSEMBLE>::type() const
-    {
-	return MultiRangeType(RangeType::ENSEMBLE);
-    }
-    
-    
-    SingleIndex<size_t,RangeType::ENSEMBLE> SingleRange<size_t,RangeType::ENSEMBLE>::begin() const
-    {
-	return SingleIndex<size_t,RangeType::ENSEMBLE>(this, 0);
-    }
-
-    
-    SingleIndex<size_t,RangeType::ENSEMBLE> SingleRange<size_t,RangeType::ENSEMBLE>::end() const
-    {
-	return SingleIndex<size_t,RangeType::ENSEMBLE>(this, size());
-    }
-
-    //
-    
-    VET SingleRange<VET,RangeType::VALUE_ERROR>::get(size_t pos) const
-    {
-	return static_cast<VET>( pos );
-    }
-
-    
-    size_t SingleRange<VET,RangeType::VALUE_ERROR>::getMeta(VET metaPos) const
-    {
-	return static_cast<size_t>( metaPos );
-    }
-
-    
-    size_t SingleRange<VET,RangeType::VALUE_ERROR>::size() const
-    {
-	return 2;
-    }
-
-    
-    MultiRangeType SingleRange<VET,RangeType::VALUE_ERROR>::type() const
-    {
-	return MultiRangeType(RangeType::VALUE_ERROR);
-    }
-    
-    
-    SingleIndex<VET,RangeType::VALUE_ERROR> SingleRange<VET,RangeType::VALUE_ERROR>::begin() const
-    {
-	return SingleIndex<VET,RangeType::VALUE_ERROR>(this, 0);
-    }
-
-    
-    SingleIndex<VET,RangeType::VALUE_ERROR> SingleRange<VET,RangeType::VALUE_ERROR>::end() const
-    {
-	return SingleIndex<VET,RangeType::VALUE_ERROR>(this, size());
-    }
-
-    // 
-
-    
-    size_t SingleRange<size_t,RangeType::LORENTZ>::get(size_t pos) const
-    {
-	return pos;
-    }
-
-    
-    size_t SingleRange<size_t,RangeType::LORENTZ>::getMeta(size_t metaPos) const
-    {
-	return metaPos;
-    }
-
-    
-    size_t SingleRange<size_t,RangeType::LORENTZ>::size() const
-    {
-	#ifdef LORENTZ_DIMENSION
-	return LORENTZ_DIMENSION;
-	#else
-	return 4; // 4
-	#endif
-    }
-
-    
-    MultiRangeType SingleRange<size_t,RangeType::LORENTZ>::type() const
-    {
-	return MultiRangeType(RangeType::LORENTZ);
-    }
-    
-    
-    SingleIndex<size_t,RangeType::LORENTZ> SingleRange<size_t,RangeType::LORENTZ>::begin() const
-    {
-	return SingleIndex<size_t,RangeType::LORENTZ>(this, 0);
-    }
-
-    
-    SingleIndex<size_t,RangeType::LORENTZ> SingleRange<size_t,RangeType::LORENTZ>::end() const
-    {
-	return SingleIndex<size_t,RangeType::LORENTZ>(this, size());
-    }
-    */    
-}
diff --git a/src/single_range.h b/src/single_range.h
index 8cc4a23..100a60a 100644
--- a/src/single_range.h
+++ b/src/single_range.h
@@ -90,133 +90,190 @@ namespace MultiArrayTools
 
 	std::vector<U> mSpace;
     };
-    /*
-    // specializaions
 
-    template <>
-    class SingleRange<int,RangeType::SPACE> : public RangeBase<SingleIndex<int,RangeType::SPACE> >
-    {
-    public:
-	typedef typename RangeBase<SingleIndex<int,RangeType::SPACE> >::IndexType IndexType;
-
-	static SingleRange<int,RangeType::SPACE> oType() { return SingleRange<int,RangeType::SPACE>(); }
-	
-	virtual size_t size() const override;
-	
-	int get(size_t pos) const;
-	size_t getMeta(int metaPos) const;
-
-	virtual MultiRangeType type() const override;
-	
-	SingleIndex<int,RangeType::SPACE> begin() const override;
-	SingleIndex<int,RangeType::SPACE> end() const override;
-	
-    protected:
-	
-	SingleRange(size_t ext);
-
-	size_t mExt;
-    };
-
-    template <>
-    class SingleRange<size_t,RangeType::DISTANCE> : public RangeBase<SingleIndex<size_t,RangeType::DISTANCE> >
-    {
-    public:
-	typedef typename RangeBase<SingleIndex<size_t,RangeType::DISTANCE> >::IndexType IndexType;
-
-	static SingleRange<size_t,RangeType::DISTANCE> oType() { return SingleRange<size_t,RangeType::DISTANCE>(); }
-	
-	virtual size_t size() const override;
-	
-	size_t get(size_t pos) const;
-	size_t getMeta(size_t metaPos) const;
-
-	virtual MultiRangeType type() const override;
-	
-	SingleIndex<size_t,RangeType::DISTANCE> begin() const override;
-	SingleIndex<size_t,RangeType::DISTANCE> end() const override;
-	
-    protected:
-	
-	SingleRange(size_t ext);
-
-	size_t mExt;
-    };
-
-    template <>
-    class SingleRange<size_t,RangeType::ENSEMBLE> : public RangeBase<SingleIndex<size_t,RangeType::ENSEMBLE> >
-    {
-    public:
-	typedef typename RangeBase<SingleIndex<size_t,RangeType::ENSEMBLE> >::IndexType IndexType;
-
-	static SingleRange<size_t,RangeType::ENSEMBLE> oType() { return SingleRange<size_t,RangeType::ENSEMBLE>(); }
-	
-	virtual size_t size() const override;
-	
-	size_t get(size_t pos) const;
-	size_t getMeta(size_t metaPos) const;
-
-	virtual MultiRangeType type() const override;
-	
-	SingleIndex<size_t,RangeType::ENSEMBLE> begin() const override;
-	SingleIndex<size_t,RangeType::ENSEMBLE> end() const override;
-	
-    protected:
-	
-	SingleRange(size_t num);
-
-	size_t mNum;
-    };
-
-    enum class VET
-    {
-	VALUE = 0,
-	ERROR = 1
-    };
-
-    std::ostream& operator<<(std::ostream& os, VET vet);
-    
-    template <>
-    class SingleRange<VET,RangeType::VALUE_ERROR> : public RangeBase<SingleIndex<VET,RangeType::VALUE_ERROR> >
-    {
-    public:
-	typedef typename RangeBase<SingleIndex<VET,RangeType::VALUE_ERROR> >::IndexType IndexType;
-
-	static SingleRange<VET,RangeType::VALUE_ERROR> oType()
-	{ return SingleRange<VET,RangeType::VALUE_ERROR>(); }
-	
-	virtual size_t size() const override;
-	
-	VET get(size_t pos) const;
-	size_t getMeta(VET metaPos) const;
-
-	virtual MultiRangeType type() const override;
-	
-	SingleIndex<VET,RangeType::VALUE_ERROR> begin() const override;
-	SingleIndex<VET,RangeType::VALUE_ERROR> end() const override;
-    };
-
-    template <>
-    class SingleRange<size_t,RangeType::LORENTZ> : public RangeBase<SingleIndex<size_t,RangeType::LORENTZ> >
-    {
-    public:
-	typedef typename RangeBase<SingleIndex<size_t,RangeType::LORENTZ> >::IndexType IndexType;
-
-	static SingleRange<size_t,RangeType::LORENTZ> oType() { return SingleRange<size_t,RangeType::LORENTZ>(); }
-	
-	virtual size_t size() const override;
-	
-	size_t get(size_t pos) const;
-	size_t getMeta(size_t metaPos) const;
-
-	virtual MultiRangeType type() const override;
-	
-	SingleIndex<size_t,RangeType::LORENTZ> begin() const override;
-	SingleIndex<size_t,RangeType::LORENTZ> end() const override;
-    };
-    */
 }
 
-#include "single_range.cc"
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    /******************
+     *  SingleIndex   *	     
+     ******************/
+
+    template <typename U, RangeType TYPE>
+    SingleIndex<U,TYPE>::SingleIndex(const std::shared_ptr<SingleRange<U,TYPE> >& range) :
+	IndexInterface<U>(range, 0) {}
+
+    template <typename U, RangeType TYPE>
+    IndexType SingleIndex<U,TYPE>::type() const
+    {
+	return IndexType::SINGLE;
+    }
+	
+    template <typename U, RangeType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	return *this;
+    }
+    
+    template <typename U, RangeType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator++()
+    {
+	++IB::mPos;
+	return *this;
+    }
+
+    template <typename U, RangeType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::operator--()
+    {
+	--IB::mPos;
+	return *this;
+    }
+
+    template <typename U, RangeType TYPE>
+    int SingleIndex<U,TYPE>::pp(std::shared_ptr<IndexBase>& idxPtr)
+    {
+	++(*this);
+	return 1;
+    }
+
+    template <typename U, RangeType TYPE>
+    int SingleIndex<U,TYPE>::mm(std::shared_ptr<IndexBase>& idxPtr)
+    {
+	--(*this);
+	return 1;
+    }
+    
+    template <typename U, RangeType TYPE>
+    U SingleIndex<U,TYPE>::meta() const
+    {
+	return std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->get( IB::pos() );
+    }
+
+    template <typename U, RangeType TYPE>
+    SingleIndex<U,TYPE>& SingleIndex<U,TYPE>::at(const U& metaPos)
+    {
+	operator=( std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->getMeta( metaPos ) );
+	return *this;
+    }
+
+    template <typename U, RangeType TYPE>
+    size_t SingleIndex<U,TYPE>::dim() const
+    {
+	return 1;
+    }
+
+    template <typename U, RangeType TYPE>
+    bool SingleIndex<U,TYPE>::last() const
+    {
+	return IB::mPos == IB::mRangePtr->size() - 1;
+    }
+
+    template <typename U, RangeType TYPE>
+    bool SingleIndex<U,TYPE>::first() const
+    {
+	return IB::mPos == 0;
+    }
+
+    template <typename U, RangeType TYPE>
+    std::shared_ptr<IndexBase> SingleIndex<U,TYPE>::getPtr(size_t n) const
+    {
+	return std::shared_ptr<IndexBase>();
+    }
+
+    template <typename U, RangeType TYPE>    
+    size_t SingleIndex<U,TYPE>::getStepSize(size_t n) const
+    {
+	return 1;
+    }
+    
+    /********************
+     *   SingleRange    *
+     ********************/
+
+    template <typename U, RangeType TYPE>
+    SingleRangeFactory<U,TYPE>::SingleRangeFactory(const std::vector<U>& space)
+    {
+	mProd = std::shared_ptr<oType>( new SingleRange<U,TYPE>( space ) );
+    }
+
+    template <typename U, RangeType TYPE>
+    std::shared_ptr<RangeBase> SingleRangeFactory<U,TYPE>::create()
+    {
+	setSelf();
+	return mProd;
+    }
+    
+    /********************
+     *   SingleRange    *
+     ********************/
+    
+    template <typename U, RangeType TYPE>
+    SingleRange<U,TYPE>::SingleRange(const std::vector<U>& space) : RangeInterface<SingleIndex<U,TYPE> >(),
+	mSpace(space) {}
+    
+    template <typename U, RangeType TYPE>
+    const U& SingleRange<U,TYPE>::get(size_t pos) const
+    {
+	return mSpace[pos];
+    }
+
+    template <typename U, RangeType TYPE>
+    size_t SingleRange<U,TYPE>::getMeta(const U& metaPos) const
+    {
+	size_t cnt = 0;
+	for(auto& x: mSpace){
+	    if(x == metaPos){
+		return cnt;
+	    }
+	    ++cnt;
+	}
+	return cnt;
+    }
+
+    template <typename U, RangeType TYPE>
+    size_t SingleRange<U,TYPE>::size() const
+    {
+	return mSpace.size();
+    }
+
+    template <typename U, RangeType TYPE>
+    size_t SingleRange<U,TYPE>::dim() const
+    {
+	return 1;
+    }
+
+    template <typename U, RangeType TYPE>
+    typename SingleRange<U,TYPE>::IndexType SingleRange<U,TYPE>::begin() const
+    {
+	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+    
+    template <typename U, RangeType TYPE>
+    typename SingleRange<U,TYPE>::IndexType SingleRange<U,TYPE>::end() const
+    {
+	SingleIndex<U,TYPE> i( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+			       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = size();
+	return i;
+    }
+
+    // put this in the interface class !!!
+    template <typename U, RangeType TYPE>
+    std::shared_ptr<IndexBase> SingleRange<U,TYPE>::index() const
+    {
+	return std::make_shared<SingleIndex<U,TYPE> >
+	    ( std::dynamic_pointer_cast<SingleRange<U,TYPE> >
+	      ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+    }
+
+}
 
 #endif