From c043008217cf4d0814934c54eba9275a379f0176 Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Thu, 19 Jul 2018 19:29:52 +0200
Subject: [PATCH] function structs corresponding to math.h functions + some
 code optimizations

---
 src/include/arith.h                  |  43 +++++++++--
 src/include/extensions/math.h        |  31 ++++++++
 src/include/multi_array.h            |  54 +++++++++++++-
 src/include/multi_array_header.h     |   1 +
 src/include/multi_array_operation.h  | 104 ++++-----------------------
 src/include/operation_def.h          |  94 ++++++++++++++++++++++++
 src/include/pack_num.h               |  14 ++++
 src/include/ranges/anonymous_range.h |   1 -
 8 files changed, 246 insertions(+), 96 deletions(-)
 create mode 100644 src/include/extensions/math.h
 create mode 100644 src/include/operation_def.h

diff --git a/src/include/arith.h b/src/include/arith.h
index a2dd583..336f71d 100644
--- a/src/include/arith.h
+++ b/src/include/arith.h
@@ -60,7 +60,7 @@ namespace MultiArrayTools
     };
     
     template <typename T>
-    struct plus
+    struct plus : public StaticFunctionBase<T, plus<T>, T, T>
     {
 	static constexpr bool FISSTATIC = true;
 	
@@ -71,7 +71,7 @@ namespace MultiArrayTools
     };
 
     template <typename T>
-    struct minus
+    struct minus : public StaticFunctionBase<T, minus<T>, T, T>
     {
 	static constexpr bool FISSTATIC = true;
 	
@@ -82,7 +82,7 @@ namespace MultiArrayTools
     };
 
     template <typename T>
-    struct multiplies
+    struct multiplies : public StaticFunctionBase<T, multiplies<T>, T, T>
     {
 	static constexpr bool FISSTATIC = true;
 	
@@ -93,7 +93,7 @@ namespace MultiArrayTools
     };
 
     template <typename T>
-    struct divides
+    struct divides : public StaticFunctionBase<T, divides<T>, T, T>
     {
 	static constexpr bool FISSTATIC = true;
 	
@@ -102,6 +102,41 @@ namespace MultiArrayTools
 	    return a1 / a2;
 	}
     };
+
+#include <cmath>
+#define regFunc1(fff) template <typename T>\
+    struct x_##fff : public StaticFunctionBase<T, x_##fff<T>, T> {\
+    static constexpr bool FISSTATIC = true;\
+    static inline T apply(T a){\
+	return fff(a); } };
+
+#include "extensions/math.h"
+#undef regFunc1
+    
+    template <size_t N>
+    struct x_ipow
+    {
+	static constexpr bool FISSTATIC = true;
+
+	template <typename T>
+	static inline T apply(T a)
+	{
+	    return a * x_ipow<N-1>::apply(a);
+	}
+    };
+
+    template <>
+    struct x_ipow<0>
+    {
+	static constexpr bool FISSTATIC = true;
+
+	template <typename T>
+	static inline T apply(T a)
+	{
+	    return a;
+	}
+    };
+	
     /*
     template <typename T, class Func>
     struct dynamic_function
diff --git a/src/include/extensions/math.h b/src/include/extensions/math.h
new file mode 100644
index 0000000..394fb0e
--- /dev/null
+++ b/src/include/extensions/math.h
@@ -0,0 +1,31 @@
+
+#ifdef regFunc1
+
+regFunc1(exp)
+regFunc1(exp2)
+regFunc1(expm1)
+regFunc1(log)
+regFunc1(log10)
+regFunc1(log2)
+regFunc1(log1p)
+regFunc1(sqrt)
+regFunc1(cbrt)
+regFunc1(sin)
+regFunc1(cos)
+regFunc1(tan)
+regFunc1(asin)
+regFunc1(acos)
+regFunc1(atan)
+regFunc1(atan2)
+regFunc1(sinh)
+regFunc1(cosh)
+regFunc1(tanh)
+regFunc1(asinh)
+regFunc1(acosh)
+regFunc1(atanh)
+regFunc1(erf)
+regFunc1(erfc)
+regFunc1(tgamma)
+regFunc1(lgamma)
+
+#endif
diff --git a/src/include/multi_array.h b/src/include/multi_array.h
index a6aefcf..9273ba9 100644
--- a/src/include/multi_array.h
+++ b/src/include/multi_array.h
@@ -3,6 +3,8 @@
 #ifndef __multi_array_h__
 #define __multi_array_h__
 
+#include <algorithm>
+
 #include "multi_array_base.h"
 #include "ranges/anonymous_range.h"
 
@@ -77,6 +79,11 @@ namespace MultiArrayTools
 	    -> decltype(ArrayCatter<T>::cat(*this));
 	
 	operator T() const;	
+
+	MultiArray& operator+=(const MultiArray& in);
+	MultiArray& operator-=(const MultiArray& in);
+	MultiArray& operator*=(const T& in);
+	MultiArray& operator/=(const T& in);
 	
 	template <typename U, class... SRanges2>
 	friend class MultiArray;
@@ -187,7 +194,8 @@ namespace MultiArrayTools
 
     template <typename T, class... SRanges>
     MultiArray<T,SRanges...>::MultiArray(MultiArray<T,AnonymousRange>& ama, SIZET<SRanges>... sizes) :
-	MutableMultiArrayBase<T,SRanges...>( ama.range()->template scast<SRanges...>(sizes...)->space() ),
+	MutableMultiArrayBase<T,SRanges...>
+	( ama.range()->template get<0>().template scast<SRanges...>(sizes...)->space() ),
 	mCont( std::move( ama.mCont ) )
     {
 	MAB::mInit = true;
@@ -297,6 +305,50 @@ namespace MultiArrayTools
 	      std::move(mCont) );
     }	
 
+    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
     {
diff --git a/src/include/multi_array_header.h b/src/include/multi_array_header.h
index ab6021d..0879676 100644
--- a/src/include/multi_array_header.h
+++ b/src/include/multi_array_header.h
@@ -16,6 +16,7 @@
 #include "multi_array.h"
 #include "functional_multi_array.h"
 #include "helper_tools.h"
+#include "operation_def.h"
 //#include "slice.h"
 //#include "manipulator.h"
 //#include "range_transformer.h"
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index b3754ae..d7fdfd6 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -33,32 +33,27 @@ namespace MultiArrayTools
 	OperationClass& THIS() { return static_cast<OperationClass&>(*this); }
 	const OperationClass& THIS() const { return static_cast<OperationClass const&>(*this); }
 
-	inline auto operator+(const T& in) const
-	    -> Operation<T,plus<T>,OperationClass,OperationValue<T> >;
-
-	// !!!
-	
 	template <class Second>
-	auto operator+(const Second& in) const
+	auto operator+(const OperationBase<T,Second>& in) const
 	    -> Operation<T,plus<T>,OperationClass,Second>;
 	
 	template <class Second>
-	auto operator-(const Second& in) const
+	auto operator-(const OperationBase<T,Second>& in) const
 	    -> Operation<T,minus<T>,OperationClass,Second>;
 	
 	template <class Second>
-	auto operator*(const Second& in) const
+	auto operator*(const OperationBase<T,Second>& in) const
 	    -> Operation<T,multiplies<T>,OperationClass,Second>;
 	
 	template <class Second>
-	auto operator/(const Second& in) const
+	auto operator/(const OperationBase<T,Second>& in) const
 	    -> Operation<T,divides<T>,OperationClass,Second>;
 
 	template <class IndexType>
 	auto c(std::shared_ptr<IndexType>& ind) const
 	    -> Contraction<T,OperationClass,IndexType>;
-		
-    private:
+
+    private:		
 	friend OperationClass;
 	friend OperationTemplate<T,OperationClass>;
 	OperationBase() = default;
@@ -306,69 +301,6 @@ namespace MultiArrayTools
 	return Operation<typename OpFunction::value_type,OpFunction,Ops...>(ops...);
     }
 
-    template <typename T, class OperationClass>
-    auto operator+(const T& a, const OperationBase<T,OperationClass>& b)
-	-> Operation<T,plus<T>,OperationValue<T>,OperationClass>
-    {
-	OperationValue<T> v(a);
-	return Operation<T,plus<T>,OperationValue<T>,OperationClass>(v, b.THIS());
-    }
-
-    template <typename T, class OperationClass>
-    auto operator-(const T& a, const OperationBase<T,OperationClass>& b)
-	-> Operation<T,minus<T>,OperationValue<T>,OperationClass>
-    {
-	OperationValue<T> v(a);
-	return Operation<T,minus<T>,OperationValue<T>,OperationClass>(v, b.THIS());
-    }
-
-    template <typename T, class OperationClass>
-    auto operator*(const T& a, const OperationBase<T,OperationClass>& b)
-	-> Operation<T,multiplies<T>,OperationValue<T>,OperationClass>
-    {
-	OperationValue<T> v(a);
-	return Operation<T,multiplies<T>,OperationValue<T>,OperationClass>(v, b.THIS());
-    }
-
-    template <typename T, class OperationClass>
-    auto operator/(const T& a, const OperationBase<T,OperationClass>& b)
-	-> Operation<T,divides<T>,OperationValue<T>,OperationClass>
-    {
-	OperationValue<T> v(a);
-	return Operation<T,divides<T>,OperationValue<T>,OperationClass>(v, b.THIS());
-    }
-
-    template <typename T, class OperationClass>
-    auto operator+(const OperationBase<T,OperationClass>& a, const T& b)
-	-> Operation<T,plus<T>,OperationClass,OperationValue<T> >
-    {
-	OperationValue<T> v(b);
-	return Operation<T,plus<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
-    }
-
-    template <typename T, class OperationClass>
-    auto operator-(const OperationBase<T,OperationClass>& a, const T& b)
-	-> Operation<T,minus<T>,OperationClass,OperationValue<T> >
-    {
-	OperationValue<T> v(b);
-	return Operation<T,minus<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
-    }
-
-    template <typename T, class OperationClass>
-    auto operator*(const OperationBase<T,OperationClass>& a, const T& b)
-	-> Operation<T,multiplies<T>,OperationClass,OperationValue<T> >
-    {
-	OperationValue<T> v(b);
-	return Operation<T,multiplies<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
-    }
-
-    template <typename T, class OperationClass>
-    auto operator/(const OperationBase<T,OperationClass>& a, const T& b)
-	-> Operation<T,divides<T>,OperationClass,OperationValue<T> >
-    {
-	OperationValue<T> v(b);
-	return Operation<T,divides<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
-    }
 
     
     template <typename T, class Op, class IndexType>
@@ -418,45 +350,37 @@ namespace MultiArrayTools
      *   OperationTemplate     *
      ***************************/
 
-
-    template <typename T, class OperationClass>
-    auto OperationBase<T,OperationClass>::operator+(const T& in) const
-	-> Operation<T,plus<T>,OperationClass,OperationValue<T> >
-    {
-	return Operation<T,plus<T>,OperationClass,OperationValue<T> >(THIS(), in);
-    }
-
     
     template <typename T, class OperationClass>
     template <class Second>
-    auto OperationBase<T,OperationClass>::operator+(const Second& in) const
+    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);
+    	return Operation<T,plus<T>,OperationClass,Second>(THIS(), in.THIS());
     }
 
     template <typename T, class OperationClass>
     template <class Second>
-    auto OperationBase<T,OperationClass>::operator-(const Second& in) const
+    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);
+    	return Operation<T,minus<T>,OperationClass,Second>(THIS(), in.THIS());
     }
     
     template <typename T, class OperationClass>
     template <class Second>
-    auto OperationBase<T,OperationClass>::operator*(const Second& in) const
+    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);
+    	return Operation<T,multiplies<T>,OperationClass,Second>(THIS(), in.THIS());
     }
 
     template <typename T, class OperationClass>
     template <class Second>
-    auto OperationBase<T,OperationClass>::operator/(const Second& in) const
+    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);
+    	return Operation<T,divides<T>,OperationClass,Second>(THIS(), in.THIS());
     }
 
     template <typename T, class OperationClass>
diff --git a/src/include/operation_def.h b/src/include/operation_def.h
new file mode 100644
index 0000000..ec4fafc
--- /dev/null
+++ b/src/include/operation_def.h
@@ -0,0 +1,94 @@
+
+#ifndef __operation_def_h__
+#define __operation_def_h__
+
+#include "multi_array_operation.h"
+
+namespace MultiArrayTools
+{
+
+    template <typename T, class OperationClass>
+    auto operator+(const T& a, const OperationBase<T,OperationClass>& b)
+	-> Operation<T,plus<T>,OperationValue<T>,OperationClass>
+    {
+	OperationValue<T> v(a);
+	return Operation<T,plus<T>,OperationValue<T>,OperationClass>(v, b.THIS());
+    }
+
+    template <typename T, class OperationClass>
+    auto operator-(const T& a, const OperationBase<T,OperationClass>& b)
+	-> Operation<T,minus<T>,OperationValue<T>,OperationClass>
+    {
+	OperationValue<T> v(a);
+	return Operation<T,minus<T>,OperationValue<T>,OperationClass>(v, b.THIS());
+    }
+
+    template <typename T, class OperationClass>
+    auto operator*(const T& a, const OperationBase<T,OperationClass>& b)
+	-> Operation<T,multiplies<T>,OperationValue<T>,OperationClass>
+    {
+	OperationValue<T> v(a);
+	return Operation<T,multiplies<T>,OperationValue<T>,OperationClass>(v, b.THIS());
+    }
+
+    template <typename T, class OperationClass>
+    auto operator/(const T& a, const OperationBase<T,OperationClass>& b)
+	-> Operation<T,divides<T>,OperationValue<T>,OperationClass>
+    {
+	OperationValue<T> v(a);
+	return Operation<T,divides<T>,OperationValue<T>,OperationClass>(v, b.THIS());
+    }
+
+    template <typename T, class OperationClass>
+    auto operator+(const OperationBase<T,OperationClass>& a, const T& b)
+	-> Operation<T,plus<T>,OperationClass,OperationValue<T> >
+    {
+	OperationValue<T> v(b);
+	return Operation<T,plus<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
+    }
+
+    template <typename T, class OperationClass>
+    auto operator-(const OperationBase<T,OperationClass>& a, const T& b)
+	-> Operation<T,minus<T>,OperationClass,OperationValue<T> >
+    {
+	OperationValue<T> v(b);
+	return Operation<T,minus<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
+    }
+
+    template <typename T, class OperationClass>
+    auto operator*(const OperationBase<T,OperationClass>& a, const T& b)
+	-> Operation<T,multiplies<T>,OperationClass,OperationValue<T> >
+    {
+	OperationValue<T> v(b);
+	return Operation<T,multiplies<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
+    }
+
+    template <typename T, class OperationClass>
+    auto operator/(const OperationBase<T,OperationClass>& a, const T& b)
+	-> Operation<T,divides<T>,OperationClass,OperationValue<T> >
+    {
+	OperationValue<T> v(b);
+	return Operation<T,divides<T>,OperationClass,OperationValue<T> >(a.THIS(), v);
+    }
+    
+#define regFunc1(fff) template <typename T, class OperationClass>	\
+    auto fff(const OperationBase<T,OperationClass>& a)	\
+	-> Operation<T,x_##fff<T>,OperationClass> {			\
+	return Operation<T,x_##fff<T>,OperationClass>(a.THIS()); }
+
+    
+#include "extensions/math.h"
+#undef regFunc1
+
+
+    template <size_t N, typename T, class OperationClass>
+    auto ipow(const OperationBase<T,OperationClass>& a)
+	-> Operation<T,x_ipow<N>,OperationClass>
+    {
+	return Operation<T,x_ipow<N>,OperationClass>(a.THIS());
+    }
+    
+
+}
+
+#endif
diff --git a/src/include/pack_num.h b/src/include/pack_num.h
index 69cf12a..bb38001 100644
--- a/src/include/pack_num.h
+++ b/src/include/pack_num.h
@@ -100,6 +100,13 @@ namespace MultiArrayHelper
 	    std::get<N+1>(blocks) = tmp;
 	    PackNum<N-1>::template mkSliceBlocks<T,Op,SRanges...>(blocks, index, op, total * tmp);
 	}
+
+	template <class... SRanges>
+	static bool checkIfSameInstance(const std::tuple<std::shared_ptr<SRanges>...>& rtp1,
+					const std::tuple<std::shared_ptr<SRanges>...>& rtp2)
+	{
+	    return std::get<N>(rtp1).get() == std::get<N>(rtp2).get() and PackNum<N-1>::checkIfSameInstance(rtp1,rtp2);
+	}
     };
     
     template<>
@@ -164,6 +171,13 @@ namespace MultiArrayHelper
 	    std::get<0>(blocks) = total * tmp; // this is not correct, but not used so far ... !!!
 	}
 
+	template <class... SRanges>
+	static bool checkIfSameInstance(const std::tuple<std::shared_ptr<SRanges>...>& rtp1,
+					const std::tuple<std::shared_ptr<SRanges>...>& rtp2)
+	{
+	    return std::get<0>(rtp1).get() == std::get<0>(rtp2).get();
+	}
+
     };
 
 
diff --git a/src/include/ranges/anonymous_range.h b/src/include/ranges/anonymous_range.h
index 7eaab27..324680e 100644
--- a/src/include/ranges/anonymous_range.h
+++ b/src/include/ranges/anonymous_range.h
@@ -100,7 +100,6 @@ namespace MultiArrayTools
 	
 	std::vector<std::shared_ptr<RangeBase> > mOrig;
     };
-
 }
 
 /* ========================= *