diff --git a/src/base_def.h b/src/base_def.h
index 0ce3c95..c8ce733 100644
--- a/src/base_def.h
+++ b/src/base_def.h
@@ -74,11 +74,15 @@ namespace MultiArrayTools
     class MultiArray;
 
     // multi_array_operation.h
-    template <typename T, class Range>
+    template <typename T>
     class MultiArrayOperationBase;
+
+    // multi_array_operation.h
+    template <typename T, class Range>
+    class MultiArrayOperationRoot;
     
     // multi_array_operation.h
-    template <typename T, class Range, class Operation, class... Ranges>
+    template <typename T, class Operation, class... MAOps>
     class MultiArrayOperation;
 
     // slice.h
diff --git a/src/multi_array.cc b/src/multi_array.cc
index 8fdcd68..cb51356 100644
--- a/src/multi_array.cc
+++ b/src/multi_array.cc
@@ -46,9 +46,9 @@ namespace MultiArrayTools
 
     template <typename T, class Range>
     template <class... NameTypes>
-    MultiArrayOperationBase<T,Range> MultiArray<T,Range>::operator()(const NameTypes&... str)
+    MultiArrayOperationRoot<T,Range> MultiArray<T,Range>::operator()(const NameTypes&... str)
     {
-	return MultiArrayOperationBase<T,Range>(*this, Name("master", str...));
+	return MultiArrayOperationRoot<T,Range>(*this, Name("master", str...));
     }
 
     template <typename T, class Range>
diff --git a/src/multi_array.h b/src/multi_array.h
index 5ab0a9f..3f6fd66 100644
--- a/src/multi_array.h
+++ b/src/multi_array.h
@@ -26,7 +26,7 @@ namespace MultiArrayTools
 	MultiArray(const Range& range, std::vector<T>&& vec);
 	
 	template <class... NameTypes>
-	MultiArrayOperationBase<T,Range> operator()(const NameTypes&... str);
+	MultiArrayOperationRoot<T,Range> operator()(const NameTypes&... str);
 	
 	T& operator[](const typename Range::IndexType& i);
 	const T& operator[](const typename Range::IndexType& i) const;
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
index 8b259ff..e59d269 100644
--- a/src/multi_array_operation.cc
+++ b/src/multi_array_operation.cc
@@ -3,140 +3,149 @@
 
 namespace MultiArrayTools
 {
+
+    /*********************************
+     *   MultiArrayOperationRoot     *
+     *********************************/
+    
+    template <typename T>
+    MultiArrayOperationBase<T>::~MultiArrayOperationBase()
+    {
+	//delete mIibPtr;
+    }
+
+    template <typename T>
+    IndefinitIndexBase* MultiArrayOperationBase<T>::index()
+    {
+	return mIibPtr;
+    }
+    
     /*********************************
      *   MultiArrayOperationBase     *
      *********************************/
 
     template <typename T, class Range>
-    MultiArrayOperationBase<T,Range>::
-    MultiArrayOperationBase(MultiArray<T,Range>& ma,
+    MultiArrayOperationRoot<T,Range>::
+    MultiArrayOperationRoot(MultiArray<T,Range>& ma,
 			    const Name& nm) : mArrayRef(ma),
 					      mIndex(mArrayRef.begin()),
-					      mIibPtr(&mIndex),
+					      MAOB::mIibPtr(&mIndex),
 					      mNm(nm)
     {
-	mIibPtr->name(nm);
+	MAOB::mIibPtr->name(nm);
 	//mIndex.name(nm);
     }
 
     template <typename T, class Range>
-    MultiArrayOperationBase<T,Range>::~MultiArrayOperationBase()
+    MultiArrayOperationRoot<T,Range>&
+    MultiArrayOperationRoot<T,Range>::operator=(const MultiArrayOperationBase<T>& in)
     {
-	//delete mIibPtr;
-    }
-
-    template <typename T, class Range>
-    MultiArrayOperationBase<T,Range>&
-    MultiArrayOperationBase<T,Range>::operator=(const MultiArrayOperationBase<T,Range>& in)
-    {
-	in.linkIndicesTo(mIibPtr);
-	IndexType& iref = dynamic_cast<IndexType&>(*mIibPtr);
+	in.linkIndicesTo(MAOB::mIibPtr);
+	IndexType& iref = dynamic_cast<IndexType&>(*MAOB::mIibPtr);
+	//if(mArrayRef.isSlice()){
+	//    linkSlice(&in.index(), MAOB::mIibPtr);
+	//    return *this;
+	//}
 	for(iref = mArrayRef.begin().pos(); iref != mArrayRef.end(); ++iref){
 	    // build in vectorization later
 	    get() = in.get();
 	}
-	mIibPtr->freeLinked();
+	MAOB::mIibPtr->freeLinked();
 	return *this;
     }
-    
+
+    /*
     template <typename T, class Range>
     template <class Range2>
-    MultiArrayOperationBase<T,Range>&
-    MultiArrayOperationBase<T,Range>::operator=(const MultiArrayOperationBase<T, Range2>& in)
+    MultiArrayOperationRoot<T,Range>&
+    MultiArrayOperationRoot<T,Range>::operator=(const MultiArrayOperationBase<T, Range2>& in)
     {
-	in.linkIndicesTo(mIibPtr);
-	IndexType& iref = dynamic_cast<IndexType&>(*mIibPtr);
+	in.linkIndicesTo(MAOB::mIibPtr);
+	IndexType& iref = dynamic_cast<IndexType&>(*MAOB::mIibPtr);
 	if(mArrayRef.isSlice()){
-	    linkSlice(&in.index(), mIibPtr);
+	    linkSlice(&in.index(), MAOB::mIibPtr);
 	    return *this;
 	}
 	for(iref = mArrayRef.begin().pos(); iref != mArrayRef.end(); ++iref){
 	    // build in vectorization later
 	    get() = in.get();
 	}
-	mIibPtr->freeLinked();
+	MAOB::mIibPtr->freeLinked();
 	return *this;
-    }
+	}*/
 
     template <typename T, class Range>
-    template <class Operation, class... Ranges>
-    MultiArrayOperation<T,Range,Operation,Ranges...>
-    MultiArrayOperationBase<T,Range>::operator()(Operation& op, const MultiArrayOperationBase<T,Ranges>&... secs)
+    template <class Operation, class... MAOps>
+    MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
+    MultiArrayOperationRoot<T,Range>::operator()(Operation& op, const MAOps&... secs)
     {
-	return MultiArrayOperation<T,Range,Operation,Ranges...>(mArrayRef, mNm, op, secs...);
+	return MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>(op, *this, secs...);
     }
 
     template <typename T, class Range>
-    template <class Operation, class... Ranges>
-    MultiArrayOperation<T,Range,Operation,Ranges...>
-    MultiArrayOperationBase<T,Range>::operator()(const Operation& op,
-						 const MultiArrayOperationBase<T,Ranges>&... secs)
+    template <class Operation, class... MAOps>
+    MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
+    MultiArrayOperationRoot<T,Range>::operator()(const Operation& op, const MAOps&... secs)
     {
-	return MultiArrayOperation<T,Range,Operation,Ranges...>(mArrayRef, mNm, op, secs...);
+	return MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>(op, *this, secs...);
     }
 
     template <typename T, class Range>
-    template <class Range2>
-    MultiArrayOperation<T,Range,std::plus<T>,Range2>
-    MultiArrayOperationBase<T,Range>::operator+(const MultiArrayOperationBase<T,Range2>& sec)
+    template <class MAOp>
+    auto MultiArrayOperationRoot<T,Range>::operator+(const MAOp& sec)
+	-> decltype(operator()(std::plus<T>(), sec))
     {
 	return operator()(std::plus<T>(), sec);
     }
 
     template <typename T, class Range>
-    template <class Range2>
-    MultiArrayOperation<T,Range,std::minus<T>,Range2>
-    MultiArrayOperationBase<T,Range>::operator-(const MultiArrayOperationBase<T,Range2>& sec)
+    template <class MAOp>
+    auto MultiArrayOperationRoot<T,Range>::operator-(const MAOp& sec)
+	-> decltype(operator()(std::minus<T>(), sec))
     {
 	return operator()(std::minus<T>(), sec);
     }
 
     template <typename T, class Range>
-    template <class Range2>
-    MultiArrayOperation<T,Range,std::multiplies<T>,Range2>
-    MultiArrayOperationBase<T,Range>::operator*(const MultiArrayOperationBase<T,Range2>& sec)
+    template <class MAOp>
+    auto MultiArrayOperationRoot<T,Range>::operator*(const MAOp& sec)
+	-> decltype(operator()(std::multiplies<T>(), sec))
     {
 	return operator()(std::multiplies<T>(), sec);
     }
 
     template <typename T, class Range>
-    template <class Range2>
-    MultiArrayOperation<T,Range,std::divides<T>,Range2>
-    MultiArrayOperationBase<T,Range>::operator/(const MultiArrayOperationBase<T,Range2>& sec)
+    template <class MAOp>
+    auto MultiArrayOperationRoot<T,Range>::operator/(const MAOp& sec)
+	-> decltype(operator()(std::divides<T>(), sec))
     {
 	return operator()(std::divides<T>(), sec);
     }
     
     template <typename T, class Range>
-    size_t MultiArrayOperationBase<T,Range>::argNum() const
+    size_t MultiArrayOperationRoot<T,Range>::argNum() const
     {
 	return 1;
     }
 
     template <typename T, class Range>
-    IndefinitIndexBase* MultiArrayOperationBase<T,Range>::index()
+    void MultiArrayOperationRoot<T,Range>::linkIndicesTo(IndefinitIndexBase* target) const
     {
-	return mIibPtr;
+	MAOB::mIibPtr->linkTo(target);
     }
 
     template <typename T, class Range>
-    void MultiArrayOperationBase<T,Range>::linkIndicesTo(IndefinitIndexBase* target) const
-    {
-	mIibPtr->linkTo(target);
-    }
-
-    template <typename T, class Range>
-    T& MultiArrayOperationBase<T,Range>::get()
+    T& MultiArrayOperationRoot<T,Range>::get()
     {
 	//return mArrayRef[mIndex];
-	return mArrayRef[*dynamic_cast<IndexType*>(mIibPtr)];
+	return mArrayRef[*dynamic_cast<IndexType*>(MAOB::mIibPtr)];
     }
 
     template <typename T, class Range>
-    const T& MultiArrayOperationBase<T,Range>::get() const
+    const T& MultiArrayOperationRoot<T,Range>::get() const
     {
 	//return mArrayRef[mIndex];
-	return mArrayRef[*dynamic_cast<IndexType*>(mIibPtr)];
+	return mArrayRef[*dynamic_cast<IndexType*>(MAOB::mIibPtr)];
     }
 
     /*****************************
@@ -167,88 +176,127 @@ namespace MultiArrayTools
     template <size_t N>
     struct OperationCall
     {
-	template <typename T, class Operation, class Tuple, class... MBases>
-	static auto callOperation(Operation& op, const Tuple& tp, const T& first, const MBases&... secs)
-	    -> decltype(OperationCall<N-1>::template callOperation(op, tp, std::get<N>(tp), secs...))
+	template <typename T, class Operation, class Tuple, class... MAOps>
+	static auto callOperation(Operation& op, const Tuple& tp, const MAOps&... args)
+	    -> decltype(OperationCall<N-1>::template callOperation(op, tp, std::get<N>(tp), args...))
 	{
-	    return OperationCall<N-1>::template callOperation(op, tp, first, std::get<N>(tp), secs...);
+	    return OperationCall<N-1>::template callOperation(op, tp, std::get<N>(tp), args...);
 	}
 
-	template <typename T, class Operation, class Tuple, class... MBases>
-	static auto callOperation(const Operation& op, const Tuple& tp, const T& first, const MBases&... secs)
-	    -> decltype(OperationCall<N-1>::template callOperation(op, tp, std::get<N>(tp), secs...))
+	template <typename T, class Operation, class Tuple, class... MAOps>
+	static auto callOperation(const Operation& op, const Tuple& tp, const MAOps&... args)
+	    -> decltype(OperationCall<N-1>::template callOperation(op, tp, std::get<N>(tp), args...))
 	{
-	    return OperationCall<N-1>::template callOperation(op, tp, first, std::get<N>(tp), secs...);
+	    return OperationCall<N-1>::template callOperation(op, tp, std::get<N>(tp), args...);
 	}
     };
 
     template <>
     struct OperationCall<0>
     {
-	template <typename T, class Operation, class Tuple, class... MBases>
-	static auto callOperation(Operation& op, const Tuple& tp, const T& first, const MBases&... secs)
-	    -> decltype(op(first, std::get<0>(tp).get(), secs.get()...))
+	template <typename T, class Operation, class Tuple, class... MAOps>
+	static auto callOperation(Operation& op, const Tuple& tp, const MAOps&... args)
+	    -> decltype(op(std::get<0>(tp).get(), args.get()...))
 	{
-	    return op(first, std::get<0>(tp).get(), secs.get()...);
+	    return op(std::get<0>(tp).get(), args.get()...);
 	}
 
-	template <typename T, class Operation, class Tuple, class... MBases>
-	static auto callOperation(const Operation& op, const Tuple& tp, const T& first, const MBases&... secs)
-	    -> decltype(op(first, std::get<0>(tp).get(), secs.get()...))
+	template <typename T, class Operation, class Tuple, class... MAOps>
+	static auto callOperation(const Operation& op, const Tuple& tp, const MAOps&... args)
+	    -> decltype(op(std::get<0>(tp).get(), args.get()...))
 	{
-	    return op(first, std::get<0>(tp).get(), secs.get()...);
+	    return op(std::get<0>(tp).get(), args.get()...);
 	}
     };
 
-    template <typename T, class Range, class Operation, class... Ranges>
-    MultiArrayOperation<T,Range,Operation,Ranges...>::
-    MultiArrayOperation(MultiArray<T,Range>& ma, const Name& nm,
-			Operation& op,
-			const MultiArrayOperationBase<T,Ranges>&... secs) :
-	MultiArrayOperationBase<T,Range>(ma, nm),
+    template <typename T, class Operation, class... MAOps>
+    MultiArrayOperation<T,Operation,MAOps...>::
+    MultiArrayOperation(Operation& op, const MAOps&... args) :
 	mOp(op),
-	mSecs(std::make_tuple(secs...)) {}
+	mArgs(std::make_tuple(args...)) {}
 
-    template <typename T, class Range, class Operation, class... Ranges>
-    MultiArrayOperation<T,Range,Operation,Ranges...>::
-    MultiArrayOperation(MultiArray<T,Range>& ma, const Name& nm,
-			const Operation& op,
-			const MultiArrayOperationBase<T,Ranges>&... secs) :
-	MultiArrayOperationBase<T,Range>(ma, nm),
+    template <typename T, class Operation, class... MAOps>
+    MultiArrayOperation<T,Operation,MAOps...>::
+    MultiArrayOperation(const Operation& op, const MAOps&... args) :
 	mOp(op),
-	mSecs(std::make_tuple(secs...)) {}
-    
-    template <typename T, class Range, class Operation, class... Ranges>
-    size_t MultiArrayOperation<T,Range,Operation,Ranges...>::argNum() const
+	mArgs(std::make_tuple(args...)) {}
+
+    template <typename T, class Operation, class... MAOps>
+    template <class Operation2, class... MAOps2>
+    MultiArrayOperation<T,Operation2,MultiArrayOperation<T,Operation,MAOps...>,MAOps2...>
+    MultiArrayOperation<T,Operation,MAOps...>::operator()(Operation2& op, const MAOps&... secs)
     {
-	return sizeof...(Ranges) + 1;
+	return MultiArrayOperation<T,Operation2,MultiArrayOperation<T,Operation,MAOps...>,
+				   MAOps2...>(op, *this, secs...);
     }
 
-    template <typename T, class Range, class Operation, class... Ranges>    
-    void MultiArrayOperation<T,Range,Operation,Ranges...>::linkIndicesTo(IndefinitIndexBase* target) const
+    template <typename T, class Operation, class... MAOps>
+    template <class Operation2, class... MAOps2>
+    MultiArrayOperation<T,Operation2,MultiArrayOperation<T,Operation,MAOps...>,MAOps2...>
+    MultiArrayOperation<T,Operation,MAOps...>::operator()(const Operation2& op, const MAOps&... secs)
     {
-	OB::mIibPtr->linkTo(target);
-	TupleIndicesLinker<sizeof...(Ranges)-1>::linkTupleIndicesTo(mSecs, target);
+	return MultiArrayOperation<T,Operation2,MultiArrayOperation<T,Operation,MAOps...>,
+				   MAOps2...>(op, *this, secs...);
+    }
+
+    template <typename T, class Operation, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayOperation<T,Operation,MAOps...>::operator+(const MAOp2& sec)
+	-> decltype(operator()(std::plus<T>(), sec))
+    {
+	return operator()(std::plus<T>(), sec);
+    }
+
+    template <typename T, class Operation, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayOperation<T,Operation,MAOps...>::operator-(const MAOp2& sec)
+	-> decltype(operator()(std::minus<T>(), sec))
+    {
+	return operator()(std::minus<T>(), sec);
+    }
+
+    template <typename T, class Operation, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayOperation<T,Operation,MAOps...>::operator*(const MAOp2& sec)
+	-> decltype(operator()(std::multiplies<T>(), sec))
+    {
+	return operator()(std::multiplies<T>(), sec);
+    }
+
+    template <typename T, class Operation, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayOperation<T,Operation,MAOps...>::operator/(const MAOp2& sec)
+	-> decltype(operator()(std::divides<T>(), sec))
+    {
+	return operator()(std::divides<T>(), sec);
     }
     
-    template <typename T, class Range, class Operation, class... Ranges>    
-    T& MultiArrayOperation<T,Range,Operation,Ranges...>::get()
+    template <typename T, class Operation, class... MAOps>
+    size_t MultiArrayOperation<T,Operation,MAOps...>::argNum() const
     {
-	mVal = OperationCall<sizeof...(Ranges)-1>::
-	    template callOperation(mOp, mSecs,
-				   //OB::mArrayRef[OB::mIndex]);
-				   OB::mArrayRef[*dynamic_cast<typename OB::IndexType*>(OB::mIibPtr)]);
+	return sizeof...(MAOps) + 1;
+    }
+
+    template <typename T, class Operation, class... MAOps>    
+    void MultiArrayOperation<T,Operation,MAOps...>::linkIndicesTo(IndefinitIndexBase* target) const
+    {
+	MAOB::mIibPtr->linkTo(target);
+	TupleIndicesLinker<sizeof...(MAOps)-1>::linkTupleIndicesTo(mArgs, target);
+    }
+    
+    template <typename T, class Operation, class... MAOps>    
+    T& MultiArrayOperation<T,Operation,MAOps...>::get()
+    {
+	mVal = OperationCall<sizeof...(MAOps)-1>::
+	    template callOperation(mOp, mArgs);
 	return mVal;
     }
 
-    template <typename T, class Range, class Operation, class... Ranges>    
-    const T& MultiArrayOperation<T,Range,Operation,Ranges...>::get() const
+    template <typename T, class Operation, class... MAOps>    
+    const T& MultiArrayOperation<T,Operation,MAOps...>::get() const
     {
-	mVal = OperationCall<sizeof...(Ranges)-1>::
-	    template callOperation(mOp, mSecs,
-				   //OB::mArrayRef[OB::mIndex]);
-				   OB::mArrayRef[*dynamic_cast<typename OB::IndexType*>(OB::mIibPtr)]);
+	mVal = OperationCall<sizeof...(MAOps)-1>::
+	    template callOperation(mOp, mArgs);
 	return mVal;
     }
-
 }
diff --git a/src/multi_array_operation.h b/src/multi_array_operation.h
index fca3f78..90253fe 100644
--- a/src/multi_array_operation.h
+++ b/src/multi_array_operation.h
@@ -12,89 +12,123 @@
 namespace MultiArrayTools
 {
 
-    template <typename T, class Range>
+    template <typename T>
     class MultiArrayOperationBase
     {
     public:
-
-	typedef decltype(MultiArray<T,Range>().begin()) IndexType;
-	
-	MultiArrayOperationBase(MultiArray<T,Range>& ma, const Name& nm);
-	MultiArrayOperationBase& operator=(const MultiArrayOperationBase& in);
-	//MultiArrayOperationBase(const MultiArrayOperationBase& in) = default;
-	
 	virtual ~MultiArrayOperationBase();
-	
-	// execute AnyOperation
-	// exception if range types are inconsitent with names
-	template <class Range2>
-	MultiArrayOperationBase& operator=(const MultiArrayOperationBase<T, Range2>& in);
-	
-	
-	template <class Operation, class... Ranges>
-	MultiArrayOperation<T,Range,Operation,Ranges...>
-	operator()(Operation& op, const MultiArrayOperationBase<T,Ranges>&... secs);
 
-	template <class Operation, class... Ranges>
-	MultiArrayOperation<T,Range,Operation,Ranges...>
-	operator()(const Operation& op, const MultiArrayOperationBase<T,Ranges>&... secs);
-	
-	template <class Range2>
-	MultiArrayOperation<T,Range,std::plus<T>,Range2> operator+(const MultiArrayOperationBase<T,Range2>& sec);
-
-	template <class Range2>
-	MultiArrayOperation<T,Range,std::minus<T>,Range2> operator-(const MultiArrayOperationBase<T,Range2>& sec);
-
-	template <class Range2>
-	MultiArrayOperation<T,Range,std::multiplies<T>,Range2> operator*(const MultiArrayOperationBase<T,Range2>& sec);
-
-	template <class Range2>
-	MultiArrayOperation<T,Range,std::divides<T>,Range2> operator/(const MultiArrayOperationBase<T,Range2>& sec);
-
-	virtual size_t argNum() const;
-
-	//IndexType& index() ;
+	virtual size_t argNum() const = 0;
 	IndefinitIndexBase* index();
-
-	virtual void linkIndicesTo(IndefinitIndexBase* target) const;
-
-	virtual T& get();
-	virtual const T& get() const;
+	virtual void linkIndicesTo(IndefinitIndexBase* target) const = 0;
+	
+	virtual T& get() = 0;
+	virtual const T& get() const = 0;
 
     protected:
-
-	MultiArray<T,Range>& mArrayRef;
-	mutable IndexType mIndex;
 	IndefinitIndexBase* mIibPtr = nullptr;
-	Name mNm;
     };
-
-    template <typename T, class Range, class Operation, class... Ranges>
-    class MultiArrayOperation : public MultiArrayOperationBase<T,Range>
+    
+    template <typename T, class Range>
+    class MultiArrayOperationRoot : public MultiArrayOperationBase<T>
     {
     public:
 
-	typedef MultiArrayOperationBase<T,Range> OB;
-	typedef std::tuple<MultiArrayOperationBase<T,Ranges>... > OBT;
+	typedef MultiArrayOperationBase<T> MAOB;
+	typedef decltype(MultiArray<T,Range>().begin()) IndexType;
 	
-	MultiArrayOperation(MultiArray<T,Range>& ma, const Name& nm,
-			    Operation& op, const MultiArrayOperationBase<T,Ranges>&... secs);
+	MultiArrayOperationRoot(MultiArray<T,Range>& ma, const Name& nm);
+	MultiArrayOperationRoot& operator=(const MultiArrayOperationBase<T>& in);
+	
+	MultiArrayOperationRoot& operator=(const MultiArrayOperationRoot& in) = delete;
+	//MultiArrayOperationRoot(const MultiArrayOperationRoot& in) = default;
+	
+	// execute AnyOperation
+	// exception if range types are inconsitent with names
+	//MultiArrayOperationRoot& operator=(const MultiArrayOperationBase<T>& in);
+	
+	
+	template <class Operation, class... MAOps>
+	MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
+	operator()(Operation& op, const MAOps&... secs);
 
-	MultiArrayOperation(MultiArray<T,Range>& ma, const Name& nm,
-			    const Operation& op, const MultiArrayOperationBase<T,Ranges>&... secs);
+	template <class Operation, class... MAOps>
+	MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
+	operator()(const Operation& op, const MAOps&... secs);
+	
+	template <class MAOp>
+	auto operator+(const MAOp& sec) -> decltype(operator()(std::plus<T>(), sec));
+
+	template <class MAOp>
+	auto operator-(const MAOp& sec) -> decltype(operator()(std::minus<T>(), sec));
+
+	template <class MAOp>
+	auto operator*(const MAOp& sec) -> decltype(operator()(std::multiplies<T>(), sec));
+
+	template <class MAOp>
+	auto operator/(const MAOp& sec) -> decltype(operator()(std::divides<T>(), sec));
 
 	virtual size_t argNum() const override;
-	
+
+	//IndexType& index() ;
+
 	virtual void linkIndicesTo(IndefinitIndexBase* target) const override;
 
 	virtual T& get() override;
 	virtual const T& get() const override;
 	
+    protected:
+	
+	MultiArray<T,Range>& mArrayRef;
+	mutable IndexType mIndex;
+	Name mNm;
+    };
+
+    template <typename T, class Operation, class... MAOps>
+    class MultiArrayOperation : public MultiArrayOperationBase<T>
+    {
+    public:
+
+	typedef MultiArrayOperationBase<T> MAOB;
+	typedef std::tuple<MAOps...> OBT;
+	
+	MultiArrayOperation(Operation& op, const MAOps&... secs);
+	MultiArrayOperation(const Operation& op, const MAOps&... secs);
+
+	template <class Operation2, class... MAOps2>
+	MultiArrayOperation<T,Operation2,MultiArrayOperation<T,Operation,MAOps...>,MAOps2...>
+	operator()(Operation2& op, const MAOps&... secs);
+	
+	template <class Operation2, class... MAOps2>
+	MultiArrayOperation<T,Operation2,MultiArrayOperation<T,Operation,MAOps...>,MAOps2...>
+	operator()(const Operation2& op, const MAOps&... secs);
+	
+	template <class MAOp2>
+	auto operator+(const MAOp2& sec) -> decltype(operator()(std::plus<T>(), sec));
+
+	template <class MAOp2>
+	auto operator-(const MAOp2& sec) -> decltype(operator()(std::minus<T>(), sec));
+
+	template <class MAOp2>
+	auto operator*(const MAOp2& sec) -> decltype(operator()(std::multiplies<T>(), sec));
+
+	template <class MAOp2>
+	auto operator/(const MAOp2& sec) -> decltype(operator()(std::divides<T>(), sec));
+	
+	virtual size_t argNum() const override;
+	
+	virtual void linkIndicesTo(IndefinitIndexBase* target) const override;
+
+	virtual T& get() override;
+	virtual const T& get() const override;
+
+	virtual void executeOp() override;
+	
     protected:
 
  	mutable T mVal;
 	Operation mOp;
-	OBT mSecs;
+	OBT mArgs; // include first arg also here !!!
     };
 
 
diff --git a/src/unit_test.cc b/src/unit_test.cc
index 2acbc9f..7dcd559 100644
--- a/src/unit_test.cc
+++ b/src/unit_test.cc
@@ -82,14 +82,18 @@ namespace {
     protected:
 	
 	typedef MAT::SingleRange<char,MAT::RangeType::ANY> Range1dAny;
+	typedef MAT::MultiRange<Range1dAny> RangeAny;
 	typedef MAT::MultiRange<Range1dAny,Range1dAny> Range2dAny;
 	typedef MAT::MultiRange<Range1dAny,Range1dAny,Range1dAny> Range3dAny;
+	typedef MAT::MultiArray<int,RangeAny> MultiArray1dAny;
 	typedef MAT::MultiArray<int,Range2dAny> MultiArray2dAny;
 	typedef MAT::MultiArray<int,Range3dAny> MultiArray3dAny;
 
 	OperationTest() : r1({'a','b','c'}), r2({'a','b','c','d'}), r3({'a','b'}),
+			  r11(r1),
 			  ra(r1,r2),
 			  r3d(r1,r2,r3),
+			  ma1d(r11, {-5,6,2}),
 			  ma(ra, {-5,6,2,1,9,54,27,-7,-13,32,90,-67}),
 			  ma3d(r3d, {-5,6,2,1,9,54,27,-7,-13,32,90,-67,
 				      -10,16,-2,101,39,-64,81,-22,14,34,95,-62}) {}
@@ -97,8 +101,10 @@ namespace {
 	Range1dAny r1;
 	Range1dAny r2;
 	Range1dAny r3;
+	RangeAny r11;
 	Range2dAny ra;
 	Range3dAny r3d;
+	MultiArray1dAny ma1d;
 	MultiArray2dAny ma;
 	MultiArray3dAny ma3d;
     };
@@ -280,6 +286,45 @@ namespace {
 	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 3, i3 = 0)], 28);
 	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 3, i3 = 1)], -129);
     }
+
+    TEST_F(OperationTest, CorrectlyAdded_3args)
+    {
+	MultiArray3dAny ma3d2(r3d);
+	//MultiArray2dAny ma2(ra);
+	auto i = ma3d2.begin();
+	auto i1 = i.template getIndex<0>();
+	auto i2 = i.template getIndex<1>();
+	auto i3 = i.template getIndex<2>();
+
+	//ma2("alpha","beta") = ma("alpha","beta") + ma1d("alpha");
+	ma3d2("alpha","beta","gamma") = ma3d("alpha","beta","gamma") + ( ma("alpha","beta") + ma1d("alpha") );
+
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 0, i3 = 0)], -15);
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 0, i3 = 1)], -4);
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 1, i3 = 0)], 3);
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 1, i3 = 1)], 2);
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 2, i3 = 0)], 6);
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 2, i3 = 1)], 51);
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 3, i3 = 0)], 23);
+	EXPECT_EQ(ma3d2[i(i1 = 0, i2 = 3, i3 = 1)], -11);
+	
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 0, i3 = 0)], 2);
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 0, i3 = 1)], 47);
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 1, i3 = 0)], 150);
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 1, i3 = 1)], -7);
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 2, i3 = 0)], 23);
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 2, i3 = 1)], 49);
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 3, i3 = 0)], -3);
+	EXPECT_EQ(ma3d2[i(i1 = 1, i2 = 3, i3 = 1)], 100);
+
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 0, i3 = 0)], 28);
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 0, i3 = 1)], -75);
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 1, i3 = 0)], 115);
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 1, i3 = 1)], 12);
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 2, i3 = 0)], 106);
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 2, i3 = 1)], 126);
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 3, i3 = 0)], 30);
+	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 3, i3 = 1)], -127);    }
     
 } // end namespace