diff --git a/src/Makefile.am b/src/Makefile.am
index 802c23f261255c7172f16e112cf6b93025b64336..dc9b7b8ff75e07ed32d9246d7f596c0937c07b43 100755
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -6,7 +6,7 @@ numerics_SOURCES = main.cpp ../../openfpm_vcluster/src/VCluster.cpp ../../openfp
 numerics_CXXFLAGS = $(AM_CXXFLAGS) $(HDF5_CPPFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
 numerics_CFLAGS = $(CUDA_CFLAGS)
 numerics_LDADD = $(LINKLIBS) -lparmetis -lmetis
-nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp
+nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp
 
 .cu.o :
 	$(NVCC) $(NVCCFLAGS) -o $@ -c $<
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index 8d20ea33ea2839d063b4ed2270d3ccfc9cd467ca..a9b554b1ee07b299e9976ce96337cbf36ceac843 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -10,12 +10,13 @@
 
 #include "Vector/vector_dist.hpp"
 
+#define PROP_POS (unsigned int)-1
+
 #define VECT_SUM 1
 #define VECT_SUB 2
 #define VECT_MUL 3
 #define VECT_DIV 4
-#define VECT_NORM 5
-#define VECT_NORM2 6
+
 #define VECT_APPLYKER_IN 7
 #define VECT_APPLYKER_OUT 8
 #define VECT_APPLYKER_REDUCE 9
@@ -29,6 +30,43 @@
 #define VECT_APPLYKER_MULTI_OUT_GEN 17
 #define VECT_APPLYKER_MULTI_REDUCE_GEN 18
 
+#define VECT_NORM 56
+#define VECT_NORM2 57
+#define VECT_ABS 58
+#define VECT_EXP 59
+#define VECT_EXP2 60
+#define VECT_EXPM1 61
+#define VECT_LOG 62
+#define VECT_LOG10 63
+#define VECT_LOG2 64
+#define VECT_LOG1P 65
+#define VECT_SQRT 67
+#define VECT_CBRT 68
+#define VECT_SIN 69
+#define VECT_COS 70
+#define VECT_TAN 71
+#define VECT_ASIN 72
+#define VECT_ACOS 73
+#define VECT_ATAN 74
+#define VECT_SINH 75
+#define VECT_COSH 76
+#define VECT_TANH 77
+#define VECT_ASINH 78
+#define VECT_ACOSH 79
+#define VECT_ATANH 80
+#define VECT_ERF 81
+#define VECT_ERFC 82
+#define VECT_TGAMMA 83
+#define VECT_LGAMMA 84
+#define VECT_CEIL 85
+#define VECT_FLOOR 86
+#define VECT_TRUNC 87
+#define VECT_ROUND 88
+#define VECT_NEARBYINT 89
+#define VECT_RINT 90
+#define VECT_PMUL 91
+#define VECT_SUB_UNI 92
+
 
 /*! \brief has_init check if a type has defined a
  * method called init
@@ -243,14 +281,30 @@ public:
 	}
 };
 
-/*! \brief norm operation
+/*! \brief selector for position or properties
  *
- * \tparam exp1 expression1
- * \tparam exp2 expression2
  *
  */
+template <typename vector, unsigned int prp>
+struct pos_or_prop
+{
+	static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(v.template getProp<prp>(k))
+	{
+		return v.template getProp<prp>(k);
+	}
+};
+
+template <typename vector>
+struct pos_or_prop<vector,PROP_POS>
+{
+	static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(getExpr(v.template getPos(k)))
+	{
+		return getExpr(v.template getPos(k));
+	}
+};
+
 template <typename exp1>
-class vector_dist_expression_op<exp1,void,VECT_NORM>
+class vector_dist_expression_op<exp1,void,VECT_SUB_UNI>
 {
 	const exp1 o1;
 
@@ -260,28 +314,17 @@ public:
 	:o1(o1)
 	{}
 
-	/*! \brief This function must be called before value
-	 *
-	 * it initialize the expression if needed
-	 *
-	 */
 	inline void init() const
 	{
 		o1.init();
 	}
 
-	/*! \brief Evaluate the expression
-	 *
-	 * \param key where to evaluate the expression
-	 *
-	 */
-	template<typename r_type=typename std::remove_reference<decltype(norm(o1.value(vect_dist_key_dx(0))))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	template<typename r_type=typename std::remove_reference<decltype(-(o1.value(vect_dist_key_dx(0))))>::type > inline r_type value(const vect_dist_key_dx & key) const
 	{
-		return norm(o1.value(key));
+		return -(o1.value(key));
 	}
 };
 
-
 /*! \brief Main class that encapsulate a vector properties
  *
  * \tparam prp property involved
@@ -312,9 +355,32 @@ public:
 	 * \param key where to evaluate the expression
 	 *
 	 */
-	inline auto value(const vect_dist_key_dx & k) const -> decltype(v.template getProp<prp>(k))
+	inline auto value(const vect_dist_key_dx & k) const -> decltype(pos_or_prop<vector,prp>::value(v,k))
 	{
-		return v.template getProp<prp>(k);
+		return pos_or_prop<vector,prp>::value(v,k);
+	}
+
+	/*! \brief Fill the vector property with the evaluated expression
+	 *
+	 * \param v_exp expression to evaluate
+	 *
+	 */
+	template<unsigned int prp2> vector & operator=(const vector_dist_expression<prp2,vector> & v_exp)
+	{
+		v_exp.init();
+
+		auto it = v.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			v.template getProp<prp>(key) = v_exp.value(key);
+
+			++it;
+		}
+
+		return v;
 	}
 
 
@@ -363,7 +429,7 @@ public:
 	}
 };
 
-/*! \Create an expression from a vector
+/*! \Create an expression from a vector property
  *
  * \tpatam prp property
  * \param v
@@ -567,6 +633,39 @@ operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
+/* \brief minus of a distributed vector expression operator
+ *
+ * \param va vector expression one
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2_, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,void,VECT_SUB_UNI>
+operator-(const vector_dist_expression_op<exp1,exp2_,op1> & va)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,void,VECT_SUB_UNI> exp_sum(va);
+
+	return exp_sum;
+}
+
+/* \brief minus of a distributed vector expression
+ *
+ * \param va vector expression one
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<unsigned int p1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,void,VECT_SUB_UNI>
+operator-(const vector_dist_expression<p1,v1> & va)
+{
+	vector_dist_expression_op<vector_dist_expression<p1,v1>,void,VECT_SUB_UNI> exp_sum(va);
+
+	return exp_sum;
+}
+
+
 /* \brief subtract two distributed vector expression
  *
  * \param va vector expression one
@@ -737,6 +836,40 @@ operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 	return exp_sum;
 }
 
+/* \brief Multiply a distributed vector expression by a number
+ *
+ * \param va vector expression
+ * \param d number
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_MUL>
+operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_MUL> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+/* \brief Multiply a distributed vector expression by a number
+ *
+ * \param d number
+ * \param vb vector expression
+ *
+ * \return an object that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL>
+operator*(double d, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL> exp_sum(vector_dist_expression<0,double>(d),vb);
+
+	return exp_sum;
+}
+
 /* \brief Divide two distributed vector expression
  *
  * \param va vector expression one
@@ -873,60 +1006,9 @@ operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 
 	return exp_sum;
 }
-/////////// NORM operator ///////////////////////
-
-/* \brief Divide two distributed vector expression
- *
- * \param va vector expression one
- * \param vb vector expression two
- *
- * \return an object that encapsulate the expression
- *
- */
-template<typename exp1, typename exp2, unsigned int op1>
-inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,void,VECT_NORM>
-norm(const vector_dist_expression_op<exp1,exp2,op1> & va)
-{
-	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,void,VECT_NORM> exp_sum(va);
-
-	return exp_sum;
-}
-
-
-/* \brief Divide two distributed vector expression
- *
- * \param va vector expression one
- * \param vb vector expression two
- *
- * \return an object that encapsulate the expression
- *
- */
-template<typename exp1, typename exp2, unsigned int op1>
-inline vector_dist_expression_op<vector_dist_expression<0,double>,void,VECT_NORM>
-norm(double d)
-{
-	vector_dist_expression_op<vector_dist_expression<0,double>,void,VECT_NORM> exp_sum( (vector_dist_expression<0,double>(d)) );
-
-	return exp_sum;
-}
-
-/* \brief Divide two distributed vector expression
- *
- * \param va vector expression one
- * \param vb vector expression two
- *
- * \return an object that encapsulate the expression
- *
- */
-template<unsigned int prp1, typename v1>
-inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,void,VECT_NORM>
-norm(const vector_dist_expression<prp1,v1> & va)
-{
-	vector_dist_expression_op<vector_dist_expression<prp1,v1>,void,VECT_NORM> exp_sum(va);
-
-	return exp_sum;
-}
 
 #include "vector_dist_operators_apply_kernel.hpp"
+#include "vector_dist_operators_functions.hpp"
+#include "vector_dist_operators_extensions.hpp"
 
 #endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_extensions.hpp b/src/Operators/Vector/vector_dist_operators_extensions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a3df229895af9d4b2a171c475c8e29733d821dda
--- /dev/null
+++ b/src/Operators/Vector/vector_dist_operators_extensions.hpp
@@ -0,0 +1,63 @@
+/*
+ * vector_dist_operators_extensions.hpp
+ *
+ *  Created on: Jul 18, 2016
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_EXTENSIONS_HPP_
+#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_EXTENSIONS_HPP_
+
+
+/*! \Create an expression from a Point
+ *
+ * \tpatam prp property
+ * \param v
+ *
+ */
+template <unsigned int dim, typename T> inline vector_dist_expression<16384,Point<dim,T> > getVExpr(Point<dim,T> & v)
+{
+	vector_dist_expression<(unsigned int)16384,Point<dim,T>> exp_v(v);
+
+	return exp_v;
+}
+
+
+/*! \brief Main class that encapsulate a vector properties
+ *
+ * \tparam prp property involved
+ * \tparam vector involved
+ *
+ */
+template<typename point>
+class vector_dist_expression<16384,point>
+{
+	point p;
+
+public:
+
+	vector_dist_expression(point p)
+	:p(p)
+	{}
+
+	/*! \brief This function must be called before value
+	 *
+	 * it initialize the expression if needed
+	 *
+	 */
+	inline void init() const
+	{}
+
+	/*! \brief Evaluate the expression
+	 *
+	 * \param key where to evaluate the expression
+	 *
+	 */
+	inline point value(const vect_dist_key_dx & k) const
+	{
+		return p;
+	}
+};
+
+
+#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_EXTENSIONS_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_functions.hpp b/src/Operators/Vector/vector_dist_operators_functions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c53919baa6a8c66e9dae2542c6bd0d322b45662e
--- /dev/null
+++ b/src/Operators/Vector/vector_dist_operators_functions.hpp
@@ -0,0 +1,215 @@
+/*
+ * vector_dist_operators_functions.hpp
+ *
+ *  Created on: Jul 17, 2016
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_FUNCTIONS_HPP_
+#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_FUNCTIONS_HPP_
+
+/*! A macro to define single value function specialization that apply the function component-wise
+ *
+ * \param fun function name
+ * \param ID function ID
+ *
+ */
+#define CREATE_VDIST_ARG_FUNC(fun_base,fun_name,OP_ID) \
+\
+\
+template <typename exp1>\
+class vector_dist_expression_op<exp1,void,OP_ID>\
+{\
+	const exp1 o1;\
+\
+public:\
+\
+	vector_dist_expression_op(const exp1 & o1)\
+	:o1(o1)\
+	{}\
+\
+	inline void init() const\
+	{\
+		o1.init();\
+	}\
+\
+	template<typename r_type=typename std::remove_reference<decltype(fun_base(o1.value(vect_dist_key_dx(0))))>::type > inline r_type value(const vect_dist_key_dx & key) const\
+	{\
+		return fun_base(o1.value(key));\
+	}\
+};\
+\
+\
+template<typename exp1, typename exp2_, unsigned int op1>\
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,void,OP_ID>\
+fun_name(const vector_dist_expression_op<exp1,exp2_,op1> & va)\
+{\
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,void,OP_ID> exp_sum(va);\
+\
+	return exp_sum;\
+}\
+\
+template<typename exp1, typename exp2_, unsigned int op1>\
+inline vector_dist_expression_op<vector_dist_expression<0,double>,void,OP_ID>\
+fun_name(double d)\
+{\
+	vector_dist_expression_op<vector_dist_expression<0,double>,void,OP_ID> exp_sum( (vector_dist_expression<0,double>(d)) );\
+\
+	return exp_sum;\
+}\
+\
+template<unsigned int prp1, typename v1>\
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,void,OP_ID>\
+fun_name(const vector_dist_expression<prp1,v1> & va)\
+{\
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,void,OP_ID> exp_sum(va);\
+\
+	return exp_sum;\
+}
+
+
+CREATE_VDIST_ARG_FUNC(norm,norm,VECT_NORM)
+CREATE_VDIST_ARG_FUNC(norm2,norm2,VECT_NORM2)
+CREATE_VDIST_ARG_FUNC(abs,abs,POINT_ABS)
+CREATE_VDIST_ARG_FUNC(exp,exp,POINT_EXP)
+CREATE_VDIST_ARG_FUNC(exp2,exp2,POINT_EXP2)
+CREATE_VDIST_ARG_FUNC(expm1,expm1,POINT_EXPM1)
+CREATE_VDIST_ARG_FUNC(log,log,POINT_LOG)
+CREATE_VDIST_ARG_FUNC(log10,log10,POINT_LOG10)
+CREATE_VDIST_ARG_FUNC(log2,log2,POINT_LOG2)
+CREATE_VDIST_ARG_FUNC(log1p,log1p,POINT_LOG1P)
+CREATE_VDIST_ARG_FUNC(sqrt,sqrt,POINT_SQRT)
+CREATE_VDIST_ARG_FUNC(cbrt,cbrt,POINT_CBRT)
+CREATE_VDIST_ARG_FUNC(sin,sin,POINT_SIN)
+CREATE_VDIST_ARG_FUNC(cos,cos,POINT_COS)
+CREATE_VDIST_ARG_FUNC(tan,tan,POINT_TAN)
+CREATE_VDIST_ARG_FUNC(asin,asin,POINT_ASIN)
+CREATE_VDIST_ARG_FUNC(acos,acos,POINT_ACOS)
+CREATE_VDIST_ARG_FUNC(atan,atan,POINT_ATAN)
+CREATE_VDIST_ARG_FUNC(sinh,sinh,POINT_SINH)
+CREATE_VDIST_ARG_FUNC(cosh,cosh,POINT_COSH)
+CREATE_VDIST_ARG_FUNC(tanh,tanh,POINT_TANH)
+CREATE_VDIST_ARG_FUNC(asinh,asinh,POINT_ASINH)
+CREATE_VDIST_ARG_FUNC(acosh,acosh,POINT_ACOSH)
+CREATE_VDIST_ARG_FUNC(atanh,atanh,POINT_ATANH)
+CREATE_VDIST_ARG_FUNC(erf,erf,POINT_ERF)
+CREATE_VDIST_ARG_FUNC(erfc,erfc,POINT_ERFC)
+CREATE_VDIST_ARG_FUNC(tgamma,tgamma,POINT_TGAMMA)
+CREATE_VDIST_ARG_FUNC(lgamma,lgamma,POINT_LGAMMA)
+CREATE_VDIST_ARG_FUNC(ceil,ceil,POINT_CEIL)
+CREATE_VDIST_ARG_FUNC(floor,floor,POINT_FLOOR)
+CREATE_VDIST_ARG_FUNC(trunc,trunc,POINT_TRUNC)
+CREATE_VDIST_ARG_FUNC(round,round,POINT_ROUND)
+CREATE_VDIST_ARG_FUNC(nearbyint,nearbyint,POINT_NEARBYINT)
+CREATE_VDIST_ARG_FUNC(rint,rint,POINT_RINT)
+
+
+/*! A macro to define single value function specialization that apply the function component-wise
+ *
+ * \param fun function name
+ * \param ID function ID
+ *
+ */
+#define CREATE_VDIST_ARG2_FUNC(fun_base,fun_name,OP_ID) \
+\
+\
+template <typename exp1,typename exp2>\
+class vector_dist_expression_op<exp1,exp2,OP_ID>\
+{\
+	const exp1 o1;\
+	const exp2 o2;\
+\
+public:\
+\
+	vector_dist_expression_op(const exp1 & o1, const exp2 & o2)\
+	:o1(o1),o2(o2)\
+	{}\
+\
+	inline void init() const\
+	{\
+		o1.init();\
+		o2.init();\
+	}\
+\
+	template<typename r_type=typename std::remove_reference<decltype(fun_base(o1.value(vect_dist_key_dx(0)),o2.value(vect_dist_key_dx(0)) ))>::type > inline r_type value(const vect_dist_key_dx & key) const\
+	{\
+		return fun_base(o1.value(key),o2.value(key));\
+	}\
+};\
+\
+\
+template<unsigned int p1, unsigned int p2, typename v1, typename v2>\
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,OP_ID>\
+fun_name(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)\
+{\
+	vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,OP_ID> exp_sum(va,vb);\
+\
+	return exp_sum;\
+}\
+\
+template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>\
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,OP_ID>\
+fun_name(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<prp1,v1> & vb)\
+{\
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,OP_ID> exp_sum(va,vb);\
+\
+	return exp_sum;\
+}\
+\
+template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>\
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,OP_ID>\
+fun_name(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)\
+{\
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,OP_ID> exp_sum(va,vb);\
+\
+	return exp_sum;\
+}\
+\
+template<typename exp1 , typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>\
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,OP_ID>\
+fun_name(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)\
+{\
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,OP_ID> exp_sum(va,vb);\
+\
+	return exp_sum;\
+}\
+\
+template<unsigned int prp1 , typename v1>\
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,OP_ID>\
+fun_name(const vector_dist_expression<prp1,v1> & va, double d)\
+{\
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,OP_ID> exp_sum(va,vector_dist_expression<0,double>(d));\
+\
+	return exp_sum;\
+}\
+\
+template<unsigned int prp1 , typename v1>\
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,OP_ID>\
+fun_name(double d, const vector_dist_expression<prp1,v1> & vb)\
+{\
+	vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,OP_ID> exp_sum(vector_dist_expression<0,double>(d),vb);\
+\
+	return exp_sum;\
+}\
+\
+template<typename exp1 , typename exp2, unsigned int op1>\
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,OP_ID>\
+fun_name(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)\
+{\
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,OP_ID> exp_sum(va,vector_dist_expression<0,double>(d));\
+\
+	return exp_sum;\
+}\
+\
+template<typename exp1 , typename exp2, unsigned int op1>\
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression_op<exp1,exp2,op1>,OP_ID>\
+fun_name(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)\
+{\
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,OP_ID> exp_sum(vector_dist_expression<0,double>(d),va);\
+\
+	return exp_sum;\
+}
+
+CREATE_VDIST_ARG2_FUNC(pmul,pmul,VECT_PMUL)
+
+#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_FUNCTIONS_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_unit_tests.hpp b/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
index b8eff9a0e7ad6baf26751e9ecce4afcc3a27ac89..208f7c8390b2d9a7dc50ed51b1a40f3a51d732c4 100644
--- a/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
+++ b/src/Operators/Vector/vector_dist_operators_unit_tests.hpp
@@ -61,6 +61,94 @@ template <typename vector> bool check_values_complex_expr(vector & vd)
 	return ret;
 }
 
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_pos_sum(vector & vd, const rtype & p)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		rtype base1 = vd.template getPos(key) + p;
+		rtype base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_pos_sub(vector & vd, const rtype & p)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		rtype base1 = vd.template getPos(key) - p;
+		rtype base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_pos_sub_minus(vector & vd, const rtype & p)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		rtype base1 = -(vd.template getPos(key) - p);
+		rtype base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_point_sub(vector & vd, const rtype & p)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		rtype base1 = -vd.template getProp<B>(key);
+		rtype base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
 template <typename rtype, typename vector, unsigned int A, unsigned int B, unsigned int C> bool check_values_sum(vector & vd, double d)
 {
 	bool ret = true;
@@ -622,6 +710,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 	auto vVB = getV<VB>(vd);
 	auto vVC = getV<VC>(vd);
 
+	auto vPOS = getV<PROP_POS>(vd);
+
 	vA = 1.0;
 	vB = 2.0f;
 	vC = 3.0;
@@ -630,6 +720,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 	check_values<B>(vd,2.0);
 	check_values<C>(vd,3.0);
 
+	vA = vB;
+	check_values<A>(vd,2.0);
+
 	fill_values(vd);
 
 	vA = vB + 2.0 + vB - 2.0*vC / 5.0;
@@ -699,30 +792,30 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 
 	// Various combination of 2 operator
 
-	vVA = vVB + 2.0f;
+	vVA = vVB + 2.0;
 	check_values_sum<VectorS<3,float>,vtype,VA,VB,VC>(vd,2.0f);
-	vVA = 2.0f + vVB;
+	vVA = 2.0 + vVB;
 	check_values_sum<VectorS<3,float>,vtype,VA,VB,VC>(vd,2.0f);
 	vVA = vVC + vVB;
 	check_values_sum<VectorS<3,float>,vtype,VA,VB,VC>(vd,vd);
 
-	vVA = vVB - 2.0f;
+	vVA = vVB - 2.0;
 	check_values_sub<VectorS<3,float>,vtype,VA,VB,VC>(vd,2.0f);
-	vVA = 2.0f - vVB;
+	vVA = 2.0 - vVB;
 	check_values_sub<VectorS<3,float>,vtype,VA,VB,VC>(2.0f,vd);
 	vVA = vVC - vVB;
 	check_values_sub<VectorS<3,float>,vtype,VA,VB,VC>(vd,vd);
 
-	vVA = vVB * 2.0f;
+	vVA = vVB * 2.0;
 	check_values_mul<VectorS<3,float>,vtype,VA,VB,VC>(vd,2.0f);
-	vVA = 2.0f * vVB;
+	vVA = 2.0 * vVB;
 	check_values_mul<VectorS<3,float>,vtype,VA,VB,VC>(vd,2.0f);
 	vVA = vVC * vVB;
 	check_values_mul<VectorS<3,float>,vtype,VA,VB,VC>(vd,vd);
 
-	vVA = vVB / 2.0f;
+	vVA = vVB / 2.0;
 	check_values_div<VectorS<3,float>,vtype,VA,VB,VC>(vd,2.0f);
-	vVA = 2.0f / vVB;
+	vVA = 2.0 / vVB;
 	check_values_div<VectorS<3,float>,vtype,VA,VB,VC>(2.0f,vd);
 	vVA = vVC / vVB;
 	check_values_div<VectorS<3,float>,vtype,VA,VB,VC>(vd,vd);
@@ -767,6 +860,61 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 
 	vA = vVB * vVC + norm(vVC + vVB) + distance(vVC,vVB);
 	check_values_scal_norm_dist(vd);
+
+	Point<3,float> p0({2.0,2.0,2.0});
+	auto p0_e = getVExpr(p0);
+
+	vVA = vPOS + p0_e;
+	check_values_pos_sum<VectorS<3,float>,vtype,VA,VB,VC>(vd,p0);
+
+	vVA = vPOS - p0_e;
+	check_values_pos_sub<Point<3,float>,vtype,VA,VB,VC>(vd,p0);
+
+	vVA = -(vPOS - p0_e);
+	check_values_pos_sub_minus<Point<3,float>,vtype,VA,VB,VC>(vd,p0);
+
+	vVA = -vVB;
+	check_values_point_sub<Point<3,float>,vtype,VA,VB,VC>(vd,p0);
+
+	// Just check it compile testing it will test the same code
+	// as the previuous one
+	vVC = exp(vVB);
+	vA = norm(vPOS);
+	vVA = vPOS + 2.0;
+	vVA = 2.0 + vPOS;
+//	vVA = vPOS + vPOS;
+	vVA = vPOS - 2.0f;
+	vVA = 2.0 - vPOS;
+//	vVA = vPOS - vPOS;
+
+/*	vVA = vPOS * 2.0;
+	vVA = 2.0 * vPOS;
+	vVA = vPOS * vPOS;
+
+	vVA = vPOS / 2.0f;
+	vVA = 2.0f / vPOS;
+	vVA = vPOS / vPOS;
+
+	// Variuos combination 3 operator
+
+	vVA = vPOS + (vPOS + vPOS);
+	vVA = (vPOS + vPOS) + vPOS;
+	vVA = (vPOS + vPOS) + (vPOS + vPOS);
+
+	vVA = vPOS - (vPOS + vPOS);
+	vVA = (vPOS + vPOS) - vPOS;
+	vVA = (vVC + vPOS) - (vPOS + vPOS);
+
+	vVA = vPOS * (vPOS + vPOS);
+	vVA = (vPOS + vPOS) * vPOS;
+	vVA = (vPOS + vPOS) * (vPOS + vPOS);
+	vA = vPOS * (vPOS + vPOS);
+	vA = (vPOS + vPOS) * vPOS;
+	vA = (vPOS + vPOS) * (vPOS + vPOS);
+
+	vVA = vPOS / (vPOS + vPOS);
+	vVA = (vPOS + vPOS) / vPOS;
+	vVA = (vPOS + vPOS) / (vPOS + vPOS);*/
 }
 
 #include "vector_dist_operators_apply_kernel_unit_tests.hpp"