diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7606be5..64c01a9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8)
 
 project(multi_array)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++14 -g -Wpedantic -Ofast -march=native -faligned-new -funroll-loops -fopenmp")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++14 -g -Wpedantic -Ofast -march=native -faligned-new -funroll-loops -fopenmp -ftemplate-backtrace-limit=0")
 #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -std=c++14 -g -Wpedantic -O0 -march=native -faligned-new -funroll-loops -fopenmp")
 
 enable_testing()
diff --git a/src/include/dynamic_operation.cc.h b/src/include/dynamic_operation.cc.h
index b2353a5..4ecd774 100644
--- a/src/include/dynamic_operation.cc.h
+++ b/src/include/dynamic_operation.cc.h
@@ -4,20 +4,23 @@
 namespace MultiArrayTools
 {
     template <typename T, class Operation>
-    const T& DynamicOperation<T,Operation>::get(const DExt& pos) const
+    const T& DynamicOperation<T,Operation>::get(const DExtT& pos) const
     {
+	VCHECK(pos.size());
+	VCHECK(sizeof(ET)/sizeof(size_t));
+	assert(0);
 	return mOp.get(pos.expl<ET>());
     }
     
     template <typename T, class Operation>
-    DynamicOperationBase& DynamicOperation<T,Operation>::set(const DExt& pos)
+    DynamicOperationBase<T>& DynamicOperation<T,Operation>::set(const DExtT& pos)
     {
 	mOp.set(pos.expl<ET>());
 	return *this;
     }
     
     template <typename T, class Operation>
-    DExt DynamicOperation<T,Operation>::rootSteps(std::intptr_t iPtrNum = 0) const
+    DExtT DynamicOperation<T,Operation>::rootSteps(std::intptr_t iPtrNum) const
     {
 	return mOp.rootSteps(iPtrNum);
     }
@@ -34,4 +37,5 @@ namespace MultiArrayTools
 	return mOp.data();
     }
 
+
 } // namespace MultiArrayTools
diff --git a/src/include/dynamic_operation.h b/src/include/dynamic_operation.h
index f8cab39..a8f2738 100644
--- a/src/include/dynamic_operation.h
+++ b/src/include/dynamic_operation.h
@@ -24,9 +24,9 @@ namespace MultiArrayTools
 	DynamicOperationBase& operator=(const DynamicOperationBase& in) = default;
 	DynamicOperationBase& operator=(DynamicOperationBase&& in) = default;
 
-	virtual const T& get(const DExt& pos) const = 0;
-	virtual DynamicOperationBase& set(const DExt& pos) = 0;
-	virtual DExt rootSteps(std::intptr_t iPtrNum = 0) const = 0;
+	virtual const T& get(const DExtT& pos) const = 0;
+	virtual DynamicOperationBase& set(const DExtT& pos) = 0;
+	virtual DExtT rootSteps(std::intptr_t iPtrNum = 0) const = 0;
 	virtual DynamicExpression loop(const DynamicExpression& exp) const = 0;
 	virtual const T* data() const = 0;
 	
@@ -48,12 +48,41 @@ namespace MultiArrayTools
 
 	DynamicOperation(const Operation& op) : mOp(op) {}
 
-	virtual const T& get(const DExt& pos) const override final;
-	virtual DynamicOperationBase& set(const DExt& pos) override final;
-	virtual DExt rootSteps(std::intptr_t iPtrNum = 0) const override final;
+	virtual const T& get(const DExtT& pos) const override final;
+	virtual DynamicOperationBase<T>& set(const DExtT& pos) override final;
+	virtual DExtT rootSteps(std::intptr_t iPtrNum = 0) const override final;
 	virtual DynamicExpression loop(const DynamicExpression& exp) const override final;
 	virtual const T* data() const override final;
     };
+
+    template <typename T>
+    class DynamicO : public OperationTemplate<T,DynamicO<T>>
+    {
+    private:
+	// NOT THREAD SAFE!!!
+	std::shared_ptr<DynamicOperationBase<T>> mOp;
+    public:
+	typedef T value_type;
+	typedef OperationBase<T,DynamicO<T>> OT;
+	
+	static constexpr size_t SIZE = 1;
+	static constexpr bool CONT = true;
+	
+	DynamicO() = default;
+	DynamicO(const DynamicO& in) = default;
+	DynamicO(DynamicO&& in) = default;
+	DynamicO& operator=(const DynamicO& in) = default;
+	DynamicO& operator=(DynamicO&& in) = default;
+
+	template <class Op>
+	DynamicO(const Op& op) : mOp(std::make_shared<DynamicOperation<T,Op>>(op)) {}
+	
+	inline const T& get(const DExtT& pos) const { return mOp->get(pos); }
+	inline DynamicO& set(const DExtT& pos) { return mOp->set(pos); }
+	inline DExtT rootSteps(std::intptr_t iPtrNum = 0) const { return mOp->rootSteps(iPtrNum); }
+	inline DynamicExpression loop(const DynamicExpression& exp) const { return mOp->loop(exp); }
+	inline const T* data() const { return mOp->data(); }
+    };
     
 } // namespace MultiArrayTools
 
diff --git a/src/include/helper_tools.cc.h b/src/include/helper_tools.cc.h
index b68f178..7a07382 100644
--- a/src/include/helper_tools.cc.h
+++ b/src/include/helper_tools.cc.h
@@ -45,6 +45,13 @@ namespace MultiArrayTools
 	return mi;
     }
 
+    template <class Index>
+    auto mkIndexW(const std::shared_ptr<Index>& ind)
+	-> std::shared_ptr<IndexW>
+    {
+	return std::make_shared<IndexWrapper<Index>>(ind);
+    }
+
     template <SpaceType STYPE, class Op, class MA, class... RangeTypes>
     auto mkGenMapR(const std::tuple<Op,MA>& f, std::shared_ptr<RangeTypes>... ranges)
 	-> std::shared_ptr<GenMapRange<Op,STYPE,RangeTypes...> >
diff --git a/src/include/helper_tools.h b/src/include/helper_tools.h
index 0de5a6a..4d99d20 100644
--- a/src/include/helper_tools.h
+++ b/src/include/helper_tools.h
@@ -61,6 +61,10 @@ namespace MultiArrayTools
     auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
 	-> decltype( getIndex( mkMulti( indices.range()... ) ) );
 
+    template <class Index>
+    auto mkIndexW(const std::shared_ptr<Index>& ind)
+	-> std::shared_ptr<IndexW>;
+
     template <class... RangeTypes>
     auto mkMulti(std::tuple<std::shared_ptr<RangeTypes>...> rangesTuple)
 	-> MultiRange<RangeTypes...>;
diff --git a/src/include/multi_array_header.cc.h b/src/include/multi_array_header.cc.h
index 1412aa9..c394ebf 100644
--- a/src/include/multi_array_header.cc.h
+++ b/src/include/multi_array_header.cc.h
@@ -6,6 +6,7 @@
 #include "multi_array_base.cc.h"
 #include "multi_array.cc.h"
 #include "slice.cc.h"
+#include "dynamic_operation.cc.h"
 //#include "expressions.cc.h"
 
 #include "ranges/dynamic_range.cc.h"
diff --git a/src/include/multi_array_header.h b/src/include/multi_array_header.h
index 4ab4151..95ffbc5 100644
--- a/src/include/multi_array_header.h
+++ b/src/include/multi_array_header.h
@@ -11,6 +11,7 @@
 #include "helper_tools.h"
 #include "operation_def.h"
 #include "map_range.h"
+#include "dynamic_operation.h"
 //#include "expressions.h"
 
 #include "multi_array_header.cc.h"
diff --git a/src/include/xfor/exttype.h b/src/include/xfor/exttype.h
index 50ecbb6..5744983 100644
--- a/src/include/xfor/exttype.h
+++ b/src/include/xfor/exttype.h
@@ -138,8 +138,6 @@ namespace MultiArrayHelper
 	template <class Z>
 	inline MExt(size_t y, const Z& z);
 
-
-
 	template <class Y, class Z>
 	inline MExt(const Y& y, const Z& z);
 
diff --git a/src/include/xfor/xfor.h b/src/include/xfor/xfor.h
index 59ec7a6..53a1d20 100644
--- a/src/include/xfor/xfor.h
+++ b/src/include/xfor/xfor.h
@@ -26,16 +26,17 @@ namespace MultiArrayHelper
     class ExtBase
     {
     public:
-        ExtBase() = default;
-        ExtBase(const ExtBase& in) = default;
-        ExtBase(ExtBase&& in) = default;
-        ExtBase& operator=(const ExtBase& in) = default;
-        ExtBase& operator=(ExtBase&& in) = default;
+	ExtBase() = default;
+	ExtBase(const ExtBase& in) = default;
+	ExtBase(ExtBase&& in) = default;
+	ExtBase& operator=(const ExtBase& in) = default;
+	ExtBase& operator=(ExtBase&& in) = default;
         
-        virtual size_t size() const = 0;
-        //virtual size_t rootSteps() const = 0;
-        virtual std::shared_ptr<ExtBase> operator+(const ExtBase& in) const = 0;
-        virtual std::shared_ptr<ExtBase> operator*(size_t in) const = 0;
+	virtual size_t size() const = 0;
+	virtual const size_t& val() const = 0;
+	//virtual size_t rootSteps() const = 0;
+	virtual std::shared_ptr<ExtBase> operator+(const ExtBase& in) const = 0;
+	virtual std::shared_ptr<ExtBase> operator*(size_t in) const = 0;
 	
 	template <class ExtType>
 	const ExtType& expl() const;
@@ -48,49 +49,82 @@ namespace MultiArrayHelper
     class ExtT : public ExtBase
     {
     private:
-        ExtType mExt;
+	ExtType mExt;
     public:
-        ExtT() = default;
-        ExtT(const ExtT& in) = default;
-        ExtT(ExtT&& in) = default;
-        ExtT& operator=(const ExtT& in) = default;
-        ExtT& operator=(ExtT&& in) = default;
+	static constexpr size_t SIZE = ExtType::SIZE;
+	static constexpr size_t NUM = ExtType::NUM;
+	
+	ExtT() = default;
+	ExtT(const ExtT& in) = default;
+	ExtT(ExtT&& in) = default;
+	ExtT& operator=(const ExtT& in) = default;
+	ExtT& operator=(ExtT&& in) = default;
 
-        ExtT(const ExtType& in) : mExt(in) {}
+	ExtT(const ExtType& in) : mExt(in) {}
         
-        virtual size_t size() const override final { return sizeof(ExtType)/sizeof(size_t); }
-        //virtual size_t size() const override final { return ExtType::MExtSize(); }
-        //virtual size_t rootSteps() const override final;
-        const ExtType& ext() const { return mExt; }
+	virtual size_t size() const override final { return sizeof(ExtType)/sizeof(size_t); }
+	//virtual size_t size() const override final { return ExtType::MExtSize(); }
+	//virtual size_t rootSteps() const override final;
+	const ExtType& ext() const { return mExt; }
+	virtual const size_t& val() const override final { return mExt.val(); }
 
-        virtual DExt operator+(const ExtBase& in) const override final
-        { return std::make_shared<ExtT>( mExt + dynamic_cast<const ExtT&>(in).mExt ); }
-        virtual DExt operator*(size_t in) const override final
-        { return std::make_shared<ExtT>( mExt * in ); }
+	virtual DExt operator+(const ExtBase& in) const override final
+	{ return std::make_shared<ExtT>( mExt + dynamic_cast<const ExtT&>(in).mExt ); }
+	virtual DExt operator*(size_t in) const override final
+	{ return std::make_shared<ExtT>( mExt * in ); }
     };
+    //class DExtT;
 
-
-    class DExtT
+    template <class X>
+    class DExtTX
     {
     private:
-        DExt mDExt = nullptr;
+	DExt mDExt = nullptr;
+	X mNext;
     public:
-        DExtT() = default;
-        DExtT(const DExtT& in) = default;
-        DExtT(DExtT&& in) = default;
-        DExtT& operator=(const DExtT& in) = default;
-        DExtT& operator=(DExtT&& in) = default;
-        DExtT(const DExt& in) : mDExt(in) {}
-        
+	static constexpr size_t NUM = X::SIZE;
+	static constexpr size_t SIZE = NUM + 1;
+	
+        DExtTX() = default;
+        DExtTX(const DExtTX& in) = default;
+        DExtTX(DExtTX&& in) = default;
+        DExtTX& operator=(const DExtTX& in) = default;
+        DExtTX& operator=(DExtTX&& in) = default;
+        DExtTX(const DExt& in) : mDExt(in) {}
+
+	template <class Y>
+	DExtTX& operator=(const Y& y) { mDExt = std::make_shared<ExtT<Y>>(y); return *this; }
+
+	template <class Y>
+	DExtTX(const Y& y) : mDExt(std::make_shared<ExtT<Y>>(y)) {}
+
+	DExtTX(const DExt& y, const X& x) : mDExt(y),
+					    mNext(x) {}
+
         virtual size_t size() const { return mDExt->size(); }
         inline const DExt& get() const { return mDExt; }
         
-        inline DExtT operator+(const DExtT& in) const
-        { if (not mDExt) return in; else return DExtT( (*mDExt) + (*in.mDExt) ) ; }
-        inline DExtT operator*(size_t in) const
-        { if (not mDExt) return *this; else return DExtT( (*mDExt) * in ) ; }
-        
+        inline DExtTX operator+(const DExtTX& in) const
+        { if (not mDExt) return in; else return DExtTX( (*mDExt) + (*in.mDExt) ) ; }
+        inline DExtTX operator*(size_t in) const
+        { if (not mDExt) return *this; else return DExtTX( (*mDExt) * in ) ; }
+
+	template <class ExtType>
+	inline const ExtType& expl() const { return mDExt->expl<ExtType>(); }
+
+	template <class Y>
+	inline auto extend(const Y& y) const -> DExtTX<decltype(mNext.extend(y))>
+	{ return DExtTX<decltype(mNext.extend(y))>(mDExt, mNext.extend(y)); }
+
+	inline const size_t& val() const { return mDExt->val(); }
+	inline const X& next() const { return mNext; }
+
+	template <size_t N>
+	inline auto nn() const -> decltype(Getter<N>::getX(*this))
+	    { return Getter<N>::getX(*this); }
     };
+
+    typedef DExtTX<None> DExtT;
     
     inline MExt<None> mkExt(size_t s) { return MExt<None>(s); }
     
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index 42c9eac..b272f7f 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -31,6 +31,11 @@ add_dependencies(op3utest multiarray)
 target_link_libraries(op3utest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
 add_test(NAME op3utest COMMAND op3utest)
 
+add_executable(op4utest op4_unit_test.cc)
+add_dependencies(op4utest multiarray)
+target_link_libraries(op4utest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
+add_test(NAME op4utest COMMAND op4utest)
+
 add_executable(opptest op_perf_test.cc)
 add_dependencies(opptest multiarray)
 target_link_libraries(opptest multiarray)
diff --git a/src/tests/op4_unit_test.cc b/src/tests/op4_unit_test.cc
new file mode 100644
index 0000000..4623943
--- /dev/null
+++ b/src/tests/op4_unit_test.cc
@@ -0,0 +1,148 @@
+
+#include "test_header.h"
+#include <ctime>
+
+namespace
+{
+
+    template <class MA>
+    void setMARandom(MA& ma, int seed)
+    {
+	std::srand(seed);
+	double n = 0;
+	double d = 0;
+	for(auto& x: ma.vdata()){
+	    double nx = 0;
+	    double dx = 0;
+	    while(nx == x or dx == d or nx == dx or dx == 0){
+		nx = static_cast<double>(std::rand() % 1000);
+		dx = static_cast<double>(std::rand() % 1000);
+	    }
+	    n = nx;
+	    d = dx;
+	    assert(n != d);
+	    x = n/d;
+	}
+    }
+
+    template <class Op>
+    DynamicO<double> mkDynOp(const Op& op)
+    {
+	return DynamicO<double>(op);
+    }
+
+    template <class Op1, class Op2>
+    DynamicO<double> mkDynMul(const Op1& op1, const Op2& op2)
+    {
+	auto op = op1 * op2;
+	return DynamicO<double>(std::make_shared<DynamicOperation<double,decltype(op)>>(op));
+    }
+    
+    class OpTest_Dyn : public ::testing::Test
+    {
+    protected:
+
+	MultiArray<double,CR,DR> ma1;
+	MultiArray<double,CR,DR> ma2;
+	MultiArray<double,DR> ma3;
+
+	MultiArray<double,CR,DR> res1;
+
+	std::map<std::string,std::shared_ptr<IndexW>> imap;
+
+	std::shared_ptr<DR> dr1;
+	std::shared_ptr<DR> dr2;
+	std::shared_ptr<DR> dr3;
+	std::shared_ptr<DR> dr4;
+	std::shared_ptr<CR> cr1;
+	
+	OpTest_Dyn()
+	{
+	    cr1 = createRangeE<CR>(5);
+	    
+	    auto cr2 = createRangeE<CR>(7);
+	    auto cr3 = createRangeE<CR>(11);
+	    auto cr4 = createRangeE<CR>(3);
+	    auto cr5 = createRangeE<CR>(13);
+
+	    dr1 = createRangeE<DR>(cr2,cr2,cr3,cr4);
+	    dr2 = createRangeE<DR>(cr3,cr3,cr4);
+	    dr3 = createRangeE<DR>(cr2,cr5);
+
+	    dr4 = createRangeE<DR>(cr2,cr3,cr4,cr4);
+
+	    ma1 = mkArray<double>(cr1,dr1);
+	    ma2 = mkArray<double>(cr1,dr2);
+	    ma3 = mkArray<double>(dr3);
+
+	    res1 = mkArray<double>(cr1,dr4);
+
+	    setMARandom(ma1, 25);
+	    setMARandom(ma2, 31);
+	    setMARandom(ma3, 47);
+
+	    imap["i2_1"] = mkIndexW(getIndex(cr2));
+	    imap["i2_2"] = mkIndexW(getIndex(cr2));
+	    imap["i3_1"] = mkIndexW(getIndex(cr3));
+	    imap["i3_2"] = mkIndexW(getIndex(cr3));
+	    imap["i4_1"] = mkIndexW(getIndex(cr4));
+	    imap["i4_2"] = mkIndexW(getIndex(cr4));
+	    imap["i5_1"] = mkIndexW(getIndex(cr5));
+	    imap["i5_2"] = mkIndexW(getIndex(cr5));
+	}
+    };
+
+    TEST_F(OpTest_Dyn, Multiply)
+    {
+	auto i1 = getIndex(cr1);
+	auto di1 = getIndex(dr1);
+	auto di2 = getIndex(dr2);
+	auto di4 = getIndex(dr4);
+
+	(*di1)({imap["i2_1"],imap["i2_2"],imap["i3_1"],imap["i4_1"]});
+	(*di2)({imap["i3_1"],imap["i3_1"],imap["i4_2"]});
+	(*di4)({imap["i2_1"],imap["i3_1"],imap["i4_1"],imap["i4_2"]});
+
+	auto resx1 = res1;
+	res1(i1,di4) = ma1(i1,di1) * ma2(i1,di2);
+
+	resx1(i1,di4) = mkDynOp(ma1(i1,di1)) * mkDynOp(ma2(i1,di2));
+	
+	auto i2_1 = imap.at("i2_1");
+	auto i2_2 = imap.at("i2_2");
+	auto i3_1 = imap.at("i3_1");
+	auto i4_1 = imap.at("i4_1");
+	auto i4_2 = imap.at("i4_2");
+	
+	for(size_t ii1 = 0; ii1 != i1->max(); ++ii1){
+	    for(size_t ii2_1 = 0; ii2_1 != i2_1->max(); ++ii2_1){
+		size_t ii2_2 = 0;
+		for(size_t ii3_1 = 0; ii3_1 != i3_1->max(); ++ii3_1){
+		    for(size_t ii4_1 = 0; ii4_1 != i4_1->max(); ++ii4_1){
+			for(size_t ii4_2 = 0; ii4_2 != i4_2->max(); ++ii4_2){
+			    const size_t jr = (((ii1*i2_1->max() + ii2_1)*i3_1->max() + ii3_1)*i4_1->max() + ii4_1)*i4_2->max() + ii4_2;
+			    const size_t j1 = (((ii1*i2_1->max() + ii2_1)*i2_2->max() + ii2_2)*i3_1->max() + ii3_1)*i4_1->max() + ii4_1;
+			    const size_t j2 = ((ii1*i3_1->max() + ii3_1)*i3_1->max() + ii3_1)*i4_2->max() + ii4_2;
+			    auto resx = xround(res1.vdata()[jr]);
+			    //std::cout << resx << "  ";
+			    auto x12 = xround(ma1.vdata()[j1]*ma2.vdata()[j2]);
+			    EXPECT_EQ( resx, x12 );
+			}
+		    }
+		}
+	    }
+	    //std::cout << std::endl;
+	}
+    }
+
+    TEST_F(OpTest_Dyn, Contract)
+    {
+
+    }
+
+    TEST_F(OpTest_Dyn, Functional)
+    {
+
+    }
+    
+}
diff --git a/src/tests/test_header.h b/src/tests/test_header.h
index 663ed8f..0baab35 100644
--- a/src/tests/test_header.h
+++ b/src/tests/test_header.h
@@ -19,6 +19,11 @@ double xround(double arg)
     
 using namespace MAT;
 
+typedef ClassicRange CR;
+typedef ClassicRF CRF;
+typedef DynamicRange DR;
+typedef DynamicRangeFactory DRF;
+
 template <class Factory>
 void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr)
 {