diff --git a/src/index_base.cc b/src/index_base.cc
index 3a425a6..ad85584 100644
--- a/src/index_base.cc
+++ b/src/index_base.cc
@@ -1,4 +1,4 @@
-
+// -*- C++ -*-
 #include "index_base.h"
 
 namespace MultiArrayTools
@@ -9,7 +9,7 @@ namespace MultiArrayTools
     
     IndefinitIndexBase::~IndefinitIndexBase()
     {
-	freeLinked();
+	//freeLinked();
 	mLinked = nullptr;
 	mMajor = nullptr;
 	mSoftLinked = nullptr;
diff --git a/src/index_base.h b/src/index_base.h
index 9ed71c6..f230d83 100644
--- a/src/index_base.h
+++ b/src/index_base.h
@@ -44,7 +44,7 @@ namespace MultiArrayTools
 	virtual bool linked() const;
 	virtual void linkTo(IndefinitIndexBase* target) = 0; 
 
-	virtual IndefinitIndexBase* getLinked(const std::string& name) = 0;
+	//virtual IndefinitIndexBase* getLinked(const std::string& name) = 0;
 	
 	virtual void setPos(size_t pos);
 
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
index f56879f..51031ac 100644
--- a/src/multi_array_operation.cc
+++ b/src/multi_array_operation.cc
@@ -154,6 +154,44 @@ namespace MultiArrayTools
     {
 	return operator()(std::divides<T>(), sec);
     }
+
+    template <typename T, class Range>
+    MultiArrayOperationRoot<T,Range> MultiArrayOperationRoot<T,Range>::copyThis()
+    {
+	return MultiArrayOperationRoot<T,Range>(mArrayRef, mNm);
+    }
+    
+    template <typename T, class Range>
+    template <class MAOp>
+    MultiArrayOperationRoot<T,Range>&
+    MultiArrayOperationRoot<T,Range>::operator+=(const MAOp& sec)
+    {
+	return (*this) = copyThis() + sec;
+    }
+
+    template <typename T, class Range>
+    template <class MAOp>
+    MultiArrayOperationRoot<T,Range>&
+    MultiArrayOperationRoot<T,Range>::operator-=(const MAOp& sec)
+    {
+	return (*this) = copyThis() - sec;
+    }
+
+    template <typename T, class Range>
+    template <class MAOp>
+    MultiArrayOperationRoot<T,Range>&
+    MultiArrayOperationRoot<T,Range>::operator*=(const MAOp& sec)
+    {
+	return (*this) = copyThis() * sec;
+    }
+
+    template <typename T, class Range>
+    template <class MAOp>
+    MultiArrayOperationRoot<T,Range>&
+    MultiArrayOperationRoot<T,Range>::operator/=(const MAOp& sec)
+    {
+	return (*this) = copyThis() / sec;
+    }
     
     template <typename T, class Range>
     size_t MultiArrayOperationRoot<T,Range>::argNum() const
diff --git a/src/multi_array_operation.h b/src/multi_array_operation.h
index a697c25..1abebb1 100644
--- a/src/multi_array_operation.h
+++ b/src/multi_array_operation.h
@@ -69,7 +69,21 @@ namespace MultiArrayTools
 
 	template <class MAOp>
 	auto operator/(const MAOp& sec) -> decltype(operator()(std::divides<T>(), sec));
+	
+	MultiArrayOperationRoot copyThis();
+	
+	template <class MAOp>
+	MultiArrayOperationRoot& operator+=(const MAOp& sec);
 
+	template <class MAOp>
+	MultiArrayOperationRoot& operator-=(const MAOp& sec);
+
+	template <class MAOp>
+	MultiArrayOperationRoot& operator*=(const MAOp& sec);
+
+	template <class MAOp>
+	MultiArrayOperationRoot& operator/=(const MAOp& sec);
+	
 	virtual size_t argNum() const override;
 
 	// set index -> implement !!!!!
diff --git a/src/unit_test.cc b/src/unit_test.cc
index 0440bea..95ac581 100644
--- a/src/unit_test.cc
+++ b/src/unit_test.cc
@@ -313,6 +313,44 @@ namespace {
 	EXPECT_EQ(ma3d2[i(i1 = 2, i2 = 3, i3 = 1)], -129);
     }
 
+    TEST_F(OperationTest, CorrectlyAdded2)
+    {
+	//MultiArray3dAny ma3d2(r3d);
+	auto i = ma3d.begin();
+	auto i1 = i.template getIndex<0>();
+	auto i2 = i.template getIndex<1>();
+	auto i3 = i.template getIndex<2>();
+
+	ma3d("alpha","beta","gamma") += ma("alpha","beta");
+
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 0, i3 = 0)], -10);
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 0, i3 = 1)], 1);
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 1, i3 = 0)], 8);
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 1, i3 = 1)], 7);
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 2, i3 = 0)], 11);
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 2, i3 = 1)], 56);
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 3, i3 = 0)], 28);
+	EXPECT_EQ(ma3d[i(i1 = 0, i2 = 3, i3 = 1)], -6);
+	
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 0, i3 = 0)], -4);
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 0, i3 = 1)], 41);
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 1, i3 = 0)], 144);
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 1, i3 = 1)], -13);
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 2, i3 = 0)], 17);
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 2, i3 = 1)], 43);
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 3, i3 = 0)], -9);
+	EXPECT_EQ(ma3d[i(i1 = 1, i2 = 3, i3 = 1)], 94);
+
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 0, i3 = 0)], 26);
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 0, i3 = 1)], -77);
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 1, i3 = 0)], 113);
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 1, i3 = 1)], 10);
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 2, i3 = 0)], 104);
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 2, i3 = 1)], 124);
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 3, i3 = 0)], 28);
+	EXPECT_EQ(ma3d[i(i1 = 2, i2 = 3, i3 = 1)], -129);
+    }
+
     TEST_F(OperationTest, CorrectlyAdded_3args)
     {
 	MultiArray3dAny ma3d2(r3d);