diff --git a/src/include/expressions.h b/src/include/expressions.h
new file mode 100644
index 0000000..4911019
--- /dev/null
+++ b/src/include/expressions.h
@@ -0,0 +1,141 @@
+
+#ifndef __expressions_h__
+#define __expressions_h__
+
+#include "ranges/dynamic_range.h"
+#include "multi_array_base.h"
+#include "multi_array_operation.h"
+#include "xfor/xfor.h"
+
+#define Y() ,
+
+namespace MultiArrayTools
+{
+
+    //class Expressions1;
+
+    template <class EC>
+    using DDMA = MultiArrayBase<double,DynamicRange<EC>>;
+
+    template <class EC>
+    using DDMMA = MutableMultiArrayBase<double,DynamicRange<EC>>;
+
+    template <class EC, class MA>
+    using oo = decltype(std::declval<MA>()(std::declval<std::shared_ptr<DynamicIndex<EC>>>()));
+
+    template <template <class> class OpF, class... oos>
+    using OO = Operation<double,OpF<double>,oos...>;
+
+    template <class EC, template <class> class OpF, class... MAs>
+    using OX = Operation<double,OpF<double>,oo<EC,MAs>...>;
+    /*
+    template <class EC>
+    class ECInterface
+    {
+    public:
+
+        EC& THIS() { return static_cast<EC&>(*this); }
+        const EC& THIS() const { return static_cast<EC const&>(*this); }
+        
+        template <class Expr>
+        auto ifor(size_t step, Expr ex) const
+            -> decltype(THIS().ifor(step, ex))
+        {
+            return THIS().ifor(step, ex);
+        }
+
+        template <class Expr>
+        auto iforh(size_t step, Expr ex) const
+            -> decltype(THIS().iforh(step, ex))
+        {
+            return THIS().iforh(step, ex);
+        }
+    };
+    */
+    template <class EC, class Index>
+    EC makeec(const std::shared_ptr<Index>& i);
+
+#define V_IFOR_X(Expr) \
+    virtual ExpressionHolder<Expr> iforx(size_t step, Expr ex) const = 0; \
+    virtual ExpressionHolder<Expr> iforhx(size_t step, Expr ex) const = 0
+
+    template <class Index>
+    class E1;
+
+    class Expressions1 //: public ECInterface<Expressions1<X>>
+    {
+    public:
+        typedef Expressions1 EX;
+        
+        template <class Index>
+        static auto make(const std::shared_ptr<Index>& i)
+            -> decltype(makeec<E1<Index>>(i))
+        {
+            return makeec<E1<Index>>(i);
+        }
+
+    private:
+        V_IFOR_X(OX<EX Y() plus Y() DDMA<EX> Y() DDMA<EX>>);
+        V_IFOR_X(OX<EX Y() minus Y() DDMA<EX> Y() DDMA<EX>>);
+        V_IFOR_X(OX<EX Y() multiplies Y() DDMA<EX> Y() DDMA<EX>>);
+        V_IFOR_X(OX<EX Y() divides Y() DDMA<EX> Y() DDMA<EX>>);
+
+    public:
+        template <class Expr>
+        auto ifor(size_t step, Expr ex) const
+            -> decltype(iforx(step, ex))
+        {
+            return iforx(step, ex);
+        }
+
+        template <class Expr>
+        auto iforh(size_t step, Expr ex) const
+            -> decltype(iforhx(step, ex))
+        {
+            return iforhx(step, ex);
+        }
+
+    };
+
+    
+#define O_IFOR_X(Expr,Ind) \
+    virtual ExpressionHolder<Expr> iforx(size_t step, Expr ex) const override final \
+    { return ExpressionHolder<Expr>(Ind->ifor(step, ex)); } \
+    virtual ExpressionHolder<Expr> iforhx(size_t step, Expr ex) const override final \
+    { return ExpressionHolder<Expr>(Ind->iforh(step, ex)); }
+    
+    template <class Index>
+    class E1 : public Expressions1
+    {
+    private:
+
+        typedef Expressions1 EX;
+        
+        E1() = default;
+        
+        std::shared_ptr<Index> mI;
+
+        O_IFOR_X(OX<EX Y() plus Y() DDMA<EX> Y() DDMA<EX>>, mI);
+        O_IFOR_X(OX<EX Y() minus Y() DDMA<EX> Y() DDMA<EX>>, mI);
+        O_IFOR_X(OX<EX Y() multiplies Y() DDMA<EX> Y() DDMA<EX>>, mI);
+        O_IFOR_X(OX<EX Y() divides Y() DDMA<EX> Y() DDMA<EX>>, mI);
+
+    public:
+        E1(const E1& in) = default;
+        E1(E1&& in) = default;
+        E1& operator=(const E1& in) = default;
+        E1& operator=(E1&& in) = default;
+
+        E1(const std::shared_ptr<Index>& i) : mI(i) {} 
+    };
+
+    template <class EC, class Index>
+    EC makeec(const std::shared_ptr<Index>& i)
+    {
+        return EC(i);
+    }
+
+
+} // namespace MultiArrayTools
+
+#endif
diff --git a/src/include/functional_multi_array.cc.h b/src/include/functional_multi_array.cc.h
new file mode 100644
index 0000000..97086a3
--- /dev/null
+++ b/src/include/functional_multi_array.cc.h
@@ -0,0 +1,124 @@
+
+#include "functional_multi_array.h"
+
+namespace MultiArrayTools
+{
+
+
+    /****************************
+     *   FunctionalMultiArray   *	     
+     ****************************/
+
+    /*
+    template <typename T, class Range, class Function>
+    FunctionalMultiArray<T,Range,Function>::FunctionalMultiArray(const Range& range) :
+	MultiArrayBase<T,SRanges...>(range), mFunc() {}
+    */
+    template <bool FISSTATIC>
+    struct Application
+    {
+	template <typename T, class Function, typename Meta>
+	static inline T apply(const std::shared_ptr<Function>& f, const Meta& m)
+	{
+	    return (*f)(m);
+	}
+    };
+
+    template <>
+    struct Application<true>
+    {
+	template <typename T, class Function, typename Meta>
+	static inline T apply(const std::shared_ptr<Function>& f, const Meta& m)
+	{
+	    return Function::apply(m);
+	}
+    };
+	
+    template <typename T, class Function, class... SRanges>
+    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges,
+								      const std::shared_ptr<Function>& func) :
+	MultiArrayBase<T,SRanges...>(ranges...), mFunc(func) {}
+
+    template <typename T, class Function, class... SRanges>
+    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges) :
+	MultiArrayBase<T,SRanges...>(ranges...) {}
+
+    template <typename T, class Function, class... SRanges>
+    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const typename CRange::Space& space) :
+	MultiArrayBase<T,SRanges...>(space) {}
+
+    template <typename T, class Function, class... SRanges>
+    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const typename CRange::Space& space,
+								      const std::shared_ptr<Function>& func) :
+	MultiArrayBase<T,SRanges...>(space), mFunc(func) {}
+
+    
+    template <typename T, class Function, class... SRanges>
+    const T& FunctionalMultiArray<T,Function,SRanges...>::operator[](const IndexType& i) const
+    {
+	mVal = Application<Function::FISSTATIC>::template apply<T,Function,typename IndexType::MetaType>(mFunc, i.meta());
+	return mVal;
+    }
+
+    template <typename T, class Function, class... SRanges>
+    const T& FunctionalMultiArray<T,Function,SRanges...>::at(const typename CRange::IndexType::MetaType& meta) const
+    {
+	mVal = Application<Function::FISSTATIC>::template apply<T,Function,typename IndexType::MetaType>(mFunc,meta);
+	return mVal;
+    }
+
+    template <typename T, class Function, class... SRanges>
+    const T* FunctionalMultiArray<T,Function,SRanges...>::data() const
+    {
+	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;
+    }
+
+    template <typename T, class Function, class... SRanges>
+    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > FunctionalMultiArray<T,Function,SRanges...>::anonymous(bool slice) const
+    {
+	assert(0); // think about it carefully
+	return nullptr;
+    }
+	
+    template <typename T, class Function, class... SRanges>
+    ConstOperationRoot<T,SRanges...> FunctionalMultiArray<T,Function,SRanges...>::
+    operator()(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	if(not mMaPtr){
+	    mMaPtr = std::make_shared<MAType>( MAB::mRange->space() );
+	    (*mMaPtr)(inds...) = exec(inds...);
+	}
+	return ConstOperationRoot<T,SRanges...>( *mMaPtr, inds... );
+    }
+
+    /*
+    template <typename T, class Function, class... SRanges>
+    auto FunctionalMultiArray<T,Function,SRanges...>::
+    exec(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
+	-> decltype( mkOperation( mFunc, ConstOperationRoot<typename SRanges::IndexType::MetaType,SRanges>( mkMAObject( inds ), inds) ... ) )
+    {
+ 	return mkOperation( mFunc, ConstOperationRoot<typename SRanges::IndexType::MetaType,SRanges>( mkMAObject( inds ), inds ) ... );
+    }
+    */
+
+    template <typename T, class Function, class... SRanges>
+    auto FunctionalMultiArray<T,Function,SRanges...>::
+    exec(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
+	-> decltype( mkOperation( mFunc, mkOpObject(inds) ... ) )
+    {
+ 	return mkOperation( mFunc, mkOpObject(inds) ... );
+    }
+
+} // namespace MultiArrayTools
diff --git a/src/include/functional_multi_array.h b/src/include/functional_multi_array.h
index daff630..cf507c4 100644
--- a/src/include/functional_multi_array.h
+++ b/src/include/functional_multi_array.h
@@ -147,126 +147,5 @@ namespace MultiArrayTools
  * ---   TEMPLATE CODE   --- *
  * ========================= */
 
-namespace MultiArrayTools
-{
-
-
-    /****************************
-     *   FunctionalMultiArray   *	     
-     ****************************/
-
-    /*
-    template <typename T, class Range, class Function>
-    FunctionalMultiArray<T,Range,Function>::FunctionalMultiArray(const Range& range) :
-	MultiArrayBase<T,SRanges...>(range), mFunc() {}
-    */
-    template <bool FISSTATIC>
-    struct Application
-    {
-	template <typename T, class Function, typename Meta>
-	static inline T apply(const std::shared_ptr<Function>& f, const Meta& m)
-	{
-	    return (*f)(m);
-	}
-    };
-
-    template <>
-    struct Application<true>
-    {
-	template <typename T, class Function, typename Meta>
-	static inline T apply(const std::shared_ptr<Function>& f, const Meta& m)
-	{
-	    return Function::apply(m);
-	}
-    };
-	
-    template <typename T, class Function, class... SRanges>
-    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges,
-								      const std::shared_ptr<Function>& func) :
-	MultiArrayBase<T,SRanges...>(ranges...), mFunc(func) {}
-
-    template <typename T, class Function, class... SRanges>
-    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const std::shared_ptr<SRanges>&... ranges) :
-	MultiArrayBase<T,SRanges...>(ranges...) {}
-
-    template <typename T, class Function, class... SRanges>
-    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const typename CRange::Space& space) :
-	MultiArrayBase<T,SRanges...>(space) {}
-
-    template <typename T, class Function, class... SRanges>
-    FunctionalMultiArray<T,Function,SRanges...>::FunctionalMultiArray(const typename CRange::Space& space,
-								      const std::shared_ptr<Function>& func) :
-	MultiArrayBase<T,SRanges...>(space), mFunc(func) {}
-
-    
-    template <typename T, class Function, class... SRanges>
-    const T& FunctionalMultiArray<T,Function,SRanges...>::operator[](const IndexType& i) const
-    {
-	mVal = Application<Function::FISSTATIC>::template apply<T,Function,typename IndexType::MetaType>(mFunc, i.meta());
-	return mVal;
-    }
-
-    template <typename T, class Function, class... SRanges>
-    const T& FunctionalMultiArray<T,Function,SRanges...>::at(const typename CRange::IndexType::MetaType& meta) const
-    {
-	mVal = Application<Function::FISSTATIC>::template apply<T,Function,typename IndexType::MetaType>(mFunc,meta);
-	return mVal;
-    }
-
-    template <typename T, class Function, class... SRanges>
-    const T* FunctionalMultiArray<T,Function,SRanges...>::data() const
-    {
-	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;
-    }
-
-    template <typename T, class Function, class... SRanges>
-    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > FunctionalMultiArray<T,Function,SRanges...>::anonymous(bool slice) const
-    {
-	assert(0); // think about it carefully
-	return nullptr;
-    }
-	
-    template <typename T, class Function, class... SRanges>
-    ConstOperationRoot<T,SRanges...> FunctionalMultiArray<T,Function,SRanges...>::
-    operator()(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
-    {
-	if(not mMaPtr){
-	    mMaPtr = std::make_shared<MAType>( MAB::mRange->space() );
-	    (*mMaPtr)(inds...) = exec(inds...);
-	}
-	return ConstOperationRoot<T,SRanges...>( *mMaPtr, inds... );
-    }
-
-    /*
-    template <typename T, class Function, class... SRanges>
-    auto FunctionalMultiArray<T,Function,SRanges...>::
-    exec(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
-	-> decltype( mkOperation( mFunc, ConstOperationRoot<typename SRanges::IndexType::MetaType,SRanges>( mkMAObject( inds ), inds) ... ) )
-    {
- 	return mkOperation( mFunc, ConstOperationRoot<typename SRanges::IndexType::MetaType,SRanges>( mkMAObject( inds ), inds ) ... );
-    }
-    */
-
-    template <typename T, class Function, class... SRanges>
-    auto FunctionalMultiArray<T,Function,SRanges...>::
-    exec(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
-	-> decltype( mkOperation( mFunc, mkOpObject(inds) ... ) )
-    {
- 	return mkOperation( mFunc, mkOpObject(inds) ... );
-    }
-
-} // namespace MultiArrayTools
 
 #endif
diff --git a/src/include/helper_tools.cc.h b/src/include/helper_tools.cc.h
new file mode 100644
index 0000000..4d11dff
--- /dev/null
+++ b/src/include/helper_tools.cc.h
@@ -0,0 +1,117 @@
+
+#include "helper_tools.h"
+
+namespace MultiArrayTools
+{
+    template <typename... T>
+    std::ostream& operator<<(std::ostream& out, const std::tuple<T...>& tp)
+    {
+	PackNum<sizeof...(T)-1>::printTuple(out, tp);
+	return out;
+    }
+    
+    template <class RangeType>
+    auto getIndex(std::shared_ptr<RangeType> range)
+	-> std::shared_ptr<typename RangeType::IndexType>
+    {
+	return std::make_shared<typename RangeType::IndexType>(range);
+    }
+
+    template <class RangeType>
+    auto getIndex()
+	-> std::shared_ptr<typename RangeType::IndexType>
+    {
+	static_assert( RangeType::defaultable,
+		       /*typeid(typename RangeType).name() + */" is not defaultable" );
+	static auto f = RangeType::factory();
+	static auto r = std::dynamic_pointer_cast<RangeType>( f.create() );
+	return std::make_shared<typename RangeType::IndexType>(r);
+    }
+    
+    template <class... RangeTypes>
+    auto mkMulti(std::shared_ptr<RangeTypes>... ranges)
+	-> std::shared_ptr<MultiRange<RangeTypes...> >
+    {
+	MultiRangeFactory<RangeTypes...> mrf( ranges... );
+	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );
+    }
+
+    template <class... IndexTypes>
+    auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
+	-> decltype( getIndex( mkMulti( indices->range()... ) ) )
+    {
+	auto mi = getIndex( mkMulti( indices->range()... ) );	    
+	(*mi)( indices... );
+	return mi;
+    }
+
+    template <class Func, class... RangeTypes>
+    auto mkMapR(const std::shared_ptr<Func>& f, std::shared_ptr<RangeTypes>... ranges)
+	-> std::shared_ptr<MapRange<FunctionalMultiArray<typename Func::value_type,Func,RangeTypes...>,
+				    RangeTypes...> >
+    {
+	typedef FunctionalMultiArray<typename Func::value_type,Func,RangeTypes...> FMA;
+	std::shared_ptr<MapRangeFactory<FMA,RangeTypes...>> mrfptr;
+	if(Func::FISSTATIC){
+	    FMA fma(ranges...);
+	    mrfptr = std::make_shared<MapRangeFactory<FMA,RangeTypes...> >( fma, ranges... );
+	}
+	else {
+	    FMA fma(ranges...,f);
+	    mrfptr = std::make_shared<MapRangeFactory<FMA,RangeTypes...> >( fma, ranges... );
+	}
+	return std::dynamic_pointer_cast<MapRange<FMA,RangeTypes...> >( mrfptr->create() );
+    }
+
+    template <class Func, class... IndexTypes>
+    auto mkMapI(const std::shared_ptr<Func>& f, std::shared_ptr<IndexTypes>... indices)
+	-> decltype( getIndex( mkMapR( f, indices->range()... ) ) )
+    {
+	auto mi = getIndex( mkMapR( f, indices->range()... ) );
+	(*mi)(indices...);
+	return mi;
+    }
+
+    template <class... RangeTypes>
+    auto mkMulti(std::tuple<std::shared_ptr<RangeTypes>...> rangesTuple)
+	-> MultiRange<RangeTypes...>
+    {
+	MultiRangeFactory<RangeTypes...> mrf( rangesTuple );
+	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );
+    }
+
+    template <class RangeFactory>
+    auto createExplicit(RangeFactory& rf)
+	-> std::shared_ptr<typename RangeFactory::oType>
+    {
+	return std::dynamic_pointer_cast<typename RangeFactory::oType>( rf.create() );
+    }
+
+    template <class RangeFactory>
+    auto createExplicit(std::shared_ptr<RangeFactory> rfp)
+	-> std::shared_ptr<typename RangeFactory::oType>
+    {
+	return std::dynamic_pointer_cast<typename RangeFactory::oType>( rfp->create() );
+    }
+
+    template <class Range>
+    auto createRange(const std::vector<char>& cvec)
+	-> std::shared_ptr<Range>
+    {
+	const char* dp = cvec.data();
+	auto ff = createRangeFactory(&dp);
+	auto rbptr = ff->create();
+	assert(rbptr->spaceType() == Range::STYPE);
+	// CATCH CAST ERROR HERE !!!
+	return std::dynamic_pointer_cast<Range>( rbptr );
+    }
+
+    template <size_t N, class MArray>
+    auto rptr(const MArray& ma)
+	-> decltype(ma.template getRangePtr<N>())
+    {
+	return ma.template getRangePtr<N>();
+    }
+
+
+}
diff --git a/src/include/helper_tools.h b/src/include/helper_tools.h
index c8663b0..b43ca02 100644
--- a/src/include/helper_tools.h
+++ b/src/include/helper_tools.h
@@ -59,132 +59,6 @@ namespace MultiArrayTools
     auto prtr(const MArray& ma)
 	-> decltype(ma.template getRangePtr<N>());
 
-    template <class IndexType>
-    inline void For(const std::shared_ptr<IndexType>& ind, const std::function<void(void)>& ll);
-
-    template <class Index>
-    inline auto mkOp(const std::shared_ptr<Index>& i)
-	-> decltype(std::declval<FunctionalMultiArray<typename Index::MetaType,
-		    identity<typename Index::MetaType>,typename Index::RangeType> >
-		    ().exec(i));
-}
-
-/* ========================= *
- * ---   TEMPLATE CODE   --- *
- * ========================= */
-
-namespace MultiArrayTools
-{
-    template <typename... T>
-    std::ostream& operator<<(std::ostream& out, const std::tuple<T...>& tp)
-    {
-	PackNum<sizeof...(T)-1>::printTuple(out, tp);
-	return out;
-    }
-    
-    template <class RangeType>
-    auto getIndex(std::shared_ptr<RangeType> range)
-	-> std::shared_ptr<typename RangeType::IndexType>
-    {
-	return std::make_shared<typename RangeType::IndexType>(range);
-    }
-
-    template <class RangeType>
-    auto getIndex()
-	-> std::shared_ptr<typename RangeType::IndexType>
-    {
-	static_assert( RangeType::defaultable,
-		       /*typeid(typename RangeType).name() + */" is not defaultable" );
-	static auto f = RangeType::factory();
-	static auto r = std::dynamic_pointer_cast<RangeType>( f.create() );
-	return std::make_shared<typename RangeType::IndexType>(r);
-    }
-    
-    template <class... RangeTypes>
-    auto mkMulti(std::shared_ptr<RangeTypes>... ranges)
-	-> std::shared_ptr<MultiRange<RangeTypes...> >
-    {
-	MultiRangeFactory<RangeTypes...> mrf( ranges... );
-	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );
-    }
-
-    template <class... IndexTypes>
-    auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
-	-> decltype( getIndex( mkMulti( indices->range()... ) ) )
-    {
-	auto mi = getIndex( mkMulti( indices->range()... ) );	    
-	(*mi)( indices... );
-	return mi;
-    }
-
-    template <class Func, class... RangeTypes>
-    auto mkMapR(const std::shared_ptr<Func>& f, std::shared_ptr<RangeTypes>... ranges)
-	-> std::shared_ptr<MapRange<FunctionalMultiArray<typename Func::value_type,Func,RangeTypes...>,
-				    RangeTypes...> >
-    {
-	typedef FunctionalMultiArray<typename Func::value_type,Func,RangeTypes...> FMA;
-	std::shared_ptr<MapRangeFactory<FMA,RangeTypes...>> mrfptr;
-	if(Func::FISSTATIC){
-	    FMA fma(ranges...);
-	    mrfptr = std::make_shared<MapRangeFactory<FMA,RangeTypes...> >( fma, ranges... );
-	}
-	else {
-	    FMA fma(ranges...,f);
-	    mrfptr = std::make_shared<MapRangeFactory<FMA,RangeTypes...> >( fma, ranges... );
-	}
-	return std::dynamic_pointer_cast<MapRange<FMA,RangeTypes...> >( mrfptr->create() );
-    }
-
-    template <class Func, class... IndexTypes>
-    auto mkMapI(const std::shared_ptr<Func>& f, std::shared_ptr<IndexTypes>... indices)
-	-> decltype( getIndex( mkMapR( f, indices->range()... ) ) )
-    {
-	auto mi = getIndex( mkMapR( f, indices->range()... ) );
-	(*mi)(indices...);
-	return mi;
-    }
-
-    template <class... RangeTypes>
-    auto mkMulti(std::tuple<std::shared_ptr<RangeTypes>...> rangesTuple)
-	-> MultiRange<RangeTypes...>
-    {
-	MultiRangeFactory<RangeTypes...> mrf( rangesTuple );
-	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );
-    }
-
-    template <class RangeFactory>
-    auto createExplicit(RangeFactory& rf)
-	-> std::shared_ptr<typename RangeFactory::oType>
-    {
-	return std::dynamic_pointer_cast<typename RangeFactory::oType>( rf.create() );
-    }
-
-    template <class RangeFactory>
-    auto createExplicit(std::shared_ptr<RangeFactory> rfp)
-	-> std::shared_ptr<typename RangeFactory::oType>
-    {
-	return std::dynamic_pointer_cast<typename RangeFactory::oType>( rfp->create() );
-    }
-
-    template <class Range>
-    auto createRange(const std::vector<char>& cvec)
-	-> std::shared_ptr<Range>
-    {
-	const char* dp = cvec.data();
-	auto ff = createRangeFactory(&dp);
-	auto rbptr = ff->create();
-	assert(rbptr->spaceType() == Range::STYPE);
-	// CATCH CAST ERROR HERE !!!
-	return std::dynamic_pointer_cast<Range>( rbptr );
-    }
-
-    template <size_t N, class MArray>
-    auto rptr(const MArray& ma)
-	-> decltype(ma.template getRangePtr<N>())
-    {
-	return ma.template getRangePtr<N>();
-    }
-
     template <class IndexType>
     inline void For(const std::shared_ptr<IndexType>& ind, const std::function<void(void)>& ll)
     {
@@ -193,6 +67,7 @@ namespace MultiArrayTools
 	}
     }
 
+
     template <class Index>
     inline auto mkOp(const std::shared_ptr<Index>& i)
 	-> decltype(std::declval<FunctionalMultiArray<typename Index::MetaType,
@@ -204,7 +79,8 @@ namespace MultiArrayTools
 			     typename Index::RangeType> fma(i->range());
 	return fma.exec(i);
     }
-
+        
 }
 
+
 #endif
diff --git a/src/include/map_range.cc.h b/src/include/map_range.cc.h
new file mode 100644
index 0000000..f7633d1
--- /dev/null
+++ b/src/include/map_range.cc.h
@@ -0,0 +1,512 @@
+
+#include "map_range.h"
+
+namespace MultiArrayTools
+{
+
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+
+    /**************
+     *   OpExpr   *
+     **************/
+
+    template <class MapF, class IndexPack, class Expr>
+    OpExpr<MapF,IndexPack,Expr>::OpExpr(const MapF& mapf, const IndexPack& ipack,
+					const std::shared_ptr<OIType>& oind, size_t step, Expr ex) :
+	mIndPtr(oind.get()), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()),
+	mStep(step), mExpr( std::forward<Expr>(ex) ),
+	mOp(mkMapOp(mapf, ipack)),
+	//mExt(ex.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )))
+	mExt( mOp.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ) ).extend
+	      ( ex.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ) ) ) )
+    {
+	assert(mIndPtr != nullptr);
+    }
+
+    template <class MapF, class IndexPack, class Expr>
+    inline void OpExpr<MapF,IndexPack,Expr>::operator()(size_t mlast,
+							 ExtType last) const
+    {
+	constexpr size_t NEXT = Op::SIZE;
+	const ExtType nxpos = last;
+	const size_t pos = mIndPtr->posAt( mOp.get( nxpos ) );
+	//VCHECK(pos);
+	const ExtType npos = last + mExt*pos;
+	//VCHECK(npos.next().next().val());
+	const size_t mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);
+	mExpr(mnpos, Getter<NEXT>::template getX<ExtType>( npos ) );
+    }
+
+    template <class MapF, class IndexPack, class Expr>
+    inline void OpExpr<MapF,IndexPack,Expr>::operator()(size_t mlast) const
+    {
+	const ExtType last;
+	constexpr size_t NEXT = Op::SIZE;
+	const ExtType nxpos = last;
+	const size_t pos = mIndPtr->posAt( mOp.get( nxpos ) );
+	const ExtType npos = last + mExt*pos;
+	const size_t mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);
+	mExpr(mnpos, Getter<NEXT>::template getX<ExtType>( npos ));
+    }
+    
+    template <class MapF, class IndexPack, class Expr>
+    auto OpExpr<MapF,IndexPack,Expr>::rootSteps(std::intptr_t iPtrNum) const
+	-> ExtType
+    {
+	return mOp.rootSteps(iPtrNum).extend( mExpr.rootSteps(iPtrNum) );
+	//return mExpr.rootSteps(iPtrNum).extend( mOp.rootSteps(iPtrNum) );
+    }
+
+    
+    // -> define in range_base.cc
+    //std::shared_ptr<RangeFactoryBase> mkMULTI(const char** dp);
+    
+    /******************
+     *   MapIndex   *
+     ******************/
+
+    /*
+    template <class MapF, class... Indices>
+    MapIndex<MapF,Indices...>::MapIndex(const MapIndex<MapF,Indices...>& in) :
+	IndexInterface<std::tuple<typename Indices::MetaType...> >(in)
+    {
+	RPackNum<sizeof...(Indices)-1>::copy(mIPack, in);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+    }
+	
+    template <class MapF, class... Indices>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator=(const MapIndex<MapF,Indices...>& in)
+    {
+	IndexI::operator=(in);
+	RPackNum<sizeof...(Indices)-1>::copy(mIPack, in);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+    */
+    
+    template <class MapF, class... Indices>
+    template <class MRange>
+    MapIndex<MapF,Indices...>::MapIndex(const std::shared_ptr<MRange>& range) :
+	IndexInterface<MapIndex<MapF,Indices...>,std::tuple<typename Indices::MetaType...> >(range, 0)
+    {
+	RPackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
+	RPackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack); // has one more element!
+	mOutIndex = std::make_shared<OIType>
+	    ( std::dynamic_pointer_cast<RangeType>( IB::mRangePtr )->outRange()->begin() );
+    }
+
+    template <class MapF, class... Indices>
+    template <size_t DIR>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::up()
+    {
+	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
+	IB::mPos += RPackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
+	RPackNum<DIR>::pp( mIPack );
+	return *this;
+    }
+
+    template <class MapF, class... Indices>
+    template <size_t DIR>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::down()
+    {
+	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
+	IB::mPos -= RPackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
+	RPackNum<DIR>::mm( mIPack );
+	return *this;
+    }
+    
+    template <class MapF, class... Indices>
+    template <size_t N>
+    auto MapIndex<MapF,Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
+    {
+	return *std::get<N>(mIPack);
+    }
+
+    template <class MapF, class... Indices>
+    template <size_t N>
+    auto MapIndex<MapF,Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>(mIPack);
+    }
+
+    template <class MapF, class... Indices>
+    auto MapIndex<MapF,Indices...>::outIndex() const -> std::shared_ptr<OIType>
+    {
+	return mOutIndex;
+    }
+    
+    template <class MapF, class... Indices>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator()(std::shared_ptr<Indices>&... indices)
+    {
+	RPackNum<sizeof...(Indices)-1>::swapIndices(mIPack, indices...);
+	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, IB::mPos);
+	return *this;
+    }
+
+    template <class MapF, class... Indices>
+    IndexType MapIndex<MapF,Indices...>::type() const
+    {
+	return IndexType::MULTI;
+    }
+
+    template <class MapF, class... Indices>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
+	return *this;
+    }
+
+    template <class MapF, class... Indices>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator++()
+    {
+	RPackNum<sizeof...(Indices)-1>::pp( mIPack );
+	++IB::mPos;
+	return *this;
+    }
+
+    template <class MapF, class... Indices>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator--()
+    {
+	RPackNum<sizeof...(Indices)-1>::mm( mIPack );
+	--IB::mPos;
+	return *this;
+    }
+
+    template <class MapF, class... Indices>
+    int MapIndex<MapF,Indices...>::pp(std::intptr_t idxPtrNum)
+    {
+	int tmp = RPackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtrNum);
+	IB::mPos += tmp;
+	return tmp;
+    }
+
+    template <class MapF, class... Indices>
+    int MapIndex<MapF,Indices...>::mm(std::intptr_t idxPtrNum)
+    {
+	int tmp = RPackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtrNum);
+	IB::mPos -= tmp;
+	return tmp;
+    }
+
+    template <class MapF, class... Indices>
+    std::string MapIndex<MapF,Indices...>::stringMeta() const
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr )->stringMeta(IB::mPos);
+    }
+    
+    template <class MapF, class... Indices>
+    typename MapIndex<MapF,Indices...>::MetaType MapIndex<MapF,Indices...>::meta() const
+    {
+	MetaType metaTuple;
+	RPackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
+	return metaTuple;
+    }
+
+    template <class MapF, class... Indices>
+    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::at(const MetaType& metaPos)
+    {
+	RPackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
+	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
+	return *this;
+    }
+
+    template <class MapF, class... Indices>
+    size_t MapIndex<MapF,Indices...>::dim() const
+    {
+	return sizeof...(Indices);
+    }
+
+    template <class MapF, class... Indices>
+    bool MapIndex<MapF,Indices...>::first() const
+    {
+	return IB::mPos == 0;
+    }
+
+    template <class MapF, class... Indices>
+    bool MapIndex<MapF,Indices...>::last() const
+    {
+	return IB::mPos == IB::mMax - 1;
+    }
+
+    template <class MapF, class... Indices>
+    std::shared_ptr<typename MapIndex<MapF,Indices...>::RangeType>
+    MapIndex<MapF,Indices...>::range() const
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
+    }
+
+    template <class MapF, class... Indices>
+    template <size_t N>
+    auto MapIndex<MapF,Indices...>::getPtr() -> decltype( std::get<N>( mIPack ) )&
+    {
+	return std::get<N>(mIPack);
+    }
+
+    template <class MapF, class... Indices>	
+    size_t MapIndex<MapF,Indices...>::getStepSize(size_t n) const
+    {
+	if(n >= sizeof...(Indices)){
+	    assert(0);
+	    // throw !!
+	}
+	return mBlockSizes[n+1];
+    }
+    
+    template <class MapF, class... Indices>
+    std::string MapIndex<MapF,Indices...>::id() const
+    {
+	return std::string("mul") + std::to_string(IB::mId);
+    }
+
+    template <class MapF, class... Indices>
+    void MapIndex<MapF,Indices...>::print(size_t offset) const
+    {
+	if(offset == 0){
+	    std::cout << " === " << std::endl;
+	}
+	for(size_t j = 0; j != offset; ++j) { std::cout << "\t"; }
+	std::cout << id() << "[" << reinterpret_cast<std::intptr_t>(this)
+		  << "]" << "(" << IB::mRangePtr << "): " << meta() << std::endl;
+	RPackNum<sizeof...(Indices)-1>::printIndex(mIPack, offset+1);
+    }
+
+    template <class MapF, class... Indices>
+    template <class Exprs>
+    auto MapIndex<MapF,Indices...>::ifor(size_t step, Exprs exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh
+		    (step, mIPack, mBlockSizes, OpExpr<MapF,IndexPack,Exprs>
+		     ( range()->map(), mIPack, mOutIndex, step, exs ) ) )
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkForh
+	    (step, mIPack, mBlockSizes, OpExpr<MapF,IndexPack,Exprs>
+	     ( range()->map(), mIPack, mOutIndex, step, exs ) );
+    }
+    /*
+    template <class MapF, class... Indices>
+    template <class Exprs>
+    auto MapIndex<MapF,Indices...>::iforh(Exprs exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs))
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs);
+    }
+    */
+    /*************************
+     *   MapRangeFactory   *
+     *************************/
+
+    template <class MapF, class... Ranges>
+    MapRangeFactory<MapF,Ranges...>::MapRangeFactory(const MapF& mapf, const std::shared_ptr<Ranges>&... rs)
+    {
+	mProd = std::shared_ptr< MapRange<MapF,Ranges...> >( new MapRange<MapF,Ranges...>( mapf, rs... ) );
+    }
+    
+    template <class MapF, class... Ranges>
+    MapRangeFactory<MapF,Ranges...>::MapRangeFactory(const MapF& mapf, const typename MapRange<MapF,Ranges...>::Space& st)
+    {
+	mProd = std::shared_ptr< MapRange<MapF,Ranges...> >( new MapRange<MapF,Ranges...>( mapf, st ) );
+    }
+    
+    template <class MapF, class... Ranges>
+    std::shared_ptr<RangeBase> MapRangeFactory<MapF,Ranges...>::create()
+    {
+	mProd = checkIfCreated( std::dynamic_pointer_cast<oType>( mProd )->mSpace );
+	setSelf();
+	return mProd;
+    }
+
+    template <class MapF, class... Ranges>
+    std::shared_ptr<RangeBase> MapRangeFactory<MapF,Ranges...>::checkIfCreated(const std::tuple<std::shared_ptr<Ranges>...>& ptp)
+    {
+	std::shared_ptr<RangeBase> out;
+	bool check = false;
+	for(auto& x: MapRangeFactoryProductMap::mAleadyCreated){
+	    if(x.second.size() == sizeof...(Ranges)){
+		check = RPackNum<sizeof...(Ranges)-1>::checkIfCreated(ptp, x.second);
+		if(check){
+		    out = x.first;
+		    break;
+		}
+	    }
+	}
+	if(not check){
+	    std::vector<std::intptr_t> pv(sizeof...(Ranges));
+	    RPackNum<sizeof...(Ranges)-1>::RangesToVec(ptp, pv);
+	    pv.push_back( reinterpret_cast<std::intptr_t>
+			  ( &std::dynamic_pointer_cast<oType>( mProd )->mMapf ) );
+	    MapRangeFactoryProductMap::mAleadyCreated[mProd] = pv;
+	    out = mProd;
+	}
+	return out;
+    }
+    
+    /******************
+     *   MapRange   *
+     ******************/
+
+    template <class MapF, class... Ranges>
+    void MapRange<MapF,Ranges...>::mkOutRange()
+    {
+	//FunctionalMultiArray<typename MapF::value_type,MapF,Ranges...> fma(mSpace, mMapf);
+	std::map<typename MapF::value_type,size_t> mult;
+	for(auto ii = mMapf.begin(); ii.max() != ii.pos(); ++ii) {
+	    mult[mMapf[ii]]++;
+	}
+	
+	std::vector<typename MapF::value_type> outmeta(mult.size());
+	std::vector<size_t> outmult(mult.size());
+
+	size_t cnt = 0;
+	for(auto& x: mult){
+	    outmeta[cnt] = x.first;
+	    outmult[cnt] = x.second;
+	    ++cnt;
+	}
+
+	ORFType orf(outmeta);
+	mOutRange = std::dynamic_pointer_cast<ORType>( orf.create() );
+	mMapMult = MultiArray<size_t,ORType>( mOutRange, outmult );
+    }
+    
+    template <class MapF, class... Ranges>
+    MapRange<MapF,Ranges...>::MapRange(const MapF& mapf, const std::shared_ptr<Ranges>&... rs) :
+	mSpace(std::make_tuple(rs...)),
+	mMapf(mapf)
+    {
+	mkOutRange();
+    }
+
+    template <class MapF, class... Ranges>
+    MapRange<MapF,Ranges...>::MapRange(const MapF& mapf, const Space& space) :
+	mSpace( space ),
+	mMapf(mapf)
+    {
+	mkOutRange();
+    }
+
+    template <class MapF, class... Ranges>
+    template <size_t N>
+    auto MapRange<MapF,Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
+    {
+	return *std::get<N>(mSpace);
+    }
+
+    template <class MapF, class... Ranges>
+    template <size_t N>
+    auto MapRange<MapF,Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
+    {
+	return std::get<N>(mSpace);
+    }
+
+    template <class MapF, class... Ranges>
+    auto MapRange<MapF,Ranges...>::outRange() const -> std::shared_ptr<ORType>
+    {
+	return mOutRange;
+    }
+
+    template <class MapF, class... Ranges>
+    const MapF& MapRange<MapF,Ranges...>::map() const
+    {
+	return mMapf;
+    }
+    
+    template <class MapF, class... Ranges>
+    size_t MapRange<MapF,Ranges...>::dim() const
+    {
+	return sdim;
+    }
+    
+    template <class MapF, class... Ranges>
+    size_t MapRange<MapF,Ranges...>::size() const
+    {
+	return mOutRange->size();
+	//return RPackNum<sizeof...(Ranges)-1>::getSize(mSpace);
+    }
+
+    template <class MapF, class... Ranges>
+    SpaceType MapRange<MapF,Ranges...>::spaceType() const
+    {
+	return SpaceType::ANY;
+    }
+    
+    template <class MapF, class... Ranges>
+    const typename MapRange<MapF,Ranges...>::Space& MapRange<MapF,Ranges...>::space() const
+    {
+	return mSpace;
+    }
+
+    template <class MapF, class... Ranges>
+    std::string MapRange<MapF,Ranges...>::stringMeta(size_t pos) const
+    {
+	auto i = begin();
+	i = pos;
+	return "[ " + RPackNum<sizeof...(Ranges)-1>::metaTupleToString(i.meta()) + " ]";
+    }
+
+    template <class MapF, class... Ranges>
+    std::vector<char> MapRange<MapF,Ranges...>::data() const
+    {
+	DataHeader h;
+	h.spaceType = static_cast<int>( SpaceType::ANY );
+	h.metaSize = sizeof...(Ranges);
+	h.multiple = 1;
+	std::vector<char> out;
+	//out.reserve(h.metaSize + sizeof(DataHeader));
+	char* hcp = reinterpret_cast<char*>(&h);
+	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
+	RPackNum<sizeof...(Ranges)-1>::fillRangeDataVec(out, mSpace);
+	return out;
+    }
+    
+    template <class MapF, class... Ranges>
+    typename MapRange<MapF,Ranges...>::IndexType MapRange<MapF,Ranges...>::begin() const
+    {
+	MapIndex<MapF,typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<MapRange<MapF,Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+	i = 0;
+	return i;
+    }
+
+    template <class MapF, class... Ranges>
+    typename MapRange<MapF,Ranges...>::IndexType MapRange<MapF,Ranges...>::end() const
+    {
+	MapIndex<MapF,typename Ranges::IndexType...>
+	    i( std::dynamic_pointer_cast<MapRange<MapF,Ranges...> >
+	       ( std::shared_ptr<RangeBase>( RB::mThis )) );
+	i = size();
+	return i;
+    }
+
+    template <class MapF, class... Ranges>
+    auto MapRange<MapF,Ranges...>::mapMultiplicity() const
+	-> const MultiArray<size_t,ORType>&
+    {
+	return mMapMult;
+    }
+
+    template <class MapF, class... Ranges>
+    auto MapRange<MapF,Ranges...>::explMapMultiplicity() const
+	-> MultiArray<size_t,MapRange>
+    {
+	auto tmp = mMapMult;
+	return tmp.format( std::dynamic_pointer_cast<MapRange<MapF,Ranges...> >
+			   ( std::shared_ptr<RangeBase>( RB::mThis )) ); 
+    }
+    
+    template <class MapF, class... Ranges>
+    template <class... ERanges>
+    auto MapRange<MapF,Ranges...>::cat(const std::shared_ptr<MapRange<ERanges...> >& erange)
+	-> std::shared_ptr<MapRange<Ranges...,ERanges...> >
+    {
+	auto crange = std::tuple_cat(mSpace, erange->space());
+	MapRangeFactory<Ranges...,ERanges...> rf(crange);
+	return std::dynamic_pointer_cast<MapRange<Ranges...,ERanges...> >(rf.create());
+    }
+
+}
diff --git a/src/include/map_range.h b/src/include/map_range.h
index 2726e86..57562d4 100644
--- a/src/include/map_range.h
+++ b/src/include/map_range.h
@@ -284,518 +284,5 @@ namespace MultiArrayTools
     
 }
 
-/* ========================= *
- * ---   TEMPLATE CODE   --- *
- * ========================= */
-
-namespace MultiArrayTools
-{
-
-    namespace
-    {
-	using namespace MultiArrayHelper;
-    }
-
-    /**************
-     *   OpExpr   *
-     **************/
-
-    template <class MapF, class IndexPack, class Expr>
-    OpExpr<MapF,IndexPack,Expr>::OpExpr(const MapF& mapf, const IndexPack& ipack,
-					const std::shared_ptr<OIType>& oind, size_t step, Expr ex) :
-	mIndPtr(oind.get()), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()),
-	mStep(step), mExpr( std::forward<Expr>(ex) ),
-	mOp(mkMapOp(mapf, ipack)),
-	//mExt(ex.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )))
-	mExt( mOp.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ) ).extend
-	      ( ex.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ) ) ) )
-    {
-	assert(mIndPtr != nullptr);
-    }
-
-    template <class MapF, class IndexPack, class Expr>
-    inline void OpExpr<MapF,IndexPack,Expr>::operator()(size_t mlast,
-							 ExtType last) const
-    {
-	constexpr size_t NEXT = Op::SIZE;
-	const ExtType nxpos = last;
-	const size_t pos = mIndPtr->posAt( mOp.get( nxpos ) );
-	//VCHECK(pos);
-	const ExtType npos = last + mExt*pos;
-	//VCHECK(npos.next().next().val());
-	const size_t mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);
-	mExpr(mnpos, Getter<NEXT>::template getX<ExtType>( npos ) );
-    }
-
-    template <class MapF, class IndexPack, class Expr>
-    inline void OpExpr<MapF,IndexPack,Expr>::operator()(size_t mlast) const
-    {
-	const ExtType last;
-	constexpr size_t NEXT = Op::SIZE;
-	const ExtType nxpos = last;
-	const size_t pos = mIndPtr->posAt( mOp.get( nxpos ) );
-	const ExtType npos = last + mExt*pos;
-	const size_t mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);
-	mExpr(mnpos, Getter<NEXT>::template getX<ExtType>( npos ));
-    }
-    
-    template <class MapF, class IndexPack, class Expr>
-    auto OpExpr<MapF,IndexPack,Expr>::rootSteps(std::intptr_t iPtrNum) const
-	-> ExtType
-    {
-	return mOp.rootSteps(iPtrNum).extend( mExpr.rootSteps(iPtrNum) );
-	//return mExpr.rootSteps(iPtrNum).extend( mOp.rootSteps(iPtrNum) );
-    }
-
-    
-    // -> define in range_base.cc
-    //std::shared_ptr<RangeFactoryBase> mkMULTI(const char** dp);
-    
-    /******************
-     *   MapIndex   *
-     ******************/
-
-    /*
-    template <class MapF, class... Indices>
-    MapIndex<MapF,Indices...>::MapIndex(const MapIndex<MapF,Indices...>& in) :
-	IndexInterface<std::tuple<typename Indices::MetaType...> >(in)
-    {
-	RPackNum<sizeof...(Indices)-1>::copy(mIPack, in);
-	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
-    }
-	
-    template <class MapF, class... Indices>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator=(const MapIndex<MapF,Indices...>& in)
-    {
-	IndexI::operator=(in);
-	RPackNum<sizeof...(Indices)-1>::copy(mIPack, in);
-	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	return *this;
-    }
-    */
-    
-    template <class MapF, class... Indices>
-    template <class MRange>
-    MapIndex<MapF,Indices...>::MapIndex(const std::shared_ptr<MRange>& range) :
-	IndexInterface<MapIndex<MapF,Indices...>,std::tuple<typename Indices::MetaType...> >(range, 0)
-    {
-	RPackNum<sizeof...(Indices)-1>::construct(mIPack, *range);
-	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	std::get<sizeof...(Indices)>(mBlockSizes) = 1;
-	RPackNum<sizeof...(Indices)-1>::initBlockSizes(mBlockSizes, mIPack); // has one more element!
-	mOutIndex = std::make_shared<OIType>
-	    ( std::dynamic_pointer_cast<RangeType>( IB::mRangePtr )->outRange()->begin() );
-    }
-
-    template <class MapF, class... Indices>
-    template <size_t DIR>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::up()
-    {
-	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
-	IB::mPos += RPackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
-	RPackNum<DIR>::pp( mIPack );
-	return *this;
-    }
-
-    template <class MapF, class... Indices>
-    template <size_t DIR>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::down()
-    {
-	static_assert(DIR < sizeof...(Indices), "DIR exceeds number of sub-indices");
-	IB::mPos -= RPackNum<sizeof...(Indices)-DIR-1>::blockSize( mIPack );
-	RPackNum<DIR>::mm( mIPack );
-	return *this;
-    }
-    
-    template <class MapF, class... Indices>
-    template <size_t N>
-    auto MapIndex<MapF,Indices...>::get() const -> decltype( *std::get<N>( mIPack ) )&
-    {
-	return *std::get<N>(mIPack);
-    }
-
-    template <class MapF, class... Indices>
-    template <size_t N>
-    auto MapIndex<MapF,Indices...>::getPtr() const -> decltype( std::get<N>( mIPack ) )&
-    {
-	return std::get<N>(mIPack);
-    }
-
-    template <class MapF, class... Indices>
-    auto MapIndex<MapF,Indices...>::outIndex() const -> std::shared_ptr<OIType>
-    {
-	return mOutIndex;
-    }
-    
-    template <class MapF, class... Indices>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator()(std::shared_ptr<Indices>&... indices)
-    {
-	RPackNum<sizeof...(Indices)-1>::swapIndices(mIPack, indices...);
-	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, IB::mPos);
-	return *this;
-    }
-
-    template <class MapF, class... Indices>
-    IndexType MapIndex<MapF,Indices...>::type() const
-    {
-	return IndexType::MULTI;
-    }
-
-    template <class MapF, class... Indices>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator=(size_t pos)
-    {
-	IB::mPos = pos;
-	RPackNum<sizeof...(Indices)-1>::setIndexPack(mIPack, pos);
-	return *this;
-    }
-
-    template <class MapF, class... Indices>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator++()
-    {
-	RPackNum<sizeof...(Indices)-1>::pp( mIPack );
-	++IB::mPos;
-	return *this;
-    }
-
-    template <class MapF, class... Indices>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::operator--()
-    {
-	RPackNum<sizeof...(Indices)-1>::mm( mIPack );
-	--IB::mPos;
-	return *this;
-    }
-
-    template <class MapF, class... Indices>
-    int MapIndex<MapF,Indices...>::pp(std::intptr_t idxPtrNum)
-    {
-	int tmp = RPackNum<sizeof...(Indices)-1>::pp(mIPack, mBlockSizes, idxPtrNum);
-	IB::mPos += tmp;
-	return tmp;
-    }
-
-    template <class MapF, class... Indices>
-    int MapIndex<MapF,Indices...>::mm(std::intptr_t idxPtrNum)
-    {
-	int tmp = RPackNum<sizeof...(Indices)-1>::mm(mIPack, mBlockSizes, idxPtrNum);
-	IB::mPos -= tmp;
-	return tmp;
-    }
-
-    template <class MapF, class... Indices>
-    std::string MapIndex<MapF,Indices...>::stringMeta() const
-    {
-	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr )->stringMeta(IB::mPos);
-    }
-    
-    template <class MapF, class... Indices>
-    typename MapIndex<MapF,Indices...>::MetaType MapIndex<MapF,Indices...>::meta() const
-    {
-	MetaType metaTuple;
-	RPackNum<sizeof...(Indices)-1>::getMetaPos(metaTuple, mIPack);
-	return metaTuple;
-    }
-
-    template <class MapF, class... Indices>
-    MapIndex<MapF,Indices...>& MapIndex<MapF,Indices...>::at(const MetaType& metaPos)
-    {
-	RPackNum<sizeof...(Indices)-1>::setMeta(mIPack, metaPos);
-	IB::mPos = RPackNum<sizeof...(Indices)-1>::makePos(mIPack);
-	return *this;
-    }
-
-    template <class MapF, class... Indices>
-    size_t MapIndex<MapF,Indices...>::dim() const
-    {
-	return sizeof...(Indices);
-    }
-
-    template <class MapF, class... Indices>
-    bool MapIndex<MapF,Indices...>::first() const
-    {
-	return IB::mPos == 0;
-    }
-
-    template <class MapF, class... Indices>
-    bool MapIndex<MapF,Indices...>::last() const
-    {
-	return IB::mPos == IB::mMax - 1;
-    }
-
-    template <class MapF, class... Indices>
-    std::shared_ptr<typename MapIndex<MapF,Indices...>::RangeType>
-    MapIndex<MapF,Indices...>::range() const
-    {
-	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
-    }
-
-    template <class MapF, class... Indices>
-    template <size_t N>
-    auto MapIndex<MapF,Indices...>::getPtr() -> decltype( std::get<N>( mIPack ) )&
-    {
-	return std::get<N>(mIPack);
-    }
-
-    template <class MapF, class... Indices>	
-    size_t MapIndex<MapF,Indices...>::getStepSize(size_t n) const
-    {
-	if(n >= sizeof...(Indices)){
-	    assert(0);
-	    // throw !!
-	}
-	return mBlockSizes[n+1];
-    }
-    
-    template <class MapF, class... Indices>
-    std::string MapIndex<MapF,Indices...>::id() const
-    {
-	return std::string("mul") + std::to_string(IB::mId);
-    }
-
-    template <class MapF, class... Indices>
-    void MapIndex<MapF,Indices...>::print(size_t offset) const
-    {
-	if(offset == 0){
-	    std::cout << " === " << std::endl;
-	}
-	for(size_t j = 0; j != offset; ++j) { std::cout << "\t"; }
-	std::cout << id() << "[" << reinterpret_cast<std::intptr_t>(this)
-		  << "]" << "(" << IB::mRangePtr << "): " << meta() << std::endl;
-	RPackNum<sizeof...(Indices)-1>::printIndex(mIPack, offset+1);
-    }
-
-    template <class MapF, class... Indices>
-    template <class Exprs>
-    auto MapIndex<MapF,Indices...>::ifor(size_t step, Exprs exs) const
-	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh
-		    (step, mIPack, mBlockSizes, OpExpr<MapF,IndexPack,Exprs>
-		     ( range()->map(), mIPack, mOutIndex, step, exs ) ) )
-    {
-	return RPackNum<sizeof...(Indices)-1>::mkForh
-	    (step, mIPack, mBlockSizes, OpExpr<MapF,IndexPack,Exprs>
-	     ( range()->map(), mIPack, mOutIndex, step, exs ) );
-    }
-    /*
-    template <class MapF, class... Indices>
-    template <class Exprs>
-    auto MapIndex<MapF,Indices...>::iforh(Exprs exs) const
-	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs))
-    {
-	return RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs);
-    }
-    */
-    /*************************
-     *   MapRangeFactory   *
-     *************************/
-
-    template <class MapF, class... Ranges>
-    MapRangeFactory<MapF,Ranges...>::MapRangeFactory(const MapF& mapf, const std::shared_ptr<Ranges>&... rs)
-    {
-	mProd = std::shared_ptr< MapRange<MapF,Ranges...> >( new MapRange<MapF,Ranges...>( mapf, rs... ) );
-    }
-    
-    template <class MapF, class... Ranges>
-    MapRangeFactory<MapF,Ranges...>::MapRangeFactory(const MapF& mapf, const typename MapRange<MapF,Ranges...>::Space& st)
-    {
-	mProd = std::shared_ptr< MapRange<MapF,Ranges...> >( new MapRange<MapF,Ranges...>( mapf, st ) );
-    }
-    
-    template <class MapF, class... Ranges>
-    std::shared_ptr<RangeBase> MapRangeFactory<MapF,Ranges...>::create()
-    {
-	mProd = checkIfCreated( std::dynamic_pointer_cast<oType>( mProd )->mSpace );
-	setSelf();
-	return mProd;
-    }
-
-    template <class MapF, class... Ranges>
-    std::shared_ptr<RangeBase> MapRangeFactory<MapF,Ranges...>::checkIfCreated(const std::tuple<std::shared_ptr<Ranges>...>& ptp)
-    {
-	std::shared_ptr<RangeBase> out;
-	bool check = false;
-	for(auto& x: MapRangeFactoryProductMap::mAleadyCreated){
-	    if(x.second.size() == sizeof...(Ranges)){
-		check = RPackNum<sizeof...(Ranges)-1>::checkIfCreated(ptp, x.second);
-		if(check){
-		    out = x.first;
-		    break;
-		}
-	    }
-	}
-	if(not check){
-	    std::vector<std::intptr_t> pv(sizeof...(Ranges));
-	    RPackNum<sizeof...(Ranges)-1>::RangesToVec(ptp, pv);
-	    pv.push_back( reinterpret_cast<std::intptr_t>
-			  ( &std::dynamic_pointer_cast<oType>( mProd )->mMapf ) );
-	    MapRangeFactoryProductMap::mAleadyCreated[mProd] = pv;
-	    out = mProd;
-	}
-	return out;
-    }
-    
-    /******************
-     *   MapRange   *
-     ******************/
-
-    template <class MapF, class... Ranges>
-    void MapRange<MapF,Ranges...>::mkOutRange()
-    {
-	//FunctionalMultiArray<typename MapF::value_type,MapF,Ranges...> fma(mSpace, mMapf);
-	std::map<typename MapF::value_type,size_t> mult;
-	for(auto ii = mMapf.begin(); ii.max() != ii.pos(); ++ii) {
-	    mult[mMapf[ii]]++;
-	}
-	
-	std::vector<typename MapF::value_type> outmeta(mult.size());
-	std::vector<size_t> outmult(mult.size());
-
-	size_t cnt = 0;
-	for(auto& x: mult){
-	    outmeta[cnt] = x.first;
-	    outmult[cnt] = x.second;
-	    ++cnt;
-	}
-
-	ORFType orf(outmeta);
-	mOutRange = std::dynamic_pointer_cast<ORType>( orf.create() );
-	mMapMult = MultiArray<size_t,ORType>( mOutRange, outmult );
-    }
-    
-    template <class MapF, class... Ranges>
-    MapRange<MapF,Ranges...>::MapRange(const MapF& mapf, const std::shared_ptr<Ranges>&... rs) :
-	mSpace(std::make_tuple(rs...)),
-	mMapf(mapf)
-    {
-	mkOutRange();
-    }
-
-    template <class MapF, class... Ranges>
-    MapRange<MapF,Ranges...>::MapRange(const MapF& mapf, const Space& space) :
-	mSpace( space ),
-	mMapf(mapf)
-    {
-	mkOutRange();
-    }
-
-    template <class MapF, class... Ranges>
-    template <size_t N>
-    auto MapRange<MapF,Ranges...>::get() const -> decltype( *std::get<N>( mSpace ) )&
-    {
-	return *std::get<N>(mSpace);
-    }
-
-    template <class MapF, class... Ranges>
-    template <size_t N>
-    auto MapRange<MapF,Ranges...>::getPtr() const -> decltype( std::get<N>( mSpace ) )&
-    {
-	return std::get<N>(mSpace);
-    }
-
-    template <class MapF, class... Ranges>
-    auto MapRange<MapF,Ranges...>::outRange() const -> std::shared_ptr<ORType>
-    {
-	return mOutRange;
-    }
-
-    template <class MapF, class... Ranges>
-    const MapF& MapRange<MapF,Ranges...>::map() const
-    {
-	return mMapf;
-    }
-    
-    template <class MapF, class... Ranges>
-    size_t MapRange<MapF,Ranges...>::dim() const
-    {
-	return sdim;
-    }
-    
-    template <class MapF, class... Ranges>
-    size_t MapRange<MapF,Ranges...>::size() const
-    {
-	return mOutRange->size();
-	//return RPackNum<sizeof...(Ranges)-1>::getSize(mSpace);
-    }
-
-    template <class MapF, class... Ranges>
-    SpaceType MapRange<MapF,Ranges...>::spaceType() const
-    {
-	return SpaceType::ANY;
-    }
-    
-    template <class MapF, class... Ranges>
-    const typename MapRange<MapF,Ranges...>::Space& MapRange<MapF,Ranges...>::space() const
-    {
-	return mSpace;
-    }
-
-    template <class MapF, class... Ranges>
-    std::string MapRange<MapF,Ranges...>::stringMeta(size_t pos) const
-    {
-	auto i = begin();
-	i = pos;
-	return "[ " + RPackNum<sizeof...(Ranges)-1>::metaTupleToString(i.meta()) + " ]";
-    }
-
-    template <class MapF, class... Ranges>
-    std::vector<char> MapRange<MapF,Ranges...>::data() const
-    {
-	DataHeader h;
-	h.spaceType = static_cast<int>( SpaceType::ANY );
-	h.metaSize = sizeof...(Ranges);
-	h.multiple = 1;
-	std::vector<char> out;
-	//out.reserve(h.metaSize + sizeof(DataHeader));
-	char* hcp = reinterpret_cast<char*>(&h);
-	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
-	RPackNum<sizeof...(Ranges)-1>::fillRangeDataVec(out, mSpace);
-	return out;
-    }
-    
-    template <class MapF, class... Ranges>
-    typename MapRange<MapF,Ranges...>::IndexType MapRange<MapF,Ranges...>::begin() const
-    {
-	MapIndex<MapF,typename Ranges::IndexType...>
-	    i( std::dynamic_pointer_cast<MapRange<MapF,Ranges...> >
-	       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
-	i = 0;
-	return i;
-    }
-
-    template <class MapF, class... Ranges>
-    typename MapRange<MapF,Ranges...>::IndexType MapRange<MapF,Ranges...>::end() const
-    {
-	MapIndex<MapF,typename Ranges::IndexType...>
-	    i( std::dynamic_pointer_cast<MapRange<MapF,Ranges...> >
-	       ( std::shared_ptr<RangeBase>( RB::mThis )) );
-	i = size();
-	return i;
-    }
-
-    template <class MapF, class... Ranges>
-    auto MapRange<MapF,Ranges...>::mapMultiplicity() const
-	-> const MultiArray<size_t,ORType>&
-    {
-	return mMapMult;
-    }
-
-    template <class MapF, class... Ranges>
-    auto MapRange<MapF,Ranges...>::explMapMultiplicity() const
-	-> MultiArray<size_t,MapRange>
-    {
-	auto tmp = mMapMult;
-	return tmp.format( std::dynamic_pointer_cast<MapRange<MapF,Ranges...> >
-			   ( std::shared_ptr<RangeBase>( RB::mThis )) ); 
-    }
-    
-    template <class MapF, class... Ranges>
-    template <class... ERanges>
-    auto MapRange<MapF,Ranges...>::cat(const std::shared_ptr<MapRange<ERanges...> >& erange)
-	-> std::shared_ptr<MapRange<Ranges...,ERanges...> >
-    {
-	auto crange = std::tuple_cat(mSpace, erange->space());
-	MapRangeFactory<Ranges...,ERanges...> rf(crange);
-	return std::dynamic_pointer_cast<MapRange<Ranges...,ERanges...> >(rf.create());
-    }
-
-}
 
 #endif
diff --git a/src/include/multi_array.cc.h b/src/include/multi_array.cc.h
new file mode 100644
index 0000000..aa4990d
--- /dev/null
+++ b/src/include/multi_array.cc.h
@@ -0,0 +1,264 @@
+
+#include "multi_array.h"
+
+namespace MultiArrayTools
+{
+    template <typename T>
+    Scalar<T> scalar(const T& in)
+    {
+	NullRF nrf;
+	return Scalar<T>( std::dynamic_pointer_cast<NullRange>( nrf.create() ), std::vector<T>( { in } ) );
+    }
+    
+    /*******************
+     *  MultiArray     *	     
+     *******************/
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const typename CRange::Space& space) :
+	MutableMultiArrayBase<T,SRanges...>(space),
+	mCont(MAB::mRange->size())
+    {
+	MAB::mInit = true;
+    }
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const typename CRange::Space& space,
+					 std::vector<T>&& vec) :
+	MutableMultiArrayBase<T,SRanges...>(space),
+	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) :
+	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 T& val) : 
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mCont(MAB::mRange->size(), val)
+    {
+	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>
+    MultiArray<T,SRanges...>::MultiArray(MultiArray<T,AnonymousRange>&& ama, SIZET<SRanges>... sizes) :
+	MutableMultiArrayBase<T,SRanges...>
+	( ama.range()->template get<0>().template scast<SRanges...>(sizes...)->space() ),
+	mCont( std::move( ama.mCont ) )
+    {
+	MAB::mInit = true;
+    }
+    
+    /*
+    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 IndexType& i)
+    {
+	return mCont[ i.pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    const T& MultiArray<T,SRanges...>::operator[](const IndexType& i) const
+    {
+	return mCont[ i.pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    T& MultiArray<T,SRanges...>::at(const typename IndexType::MetaType& meta)
+    {
+	return mCont[ MAB::beginIndex().at(meta).pos() ];
+    }
+
+    template <typename T, class... SRanges>
+    const T& MultiArray<T,SRanges...>::at(const typename 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>
+    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > MultiArray<T,SRanges...>::anonymous(bool slice) const
+    {
+	AnonymousRangeFactory arf(MAB::mRange->space());
+	if(slice){
+	    return std::make_shared<ConstSlice<T,AnonymousRange> >
+		( std::dynamic_pointer_cast<AnonymousRange>( arf.create() ),
+		  data() );
+	}
+	else {
+	    return std::make_shared<MultiArray<T,AnonymousRange> >
+		( std::dynamic_pointer_cast<AnonymousRange>( arf.create() ),
+		  mCont );
+	}
+    }	
+
+    template <typename T, class... SRanges>
+    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > MultiArray<T,SRanges...>::anonymousMove()
+    {
+	AnonymousRangeFactory arf(MAB::mRange->space());
+	MAB::mInit = false;
+	return std::make_shared<MultiArray<T,AnonymousRange> >
+	    ( std::dynamic_pointer_cast<AnonymousRange>( arf.create() ),
+	      std::move(mCont) );
+    }	
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator=(const T& in)
+    {
+	for(auto& x: mCont){
+	    x = in;
+	}
+	return *this;
+    }
+    
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator+=(const MultiArray& in)
+    {
+	if(not MAB::mInit){ // not initialized by default constructor !!
+	    (*this) = in;
+	}
+	else {
+	    assert( PackNum<sizeof...(SRanges)-1>::checkIfSameInstance( MAB::mRange->space(), in.mRange->space() ) );
+	    std::transform(mCont.begin(), mCont.end(), in.mCont.begin(), mCont.begin(), std::plus<T>());
+	}
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator-=(const MultiArray& in)
+    {
+	if(not MAB::mInit){ // not initialized by default constructor !!
+	    (*this) = in;
+	}
+	else {
+	    assert( PackNum<sizeof...(SRanges)-1>::checkIfSameInstance( MAB::mRange->space(), in.mRange->space() ) );
+	    std::transform(mCont.begin(), mCont.end(), in.mCont.begin(), mCont.begin(), std::minus<T>());
+	}
+	return *this;
+    }
+    
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator*=(const T& in)
+    {
+	for(auto& x: mCont){
+	    x *= in;
+	}
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator/=(const T& in)
+    {
+	for(auto& x: mCont){
+	    x /= in;
+	}
+	return *this;
+    }
+    
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::operator T() const
+    {
+	static_assert( sizeof...(SRanges) == 1, "try to cast non-scalar type into scalar" );
+	// TODO: check that SIZE is statically = 1 !!!
+	return mCont[0];
+    }
+
+    template <typename T, class... SRanges>
+    auto MultiArray<T,SRanges...>::cat() const
+	-> decltype(ArrayCatter<T>::cat(*this))
+    {
+	return ArrayCatter<T>::cat(*this);
+    }
+}
+
diff --git a/src/include/multi_array.h b/src/include/multi_array.h
index 6ab76db..766f1c1 100644
--- a/src/include/multi_array.h
+++ b/src/include/multi_array.h
@@ -131,265 +131,5 @@ namespace MultiArrayTools
  * ---   TEMPLATE CODE   --- *
  * ========================= */
 
-namespace MultiArrayTools
-{
-    template <typename T>
-    Scalar<T> scalar(const T& in)
-    {
-	NullRF nrf;
-	return Scalar<T>( std::dynamic_pointer_cast<NullRange>( nrf.create() ), std::vector<T>( { in } ) );
-    }
-    
-    /*******************
-     *  MultiArray     *	     
-     *******************/
-
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::MultiArray(const typename CRange::Space& space) :
-	MutableMultiArrayBase<T,SRanges...>(space),
-	mCont(MAB::mRange->size())
-    {
-	MAB::mInit = true;
-    }
-
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::MultiArray(const typename CRange::Space& space,
-					 std::vector<T>&& vec) :
-	MutableMultiArrayBase<T,SRanges...>(space),
-	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) :
-	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 T& val) : 
-	MutableMultiArrayBase<T,SRanges...>(ranges...),
-	mCont(MAB::mRange->size(), val)
-    {
-	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>
-    MultiArray<T,SRanges...>::MultiArray(MultiArray<T,AnonymousRange>&& ama, SIZET<SRanges>... sizes) :
-	MutableMultiArrayBase<T,SRanges...>
-	( ama.range()->template get<0>().template scast<SRanges...>(sizes...)->space() ),
-	mCont( std::move( ama.mCont ) )
-    {
-	MAB::mInit = true;
-    }
-    
-    /*
-    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 IndexType& i)
-    {
-	return mCont[ i.pos() ];
-    }
-
-    template <typename T, class... SRanges>
-    const T& MultiArray<T,SRanges...>::operator[](const IndexType& i) const
-    {
-	return mCont[ i.pos() ];
-    }
-
-    template <typename T, class... SRanges>
-    T& MultiArray<T,SRanges...>::at(const typename IndexType::MetaType& meta)
-    {
-	return mCont[ MAB::beginIndex().at(meta).pos() ];
-    }
-
-    template <typename T, class... SRanges>
-    const T& MultiArray<T,SRanges...>::at(const typename 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>
-    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > MultiArray<T,SRanges...>::anonymous(bool slice) const
-    {
-	AnonymousRangeFactory arf(MAB::mRange->space());
-	if(slice){
-	    return std::make_shared<ConstSlice<T,AnonymousRange> >
-		( std::dynamic_pointer_cast<AnonymousRange>( arf.create() ),
-		  data() );
-	}
-	else {
-	    return std::make_shared<MultiArray<T,AnonymousRange> >
-		( std::dynamic_pointer_cast<AnonymousRange>( arf.create() ),
-		  mCont );
-	}
-    }	
-
-    template <typename T, class... SRanges>
-    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > MultiArray<T,SRanges...>::anonymousMove()
-    {
-	AnonymousRangeFactory arf(MAB::mRange->space());
-	MAB::mInit = false;
-	return std::make_shared<MultiArray<T,AnonymousRange> >
-	    ( std::dynamic_pointer_cast<AnonymousRange>( arf.create() ),
-	      std::move(mCont) );
-    }	
-
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator=(const T& in)
-    {
-	for(auto& x: mCont){
-	    x = in;
-	}
-	return *this;
-    }
-    
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator+=(const MultiArray& in)
-    {
-	if(not MAB::mInit){ // not initialized by default constructor !!
-	    (*this) = in;
-	}
-	else {
-	    assert( PackNum<sizeof...(SRanges)-1>::checkIfSameInstance( MAB::mRange->space(), in.mRange->space() ) );
-	    std::transform(mCont.begin(), mCont.end(), in.mCont.begin(), mCont.begin(), std::plus<T>());
-	}
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator-=(const MultiArray& in)
-    {
-	if(not MAB::mInit){ // not initialized by default constructor !!
-	    (*this) = in;
-	}
-	else {
-	    assert( PackNum<sizeof...(SRanges)-1>::checkIfSameInstance( MAB::mRange->space(), in.mRange->space() ) );
-	    std::transform(mCont.begin(), mCont.end(), in.mCont.begin(), mCont.begin(), std::minus<T>());
-	}
-	return *this;
-    }
-    
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator*=(const T& in)
-    {
-	for(auto& x: mCont){
-	    x *= in;
-	}
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>& MultiArray<T,SRanges...>::operator/=(const T& in)
-    {
-	for(auto& x: mCont){
-	    x /= in;
-	}
-	return *this;
-    }
-    
-    template <typename T, class... SRanges>
-    MultiArray<T,SRanges...>::operator T() const
-    {
-	static_assert( sizeof...(SRanges) == 1, "try to cast non-scalar type into scalar" );
-	// TODO: check that SIZE is statically = 1 !!!
-	return mCont[0];
-    }
-
-    template <typename T, class... SRanges>
-    auto MultiArray<T,SRanges...>::cat() const
-	-> decltype(ArrayCatter<T>::cat(*this))
-    {
-	return ArrayCatter<T>::cat(*this);
-    }
-}
 
 #endif
diff --git a/src/include/multi_array_base.cc.h b/src/include/multi_array_base.cc.h
new file mode 100644
index 0000000..9d9dbd9
--- /dev/null
+++ b/src/include/multi_array_base.cc.h
@@ -0,0 +1,223 @@
+
+#include "multi_array_base.h"
+
+namespace MultiArrayTools
+{
+        
+    /**********************
+     *  MultiArrayBase    *	     
+     **********************/
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::MultiArrayBase(const std::shared_ptr<SRanges>&... ranges)
+    {
+	ContainerRangeFactory<T,SRanges...> crf(ranges...);
+	mRange = std::dynamic_pointer_cast<ContainerRange<T,SRanges...> >( crf.create() );
+	mProtoI = std::make_shared<IndexType>( mRange, reinterpret_cast<std::intptr_t>(this) );
+    }
+
+    template <typename T, class... SRanges>
+    MultiArrayBase<T,SRanges...>::MultiArrayBase(const typename CRange::Space& space)
+    {
+	ContainerRangeFactory<T,SRanges...> crf(space);
+	mRange = std::dynamic_pointer_cast<ContainerRange<T,SRanges...> >( crf.create() );
+	mProtoI = std::make_shared<IndexType>( mRange, reinterpret_cast<std::intptr_t>(this) );
+    }
+
+    template <typename T, class... SRanges>
+    template <typename X>
+    const T& MultiArrayBase<T,SRanges...>::operator[](const ContainerIndex<X,typename SRanges::IndexType...>& i)
+    {
+	IndexType ii(*mProtoI);
+	ii = i;
+	return (*this)[ii];
+    }
+
+    
+    template <typename T, class... SRanges>
+    size_t MultiArrayBase<T,SRanges...>::size() const
+    {
+	return mRange->size();
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType MultiArrayBase<T,SRanges...>::begin() const
+    {
+	IndexType i(*mProtoI);
+	i = 0;
+	return i.setData(data());
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType MultiArrayBase<T,SRanges...>::end() const
+    {
+	IndexType i(*mProtoI);
+	i = i.max();
+	//i = mRange->size();
+	return i.setData(data());
+    }
+    
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::beginIndex() const
+    {
+	IndexType i(*mProtoI);
+	i = 0;
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    typename MultiArrayBase<T,SRanges...>::IndexType
+    MultiArrayBase<T,SRanges...>::endIndex() const
+    {
+	IndexType i(*mProtoI);
+	i = i.max();
+	//i = mRange->size();
+	return i.setData(data());
+    }
+
+    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()(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	return ConstOperationRoot<T,SRanges...>(*this, inds...);
+    }
+
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MultiArrayBase<T,SRanges...>::op(const std::shared_ptr<IndexType>& ind) const
+    {
+	return ConstOperationRoot<T,SRanges...>(data(), *ind);
+    }
+
+    template <typename T, class... SRanges>
+    template <class... MappedRanges>
+    ConstOperationRoot<T,MappedRanges...>
+    MultiArrayBase<T,SRanges...>::m(const std::shared_ptr<typename MappedRanges::IndexType>&... inds) const
+    {
+	static_assert(sizeof...(SRanges) == sizeof...(MappedRanges),
+		      "number of mapped ranges must be equal to number of original ranges");
+	return ConstOperationRoot<T,MappedRanges...>(*this, inds...);
+    }
+    
+    template <typename T, class... SRanges>
+    bool MultiArrayBase<T,SRanges...>::isInit() const
+    {
+	return mInit;
+    }
+    
+    template <typename T, class... SRanges>
+    template <size_t N>
+    auto MultiArrayBase<T,SRanges...>::getRangePtr() const
+	-> decltype(mRange->template getPtr<N>())
+    {
+	return mRange->template getPtr<N>();
+    }
+
+    
+    /******************************
+     *  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>
+    MutableMultiArrayBase<T,SRanges...>::MutableMultiArrayBase(const typename CRange::Space& space) :
+	MultiArrayBase<T,SRanges...>(space) {}
+    /*
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::IndexType MutableMultiArrayBase<T,SRanges...>::begin()
+    {
+	auto i = mRange->begin();
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    typename MutableMultiArrayBase<T,SRanges...>::IndexType MutableMultiArrayBase<T,SRanges...>::end()
+    {
+	auto i = mRange->end();
+	return i.setData(data());
+    }
+    */
+
+    template <typename T, class... SRanges>
+    template <typename X>
+    T& MutableMultiArrayBase<T,SRanges...>::operator[](const ContainerIndex<X,typename SRanges::IndexType...>& i)
+    {
+	IndexType ii(*MAB::mProtoI);
+	ii = i;
+	return (*this)[ii];
+    }
+
+    
+    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()(const std::shared_ptr<typename SRanges::IndexType>&... inds)
+    {
+	return OperationRoot<T,SRanges...>(*this, inds...);
+    }
+
+    template <typename T, class... SRanges>
+    OperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::op(const std::shared_ptr<IndexType>& ind)
+    {
+	return OperationRoot<T,SRanges...>(data(), *ind);
+    }
+    
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::operator()(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	return ConstOperationRoot<T,SRanges...>(*this, inds...);
+    }
+
+    template <typename T, class... SRanges>
+    ConstOperationRoot<T,SRanges...>
+    MutableMultiArrayBase<T,SRanges...>::op(const std::shared_ptr<IndexType>& ind) const
+    {
+	return ConstOperationRoot<T,SRanges...>(data(), *ind);
+    }
+
+    template <typename T, class... SRanges>
+    template <class... MappedRanges>
+    OperationRoot<T,MappedRanges...>
+    MutableMultiArrayBase<T,SRanges...>::m(const std::shared_ptr<typename MappedRanges::IndexType>&... inds)
+    {
+	static_assert(sizeof...(SRanges) == sizeof...(MappedRanges),
+		      "number of mapped ranges must be equal to number of original ranges");
+	return OperationRoot<T,MappedRanges...>(*this, inds...);
+    }
+
+    template <typename T, class... SRanges>
+    template <class... MappedRanges>
+    ConstOperationRoot<T,MappedRanges...>
+    MutableMultiArrayBase<T,SRanges...>::m(const std::shared_ptr<typename MappedRanges::IndexType>&... inds) const
+    {
+	static_assert(sizeof...(SRanges) == sizeof...(MappedRanges),
+		      "number of mapped ranges must be equal to number of original ranges");
+	return ConstOperationRoot<T,MappedRanges...>(*this, inds...);
+    }
+
+} // end namespace MultiArrayTools
+
diff --git a/src/include/multi_array_base.h b/src/include/multi_array_base.h
index a01343d..e771d6c 100644
--- a/src/include/multi_array_base.h
+++ b/src/include/multi_array_base.h
@@ -142,224 +142,5 @@ namespace MultiArrayTools
  * ---   TEMPLATE CODE   --- *
  * ========================= */
 
-namespace MultiArrayTools
-{
-        
-    /**********************
-     *  MultiArrayBase    *	     
-     **********************/
-
-    template <typename T, class... SRanges>
-    MultiArrayBase<T,SRanges...>::MultiArrayBase(const std::shared_ptr<SRanges>&... ranges)
-    {
-	ContainerRangeFactory<T,SRanges...> crf(ranges...);
-	mRange = std::dynamic_pointer_cast<ContainerRange<T,SRanges...> >( crf.create() );
-	mProtoI = std::make_shared<IndexType>( mRange, reinterpret_cast<std::intptr_t>(this) );
-    }
-
-    template <typename T, class... SRanges>
-    MultiArrayBase<T,SRanges...>::MultiArrayBase(const typename CRange::Space& space)
-    {
-	ContainerRangeFactory<T,SRanges...> crf(space);
-	mRange = std::dynamic_pointer_cast<ContainerRange<T,SRanges...> >( crf.create() );
-	mProtoI = std::make_shared<IndexType>( mRange, reinterpret_cast<std::intptr_t>(this) );
-    }
-
-    template <typename T, class... SRanges>
-    template <typename X>
-    const T& MultiArrayBase<T,SRanges...>::operator[](const ContainerIndex<X,typename SRanges::IndexType...>& i)
-    {
-	IndexType ii(*mProtoI);
-	ii = i;
-	return (*this)[ii];
-    }
-
-    
-    template <typename T, class... SRanges>
-    size_t MultiArrayBase<T,SRanges...>::size() const
-    {
-	return mRange->size();
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::IndexType MultiArrayBase<T,SRanges...>::begin() const
-    {
-	IndexType i(*mProtoI);
-	i = 0;
-	return i.setData(data());
-    }
-    
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::IndexType MultiArrayBase<T,SRanges...>::end() const
-    {
-	IndexType i(*mProtoI);
-	i = i.max();
-	//i = mRange->size();
-	return i.setData(data());
-    }
-    
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::IndexType
-    MultiArrayBase<T,SRanges...>::beginIndex() const
-    {
-	IndexType i(*mProtoI);
-	i = 0;
-	return i.setData(data());
-    }
-
-    template <typename T, class... SRanges>
-    typename MultiArrayBase<T,SRanges...>::IndexType
-    MultiArrayBase<T,SRanges...>::endIndex() const
-    {
-	IndexType i(*mProtoI);
-	i = i.max();
-	//i = mRange->size();
-	return i.setData(data());
-    }
-
-    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()(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
-    {
-	return ConstOperationRoot<T,SRanges...>(*this, inds...);
-    }
-
-    template <typename T, class... SRanges>
-    ConstOperationRoot<T,SRanges...>
-    MultiArrayBase<T,SRanges...>::op(const std::shared_ptr<IndexType>& ind) const
-    {
-	return ConstOperationRoot<T,SRanges...>(data(), *ind);
-    }
-
-    template <typename T, class... SRanges>
-    template <class... MappedRanges>
-    ConstOperationRoot<T,MappedRanges...>
-    MultiArrayBase<T,SRanges...>::m(const std::shared_ptr<typename MappedRanges::IndexType>&... inds) const
-    {
-	static_assert(sizeof...(SRanges) == sizeof...(MappedRanges),
-		      "number of mapped ranges must be equal to number of original ranges");
-	return ConstOperationRoot<T,MappedRanges...>(*this, inds...);
-    }
-    
-    template <typename T, class... SRanges>
-    bool MultiArrayBase<T,SRanges...>::isInit() const
-    {
-	return mInit;
-    }
-    
-    template <typename T, class... SRanges>
-    template <size_t N>
-    auto MultiArrayBase<T,SRanges...>::getRangePtr() const
-	-> decltype(mRange->template getPtr<N>())
-    {
-	return mRange->template getPtr<N>();
-    }
-
-    
-    /******************************
-     *  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>
-    MutableMultiArrayBase<T,SRanges...>::MutableMultiArrayBase(const typename CRange::Space& space) :
-	MultiArrayBase<T,SRanges...>(space) {}
-    /*
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::IndexType MutableMultiArrayBase<T,SRanges...>::begin()
-    {
-	auto i = mRange->begin();
-	return i.setData(data());
-    }
-
-    template <typename T, class... SRanges>
-    typename MutableMultiArrayBase<T,SRanges...>::IndexType MutableMultiArrayBase<T,SRanges...>::end()
-    {
-	auto i = mRange->end();
-	return i.setData(data());
-    }
-    */
-
-    template <typename T, class... SRanges>
-    template <typename X>
-    T& MutableMultiArrayBase<T,SRanges...>::operator[](const ContainerIndex<X,typename SRanges::IndexType...>& i)
-    {
-	IndexType ii(*MAB::mProtoI);
-	ii = i;
-	return (*this)[ii];
-    }
-
-    
-    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()(const std::shared_ptr<typename SRanges::IndexType>&... inds)
-    {
-	return OperationRoot<T,SRanges...>(*this, inds...);
-    }
-
-    template <typename T, class... SRanges>
-    OperationRoot<T,SRanges...>
-    MutableMultiArrayBase<T,SRanges...>::op(const std::shared_ptr<IndexType>& ind)
-    {
-	return OperationRoot<T,SRanges...>(data(), *ind);
-    }
-    
-    template <typename T, class... SRanges>
-    ConstOperationRoot<T,SRanges...>
-    MutableMultiArrayBase<T,SRanges...>::operator()(const std::shared_ptr<typename SRanges::IndexType>&... inds) const
-    {
-	return ConstOperationRoot<T,SRanges...>(*this, inds...);
-    }
-
-    template <typename T, class... SRanges>
-    ConstOperationRoot<T,SRanges...>
-    MutableMultiArrayBase<T,SRanges...>::op(const std::shared_ptr<IndexType>& ind) const
-    {
-	return ConstOperationRoot<T,SRanges...>(data(), *ind);
-    }
-
-    template <typename T, class... SRanges>
-    template <class... MappedRanges>
-    OperationRoot<T,MappedRanges...>
-    MutableMultiArrayBase<T,SRanges...>::m(const std::shared_ptr<typename MappedRanges::IndexType>&... inds)
-    {
-	static_assert(sizeof...(SRanges) == sizeof...(MappedRanges),
-		      "number of mapped ranges must be equal to number of original ranges");
-	return OperationRoot<T,MappedRanges...>(*this, inds...);
-    }
-
-    template <typename T, class... SRanges>
-    template <class... MappedRanges>
-    ConstOperationRoot<T,MappedRanges...>
-    MutableMultiArrayBase<T,SRanges...>::m(const std::shared_ptr<typename MappedRanges::IndexType>&... inds) const
-    {
-	static_assert(sizeof...(SRanges) == sizeof...(MappedRanges),
-		      "number of mapped ranges must be equal to number of original ranges");
-	return ConstOperationRoot<T,MappedRanges...>(*this, inds...);
-    }
-
-} // end namespace MultiArrayTools
     
 #endif
diff --git a/src/include/multi_array_header.cc.h b/src/include/multi_array_header.cc.h
new file mode 100644
index 0000000..dfca6df
--- /dev/null
+++ b/src/include/multi_array_header.cc.h
@@ -0,0 +1,10 @@
+
+#include "multi_array_operation.cc.h"
+#include "functional_multi_array.cc.h"
+#include "helper_tools.cc.h"
+#include "map_range.cc.h"
+#include "multi_array_base.cc.h"
+#include "multi_array.cc.h"
+#include "slice.cc.h"
+
+
diff --git a/src/include/multi_array_header.h b/src/include/multi_array_header.h
index 61effc1..c78406f 100644
--- a/src/include/multi_array_header.h
+++ b/src/include/multi_array_header.h
@@ -18,6 +18,7 @@
 #include "helper_tools.h"
 #include "operation_def.h"
 #include "map_range.h"
+#include "expressions.h"
 //#include "slice.h"
 //#include "manipulator.h"
 //#include "range_transformer.h"
@@ -61,4 +62,6 @@ namespace MultiArrayTools
 
 }
 
+#include "multi_array_header.cc.h"
+
 #endif
diff --git a/src/include/multi_array_operation.cc.h b/src/include/multi_array_operation.cc.h
new file mode 100644
index 0000000..769344a
--- /dev/null
+++ b/src/include/multi_array_operation.cc.h
@@ -0,0 +1,548 @@
+
+#include "multi_array_operation.h"
+
+/* ========================= *
+ * ---   TEMPLATE CODE   --- *
+ * ========================= */
+
+namespace MultiArrayTools
+{
+    namespace
+    {
+	using namespace MultiArrayHelper;
+    }
+    
+    /***************************
+     *   OperationTemplate     *
+     ***************************/
+
+    
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationBase<T,OperationClass>::operator+(const OperationBase<T,Second>& in) const
+	-> Operation<T,plus<T>,OperationClass,Second>
+    {
+    	return Operation<T,plus<T>,OperationClass,Second>(THIS(), in.THIS());
+    }
+
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationBase<T,OperationClass>::operator-(const OperationBase<T,Second>& in) const
+	-> Operation<T,minus<T>,OperationClass,Second>
+    {
+    	return Operation<T,minus<T>,OperationClass,Second>(THIS(), in.THIS());
+    }
+    
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationBase<T,OperationClass>::operator*(const OperationBase<T,Second>& in) const
+	-> Operation<T,multiplies<T>,OperationClass,Second>
+    {
+    	return Operation<T,multiplies<T>,OperationClass,Second>(THIS(), in.THIS());
+    }
+
+    template <typename T, class OperationClass>
+    template <class Second>
+    auto OperationBase<T,OperationClass>::operator/(const OperationBase<T,Second>& in) const
+	-> Operation<T,divides<T>,OperationClass,Second>
+    {
+    	return Operation<T,divides<T>,OperationClass,Second>(THIS(), in.THIS());
+    }
+
+    template <typename T, class OperationClass>
+    template <class IndexType>
+    auto OperationBase<T,OperationClass>::c(const std::shared_ptr<IndexType>& ind) const
+	-> Contraction<T,OperationClass,IndexType>
+    {
+	return Contraction<T,OperationClass,IndexType>(THIS(), ind);
+    }
+
+    template <typename T, class OperationClass>
+    template <class... Indices>
+    auto OperationBase<T,OperationClass>::sl(const std::shared_ptr<Indices>&... inds) const
+	-> ConstSlice<T,typename Indices::RangeType...>
+    {
+	ConstSlice<T,typename Indices::RangeType...> out(inds->range()...);
+	out.define(inds...) = THIS();
+	return out;
+    }
+
+    template <typename T, class OperationClass>
+    template <class... Indices>
+    auto OperationBase<T,OperationClass>::slc(const std::shared_ptr<Indices>&... inds) const
+	-> SliceContraction<T,OperationClass,Indices...>
+    {
+	return SliceContraction<T,OperationClass,Indices...>
+	    (THIS(), inds...);
+    }
+
+    /*****************************************
+     *   OperationMaster::AssignmentExpr     *
+     *****************************************/
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
+    AssignmentExpr(OperationMaster& m, const OpClass& sec) :
+	mM(m), mSec(sec) {}
+
+    template <typename T, class OpClass, class... Ranges>
+    inline void OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
+    operator()(size_t start, ExtType last) const
+    {
+	//VCHECK(mSec.template get<ExtType>(last));
+	mM.add(start, mSec.template get<ExtType>(last) );
+    }
+    
+    template <typename T, class OpClass, class... Ranges>
+    typename OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::ExtType
+    OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
+    rootSteps(std::intptr_t iPtrNum) const
+    {
+	return mSec.rootSteps(iPtrNum);
+    }
+    
+    
+    /*************************
+     *   OperationMaster     *
+     *************************/
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::
+    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+		    IndexType& index) :
+	mSecond(second), mDataPtr(ma.data()),
+	mIndex(index)
+    {
+	performAssignment(0);
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::
+    OperationMaster(T* data, const OpClass& second,
+		    IndexType& index) :
+	mSecond(second), mDataPtr(data),
+	mIndex(index)
+    {
+	performAssignment(0);
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    void OperationMaster<T,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
+    {
+	AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop
+	const auto loop = mSecond.template loop<decltype(mIndex.ifor(1,ae))>( mIndex.ifor(1,ae) );
+	// hidden Loops outside ! -> auto vectorizable
+	loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    inline T OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
+    {
+	return mDataPtr[pos];
+    }
+
+   
+    /****************************
+     *   ConstOperationRoot     *
+     ****************************/
+
+    template <typename T, class... Ranges>
+    ConstOperationRoot<T,Ranges...>::
+    ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
+		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	mDataPtr(ma.data()),
+	mIndex( ma.begin() )
+    {
+	//VCHECK(ma.data());
+	mIndex(indices...);
+    }
+
+    template <typename T, class... Ranges>
+    ConstOperationRoot<T,Ranges...>::
+    ConstOperationRoot(std::shared_ptr<MultiArrayBase<T,Ranges...> > maptr,
+		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	mDataPtr(maptr->data()),
+	mIndex(maptr->begin()),
+	mMaPtr(maptr)
+    {
+	mIndex(indices...);
+    }
+
+    template <typename T, class... Ranges>
+    ConstOperationRoot<T,Ranges...>::
+    ConstOperationRoot(const T* data, const IndexType& ind) :
+	mDataPtr(data),
+	mIndex( ind ) { }
+
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline T ConstOperationRoot<T,Ranges...>::get(ET pos) const
+    {
+	return mDataPtr[pos.val()+mOff];
+    }
+
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline const ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::set(ET pos) const
+    {
+	mIndex = pos.val();
+	mOff = mIndex.pos();
+	return *this;
+    }
+
+    template <typename T, class... Ranges>
+    const T* ConstOperationRoot<T,Ranges...>::data() const
+    {
+	return mDataPtr + mIndex().pos();
+    }
+    
+    template <typename T, class... Ranges>
+    MExt<void> ConstOperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
+    {
+	return MExt<void>(getStepSize( mIndex, iPtrNum ));
+	//return MExt<void>(getStepSize( getRootIndices( mIndex->info() ), iPtrNum ));
+    }
+
+
+    template <typename T, class... Ranges>
+    template <class Expr>
+    Expr ConstOperationRoot<T,Ranges...>::loop(Expr exp) const
+    {
+	return exp;
+    }
+
+    /********************
+     *   StaticCast     *
+     ********************/
+
+    template <typename T, class Op>
+    StaticCast<T,Op>::StaticCast(const Op& op) : mOp(op) {}
+
+    template <typename T, class Op>
+    template <class ET>
+    inline T StaticCast<T,Op>::get(ET pos) const
+    {
+	return static_cast<T>( mOp.get(pos) );
+    }
+
+    template <typename T, class Op>
+    template <class ET>
+    inline const StaticCast<T,Op>& StaticCast<T,Op>::set(ET pos) const
+    {
+	mOp.set(pos);
+	return *this;
+    }
+    
+    template <typename T, class Op>
+    auto StaticCast<T,Op>::rootSteps(std::intptr_t iPtrNum) const
+	-> decltype(mOp.rootSteps(iPtrNum))
+    {
+	return mOp.rootSteps(iPtrNum);
+    }
+
+    template <typename T, class Op>
+    template <class Expr>
+    Expr StaticCast<T,Op>::loop(Expr exp) const
+    {
+	return mOp.loop(exp);
+    }
+    
+    
+    /****************************
+     *   MetaOperationRoot     *
+     ****************************/
+
+    template <class... Ranges>
+    MetaOperationRoot<Ranges...>::
+    MetaOperationRoot(const IndexType& ind) :
+	mIndex( ind ) { }
+
+    template <class... Ranges>
+    template <class ET>
+    inline typename MetaOperationRoot<Ranges...>::value_type
+    MetaOperationRoot<Ranges...>::get(ET pos) const
+    {
+	//VCHECK(pos.val());
+	//VCHECK(mDataPtr);
+	//VCHECK(mDataPtr[pos.val()])
+	return mIndex.meta(pos.val());
+    }
+
+    template <class... Ranges>
+    template <class ET>
+    inline const MetaOperationRoot<Ranges...>& MetaOperationRoot<Ranges...>::set(ET pos) const
+    {
+	mIndex = pos.val();
+	return *this;
+    }
+
+    template <class... Ranges>
+    MExt<void> MetaOperationRoot<Ranges...>::rootSteps(std::intptr_t iPtrNum) const
+    {
+	return MExt<void>(getStepSize( mIndex, iPtrNum ));
+	//return MExt<void>(getStepSize( getRootIndices( mIndex->info() ), iPtrNum ));
+    }
+
+
+    template <class... Ranges>
+    template <class Expr>
+    Expr MetaOperationRoot<Ranges...>::loop(Expr exp) const
+    {
+	return exp;
+    }
+
+    /***********************
+     *   OperationRoot     *
+     ***********************/
+    
+    template <typename T, class... Ranges>
+    OperationRoot<T,Ranges...>::
+    OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
+		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
+	mDataPtr(ma.data()),
+	mIndex( ma.begin() )
+    {
+	mIndex(indices...);
+    }
+
+    template <typename T, class... Ranges>
+    OperationRoot<T,Ranges...>::
+    OperationRoot(T* data, const IndexType& ind) :
+	mDataPtr(data),
+	mIndex( ind ) { }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
+    {
+	return OperationMaster<T,OpClass,Ranges...>(mDataPtr, in, mIndex);
+    }
+
+    template <typename T, class... Ranges>
+    OperationMaster<T,OperationRoot<T,Ranges...>,Ranges...>
+    OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
+    {
+	return operator=<OperationRoot<T,Ranges...> >(in);
+    }
+    
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline T OperationRoot<T,Ranges...>::get(ET pos) const
+    {
+	return mDataPtr[pos.val()+mOff];
+    }
+
+    template <typename T, class... Ranges>
+    template <class ET>
+    inline const OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::set(ET pos) const
+    {
+	mIndex = pos.val();
+	mOff = mIndex.pos();
+	return *this;
+    }
+
+    template <typename T, class... Ranges>
+    MExt<void> OperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
+    {
+	return MExt<void>(getStepSize( mIndex, iPtrNum ));
+	//return MExt<void>(getStepSize( mIndex.info(), iPtrNum ));
+    }
+
+    template <typename T, class... Ranges>
+    template <class Expr>
+    Expr OperationRoot<T,Ranges...>::loop(Expr exp) const
+    {
+	return exp;
+    }
+
+    template <typename T, class... Ranges>
+    T* OperationRoot<T,Ranges...>::data() const
+    {
+	return mDataPtr + mIndex().pos();
+    }
+
+    template <typename T, class... Ranges>
+    template <class... Indices>
+    auto OperationRoot<T,Ranges...>::sl(const std::shared_ptr<Indices>&... inds)
+	-> Slice<T,typename Indices::RangeType...>
+    {
+	Slice<T,typename Indices::RangeType...> out(inds->range()...);
+	out.define(inds...) = *this;
+	return out;
+    }
+
+
+    /************************
+     *   OperationValue     *
+     ************************/
+
+    template <typename T>
+    OperationValue<T>::OperationValue(const T& val) : mVal(val) {}
+
+    template <typename T>
+    template <class ET>
+    inline T OperationValue<T>::get(ET pos) const
+    {
+	return mVal;
+    }
+
+    template <typename T>
+    template <class ET>
+    inline const OperationValue<T>& OperationValue<T>::set(ET pos) const
+    {
+	return *this;
+    }
+
+    template <typename T>
+    MExt<void> OperationValue<T>::rootSteps(std::intptr_t iPtrNum) const
+    {
+	return MExt<void>(0);
+    }
+
+    template <typename T>
+    template <class Expr>
+    Expr OperationValue<T>::loop(Expr exp) const
+    {
+	return exp;
+    }
+
+    
+    /*******************
+     *   Operation     *
+     *******************/
+    
+    template <typename T, class OpFunction, class... Ops>
+    Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
+	mOps(ops...)
+    {
+	static_assert( FISSTATIC, "need function instance for non-static function" );
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    Operation<T,OpFunction,Ops...>::Operation(std::shared_ptr<OpFunction> ff,
+					      const Ops&... ops) :
+	mOps(ops...),
+	mF(ff)
+    {
+	static_assert( not FISSTATIC, "using instance of function supposed to be static" );
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    template <class ET>
+    inline T Operation<T,OpFunction,Ops...>::get(ET pos) const
+    {
+	typedef std::tuple<Ops...> OpTuple;
+	return PackNum<sizeof...(Ops)-1>::
+	    template mkOpExpr<SIZE,T,ET,OpTuple,OpFunction>(mF, pos, mOps);
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    template <class ET>
+    inline const Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::set(ET pos) const
+    {
+	PackNum<sizeof...(Ops)-1>::setOpPos(mOps,pos);
+	return *this;
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    auto Operation<T,OpFunction,Ops...>::rootSteps(std::intptr_t iPtrNum) const
+	-> decltype(PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps))
+    {
+	return PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps);
+    }
+
+    template <typename T, class OpFunction, class... Ops>
+    template <class Expr>
+    auto Operation<T,OpFunction,Ops...>::loop(Expr exp) const
+	-> decltype(PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp ))
+    {
+	return PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp );
+    }
+
+    
+    /*********************
+     *   Contraction     *
+     *********************/
+
+    template <typename T, class Op, class IndexType>
+    Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
+	mOp(op),
+	mInd(ind) {}
+
+    // forward loop !!!!
+    template <typename T, class Op, class IndexType>
+    template <class ET>
+    inline T Contraction<T,Op,IndexType>::get(ET pos) const
+    {
+	return mOp.template get<ET>(pos);
+    }
+
+    template <typename T, class Op, class IndexType>
+    template <class ET>
+    inline const Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::set(ET pos) const
+    {
+	mOp.set(pos);
+	return *this;
+    }
+
+    template <typename T, class Op, class IndexType>
+    auto Contraction<T,Op,IndexType>::rootSteps(std::intptr_t iPtrNum) const
+	-> decltype(mOp.rootSteps(iPtrNum))
+    {
+	return mOp.rootSteps(iPtrNum);
+    }
+
+    template <typename T, class Op, class IndexType>
+    template <class Expr>
+    auto Contraction<T,Op,IndexType>::loop(Expr exp) const -> decltype(mInd->iforh(1,exp))
+    {
+	return mInd->iforh(1,exp);
+    }
+
+    /**************************
+     *   SliceContraction     *
+     **************************/
+
+    template <typename T, class Op, class... Indices>
+    SliceContraction<T,Op,Indices...>::SliceContraction(const Op& op,
+							std::shared_ptr<Indices>... ind) :
+	mOp(op),
+	mCont(std::make_shared<MultiArray<T,typename Indices::RangeType...> >(ind->range()...)),
+	mTarOp(*mCont,ind...)
+    { }
+
+    // forward loop !!!!
+    template <typename T, class Op, class... Indices>
+    template <class ET>
+    inline const MultiArray<T,typename Indices::RangeType...>&
+    SliceContraction<T,Op,Indices...>::get(ET pos) const
+    {
+	*mCont = 0;
+	mOp.set(pos);
+	mTarOp = mOp;
+	return *mCont;
+    }
+
+    template <typename T, class Op, class... Indices>
+    template <class ET>
+    inline const SliceContraction<T,Op,Indices...>& SliceContraction<T,Op,Indices...>::set(ET pos) const
+    {
+	mOp.set(pos);
+	return *this;
+    }
+
+    template <typename T, class Op, class... Indices>
+    auto SliceContraction<T,Op,Indices...>::rootSteps(std::intptr_t iPtrNum) const
+	-> decltype(mOp.rootSteps(iPtrNum))
+    {
+	return mOp.rootSteps(iPtrNum);
+    }
+
+    template <typename T, class Op, class... Indices>
+    template <class Expr>
+    auto SliceContraction<T,Op,Indices...>::loop(Expr exp) const -> decltype(mOp.loop(exp))
+    {
+	return mOp.loop(exp);
+    }
+
+}
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index 7d1b602..b6aaaa4 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -505,551 +505,6 @@ namespace MultiArrayTools
     };	
 }
 
-/* ========================= *
- * ---   TEMPLATE CODE   --- *
- * ========================= */
-
-namespace MultiArrayTools
-{
-    namespace
-    {
-	using namespace MultiArrayHelper;
-    }
-    
-    /***************************
-     *   OperationTemplate     *
-     ***************************/
-
-    
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationBase<T,OperationClass>::operator+(const OperationBase<T,Second>& in) const
-	-> Operation<T,plus<T>,OperationClass,Second>
-    {
-    	return Operation<T,plus<T>,OperationClass,Second>(THIS(), in.THIS());
-    }
-
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationBase<T,OperationClass>::operator-(const OperationBase<T,Second>& in) const
-	-> Operation<T,minus<T>,OperationClass,Second>
-    {
-    	return Operation<T,minus<T>,OperationClass,Second>(THIS(), in.THIS());
-    }
-    
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationBase<T,OperationClass>::operator*(const OperationBase<T,Second>& in) const
-	-> Operation<T,multiplies<T>,OperationClass,Second>
-    {
-    	return Operation<T,multiplies<T>,OperationClass,Second>(THIS(), in.THIS());
-    }
-
-    template <typename T, class OperationClass>
-    template <class Second>
-    auto OperationBase<T,OperationClass>::operator/(const OperationBase<T,Second>& in) const
-	-> Operation<T,divides<T>,OperationClass,Second>
-    {
-    	return Operation<T,divides<T>,OperationClass,Second>(THIS(), in.THIS());
-    }
-
-    template <typename T, class OperationClass>
-    template <class IndexType>
-    auto OperationBase<T,OperationClass>::c(const std::shared_ptr<IndexType>& ind) const
-	-> Contraction<T,OperationClass,IndexType>
-    {
-	return Contraction<T,OperationClass,IndexType>(THIS(), ind);
-    }
-
-    template <typename T, class OperationClass>
-    template <class... Indices>
-    auto OperationBase<T,OperationClass>::sl(const std::shared_ptr<Indices>&... inds) const
-	-> ConstSlice<T,typename Indices::RangeType...>
-    {
-	ConstSlice<T,typename Indices::RangeType...> out(inds->range()...);
-	out.define(inds...) = THIS();
-	return out;
-    }
-
-    template <typename T, class OperationClass>
-    template <class... Indices>
-    auto OperationBase<T,OperationClass>::slc(const std::shared_ptr<Indices>&... inds) const
-	-> SliceContraction<T,OperationClass,Indices...>
-    {
-	return SliceContraction<T,OperationClass,Indices...>
-	    (THIS(), inds...);
-    }
-
-    /*****************************************
-     *   OperationMaster::AssignmentExpr     *
-     *****************************************/
-
-    template <typename T, class OpClass, class... Ranges>
-    OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
-    AssignmentExpr(OperationMaster& m, const OpClass& sec) :
-	mM(m), mSec(sec) {}
-
-    template <typename T, class OpClass, class... Ranges>
-    inline void OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
-    operator()(size_t start, ExtType last) const
-    {
-	//VCHECK(mSec.template get<ExtType>(last));
-	mM.add(start, mSec.template get<ExtType>(last) );
-    }
-    
-    template <typename T, class OpClass, class... Ranges>
-    typename OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::ExtType
-    OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
-    rootSteps(std::intptr_t iPtrNum) const
-    {
-	return mSec.rootSteps(iPtrNum);
-    }
-    
-    
-    /*************************
-     *   OperationMaster     *
-     *************************/
-
-    template <typename T, class OpClass, class... Ranges>
-    OperationMaster<T,OpClass,Ranges...>::
-    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
-		    IndexType& index) :
-	mSecond(second), mDataPtr(ma.data()),
-	mIndex(index)
-    {
-	performAssignment(0);
-    }
-
-    template <typename T, class OpClass, class... Ranges>
-    OperationMaster<T,OpClass,Ranges...>::
-    OperationMaster(T* data, const OpClass& second,
-		    IndexType& index) :
-	mSecond(second), mDataPtr(data),
-	mIndex(index)
-    {
-	performAssignment(0);
-    }
-
-    template <typename T, class OpClass, class... Ranges>
-    void OperationMaster<T,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
-    {
-	AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop
-	const auto loop = mSecond.template loop<decltype(mIndex.ifor(1,ae))>( mIndex.ifor(1,ae) );
-	// hidden Loops outside ! -> auto vectorizable
-	loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions
-    }
-
-    template <typename T, class OpClass, class... Ranges>
-    inline T OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
-    {
-	return mDataPtr[pos];
-    }
-
-   
-    /****************************
-     *   ConstOperationRoot     *
-     ****************************/
-
-    template <typename T, class... Ranges>
-    ConstOperationRoot<T,Ranges...>::
-    ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
-		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	mDataPtr(ma.data()),
-	mIndex( ma.begin() )
-    {
-	//VCHECK(ma.data());
-	mIndex(indices...);
-    }
-
-    template <typename T, class... Ranges>
-    ConstOperationRoot<T,Ranges...>::
-    ConstOperationRoot(std::shared_ptr<MultiArrayBase<T,Ranges...> > maptr,
-		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	mDataPtr(maptr->data()),
-	mIndex(maptr->begin()),
-	mMaPtr(maptr)
-    {
-	mIndex(indices...);
-    }
-
-    template <typename T, class... Ranges>
-    ConstOperationRoot<T,Ranges...>::
-    ConstOperationRoot(const T* data, const IndexType& ind) :
-	mDataPtr(data),
-	mIndex( ind ) { }
-
-    template <typename T, class... Ranges>
-    template <class ET>
-    inline T ConstOperationRoot<T,Ranges...>::get(ET pos) const
-    {
-	return mDataPtr[pos.val()+mOff];
-    }
-
-    template <typename T, class... Ranges>
-    template <class ET>
-    inline const ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::set(ET pos) const
-    {
-	mIndex = pos.val();
-	mOff = mIndex.pos();
-	return *this;
-    }
-
-    template <typename T, class... Ranges>
-    const T* ConstOperationRoot<T,Ranges...>::data() const
-    {
-	return mDataPtr + mIndex().pos();
-    }
-    
-    template <typename T, class... Ranges>
-    MExt<void> ConstOperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
-    {
-	return MExt<void>(getStepSize( mIndex, iPtrNum ));
-	//return MExt<void>(getStepSize( getRootIndices( mIndex->info() ), iPtrNum ));
-    }
-
-
-    template <typename T, class... Ranges>
-    template <class Expr>
-    Expr ConstOperationRoot<T,Ranges...>::loop(Expr exp) const
-    {
-	return exp;
-    }
-
-    /********************
-     *   StaticCast     *
-     ********************/
-
-    template <typename T, class Op>
-    StaticCast<T,Op>::StaticCast(const Op& op) : mOp(op) {}
-
-    template <typename T, class Op>
-    template <class ET>
-    inline T StaticCast<T,Op>::get(ET pos) const
-    {
-	return static_cast<T>( mOp.get(pos) );
-    }
-
-    template <typename T, class Op>
-    template <class ET>
-    inline const StaticCast<T,Op>& StaticCast<T,Op>::set(ET pos) const
-    {
-	mOp.set(pos);
-	return *this;
-    }
-    
-    template <typename T, class Op>
-    auto StaticCast<T,Op>::rootSteps(std::intptr_t iPtrNum) const
-	-> decltype(mOp.rootSteps(iPtrNum))
-    {
-	return mOp.rootSteps(iPtrNum);
-    }
-
-    template <typename T, class Op>
-    template <class Expr>
-    Expr StaticCast<T,Op>::loop(Expr exp) const
-    {
-	return mOp.loop(exp);
-    }
-    
-    
-    /****************************
-     *   MetaOperationRoot     *
-     ****************************/
-
-    template <class... Ranges>
-    MetaOperationRoot<Ranges...>::
-    MetaOperationRoot(const IndexType& ind) :
-	mIndex( ind ) { }
-
-    template <class... Ranges>
-    template <class ET>
-    inline typename MetaOperationRoot<Ranges...>::value_type
-    MetaOperationRoot<Ranges...>::get(ET pos) const
-    {
-	//VCHECK(pos.val());
-	//VCHECK(mDataPtr);
-	//VCHECK(mDataPtr[pos.val()])
-	return mIndex.meta(pos.val());
-    }
-
-    template <class... Ranges>
-    template <class ET>
-    inline const MetaOperationRoot<Ranges...>& MetaOperationRoot<Ranges...>::set(ET pos) const
-    {
-	mIndex = pos.val();
-	return *this;
-    }
-
-    template <class... Ranges>
-    MExt<void> MetaOperationRoot<Ranges...>::rootSteps(std::intptr_t iPtrNum) const
-    {
-	return MExt<void>(getStepSize( mIndex, iPtrNum ));
-	//return MExt<void>(getStepSize( getRootIndices( mIndex->info() ), iPtrNum ));
-    }
-
-
-    template <class... Ranges>
-    template <class Expr>
-    Expr MetaOperationRoot<Ranges...>::loop(Expr exp) const
-    {
-	return exp;
-    }
-
-    /***********************
-     *   OperationRoot     *
-     ***********************/
-    
-    template <typename T, class... Ranges>
-    OperationRoot<T,Ranges...>::
-    OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
-		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
-	mDataPtr(ma.data()),
-	mIndex( ma.begin() )
-    {
-	mIndex(indices...);
-    }
-
-    template <typename T, class... Ranges>
-    OperationRoot<T,Ranges...>::
-    OperationRoot(T* data, const IndexType& ind) :
-	mDataPtr(data),
-	mIndex( ind ) { }
-
-    template <typename T, class... Ranges>
-    template <class OpClass>
-    OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
-    {
-	return OperationMaster<T,OpClass,Ranges...>(mDataPtr, in, mIndex);
-    }
-
-    template <typename T, class... Ranges>
-    OperationMaster<T,OperationRoot<T,Ranges...>,Ranges...>
-    OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
-    {
-	return operator=<OperationRoot<T,Ranges...> >(in);
-    }
-    
-    template <typename T, class... Ranges>
-    template <class ET>
-    inline T OperationRoot<T,Ranges...>::get(ET pos) const
-    {
-	return mDataPtr[pos.val()+mOff];
-    }
-
-    template <typename T, class... Ranges>
-    template <class ET>
-    inline const OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::set(ET pos) const
-    {
-	mIndex = pos.val();
-	mOff = mIndex.pos();
-	return *this;
-    }
-
-    template <typename T, class... Ranges>
-    MExt<void> OperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
-    {
-	return MExt<void>(getStepSize( mIndex, iPtrNum ));
-	//return MExt<void>(getStepSize( mIndex.info(), iPtrNum ));
-    }
-
-    template <typename T, class... Ranges>
-    template <class Expr>
-    Expr OperationRoot<T,Ranges...>::loop(Expr exp) const
-    {
-	return exp;
-    }
-
-    template <typename T, class... Ranges>
-    T* OperationRoot<T,Ranges...>::data() const
-    {
-	return mDataPtr + mIndex().pos();
-    }
-
-    template <typename T, class... Ranges>
-    template <class... Indices>
-    auto OperationRoot<T,Ranges...>::sl(const std::shared_ptr<Indices>&... inds)
-	-> Slice<T,typename Indices::RangeType...>
-    {
-	Slice<T,typename Indices::RangeType...> out(inds->range()...);
-	out.define(inds...) = *this;
-	return out;
-    }
-
-
-    /************************
-     *   OperationValue     *
-     ************************/
-
-    template <typename T>
-    OperationValue<T>::OperationValue(const T& val) : mVal(val) {}
-
-    template <typename T>
-    template <class ET>
-    inline T OperationValue<T>::get(ET pos) const
-    {
-	return mVal;
-    }
-
-    template <typename T>
-    template <class ET>
-    inline const OperationValue<T>& OperationValue<T>::set(ET pos) const
-    {
-	return *this;
-    }
-
-    template <typename T>
-    MExt<void> OperationValue<T>::rootSteps(std::intptr_t iPtrNum) const
-    {
-	return MExt<void>(0);
-    }
-
-    template <typename T>
-    template <class Expr>
-    Expr OperationValue<T>::loop(Expr exp) const
-    {
-	return exp;
-    }
-
-    
-    /*******************
-     *   Operation     *
-     *******************/
-    
-    template <typename T, class OpFunction, class... Ops>
-    Operation<T,OpFunction,Ops...>::Operation(const Ops&... ops) :
-	mOps(ops...)
-    {
-	static_assert( FISSTATIC, "need function instance for non-static function" );
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    Operation<T,OpFunction,Ops...>::Operation(std::shared_ptr<OpFunction> ff,
-					      const Ops&... ops) :
-	mOps(ops...),
-	mF(ff)
-    {
-	static_assert( not FISSTATIC, "using instance of function supposed to be static" );
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    template <class ET>
-    inline T Operation<T,OpFunction,Ops...>::get(ET pos) const
-    {
-	typedef std::tuple<Ops...> OpTuple;
-	return PackNum<sizeof...(Ops)-1>::
-	    template mkOpExpr<SIZE,T,ET,OpTuple,OpFunction>(mF, pos, mOps);
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    template <class ET>
-    inline const Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::set(ET pos) const
-    {
-	PackNum<sizeof...(Ops)-1>::setOpPos(mOps,pos);
-	return *this;
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    auto Operation<T,OpFunction,Ops...>::rootSteps(std::intptr_t iPtrNum) const
-	-> decltype(PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps))
-    {
-	return PackNum<sizeof...(Ops)-1>::mkSteps(iPtrNum, mOps);
-    }
-
-    template <typename T, class OpFunction, class... Ops>
-    template <class Expr>
-    auto Operation<T,OpFunction,Ops...>::loop(Expr exp) const
-	-> decltype(PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp ))
-    {
-	return PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp );
-    }
-
-    
-    /*********************
-     *   Contraction     *
-     *********************/
-
-    template <typename T, class Op, class IndexType>
-    Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
-	mOp(op),
-	mInd(ind) {}
-
-    // forward loop !!!!
-    template <typename T, class Op, class IndexType>
-    template <class ET>
-    inline T Contraction<T,Op,IndexType>::get(ET pos) const
-    {
-	return mOp.template get<ET>(pos);
-    }
-
-    template <typename T, class Op, class IndexType>
-    template <class ET>
-    inline const Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::set(ET pos) const
-    {
-	mOp.set(pos);
-	return *this;
-    }
-
-    template <typename T, class Op, class IndexType>
-    auto Contraction<T,Op,IndexType>::rootSteps(std::intptr_t iPtrNum) const
-	-> decltype(mOp.rootSteps(iPtrNum))
-    {
-	return mOp.rootSteps(iPtrNum);
-    }
-
-    template <typename T, class Op, class IndexType>
-    template <class Expr>
-    auto Contraction<T,Op,IndexType>::loop(Expr exp) const -> decltype(mInd->iforh(1,exp))
-    {
-	return mInd->iforh(1,exp);
-    }
-
-    /**************************
-     *   SliceContraction     *
-     **************************/
-
-    template <typename T, class Op, class... Indices>
-    SliceContraction<T,Op,Indices...>::SliceContraction(const Op& op,
-							std::shared_ptr<Indices>... ind) :
-	mOp(op),
-	mCont(std::make_shared<MultiArray<T,typename Indices::RangeType...> >(ind->range()...)),
-	mTarOp(*mCont,ind...)
-    { }
-
-    // forward loop !!!!
-    template <typename T, class Op, class... Indices>
-    template <class ET>
-    inline const MultiArray<T,typename Indices::RangeType...>&
-    SliceContraction<T,Op,Indices...>::get(ET pos) const
-    {
-	*mCont = 0;
-	mOp.set(pos);
-	mTarOp = mOp;
-	return *mCont;
-    }
-
-    template <typename T, class Op, class... Indices>
-    template <class ET>
-    inline const SliceContraction<T,Op,Indices...>& SliceContraction<T,Op,Indices...>::set(ET pos) const
-    {
-	mOp.set(pos);
-	return *this;
-    }
-
-    template <typename T, class Op, class... Indices>
-    auto SliceContraction<T,Op,Indices...>::rootSteps(std::intptr_t iPtrNum) const
-	-> decltype(mOp.rootSteps(iPtrNum))
-    {
-	return mOp.rootSteps(iPtrNum);
-    }
-
-    template <typename T, class Op, class... Indices>
-    template <class Expr>
-    auto SliceContraction<T,Op,Indices...>::loop(Expr exp) const -> decltype(mOp.loop(exp))
-    {
-	return mOp.loop(exp);
-    }
-
-}
 
 #include "type_operations.h"
 
diff --git a/src/include/ranges/dynamic_range.cc.h b/src/include/ranges/dynamic_range.cc.h
new file mode 100644
index 0000000..6498268
--- /dev/null
+++ b/src/include/ranges/dynamic_range.cc.h
@@ -0,0 +1,448 @@
+
+//#include "ranges/dynamic_range.h"
+//#include "ranges/dynamic_meta.h"
+
+#include "rpack_num.h"
+
+namespace MultiArrayTools
+{
+    namespace
+    {
+        using namespace MultiArrayHelper;
+    }
+    
+    /****************************
+     *   DynamicRangeFactory    *
+     ****************************/
+
+    template <class EC>
+    template <class... RangeTypes>
+    DynamicRangeFactory<EC>::DynamicRangeFactory(const std::tuple<std::shared_ptr<RangeTypes>...>& origs)
+    {
+	mProd = std::shared_ptr<oType>( new DynamicRange<EC>( origs ) );
+    }
+
+    template <class EC>
+    template <class... RangeTypes>
+    DynamicRangeFactory<EC>::DynamicRangeFactory(std::shared_ptr<RangeTypes>... origs)
+    {
+	mProd = std::shared_ptr<oType>( new DynamicRange<EC>( origs... ) );
+    }
+
+    template <class EC>
+    template <class Range>
+    void DynamicRangeFactory<EC>::append(std::shared_ptr<Range> r)
+    {
+	if(mProductCreated){
+	    
+	    mProd = std::shared_ptr<oType>( new DynamicRange<EC>( *std::dynamic_pointer_cast<oType>(mProd) ) );
+	    mProductCreated = false;
+	}
+	std::dynamic_pointer_cast<oType>(mProd)->mOrig.push_back(r);
+	std::dynamic_pointer_cast<oType>(mProd)->mSize *= r->size();
+	std::dynamic_pointer_cast<oType>(mProd)->mEmpty = false;
+    }
+
+    template <class EC>
+    DynamicRangeFactory<EC>::DynamicRangeFactory()
+    {
+	mProd = std::shared_ptr<oType>( new DynamicRange<EC>() );
+    }
+
+    // INSTANCIATE IF NEEDED!!
+    template <class EC>
+    std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > DynamicRangeFactory<EC>::mAleadyCreated;
+
+    template <class EC>
+    std::shared_ptr<RangeBase> DynamicRangeFactory<EC>::checkIfCreated(const std::vector<std::shared_ptr<RangeBase> >& pvec)
+    {
+	std::shared_ptr<RangeBase> out;
+	bool check = false;
+	for(auto& x: mAleadyCreated){
+	    if(x.second.size() == pvec.size()){
+		check = true;
+		for(size_t i = 0; i != x.second.size(); ++i){
+		    if(x.second[i] != reinterpret_cast<std::intptr_t>( pvec[i].get() ) ){
+			check = false;
+			break;
+		    }
+		}
+		if(check == true){
+		    out = x.first;
+		    break;
+		}
+	    }
+	}
+	if(not check){
+	    std::vector<std::intptr_t> app(pvec.size());
+	    for(size_t i = 0; i != app.size(); ++i){
+		app[i] = reinterpret_cast<std::intptr_t>( pvec[i].get() );
+	    }
+	    mAleadyCreated[mProd] = app;
+	    out = mProd;
+	}
+	return out;
+    }
+
+
+    template <class EC>
+    std::shared_ptr<RangeBase> DynamicRangeFactory<EC>::create()
+    {
+	mProd = checkIfCreated(std::dynamic_pointer_cast<DynamicRange<EC>>(mProd)->mOrig); 
+	setSelf();
+	mProductCreated = true;
+	return mProd;
+    }
+
+    
+    /*********************
+     *   DynamicIndex    *
+     *********************/
+
+    template <class EC>
+    DynamicIndex<EC>::DynamicIndex(const std::shared_ptr<DynamicRange<EC> >& range) :
+	IndexInterface<DynamicIndex,MetaType>(range, 0)//,
+	//mExplicitRangePtr(std::dynamic_pointer_cast<RangeType>(IB::mRangePtr)),
+    {}
+
+    template <class EC>
+    IndexType DynamicIndex<EC>::type() const
+    {
+	return IndexType::SINGLE;
+    }
+
+    template <class EC>
+    DynamicIndex<EC>& DynamicIndex<EC>::operator=(size_t pos)
+    {
+	IB::mPos = pos;
+	return *this;
+    }
+
+    template <class EC>
+    DynamicIndex<EC>& DynamicIndex<EC>::operator++()
+    {
+        size_t ipos = mIVec.size()-1;
+        auto& ii = mIVec[ipos].first;
+        auto& jj = mIVec[ipos-1].first;
+        while(ii->pos() == ii->max()-1 and ipos != 0) {
+            (*ii) = 0;
+            ++(*jj);
+            --ipos;
+        }
+	return *this;
+    }
+
+    template <class EC>
+    DynamicIndex<EC>& DynamicIndex<EC>::operator--()
+    {
+        size_t ipos = mIVec.size()-1;
+        auto& ii = mIVec[ipos].first;
+        auto& jj = mIVec[ipos-1].first;
+        while(ii->pos() == 0 and ipos != 0) {
+            (*ii) = ii->max()-1;
+            --(*jj);
+            --ipos;
+        }
+	return *this;
+    }
+
+    template <class EC>
+    DynamicIndex<EC>& DynamicIndex<EC>::operator()(const IVecT& ivec)
+    {
+        mIVec = ivec;
+        return *this;
+    }
+    
+    template <class EC>
+    int DynamicIndex<EC>::pp(std::intptr_t idxPtrNum)
+    {
+	++(*this);
+	return 1;
+    }
+
+    template <class EC>
+    int DynamicIndex<EC>::mm(std::intptr_t idxPtrNum)
+    {
+	--(*this);
+	return 1;
+    }
+
+    template <class EC>
+    std::string DynamicIndex<EC>::stringMeta() const
+    {
+	return std::dynamic_pointer_cast<DynamicRange<EC> const>( IB::mRangePtr )->stringMeta(IB::mPos);
+    }
+    
+    template <class EC>
+    typename DynamicIndex<EC>::MetaType DynamicIndex<EC>::meta() const
+    {
+	return std::dynamic_pointer_cast<DynamicRange<EC> const>( IB::mRangePtr )->get(IB::mPos);
+    }
+    
+    template <class EC>
+    const typename DynamicIndex<EC>::MetaType* DynamicIndex<EC>::metaPtr() const
+    {
+	return nullptr;
+    }
+
+    /*
+    bool DynamicIndex<EC>::isMeta(const MetaType& metaPos) const
+    {
+	return mExplicitRangePtr->isMeta(metaPos);
+        }*/
+    
+    template <class EC>
+    DynamicIndex<EC>& DynamicIndex<EC>::at(const MetaType& metaPos)
+    {
+	(*this) = std::dynamic_pointer_cast<DynamicRange<EC> const>( IB::mRangePtr )->getMeta( metaPos );
+	return *this;
+    }
+
+    template <class EC>
+    size_t DynamicIndex<EC>::posAt(const MetaType& metaPos) const
+    {
+	return std::dynamic_pointer_cast<DynamicRange<EC> const>( IB::mRangePtr )->getMeta( metaPos );
+    }
+
+    template <class EC>
+    size_t DynamicIndex<EC>::dim() const // = 1
+    {
+	return 1;
+    }
+
+    template <class EC>
+    bool DynamicIndex<EC>::last() const
+    {
+	return IB::mPos == IB::mMax - 1;
+    }
+
+    template <class EC>
+    bool DynamicIndex<EC>::first() const
+    {
+	return IB::mPos == 0;
+    }
+
+    template <class EC>
+    const IndexW<EC>& DynamicIndex<EC>::get(size_t n) const
+    {
+        return *mIVec[n].first;
+    }
+
+    template <class EC>
+    std::shared_ptr<typename DynamicIndex<EC>::RangeType> DynamicIndex<EC>::range()
+    {
+	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
+    }
+
+
+    template <class EC>
+    template <size_t N>
+    void DynamicIndex<EC>::getPtr() {}
+    
+    template <class EC>
+    size_t DynamicIndex<EC>::getStepSize(size_t n) const
+    {
+	return mIVec[n].second;
+    }
+
+    template <class EC>
+    std::string DynamicIndex<EC>::id() const
+    {
+	return std::string("dyn") + std::to_string(IB::mId);
+    }
+
+    template <class EC>    
+    void DynamicIndex<EC>::print(size_t offset)
+    {
+	if(offset == 0){
+	    std::cout << " === " << std::endl;
+	}
+	for(size_t j = 0; j != offset; ++j) { std::cout << "\t"; }
+	std::cout << id() << "[" << reinterpret_cast<std::intptr_t>(this)
+		  << "](" << IB::mRangePtr << "): " /*<< meta().first*/ << std::endl;
+    }
+    
+    template <class EC>
+    inline DynamicExpression DynamicIndex<EC>::mkFor(size_t i, size_t step,
+						   DynamicExpression ex, bool hidden) const
+    {
+	if(i != 0){
+	    auto& ii = *mIVec[i].first;
+	    return mkFor(i-1, step,  hidden ? ii.iforh(step, ex) : ii.ifor(step, ex));
+	}
+	else {
+	    auto& ii = *mIVec[0].first;
+	    return hidden ? ii.iforh(step, ex) : ii.ifor(step, ex);
+	}		
+    }
+
+    template <class EC>
+    template <class Expr>
+    auto DynamicIndex<EC>::ifor(size_t step, Expr ex) const
+        -> DynamicExpression
+    {
+	DynamicExpression expr(std::make_shared<Expr>(ex));
+	return mkFor(mIVec.size()-1, step, expr);
+    }
+
+    template <class EC>
+    template <class Expr>
+    auto DynamicIndex<EC>::iforh(size_t step, Expr ex) const
+        -> DynamicExpression
+    {
+	DynamicExpression expr(std::make_shared<Expr>(ex));
+	return mkFor(mIVec.size()-1, step, expr, true);
+    }
+
+
+
+    /***********************
+     *   DynamicRange    *
+     ***********************/
+
+    template <class EC>
+    typename DynamicRange<EC>::MetaType DynamicRange<EC>::get(size_t pos) const
+    {
+	return MetaType(); // !!!
+    }
+
+    template <class EC>
+    size_t DynamicRange<EC>::getMeta(const MetaType& metaPos) const
+    {
+	return 0; // !!!
+    }
+
+    template <class EC>
+    size_t DynamicRange<EC>::size() const
+    {
+	return mSize;
+    }
+
+    template <class EC>
+    size_t DynamicRange<EC>::dim() const
+    {
+        return mOrig.size();
+    }
+
+    template <class EC>
+    SpaceType DynamicRange<EC>::spaceType() const
+    {
+	return SpaceType::DYN;
+    }
+
+    template <class EC>
+    bool DynamicRange<EC>::isEmpty() const
+    {
+	return mEmpty;
+    }
+
+    template <class EC>
+    std::string DynamicRange<EC>::stringMeta(size_t pos) const
+    {
+	std::string out = "[ ";
+	//size_t xpos = pos;
+	for(size_t i = mOrig.size(); i != 0; --i) {
+	    auto& x = mOrig[i-1];
+	    const size_t redpos = pos % x->size();
+	    out = ( (i == mOrig.size()) ? out : out + " , " ) + x->stringMeta(redpos);
+	    pos -= redpos;
+	    pos /= x->size();
+	}
+	out += " ]";
+	return out;
+    }
+
+    template <class EC>
+    std::vector<char> DynamicRange<EC>::data() const
+    {
+	DataHeader h;
+	h.spaceType = static_cast<int>( SpaceType::DYN );
+	h.metaSize = mOrig.size();
+	h.multiple = 1;
+	std::vector<char> out;
+	char* hcp = reinterpret_cast<char*>(&h);
+	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
+	for(auto& x: mOrig){
+	    auto part = x->data();
+	    out.insert(out.end(), part.begin(), part.end());
+	}
+	return out;
+    }
+
+    template <class EC>
+    typename DynamicRange<EC>::IndexType DynamicRange<EC>::begin() const
+    {
+	DynamicIndex<EC> i
+	    (std::dynamic_pointer_cast<DynamicRange>
+	     ( std::shared_ptr<RangeBase>(RB::mThis) ) );
+	i = 0;
+	return i;
+    }
+
+    template <class EC>
+    typename DynamicRange<EC>::IndexType DynamicRange<EC>::end() const
+    {
+	DynamicIndex<EC> i
+	    (std::dynamic_pointer_cast<DynamicRange>
+	     ( std::shared_ptr<RangeBase>(RB::mThis) ) );
+	i = size();
+	return i;
+    }
+
+    template <class EC>
+    std::shared_ptr<RangeBase> DynamicRange<EC>::sub(size_t num) const
+    {
+	return mOrig.at(num);
+    }
+
+    template <class EC>
+    void DynamicRange<EC>::sreplace(const std::shared_ptr<RangeBase> in, size_t num)
+    {
+	assert(mOrig[num]->size() == in->size());
+	mOrig[num] = in;
+    }
+
+    template <class EC>
+    template <class... RangeTypes>
+    DynamicRange<EC>::DynamicRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs) :
+	RangeInterface<DynamicIndex<EC>>()
+    {
+	RPackNum<sizeof...(RangeTypes)-1>::RangesToVec( origs, mOrig );
+	mSize = RPackNum<sizeof...(RangeTypes)-1>::getSize( origs );
+	if(sizeof...(RangeTypes)){
+	    mEmpty = false;
+	}
+    }
+
+    template <class EC>
+    template <class... RangeTypes>
+    DynamicRange<EC>::DynamicRange(std::shared_ptr<RangeTypes>... origs) :
+	RangeInterface<DynamicIndex<EC>>()
+    {
+	auto rst = std::make_tuple(origs...);
+	RPackNum<sizeof...(RangeTypes)-1>::RangesToVec( rst, mOrig );
+	mSize = RPackNum<sizeof...(RangeTypes)-1>::getSize( rst );
+	if(sizeof...(RangeTypes)){
+	    mEmpty = false;
+	}
+    }
+
+    template <class EC>
+    template <class Range>
+    std::shared_ptr<Range> DynamicRange<EC>::fullsub(size_t num) const
+    {
+	return std::dynamic_pointer_cast<Range>( mOrig.at(num) );
+    }
+
+    template <class EC>
+    template <class... Ranges>
+    std::shared_ptr<MultiRange<Ranges...> > DynamicRange<EC>::scast(SIZET<Ranges>... sizes) const
+    {
+	std::tuple<std::shared_ptr<Ranges>...> rtp;
+	RPackNum<sizeof...(Ranges)-1>::resolveRangeType(mOrig, rtp, 0, sizes...);
+	MultiRangeFactory<Ranges...> mrf(rtp);
+	return std::dynamic_pointer_cast<MultiRange<Ranges...> >( mrf.create() );
+    }
+
+    
+} // end namespace MultiArrayTools
diff --git a/src/include/ranges/dynamic_range.h b/src/include/ranges/dynamic_range.h
index 9995fc7..ce6f9cc 100644
--- a/src/include/ranges/dynamic_range.h
+++ b/src/include/ranges/dynamic_range.h
@@ -2,6 +2,8 @@
 #ifndef __dynamic_range_h__
 #define __dynamic_range_h__
 
+#include <cassert>
+
 #include "ranges/rbase_def.h"
 #include "ranges/range_base.h"
 #include "ranges/index_base.h"
@@ -9,7 +11,7 @@
 #include "xfor/xfor.h"
 
 #include <map>
-#include "ranges/rpheader.h"
+//#include "ranges/rpheader.h"
 #include "ranges/x_to_string.h"
 #include "ranges/type_map.h"
 
@@ -18,10 +20,13 @@
 namespace MultiArrayTools
 {
 
-    using MultiArrayHelper::DynamicalExpression;
-    
+    using MultiArrayHelper::DynamicExpression;
+
+    template <class ExpressionCollection>
     class IndexWrapperBase
     {
+    protected:
+        std::shared_ptr<ExpressionCollection> mEc;
     public:
 
         IndexWrapperBase() = default;
@@ -60,15 +65,26 @@ namespace MultiArrayTools
         
         virtual std::intptr_t get() const = 0;
 
-	virtual DynamicalExpression ifor(size_t step, DynamicalExpression ex) const = 0;
-	virtual DynamicalExpression iforh(size_t step, DynamicalExpression ex) const = 0;
+        template <class Expr>
+        auto ifor(size_t step, Expr ex) const
+            -> decltype(mEc->ifor(step, ex))
+        { return mEc->ifor(step, ex); }
+
+        template <class Expr>
+        auto iforh(size_t step, Expr ex) const
+            -> decltype(mEc->iforh(step, ex))
+        { return mEc->iforh(step, ex); }
+
     };
 
-    typedef IndexWrapperBase IndexW;
+    template <class EC>
+    using IndexW = IndexWrapperBase<EC>;
     
-    template <class Index, class Expr>
-    class IndexWrapper : public IndexWrapperBase
+    template <class Index, class ExpressionCollection>
+    class IndexWrapper : public IndexWrapperBase<ExpressionCollection>
     {
+    public:
+        typedef IndexWrapperBase<ExpressionCollection> IWB;
     protected:
         IndexWrapper() = default;
 
@@ -81,13 +97,14 @@ namespace MultiArrayTools
         IndexWrapper& operator=(const IndexWrapper& in) = default;
         IndexWrapper& operator=(IndexWrapper&& in) = default;
 
-        IndexWrapper(const std::shared_ptr<Index>& i) : mI(i) {}
+        IndexWrapper(const std::shared_ptr<Index>& i) : mI(i)
+        { IWB::mEc = ExpressionCollection::make(mI); }
         
         virtual IndexType type() const final { return mI->type(); }
 	
-	virtual IndexWrapperBase& operator=(size_t pos) final { (*mI) = pos; return *this; }
-	virtual IndexWrapperBase& operator++() final { ++(*mI); return *this; }
-	virtual IndexWrapperBase& operator--() final { --(*mI); return *this; }
+	virtual IndexWrapper& operator=(size_t pos) final { (*mI) = pos; return *this; }
+	virtual IndexWrapper& operator++() final { ++(*mI); return *this; }
+	virtual IndexWrapper& operator--() final { --(*mI); return *this; }
 
         virtual size_t pos() const final { return mI->pos(); }
         virtual size_t max() const final { return mI->max(); }
@@ -98,7 +115,7 @@ namespace MultiArrayTools
 	virtual std::string stringMeta() const final { return mI->stringMeta(); }
 	virtual DynamicMetaT meta() const final { return DynamicMetaT(mI->meta()); }
 	virtual const DynamicMetaT* metaPtr() const final { return nullptr; }
-	IndexWrapperBase& at(const typename Index::MetaType& metaPos) { mI->at(metaPos); return *this; }
+	IndexWrapper& at(const typename Index::MetaType& metaPos) { mI->at(metaPos); return *this; }
 	size_t posAt(const typename Index::MetaType& metaPos) const { return mI->posAt(metaPos); }
 
 	//virtual bool isMeta(const U& metaPos) const final { return mI->isMeta(); }
@@ -114,30 +131,28 @@ namespace MultiArrayTools
         
         virtual std::intptr_t get() const final { return reinterpret_cast<std::intptr_t>(mI.get()); }
 	
-	virtual DynamicalExpression ifor(size_t step, DynamicalExpression ex) const final { return mI->ifor(step, ex); }
-	virtual DynamicalExpression iforh(size_t step, DynamicalExpression ex) const final { return mI->iforh(step, ex); }
-	
     };
     
-    typedef SingleRange<size_t,SpaceType::DYN> DynamicRange;
-    
-    class DynamicIndex : public IndexInterface<DynamicIndex,DynamicMetaT>
+    //typedef SingleRange<size_t,SpaceType::DYN> DynamicRange;
+
+    template <class EC>
+    class DynamicIndex : public IndexInterface<DynamicIndex<EC>,DynamicMetaT>
     {
     private:
-        typedef std::vector<std::pair<std::shared_ptr<IndexW>,size_t>> IVecT;
+        typedef std::vector<std::pair<std::shared_ptr<IndexW<EC>>,size_t>> IVecT;
         
         IVecT mIVec;
 
-	inline DynamicalExpression mkFor(size_t i, size_t step,
-					 DynamicalExpression ex, bool hidden = false) const;	
+	inline DynamicExpression mkFor(size_t i, size_t step,
+					 DynamicExpression ex, bool hidden = false) const;	
 
     public:
-        typedef IndexInterface<DynamicIndex,DynamicMetaT> IB;
+        typedef IndexInterface<DynamicIndex<EC>,DynamicMetaT> IB;
 	typedef DynamicMetaT MetaType;
-	typedef DynamicRange RangeType;
+	typedef DynamicRange<EC> RangeType;
 	typedef DynamicIndex IType;
 
-        DynamicIndex(const std::shared_ptr<DynamicRange>& range);
+        DynamicIndex(const std::shared_ptr<RangeType>& range);
 
         static constexpr IndexType sType() { return IndexType::SINGLE; }
 	static constexpr size_t totalDim() { return 1; }
@@ -168,7 +183,7 @@ namespace MultiArrayTools
 	bool last() const;
 	bool first() const;
 
-        const IndexW& get(size_t n) const;
+        const IndexW<EC>& get(size_t n) const;
 
 	std::shared_ptr<RangeType> range();
 	
@@ -182,21 +197,22 @@ namespace MultiArrayTools
 
 	template <class Expr>
 	auto ifor(size_t step, Expr ex) const
-	    -> DynamicalExpression;
+	    -> DynamicExpression;
 
 	template <class Expr>
 	auto iforh(size_t step, Expr ex) const
-	    -> DynamicalExpression;
+	    -> DynamicExpression;
 
     };    
 
     
     // NOT THREAD SAVE!!
+    template <class EC>
     class DynamicRangeFactory : public RangeFactoryBase
     {
     public:
 	
-	typedef DynamicRange oType;
+	typedef DynamicRange<EC> oType;
 
 	DynamicRangeFactory();
 
@@ -220,8 +236,8 @@ namespace MultiArrayTools
 	bool mProductCreated = false;
     };
 
-    template <>
-    class SingleRange<size_t,SpaceType::DYN> : public RangeInterface<DynamicIndex>
+    template <class EC>
+    class DynamicRange : public RangeInterface<DynamicIndex<EC>>
     {
     public:
         static constexpr bool defaultable = true;
@@ -230,19 +246,19 @@ namespace MultiArrayTools
 	static constexpr bool HASMETACONT = false;
 	
 	typedef RangeBase RB;
-	typedef DynamicIndex IndexType;
+	typedef DynamicIndex<EC> IndexType;
 	typedef DynamicRange RangeType;
 	typedef DynamicMetaT MetaType;
 
     private:
-	SingleRange() = default;
-	SingleRange(const SingleRange& in) = default;
+	DynamicRange() = default;
+	DynamicRange(const DynamicRange& in) = default;
 
 	template <class... RangeTypes>
-	SingleRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs);
+	DynamicRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs);
 	
 	template <class... RangeTypes>
-	SingleRange(std::shared_ptr<RangeTypes>... origs);
+	DynamicRange(std::shared_ptr<RangeTypes>... origs);
 
 	size_t mSize = 1;
 	bool mEmpty = true;
@@ -276,79 +292,36 @@ namespace MultiArrayTools
 	
 	bool isEmpty() const;
 	
-	friend DynamicRangeFactory;
+	friend DynamicRangeFactory<EC>;
 	
-	static DynamicRangeFactory factory()
-	{ return DynamicRangeFactory(); }
+	static DynamicRangeFactory<EC> factory()
+	{ return DynamicRangeFactory<EC>(); }
 
     };
 
 } // namespace MultiArrayTools
 
 
-/* ========================= *
- * ---   TEMPLATE CODE   --- *
- * ========================= */
-
-namespace MultiArrayTools
-{
-       
-    /***********************
-     *   DynamicRange    *
-     ***********************/
-    
-    template <class... RangeTypes>
-    DynamicRangeFactory::DynamicRangeFactory(const std::tuple<std::shared_ptr<RangeTypes>...>& origs)
-    {
-	mProd = std::shared_ptr<oType>( new DynamicRange( origs ) );
-    }
-
-    template <class... RangeTypes>
-    DynamicRangeFactory::DynamicRangeFactory(std::shared_ptr<RangeTypes>... origs)
-    {
-	mProd = std::shared_ptr<oType>( new DynamicRange( origs... ) );
-    }
-
-    template <class Range>
-    void DynamicRangeFactory::append(std::shared_ptr<Range> r)
-    {
-	if(mProductCreated){
-	    
-	    mProd = std::shared_ptr<oType>( new DynamicRange( *std::dynamic_pointer_cast<oType>(mProd) ) );
-	    mProductCreated = false;
-	}
-	std::dynamic_pointer_cast<oType>(mProd)->mOrig.push_back(r);
-	std::dynamic_pointer_cast<oType>(mProd)->mSize *= r->size();
-	std::dynamic_pointer_cast<oType>(mProd)->mEmpty = false;
-    }
-
-    /*****************
-     *   Functions   *
-     *****************/
-
-    //std::shared_ptr<DynamicRange> defaultRange(size_t size = 0);
-}
-
 namespace MultiArrayHelper
 {
     using namespace MultiArrayTools;
 
-    template <>
-    inline void resolveSetRange<DynamicRange>(std::shared_ptr<DynamicRange>& rp,
-						const std::vector<std::shared_ptr<RangeBase> >& orig,
-						size_t origpos, size_t size)
+    template <class EC>
+    inline void resolveSetRange(std::shared_ptr<DynamicRange<EC>>& rp,
+                                const std::vector<std::shared_ptr<RangeBase> >& orig,
+                                size_t origpos, size_t size)
     {
-    	DynamicRangeFactory arf;
+    	DynamicRangeFactory<EC> arf;
 	for(size_t op = origpos; op != origpos + size; ++op){
 	    //VCHECK(op);
 	    arf.append(orig[op]);
 	}
-    	rp = std::dynamic_pointer_cast<DynamicRange>( arf.create() );
+    	rp = std::dynamic_pointer_cast<DynamicRange<EC>>( arf.create() );
     }
 
-    template <>
-    inline void setRangeToVec<DynamicRange>(std::vector<std::shared_ptr<RangeBase> >& v,
-					      std::shared_ptr<DynamicRange> r)
+    template <class EC>
+    inline void setRangeToVec(std::vector<std::shared_ptr<RangeBase> >& v,
+                              std::shared_ptr<DynamicRange<EC>> r)
     {
 	if(not r->isEmpty()){
 	    for(size_t i = r->dim(); i != 0; --i){
@@ -357,8 +330,8 @@ namespace MultiArrayHelper
 	}
     }
 
-    template <>
-    inline size_t getStepSize<DynamicIndex>(const DynamicIndex& ii, std::intptr_t j)
+    template <class EC>
+    inline size_t getStepSize(const DynamicIndex<EC>& ii, std::intptr_t j)
     {
         size_t ss = 0;
         size_t sx = 1;
@@ -374,79 +347,6 @@ namespace MultiArrayHelper
 
 }
 
-namespace MultiArrayTools
-{
-    /***********************
-     *   DynamicRange    *
-     ***********************/
-
-    template <class... RangeTypes>
-    SingleRange<size_t,SpaceType::DYN>::SingleRange(const std::tuple<std::shared_ptr<RangeTypes>...>& origs) :
-	RangeInterface<DynamicIndex>()
-    {
-	RPackNum<sizeof...(RangeTypes)-1>::RangesToVec( origs, mOrig );
-	mSize = RPackNum<sizeof...(RangeTypes)-1>::getSize( origs );
-	if(sizeof...(RangeTypes)){
-	    mEmpty = false;
-	}
-    }
-    
-    template <class... RangeTypes>
-    SingleRange<size_t,SpaceType::DYN>::SingleRange(std::shared_ptr<RangeTypes>... origs) :
-	RangeInterface<DynamicIndex>()
-    {
-	auto rst = std::make_tuple(origs...);
-	RPackNum<sizeof...(RangeTypes)-1>::RangesToVec( rst, mOrig );
-	mSize = RPackNum<sizeof...(RangeTypes)-1>::getSize( rst );
-	if(sizeof...(RangeTypes)){
-	    mEmpty = false;
-	}
-    }
-
-    template <class Range>
-    std::shared_ptr<Range> SingleRange<size_t,SpaceType::DYN>::fullsub(size_t num) const
-    {
-	return std::dynamic_pointer_cast<Range>( mOrig.at(num) );
-    }
-
-    template <class... Ranges>
-    std::shared_ptr<MultiRange<Ranges...> > SingleRange<size_t,SpaceType::DYN>::scast(SIZET<Ranges>... sizes) const
-    {
-	std::tuple<std::shared_ptr<Ranges>...> rtp;
-	RPackNum<sizeof...(Ranges)-1>::resolveRangeType(mOrig, rtp, 0, sizes...);
-	MultiRangeFactory<Ranges...> mrf(rtp);
-	return std::dynamic_pointer_cast<MultiRange<Ranges...> >( mrf.create() );
-    }
-
-    inline DynamicalExpression DynamicIndex::mkFor(size_t i, size_t step,
-						   DynamicalExpression ex, bool hidden) const
-    {
-	if(i != 0){
-	    auto& ii = *mIVec[i].first;
-	    return mkFor(i-1, step,  hidden ? ii.iforh(step, ex) : ii.ifor(step, ex));
-	}
-	else {
-	    auto& ii = *mIVec[0].first;
-	    return hidden ? ii.iforh(step, ex) : ii.ifor(step, ex);
-	}		
-    }
-    
-    template <class Expr>
-    auto DynamicIndex::ifor(size_t step, Expr ex) const
-        -> DynamicalExpression
-    {
-	DynamicalExpression expr(std::make_shared<Expr>(ex));
-	return mkFor(mIVec.size()-1, step, expr);
-    }
-    
-    template <class Expr>
-    auto DynamicIndex::iforh(size_t step, Expr ex) const
-        -> DynamicalExpression
-    {
-	DynamicalExpression expr(std::make_shared<Expr>(ex));
-	return mkFor(mIVec.size()-1, step, expr, true);
-    }
-    
-}
+#include "dynamic_range.cc.h"
 
 #endif
diff --git a/src/include/ranges/rbase_def.h b/src/include/ranges/rbase_def.h
index 80284c5..954f310 100644
--- a/src/include/ranges/rbase_def.h
+++ b/src/include/ranges/rbase_def.h
@@ -79,13 +79,16 @@ namespace MultiArrayTools
     //class AnonymousRange;
 
     // dynamic_range.h
+    template <class EC>
     class DynamicIndex;
     
     // dynamic_range.h
+    template <class EC>
     class DynamicRangeFactory;
     
     // dynamic_range.h
-    //class DynamicRange;
+    template <class EC>
+    class DynamicRange;
 
     // value_range.h
     template <typename U>
diff --git a/src/include/ranges/rheader.h b/src/include/ranges/rheader.h
index 2b68d82..6c88ee0 100644
--- a/src/include/ranges/rheader.h
+++ b/src/include/ranges/rheader.h
@@ -2,9 +2,9 @@
 //#ifndef __rheader_h__
 //#define __rheader_h__
 
+#include "dynamic_range.h"
 #include "rpheader.h"
 #include "anonymous_range.h"
-#include "dynamic_range.h"
 #include "value_range.h"
 
 //#endif
diff --git a/src/include/slice.cc.h b/src/include/slice.cc.h
new file mode 100644
index 0000000..72bbf94
--- /dev/null
+++ b/src/include/slice.cc.h
@@ -0,0 +1,256 @@
+
+#include "slice.h"
+
+namespace MultiArrayTools
+{
+
+    /*******************
+     *   ConstSlice    *
+     *******************/
+
+    template <typename T, class... SRanges>
+    void ConstSlice<T,SRanges...>::format(const std::array<size_t,sizeof...(SRanges)+1>& blocks)
+    {
+	MAB::mProtoI->format(blocks);
+    }
+    
+    template <typename T, class... SRanges>
+    ConstSlice<T,SRanges...>::ConstSlice(const std::shared_ptr<SRanges>&... ranges, const T* data) :
+	MultiArrayBase<T,SRanges...>(ranges...),
+	mData(data)
+    {
+	MAB::mInit = true;
+    }
+
+    template <typename T, class... SRanges>
+    ConstSlice<T,SRanges...>::ConstSlice(const MultiArrayBase<T,AnonymousRange>& ma, SIZET<SRanges>... sizes) :
+	MultiArrayBase<T,SRanges...>
+	( ma.range()->template get<0>().template scast<SRanges...>(sizes...)->space() ),
+	mData( ma.data() )
+    {
+	MAB::mInit = true;
+    }
+    
+    template <typename T, class... SRanges>
+    const T& ConstSlice<T,SRanges...>::operator[](const IType& i) const
+    {
+	//assert(i.sliceMode()); // -> compare objects !!!!!
+	assert(i.container() == reinterpret_cast<std::intptr_t>(this));
+	return mData[ i.pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    const T& ConstSlice<T,SRanges...>::at(const typename IType::MetaType& meta) const
+    {
+	//auto x = begin().at(meta);
+	//VCHECK(x.pos());
+	return mData[ begin().at(meta).pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    const T* ConstSlice<T,SRanges...>::data() const
+    {
+	return mData;
+    }
+
+    template <typename T, class... SRanges>
+    bool ConstSlice<T,SRanges...>::isSlice() const
+    {
+	return true;
+    }
+
+    template <typename T, class... SRanges>
+    auto ConstSlice<T,SRanges...>::begin() const -> ConstSlice<T,SRanges...>::IType 
+    {
+	IType i(*MAB::mProtoI);
+	i = 0;
+	//i = mStartPos;
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    auto ConstSlice<T,SRanges...>::end() const -> ConstSlice<T,SRanges...>::IType 
+    {
+	IType i(*MAB::mProtoI);
+	i = i.max(); // CHECK !!!
+	//i = std::get<sizeof...(SRanges)>(mBlockSizes);
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > ConstSlice<T,SRanges...>::anonymous(bool slice) const
+    {
+	assert(0); // think about carefully!!!!
+	return nullptr;
+    }
+    
+    template <typename T, class... SRanges>
+    auto ConstSlice<T,SRanges...>::define(const std::shared_ptr<typename SRanges::IndexType>&... inds)
+	-> ConstSliceDef<T,SRanges...>
+    {
+	return ConstSliceDef<T,SRanges...>(*this, inds...);
+    }
+
+    
+    /**************
+     *   Slice    *
+     **************/
+
+    template <typename T, class... SRanges>
+    void Slice<T,SRanges...>::format(const std::array<size_t,sizeof...(SRanges)+1>& blocks)
+    {
+	MAB::mProtoI->format(blocks);
+    }
+    
+    template <typename T, class... SRanges>
+    Slice<T,SRanges...>::Slice(const std::shared_ptr<SRanges>&... ranges, T* data) :
+	MutableMultiArrayBase<T,SRanges...>(ranges...),
+	mData(data) {}
+    
+    template <typename T, class... SRanges>
+    const T& Slice<T,SRanges...>::operator[](const IType& i) const
+    {
+	//assert(i.sliceMode()); // -> compare objects !!!!!
+	assert(i.container() == reinterpret_cast<std::intptr_t>(this));
+	return mData[ i.pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    T& Slice<T,SRanges...>::operator[](const IType& i)
+    {
+	//assert(i.sliceMode());
+	assert(i.container() == reinterpret_cast<std::intptr_t>(this));
+	return mData[ i.pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    const T& Slice<T,SRanges...>::at(const typename IType::MetaType& meta) const
+    {
+	//auto x = begin().at(meta);
+	//VCHECK(x.pos());
+	return mData[ begin().at(meta).pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    T& Slice<T,SRanges...>::at(const typename IType::MetaType& meta)
+    {
+	//auto x = begin().at(meta);
+	//VCHECK(x.pos());
+	return mData[ begin().at(meta).pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    const T* Slice<T,SRanges...>::data() const
+    {
+	return mData;
+    }
+
+    template <typename T, class... SRanges>
+    T* Slice<T,SRanges...>::data()
+    {
+	return mData;
+    }
+
+    template <typename T, class... SRanges>
+    bool Slice<T,SRanges...>::isSlice() const
+    {
+	return true;
+    }
+
+    template <typename T, class... SRanges>
+    auto Slice<T,SRanges...>::begin() const -> Slice<T,SRanges...>::IType 
+    {
+	IType i(*MAB::mProtoI);
+	i = 0;
+	//i = mStartPos;
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    auto Slice<T,SRanges...>::end() const -> Slice<T,SRanges...>::IType 
+    {
+	IType i(*MAB::mProtoI);
+	i = i.max(); // CHECK !!!
+	//i = std::get<sizeof...(SRanges)>(mBlockSizes);
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > Slice<T,SRanges...>::anonymous(bool slice) const
+    {
+	assert(0); // think about carefully!!!!
+	return nullptr;
+    }
+    
+    template <typename T, class... SRanges>
+    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > Slice<T,SRanges...>::anonymousMove()
+    {
+	assert(0); // think about carefully!!!!
+	return nullptr;
+    }
+    
+    template <typename T, class... SRanges>
+    auto Slice<T,SRanges...>::define(const std::shared_ptr<typename SRanges::IndexType>&... inds)
+	-> SliceDef<T,SRanges...>
+    {
+	return SliceDef<T,SRanges...>(*this, inds...);
+    }
+    
+
+    template <typename T, class... SRanges>
+    SliceDef<T,SRanges...>::SliceDef(Slice<T,SRanges...>& sl,
+				     const std::shared_ptr<typename SRanges::IndexType>&... inds) :
+	mIndex(sl.begin()),
+	mSl(sl)
+    {
+	mIndex(inds...);
+    }
+
+    template <typename T, class... SRanges>
+    template <class... ORanges>
+    SliceDef<T,SRanges...>& SliceDef<T,SRanges...>::operator=(const OperationRoot<T,ORanges...>& op)
+    {
+	std::array<size_t,sizeof...(SRanges)+1> blocks;
+	PackNum<sizeof...(SRanges)-1>::
+	    template mkSliceBlocks<T,OperationRoot<T,ORanges...>,SRanges...>(blocks, mIndex, op);
+	mSl.format(blocks);
+	mSl.mData = op.data();
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    ConstSliceDef<T,SRanges...>::ConstSliceDef(ConstSlice<T,SRanges...>& sl,
+					       const std::shared_ptr<typename SRanges::IndexType>&... inds) :
+	mIndex(sl.begin()),
+	mSl(sl)
+    {
+	mIndex(inds...);
+    }
+
+    template <typename T, class... SRanges>
+    template <class... ORanges>
+    ConstSliceDef<T,SRanges...>& ConstSliceDef<T,SRanges...>::operator=(const ConstOperationRoot<T,ORanges...>& op)
+    {
+	std::array<size_t,sizeof...(SRanges)+1> blocks;
+	PackNum<sizeof...(SRanges)-1>::
+	    template mkSliceBlocks<T,ConstOperationRoot<T,ORanges...>,SRanges...>(blocks, mIndex, op);
+	mSl.format(blocks);
+	mSl.mData = op.data();
+	return *this;
+    }
+
+    template <typename T, class... SRanges>
+    template <class... ORanges>
+    ConstSliceDef<T,SRanges...>& ConstSliceDef<T,SRanges...>::operator=(const OperationRoot<T,ORanges...>& op)
+    {
+	std::array<size_t,sizeof...(SRanges)+1> blocks;
+	PackNum<sizeof...(SRanges)-1>::
+	    template mkSliceBlocks<T,OperationRoot<T,ORanges...>,SRanges...>(blocks, mIndex, op);
+	mSl.format(blocks);
+	mSl.mData = op.data();
+	return *this;
+    }
+
+
+} // end namespace MultiArrayTools
+
diff --git a/src/include/slice.h b/src/include/slice.h
index 2937fa4..7b976fb 100644
--- a/src/include/slice.h
+++ b/src/include/slice.h
@@ -179,257 +179,5 @@ namespace MultiArrayTools
  * ---   TEMPLATE CODE   --- *
  * ========================= */
 
-namespace MultiArrayTools
-{
-
-    /*******************
-     *   ConstSlice    *
-     *******************/
-
-    template <typename T, class... SRanges>
-    void ConstSlice<T,SRanges...>::format(const std::array<size_t,sizeof...(SRanges)+1>& blocks)
-    {
-	MAB::mProtoI->format(blocks);
-    }
-    
-    template <typename T, class... SRanges>
-    ConstSlice<T,SRanges...>::ConstSlice(const std::shared_ptr<SRanges>&... ranges, const T* data) :
-	MultiArrayBase<T,SRanges...>(ranges...),
-	mData(data)
-    {
-	MAB::mInit = true;
-    }
-
-    template <typename T, class... SRanges>
-    ConstSlice<T,SRanges...>::ConstSlice(const MultiArrayBase<T,AnonymousRange>& ma, SIZET<SRanges>... sizes) :
-	MultiArrayBase<T,SRanges...>
-	( ma.range()->template get<0>().template scast<SRanges...>(sizes...)->space() ),
-	mData( ma.data() )
-    {
-	MAB::mInit = true;
-    }
-    
-    template <typename T, class... SRanges>
-    const T& ConstSlice<T,SRanges...>::operator[](const IType& i) const
-    {
-	//assert(i.sliceMode()); // -> compare objects !!!!!
-	assert(i.container() == reinterpret_cast<std::intptr_t>(this));
-	return mData[ i.pos() ];
-    }
-    
-    template <typename T, class... SRanges>
-    const T& ConstSlice<T,SRanges...>::at(const typename IType::MetaType& meta) const
-    {
-	//auto x = begin().at(meta);
-	//VCHECK(x.pos());
-	return mData[ begin().at(meta).pos() ];
-    }
-    
-    template <typename T, class... SRanges>
-    const T* ConstSlice<T,SRanges...>::data() const
-    {
-	return mData;
-    }
-
-    template <typename T, class... SRanges>
-    bool ConstSlice<T,SRanges...>::isSlice() const
-    {
-	return true;
-    }
-
-    template <typename T, class... SRanges>
-    auto ConstSlice<T,SRanges...>::begin() const -> ConstSlice<T,SRanges...>::IType 
-    {
-	IType i(*MAB::mProtoI);
-	i = 0;
-	//i = mStartPos;
-	return i.setData(data());
-    }
-
-    template <typename T, class... SRanges>
-    auto ConstSlice<T,SRanges...>::end() const -> ConstSlice<T,SRanges...>::IType 
-    {
-	IType i(*MAB::mProtoI);
-	i = i.max(); // CHECK !!!
-	//i = std::get<sizeof...(SRanges)>(mBlockSizes);
-	return i.setData(data());
-    }
-
-    template <typename T, class... SRanges>
-    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > ConstSlice<T,SRanges...>::anonymous(bool slice) const
-    {
-	assert(0); // think about carefully!!!!
-	return nullptr;
-    }
-    
-    template <typename T, class... SRanges>
-    auto ConstSlice<T,SRanges...>::define(const std::shared_ptr<typename SRanges::IndexType>&... inds)
-	-> ConstSliceDef<T,SRanges...>
-    {
-	return ConstSliceDef<T,SRanges...>(*this, inds...);
-    }
-
-    
-    /**************
-     *   Slice    *
-     **************/
-
-    template <typename T, class... SRanges>
-    void Slice<T,SRanges...>::format(const std::array<size_t,sizeof...(SRanges)+1>& blocks)
-    {
-	MAB::mProtoI->format(blocks);
-    }
-    
-    template <typename T, class... SRanges>
-    Slice<T,SRanges...>::Slice(const std::shared_ptr<SRanges>&... ranges, T* data) :
-	MutableMultiArrayBase<T,SRanges...>(ranges...),
-	mData(data) {}
-    
-    template <typename T, class... SRanges>
-    const T& Slice<T,SRanges...>::operator[](const IType& i) const
-    {
-	//assert(i.sliceMode()); // -> compare objects !!!!!
-	assert(i.container() == reinterpret_cast<std::intptr_t>(this));
-	return mData[ i.pos() ];
-    }
-    
-    template <typename T, class... SRanges>
-    T& Slice<T,SRanges...>::operator[](const IType& i)
-    {
-	//assert(i.sliceMode());
-	assert(i.container() == reinterpret_cast<std::intptr_t>(this));
-	return mData[ i.pos() ];
-    }
-    
-    template <typename T, class... SRanges>
-    const T& Slice<T,SRanges...>::at(const typename IType::MetaType& meta) const
-    {
-	//auto x = begin().at(meta);
-	//VCHECK(x.pos());
-	return mData[ begin().at(meta).pos() ];
-    }
-    
-    template <typename T, class... SRanges>
-    T& Slice<T,SRanges...>::at(const typename IType::MetaType& meta)
-    {
-	//auto x = begin().at(meta);
-	//VCHECK(x.pos());
-	return mData[ begin().at(meta).pos() ];
-    }
-    
-    template <typename T, class... SRanges>
-    const T* Slice<T,SRanges...>::data() const
-    {
-	return mData;
-    }
-
-    template <typename T, class... SRanges>
-    T* Slice<T,SRanges...>::data()
-    {
-	return mData;
-    }
-
-    template <typename T, class... SRanges>
-    bool Slice<T,SRanges...>::isSlice() const
-    {
-	return true;
-    }
-
-    template <typename T, class... SRanges>
-    auto Slice<T,SRanges...>::begin() const -> Slice<T,SRanges...>::IType 
-    {
-	IType i(*MAB::mProtoI);
-	i = 0;
-	//i = mStartPos;
-	return i.setData(data());
-    }
-
-    template <typename T, class... SRanges>
-    auto Slice<T,SRanges...>::end() const -> Slice<T,SRanges...>::IType 
-    {
-	IType i(*MAB::mProtoI);
-	i = i.max(); // CHECK !!!
-	//i = std::get<sizeof...(SRanges)>(mBlockSizes);
-	return i.setData(data());
-    }
-
-    template <typename T, class... SRanges>
-    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > Slice<T,SRanges...>::anonymous(bool slice) const
-    {
-	assert(0); // think about carefully!!!!
-	return nullptr;
-    }
-    
-    template <typename T, class... SRanges>
-    std::shared_ptr<MultiArrayBase<T,AnonymousRange> > Slice<T,SRanges...>::anonymousMove()
-    {
-	assert(0); // think about carefully!!!!
-	return nullptr;
-    }
-    
-    template <typename T, class... SRanges>
-    auto Slice<T,SRanges...>::define(const std::shared_ptr<typename SRanges::IndexType>&... inds)
-	-> SliceDef<T,SRanges...>
-    {
-	return SliceDef<T,SRanges...>(*this, inds...);
-    }
-    
-
-    template <typename T, class... SRanges>
-    SliceDef<T,SRanges...>::SliceDef(Slice<T,SRanges...>& sl,
-				     const std::shared_ptr<typename SRanges::IndexType>&... inds) :
-	mIndex(sl.begin()),
-	mSl(sl)
-    {
-	mIndex(inds...);
-    }
-
-    template <typename T, class... SRanges>
-    template <class... ORanges>
-    SliceDef<T,SRanges...>& SliceDef<T,SRanges...>::operator=(const OperationRoot<T,ORanges...>& op)
-    {
-	std::array<size_t,sizeof...(SRanges)+1> blocks;
-	PackNum<sizeof...(SRanges)-1>::
-	    template mkSliceBlocks<T,OperationRoot<T,ORanges...>,SRanges...>(blocks, mIndex, op);
-	mSl.format(blocks);
-	mSl.mData = op.data();
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    ConstSliceDef<T,SRanges...>::ConstSliceDef(ConstSlice<T,SRanges...>& sl,
-					       const std::shared_ptr<typename SRanges::IndexType>&... inds) :
-	mIndex(sl.begin()),
-	mSl(sl)
-    {
-	mIndex(inds...);
-    }
-
-    template <typename T, class... SRanges>
-    template <class... ORanges>
-    ConstSliceDef<T,SRanges...>& ConstSliceDef<T,SRanges...>::operator=(const ConstOperationRoot<T,ORanges...>& op)
-    {
-	std::array<size_t,sizeof...(SRanges)+1> blocks;
-	PackNum<sizeof...(SRanges)-1>::
-	    template mkSliceBlocks<T,ConstOperationRoot<T,ORanges...>,SRanges...>(blocks, mIndex, op);
-	mSl.format(blocks);
-	mSl.mData = op.data();
-	return *this;
-    }
-
-    template <typename T, class... SRanges>
-    template <class... ORanges>
-    ConstSliceDef<T,SRanges...>& ConstSliceDef<T,SRanges...>::operator=(const OperationRoot<T,ORanges...>& op)
-    {
-	std::array<size_t,sizeof...(SRanges)+1> blocks;
-	PackNum<sizeof...(SRanges)-1>::
-	    template mkSliceBlocks<T,OperationRoot<T,ORanges...>,SRanges...>(blocks, mIndex, op);
-	mSl.format(blocks);
-	mSl.mData = op.data();
-	return *this;
-    }
-
-
-} // end namespace MultiArrayTools
 
 #endif
diff --git a/src/include/xfor/xfor.h b/src/include/xfor/xfor.h
index 38b5db9..9382f49 100644
--- a/src/include/xfor/xfor.h
+++ b/src/include/xfor/xfor.h
@@ -84,42 +84,6 @@ namespace MultiArrayHelper
 	}
     };
 
-    template <class Expr>
-    class ExpressionHolder : public ExpressionBase
-    {
-    private:
-	ExpressionHolder() = default;
-	
-        Expr mExpr;
-	typedef decltype(mExpr.rootSteps()) ExtType;
-	ExtType mExt;
-
-        mutable ExtType mRootSteps;
-
-    public:
-	typedef ExpressionBase EB;
-	
-	static constexpr size_t LAYER = Expr::LAYER + 1;
-	static constexpr size_t SIZE = Expr::SIZE;
-
-	ExpressionHolder(const ExpressionHolder& in) = default;
-	ExpressionHolder(ExpressionHolder&& in) = default;
-	ExpressionHolder& operator=(const ExpressionHolder& in) = default;
-	ExpressionHolder& operator=(ExpressionHolder&& in) = default;
-
-	ExpressionHolder(Expr expr);
-	
-	inline void operator()(size_t mlast, DExt last) const override final;
-	inline void operator()(size_t mlast, ExtType last) const;
-	inline void operator()(size_t mlast = 0) const override final;
-
-        DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
-        DExt dExtension() const override final;
-
-        auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
-	auto extension() const -> ExtType;
-
-    };
     
     template <class IndexClass, class Expr>
     class SingleExpression : public ExpressionBase
@@ -218,24 +182,27 @@ namespace MultiArrayHelper
     template <>
     inline size_t exceptMax<1>(size_t max) { return 1; }
 
-    class DynamicalExpression : public ExpressionBase
+    class DynamicExpression : public ExpressionBase
     {
     private:
-	DynamicalExpression() = default;
+	DynamicExpression() = default;
 
 	std::shared_ptr<ExpressionBase> mNext;
 
     public:
 	
-	DynamicalExpression(const DynamicalExpression& in) = default;
-	DynamicalExpression(DynamicalExpression&& in) = default;
-	DynamicalExpression& operator=(const DynamicalExpression& in) = default;
-	DynamicalExpression& operator=(DynamicalExpression&& in) = default;
+	DynamicExpression(const DynamicExpression& in) = default;
+	DynamicExpression(DynamicExpression&& in) = default;
+	DynamicExpression& operator=(const DynamicExpression& in) = default;
+	DynamicExpression& operator=(DynamicExpression&& in) = default;
 	
-	DynamicalExpression(const std::shared_ptr<ExpressionBase>& next) :
+	DynamicExpression(const std::shared_ptr<ExpressionBase>& next) :
 	    mNext(next)
 	{}
 
+        template <class Expr>
+        DynamicExpression(Expr ex) : mNext( std::make_shared<Expr>(ex) ) {}
+        
 	inline void operator()(size_t mlast, DExt last) const override final;
 	inline void operator()(size_t mlast = 0) const override final;
 
@@ -243,7 +210,45 @@ namespace MultiArrayHelper
         inline DExt dExtension() const override final;
 
     };
-    
+
+
+    template <class Expr>
+    class ExpressionHolder : public ExpressionBase
+    {
+    private:
+	ExpressionHolder() = default;
+	
+        DynamicExpression mExpr;
+	typedef decltype(std::declval<Expr>().rootSteps()) ExtType;
+	ExtType mExt;
+
+        mutable ExtType mRootSteps;
+
+    public:
+	typedef ExpressionBase EB;
+	
+	static constexpr size_t LAYER = Expr::LAYER + 1;
+	static constexpr size_t SIZE = Expr::SIZE;
+
+	ExpressionHolder(const ExpressionHolder& in) = default;
+	ExpressionHolder(ExpressionHolder&& in) = default;
+	ExpressionHolder& operator=(const ExpressionHolder& in) = default;
+	ExpressionHolder& operator=(ExpressionHolder&& in) = default;
+
+	ExpressionHolder(DynamicExpression expr);
+	
+	inline void operator()(size_t mlast, DExt last) const override final;
+	inline void operator()(size_t mlast, ExtType last) const;
+	inline void operator()(size_t mlast = 0) const override final;
+
+        DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
+        DExt dExtension() const override final;
+
+        auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
+	auto extension() const -> ExtType;
+
+    };
+
 } // namespace MultiArrayHelper
 
 /* ========================= *
@@ -429,25 +434,25 @@ namespace MultiArrayHelper
     }
     
     /***************************
-     *   DynamicalExpression   *
+     *   DynamicExpression   *
      ***************************/
 
-    inline void DynamicalExpression::operator()(size_t mlast, DExt last) const
+    inline void DynamicExpression::operator()(size_t mlast, DExt last) const
     {
 	(*mNext)(mlast,last);
     }
 
-    inline void DynamicalExpression::operator()(size_t mlast) const
+    inline void DynamicExpression::operator()(size_t mlast) const
     {
 	(*mNext)(mlast);
     }
 
-    inline DExt DynamicalExpression::dRootSteps(std::intptr_t iPtrNum) const
+    inline DExt DynamicExpression::dRootSteps(std::intptr_t iPtrNum) const
     {
         return mNext->dRootSteps(iPtrNum);
     }
 
-    inline DExt DynamicalExpression::dExtension() const
+    inline DExt DynamicExpression::dExtension() const
     {
         return mNext->dExtension();
     }
@@ -457,7 +462,7 @@ namespace MultiArrayHelper
      ************************/
 
     template <class Expr>
-    ExpressionHolder<Expr>::ExpressionHolder(Expr expr) : mExpr(expr) {}
+    ExpressionHolder<Expr>::ExpressionHolder(DynamicExpression expr) : mExpr(expr) {}
 
     template <class Expr>
     inline void ExpressionHolder<Expr>::operator()(size_t mlast, DExt last) const
@@ -468,7 +473,10 @@ namespace MultiArrayHelper
     template <class Expr>
     inline void ExpressionHolder<Expr>::operator()(size_t mlast, ExtType last) const
     {
-	mExpr(mlast,last);
+	mExpr(mlast,
+              std::make_pair<size_t const*,size_t>
+              (reinterpret_cast<size_t const*>(&last),
+               sizeof(ExtType)/sizeof(size_t)));
     }
 
     template <class Expr>
@@ -493,14 +501,14 @@ namespace MultiArrayHelper
     auto ExpressionHolder<Expr>::rootSteps(std::intptr_t iPtrNum) const
 	-> ExtType
     {
-	return mExpr.rootSteps(iPtrNum);
+	return *reinterpret_cast<ExtType*>( mExpr.dRootSteps(iPtrNum).first );
     }
 
     template <class Expr>
     auto ExpressionHolder<Expr>::extension() const
 	-> ExtType
     {
-	return mExpr.extension();
+	return *reinterpret_cast<ExtType*>( mExpr.dExtension().first );
     }
 
 
diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt
index 0bda828..3c86dfb 100644
--- a/src/lib/CMakeLists.txt
+++ b/src/lib/CMakeLists.txt
@@ -2,7 +2,7 @@
 set(libmultiarray_a_SOURCES
   ${CMAKE_SOURCE_DIR}/src/lib/ranges/range_base.cc
   ${CMAKE_SOURCE_DIR}/src/lib/ranges/anonymous_range.cc
-  ${CMAKE_SOURCE_DIR}/src/lib/ranges/dynamic_range.cc
+  ${CMAKE_SOURCE_DIR}/src/lib/ranges/dynamic_meta.cc
   ${CMAKE_SOURCE_DIR}/src/lib/ranges/multi_range_factory_product_map.cc
   ${CMAKE_SOURCE_DIR}/src/lib/map_range_factory_product_map.cc
   )
diff --git a/src/lib/ranges/dynamic_meta.cc b/src/lib/ranges/dynamic_meta.cc
new file mode 100644
index 0000000..8b2c6b8
--- /dev/null
+++ b/src/lib/ranges/dynamic_meta.cc
@@ -0,0 +1,47 @@
+
+#include <cstdlib>
+#include <vector>
+#include "ranges/dynamic_meta.h"
+
+namespace MultiArrayTools
+{
+
+    /*********************
+     *   DynamicMetaT    *
+     *********************/
+
+    bool DynamicMetaT::operator==(const DynamicMetaT& in) const
+    {
+        if(in.size() != mMeta.size()) { return false; }
+        for(size_t i = 0; i != mMeta.size(); ++i){
+            if(in[i].second != mMeta[i].second) { return false; }
+            for(size_t j = 0; j != mMeta[i].second; ++j){
+                if(in[i].first[j] != mMeta[i].first[j]) { return false; }
+            }
+        }
+        return true;
+    }
+    
+    bool DynamicMetaT::operator!=(const DynamicMetaT& in) const
+    {
+        return not operator==(in);
+    }
+    
+    size_t DynamicMetaT::size() const
+    {
+        return mMeta.size();
+    }
+    
+    DynamicMetaElem& DynamicMetaT::operator[](size_t pos)
+    {
+        return mMeta[pos];
+    }
+    
+    const DynamicMetaElem& DynamicMetaT::operator[](size_t pos) const
+    {
+        return mMeta[pos];
+    }
+
+
+    
+} // namespace MultiArrayTools
diff --git a/src/lib/ranges/dynamic_range.cc b/src/lib/ranges/dynamic_range.cc
deleted file mode 100644
index f3d31e8..0000000
--- a/src/lib/ranges/dynamic_range.cc
+++ /dev/null
@@ -1,363 +0,0 @@
-
-#include "ranges/dynamic_range.h"
-#include "ranges/dynamic_meta.h"
-
-namespace MultiArrayTools
-{
-    /*********************
-     *   DynamicMetaT    *
-     *********************/
-
-    bool DynamicMetaT::operator==(const DynamicMetaT& in) const
-    {
-        if(in.size() != mMeta.size()) { return false; }
-        for(size_t i = 0; i != mMeta.size(); ++i){
-            if(in[i].second != mMeta[i].second) { return false; }
-            for(size_t j = 0; j != mMeta[i].second; ++j){
-                if(in[i].first[j] != mMeta[i].first[j]) { return false; }
-            }
-        }
-        return true;
-    }
-    
-    bool DynamicMetaT::operator!=(const DynamicMetaT& in) const
-    {
-        return not operator==(in);
-    }
-    
-    size_t DynamicMetaT::size() const
-    {
-        return mMeta.size();
-    }
-    
-    DynamicMetaElem& DynamicMetaT::operator[](size_t pos)
-    {
-        return mMeta[pos];
-    }
-    
-    const DynamicMetaElem& DynamicMetaT::operator[](size_t pos) const
-    {
-        return mMeta[pos];
-    }
-
-    
-    /*********************
-     *   DynamicIndex    *
-     *********************/
-
-    
-    DynamicIndex::DynamicIndex(const std::shared_ptr<DynamicRange >& range) :
-	IndexInterface<DynamicIndex,MetaType>(range, 0)//,
-	//mExplicitRangePtr(std::dynamic_pointer_cast<RangeType>(IB::mRangePtr)),
-    {}
-
-    
-    IndexType DynamicIndex::type() const
-    {
-	return IndexType::SINGLE;
-    }
-
-    
-    DynamicIndex& DynamicIndex::operator=(size_t pos)
-    {
-	IB::mPos = pos;
-	return *this;
-    }
-    
-    DynamicIndex& DynamicIndex::operator++()
-    {
-        size_t ipos = mIVec.size()-1;
-        auto& ii = mIVec[ipos].first;
-        auto& jj = mIVec[ipos-1].first;
-        while(ii->pos() == ii->max()-1 and ipos != 0) {
-            (*ii) = 0;
-            ++(*jj);
-            --ipos;
-        }
-	return *this;
-    }
-    
-    DynamicIndex& DynamicIndex::operator--()
-    {
-        size_t ipos = mIVec.size()-1;
-        auto& ii = mIVec[ipos].first;
-        auto& jj = mIVec[ipos-1].first;
-        while(ii->pos() == 0 and ipos != 0) {
-            (*ii) = ii->max()-1;
-            --(*jj);
-            --ipos;
-        }
-	return *this;
-    }
-
-    DynamicIndex& DynamicIndex::operator()(const IVecT& ivec)
-    {
-        mIVec = ivec;
-        return *this;
-    }
-    
-    int DynamicIndex::pp(std::intptr_t idxPtrNum)
-    {
-	++(*this);
-	return 1;
-    }
-
-    
-    int DynamicIndex::mm(std::intptr_t idxPtrNum)
-    {
-	--(*this);
-	return 1;
-    }
-
-    
-    std::string DynamicIndex::stringMeta() const
-    {
-	return std::dynamic_pointer_cast<DynamicRange const>( IB::mRangePtr )->stringMeta(IB::mPos);
-    }
-    
-    
-    DynamicIndex::MetaType DynamicIndex::meta() const
-    {
-	return std::dynamic_pointer_cast<DynamicRange const>( IB::mRangePtr )->get(IB::mPos);
-    }
-    
-    
-    const DynamicIndex::MetaType* DynamicIndex::metaPtr() const
-    {
-	return nullptr;
-    }
-
-    /*
-    bool DynamicIndex::isMeta(const MetaType& metaPos) const
-    {
-	return mExplicitRangePtr->isMeta(metaPos);
-        }*/
-    
-    
-    DynamicIndex& DynamicIndex::at(const MetaType& metaPos)
-    {
-	(*this) = std::dynamic_pointer_cast<DynamicRange const>( IB::mRangePtr )->getMeta( metaPos );
-	return *this;
-    }
-
-    
-    size_t DynamicIndex::posAt(const MetaType& metaPos) const
-    {
-	return std::dynamic_pointer_cast<DynamicRange const>( IB::mRangePtr )->getMeta( metaPos );
-    }
-
-    
-    size_t DynamicIndex::dim() const // = 1
-    {
-	return 1;
-    }
-
-    
-    bool DynamicIndex::last() const
-    {
-	return IB::mPos == IB::mMax - 1;
-    }
-
-    
-    bool DynamicIndex::first() const
-    {
-	return IB::mPos == 0;
-    }
-
-    const IndexW& DynamicIndex::get(size_t n) const
-    {
-        return *mIVec[n].first;
-    }
-    
-    std::shared_ptr<typename DynamicIndex::RangeType> DynamicIndex::range()
-    {
-	return std::dynamic_pointer_cast<RangeType>( IB::mRangePtr );
-    }
-
-    
-    template <size_t N>
-    void DynamicIndex::getPtr() {}
-    
-    
-    size_t DynamicIndex::getStepSize(size_t n) const
-    {
-	return mIVec[n].second;
-    }
-
-    std::string DynamicIndex::id() const
-    {
-	return std::string("dyn") + std::to_string(IB::mId);
-    }
-
-    
-    void DynamicIndex::print(size_t offset)
-    {
-	if(offset == 0){
-	    std::cout << " === " << std::endl;
-	}
-	for(size_t j = 0; j != offset; ++j) { std::cout << "\t"; }
-	std::cout << id() << "[" << reinterpret_cast<std::intptr_t>(this)
-		  << "](" << IB::mRangePtr << "): " /*<< meta().first*/ << std::endl;
-    }
-    
-    
-    /******************************
-     *   DynamicRangeFactory    *
-     ******************************/
-
-    DynamicRangeFactory::DynamicRangeFactory()
-    {
-	mProd = std::shared_ptr<oType>( new DynamicRange() );
-    }
-    
-    std::map<std::shared_ptr<RangeBase>,std::vector<std::intptr_t> > DynamicRangeFactory::mAleadyCreated;
-    
-    std::shared_ptr<RangeBase> DynamicRangeFactory::checkIfCreated(const std::vector<std::shared_ptr<RangeBase> >& pvec)
-    {
-	std::shared_ptr<RangeBase> out;
-	bool check = false;
-	for(auto& x: mAleadyCreated){
-	    if(x.second.size() == pvec.size()){
-		check = true;
-		for(size_t i = 0; i != x.second.size(); ++i){
-		    if(x.second[i] != reinterpret_cast<std::intptr_t>( pvec[i].get() ) ){
-			check = false;
-			break;
-		    }
-		}
-		if(check == true){
-		    out = x.first;
-		    break;
-		}
-	    }
-	}
-	if(not check){
-	    std::vector<std::intptr_t> app(pvec.size());
-	    for(size_t i = 0; i != app.size(); ++i){
-		app[i] = reinterpret_cast<std::intptr_t>( pvec[i].get() );
-	    }
-	    mAleadyCreated[mProd] = app;
-	    out = mProd;
-	}
-	return out;
-    }
-
-    
-    std::shared_ptr<RangeBase> DynamicRangeFactory::create()
-    {
-	mProd = checkIfCreated(std::dynamic_pointer_cast<DynamicRange>(mProd)->mOrig); 
-	setSelf();
-	mProductCreated = true;
-	return mProd;
-    }
-
-    /***********************
-     *   DynamicRange    *
-     ***********************/
-    
-    DynamicRange::MetaType DynamicRange::get(size_t pos) const
-    {
-	return MetaType(); // !!!
-    }
-        
-    size_t DynamicRange::getMeta(const MetaType& metaPos) const
-    {
-	return 0; // !!!
-    }
-        
-    size_t DynamicRange::size() const
-    {
-	return mSize;
-    }
-        
-    size_t DynamicRange::dim() const
-    {
-        return mOrig.size();
-    }
-
-    SpaceType DynamicRange::spaceType() const
-    {
-	return SpaceType::DYN;
-    }
-
-    bool DynamicRange::isEmpty() const
-    {
-	return mEmpty;
-    }
-    
-    std::string DynamicRange::stringMeta(size_t pos) const
-    {
-	std::string out = "[ ";
-	//size_t xpos = pos;
-	for(size_t i = mOrig.size(); i != 0; --i) {
-	    auto& x = mOrig[i-1];
-	    const size_t redpos = pos % x->size();
-	    out = ( (i == mOrig.size()) ? out : out + " , " ) + x->stringMeta(redpos);
-	    pos -= redpos;
-	    pos /= x->size();
-	}
-	out += " ]";
-	return out;
-    }
-
-    std::vector<char> DynamicRange::data() const
-    {
-	DataHeader h;
-	h.spaceType = static_cast<int>( SpaceType::DYN );
-	h.metaSize = mOrig.size();
-	h.multiple = 1;
-	std::vector<char> out;
-	char* hcp = reinterpret_cast<char*>(&h);
-	out.insert(out.end(), hcp, hcp + sizeof(DataHeader));
-	for(auto& x: mOrig){
-	    auto part = x->data();
-	    out.insert(out.end(), part.begin(), part.end());
-	}
-	return out;
-    }
-    
-    typename DynamicRange::IndexType DynamicRange::begin() const
-    {
-	DynamicIndex i
-	    (std::dynamic_pointer_cast<DynamicRange>
-	     ( std::shared_ptr<RangeBase>(RB::mThis) ) );
-	i = 0;
-	return i;
-    }
-    
-    typename DynamicRange::IndexType DynamicRange::end() const
-    {
-	DynamicIndex i
-	    (std::dynamic_pointer_cast<DynamicRange>
-	     ( std::shared_ptr<RangeBase>(RB::mThis) ) );
-	i = size();
-	return i;
-    }
-
-    
-    std::shared_ptr<RangeBase> DynamicRange::sub(size_t num) const
-    {
-	return mOrig.at(num);
-    }
-
-    void DynamicRange::sreplace(const std::shared_ptr<RangeBase> in, size_t num)
-    {
-	assert(mOrig[num]->size() == in->size());
-	mOrig[num] = in;
-    }
-    
-    /*****************
-     *   Functions   *
-     *****************/
-    /*
-    std::shared_ptr<DynamicRange> defaultRange(size_t size )
-    {
-	DynamicRangeFactory arf
-	    ( std::dynamic_pointer_cast<DynamicRange>
-	      (DynamicRange::factory().create() ) );
-	return std::dynamic_pointer_cast<DynamicRange>( arf.create() );
-    }
-    */
-
-    
-} // end namespace MultiArrayTools
diff --git a/src/tests/ranges/anonymous_unit_test.cc b/src/tests/ranges/anonymous_unit_test.cc
index dd77c81..b12065e 100644
--- a/src/tests/ranges/anonymous_unit_test.cc
+++ b/src/tests/ranges/anonymous_unit_test.cc
@@ -5,8 +5,7 @@
 #include <iostream>
 
 #include "ranges/rheader.h"
-
-//#include "multi_array_header.h"
+#include "multi_array_header.h"
 
 namespace MAT = MultiArrayTools;
 
@@ -14,6 +13,8 @@ namespace {
 
     using namespace MAT;
 
+    typedef Expressions1 EC1;
+    
     template <class Factory, typename T>
     void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
     {
@@ -125,15 +126,15 @@ namespace {
 
     TEST_F(AnonymousTest, DCast1)
     {
-	DynamicRangeFactory arf1(sr1ptr,m3rptr);
-	auto ar1a = std::dynamic_pointer_cast<DynamicRange>( arf1.create() );
+	DynamicRangeFactory<EC1> arf1(sr1ptr,m3rptr);
+	auto ar1a = std::dynamic_pointer_cast<DynamicRange<EC1>>( arf1.create() );
 
 	arf1.append(sr2ptr);
 	
-	auto ar1b = std::dynamic_pointer_cast<DynamicRange>( arf1.create() );
+	auto ar1b = std::dynamic_pointer_cast<DynamicRange<EC1>>( arf1.create() );
 	
-	DynamicRangeFactory arf2(sr1ptr,m3rptr,sr2ptr);
-	auto ar2 = std::dynamic_pointer_cast<DynamicRange>( arf2.create() );
+	DynamicRangeFactory<EC1> arf2(sr1ptr,m3rptr,sr2ptr);
+	auto ar2 = std::dynamic_pointer_cast<DynamicRange<EC1>>( arf2.create() );
 
 	EXPECT_EQ(ar1b.get(), ar2.get());
 	EXPECT_EQ(ar1a->size(), sr1ptr->size() * m3rptr->size());
@@ -149,9 +150,9 @@ namespace {
 
     TEST_F(AnonymousTest, DCast2)
     {
-	DynamicRangeFactory arf2(sr1ptr,m3rptr,sr2ptr);
-	auto ar = std::dynamic_pointer_cast<DynamicRange>( arf2.create() );
-	auto mr = ar->template scast<SRange,DynamicRange>(1,2);
+	DynamicRangeFactory<EC1> arf2(sr1ptr,m3rptr,sr2ptr);
+	auto ar = std::dynamic_pointer_cast<DynamicRange<EC1>>( arf2.create() );
+	auto mr = ar->template scast<SRange,DynamicRange<EC1>>(1,2);
 
 	EXPECT_EQ(mr->template getPtr<0>()->size(), sr1ptr->size());
 	EXPECT_EQ(mr->template getPtr<1>()->size(), m3rptr->size() * sr2ptr->size());