Point_operators.hpp

/*
 * Point_operators.hpp
 *
 *  Created on: Jun 14, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_DATA_SRC_SPACE_SHAPE_POINT_OPERATORS_HPP_
#define OPENFPM_DATA_SRC_SPACE_SHAPE_POINT_OPERATORS_HPP_

#include "util/multi_array_openfpm/array_openfpm.hpp"

template<unsigned int dim ,typename T> class Point;

#define POINT_SUM 1
#define POINT_SUB 2
#define POINT_MUL 3
#define POINT_DIV 4
#define POINT_MUL_POINT 5
#define POINT_NORM 6
#define POINT_NORM2 7

// Functions

#define POINT_ABS 8
#define POINT_EXP 9
#define POINT_EXP2 10
#define POINT_EXPM1 11
#define POINT_LOG 12
#define POINT_LOG10 13
#define POINT_LOG2 14
#define POINT_LOG1P 15
#define POINT_SQRT 17
#define POINT_CBRT 18
#define POINT_SIN 19
#define POINT_COS 20
#define POINT_TAN 21
#define POINT_ASIN 22
#define POINT_ACOS 23
#define POINT_ATAN 24
#define POINT_SINH 25
#define POINT_COSH 26
#define POINT_TANH 27
#define POINT_ASINH 28
#define POINT_ACOSH 29
#define POINT_ATANH 30
#define POINT_ERF 31
#define POINT_ERFC 32
#define POINT_TGAMMA 33
#define POINT_LGAMMA 34
#define POINT_CEIL 35
#define POINT_FLOOR 36
#define POINT_TRUNC 37
#define POINT_ROUND 38
#define POINT_NEARBYINT 39
#define POINT_RINT 40
#define POINT_SUB_UNI 41
#define POINT_NORM_INF 42



template<bool cond, typename exp1, typename exp2>
struct first_or_second_pt
{
    typedef typename exp2::coord_type coord_type;
};

template<typename exp1, typename exp2>
struct first_or_second_pt<true,exp1,exp2>
{
    typedef typename exp1::coord_type coord_type;
};

template<typename T, typename Sfinae = void>
struct has_coordtype: std::false_type {};

template<typename T>
struct has_coordtype<T, typename Void<typename T::coord_type>::type> : std::true_type
{};

template<typename source1, typename source2>
struct best_conv_impl
{
    typedef source1 type;
};

template<typename source2>
struct best_conv_impl<int,source2>
{
    typedef source2 type;
};

template<typename source2>
struct best_conv_impl<long int,source2>
{
    typedef source2 type;
};

template<typename source2>
struct best_conv_impl<unsigned int,source2>
{
    typedef source2 type;
};

template<typename source2>
struct best_conv_impl<unsigned long int,source2>
{
    typedef source2 type;
};

template<typename source1>
struct best_conv_impl<source1,int>
{
    typedef source1 type;
};

template<typename source1>
struct best_conv_impl<source1,long int>
{
    typedef source1 type;
};

template<typename source1>
struct best_conv_impl<source1,unsigned int>
{
    typedef source1 type;
};

template<typename source1>
struct best_conv_impl<source1,unsigned long int>
{
    typedef source1 type;
};

template<typename source1, typename source2>
struct best_conv
{
    typedef typename best_conv_impl<typename std::remove_const<source1>::type,typename std::remove_const<source2>::type>::type type;
};


constexpr unsigned int max_expr(unsigned int dim1, unsigned int dim2)
{
    return (dim1 > dim2)?dim1:dim2;
}

template <unsigned int op1_dim, unsigned int op2_dim, unsigned int op>
struct r_type_dim
{
    enum
    {
        value = max_expr(op1_dim,op2_dim),
    };
};

template <>
struct r_type_dim<1,1,POINT_SUM>
{
    enum
    {
        value = 1,
    };
};

template <>
struct r_type_dim<1,1,POINT_SUB>
{
    enum
    {
        value = 1,
    };
};

template <unsigned int op1_dim, unsigned int op2_dim>
struct r_type_dim<op1_dim,op2_dim,POINT_MUL_POINT>
{
    enum
    {
        value = 1,
    };
};

template <>
struct r_type_dim<1,1,POINT_MUL>
{
    enum
    {
        value = 1,
    };
};

template <>
struct r_type_dim<1,1,POINT_DIV>
{
    enum
    {
        value = 1,
    };
};

template <unsigned int r, typename orig>
struct r_type_p
{
    typedef orig type;
};

template <typename orig>
struct r_type_p<1,orig>
{
    typedef typename orig::coord_type type;
};


template<typename T>
class point_expression
{
    T d;

public:

    typedef int has_init;

    typedef int is_expression;

    static const unsigned int nvals = 1;

    typedef T coord_type;

    __device__ __host__ inline point_expression(T & d)
    :d(d)
    {}

    __device__ __host__  inline void init() const
    {
    }

    __device__ __host__  inline const T & value(const int k) const
    {
        return d;
    }

    __device__ __host__  inline T & value(const int k)
    {
        return d;
    }
};



template <typename orig, typename exp1, typename exp2, unsigned int op>
class point_expression_op
{

};

template <typename orig, typename exp1, typename exp2>
class point_expression_op<orig,exp1,exp2,POINT_SUM>
{
    const exp1 o1;
    const exp2 o2;

public:

    typedef typename first_or_second_pt<has_coordtype<exp1>::value,exp1,exp2>::coord_type coord_type;

    typedef orig orig_type;

    typedef int is_expression;

    typedef int has_init;

    typedef typename r_type_p<r_type_dim<exp1::nvals,exp2::nvals,POINT_SUM>::value,orig >::type return_type;

    static const unsigned int nvals = r_type_dim<exp1::nvals,exp2::nvals,POINT_SUM>::value;

    __device__ __host__ inline point_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    __device__ __host__  inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename best_conv<typename std::remove_reference<decltype(o1.value(0))>::type,
                                                typename std::remove_reference<decltype(o2.value(0))>::type>::type >
    __device__ __host__  inline r_type value(size_t k) const
    {
        return o1.value(k) + o2.value(k);
    }

    template<typename T, typename test=typename boost::disable_if_c< std::is_same<T,orig>::value || exp1::nvals != 1 || exp2::nvals != 1 >::type  >
    __device__ __host__  inline operator T() const
    {
        init();
        return o1.value(0) + o2.value(0);
    }

    template<typename r_type=typename best_conv<typename std::remove_reference<decltype(o1.value(0))>::type,
            typename std::remove_reference<decltype(o2.value(0))>::type>::type >
    __device__ __host__  inline r_type operator[](size_t i)
    {
        return o1.value(i) + o2.value(i);
    }
};

template <typename orig,typename exp1, typename exp2>
class point_expression_op<orig, exp1,exp2,POINT_SUB>
{
    const exp1 o1;
    const exp2 o2;

public:

    typedef typename first_or_second_pt<has_coordtype<exp1>::value,exp1,exp2>::coord_type coord_type;

    typedef orig orig_type;

    typedef int has_init;

    typedef int is_expression;

    typedef orig return_type;

    static const unsigned int nvals = r_type_dim<exp1::nvals,exp2::nvals,POINT_SUB>::value;

    __device__ __host__  inline point_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    __device__ __host__  inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename best_conv<typename std::remove_reference<decltype(o1.value(0))>::type,
                                                typename std::remove_reference<decltype(o2.value(0))>::type>::type >
    __device__ __host__  inline r_type value(size_t k) const
    {
        return o1.value(k) - o2.value(k);
    }

    template<typename T, typename test=typename boost::disable_if_c< std::is_same<T,orig>::value || exp1::nvals != 1 || exp2::nvals != 1 >::type  >
    __device__ __host__  operator T() const
    {
        init();
        return o1.value(0) - o2.value(0);
    }


};

template <typename orig, typename exp1, typename exp2>
class point_expression_op<orig,exp1,exp2, POINT_SUB_UNI >
{
    const exp1 o1;

    mutable typename orig::coord_type scal;

public:

    typedef orig orig_type;

    typedef int is_expression;

    typedef int has_init;

    typedef typename orig::coord_type return_type;

    static const unsigned int nvals = exp1::nvals;

    typedef typename first_or_second_pt<has_coordtype<exp1>::value,exp1,exp2>::coord_type coord_type;

    __device__ __host__  inline point_expression_op(const exp1 & o1)
    :o1(o1),scal(0.0)
    {}

    __device__ __host__  inline void init() const
    {
        o1.init();
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(0))>::type >
    __device__ __host__  inline r_type value(size_t k) const
    {
        return -(o1.value(k));
    }

    template <typename T, typename check = typename std::enable_if<!std::is_same<T,orig>::value >::type  >
    __device__ __host__  operator T() const
    {
        init();
        return -(o1.value(0));
    }
};


template <typename orig, typename exp1, typename exp2>
class point_expression_op<orig,exp1,exp2,POINT_MUL_POINT>
{
    const exp1 o1;
    const exp2 o2;

    mutable typename std::remove_const<typename orig::coord_type>::type scal;

public:

    typedef typename first_or_second_pt<has_coordtype<exp1>::value,exp1,exp2>::coord_type coord_type;

    typedef orig orig_type;

    typedef int has_init;

    typedef int is_expression;

    typedef typename orig::coord_type return_type;

    static const unsigned int nvals = 1;

    __device__ __host__  inline point_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2),scal(0.0)
    {}

    __device__ __host__  inline void init() const
    {
        o1.init();
        o2.init();

        for (size_t i = 0 ; i < orig::dims ; i++)
            scal += o1.value(i) * o2.value(i);
    }

    template<typename r_type=typename best_conv<typename std::remove_reference<decltype(o1.value(0))>::type,
                                       typename std::remove_reference<decltype(o2.value(0))>::type>::type >
    __device__ __host__  inline r_type value(size_t k) const
    {
        return scal;
    }

    template<typename T, typename test=typename boost::disable_if_c< std::is_same<T,orig>::value >::type >
    __device__ __host__  operator T() const
    {
        init();
        return scal;
    }
};


template <typename orig, typename exp1, typename exp2>
class point_expression_op<orig,exp1,exp2,POINT_MUL>
{
    const exp1 o1;
    const exp2 o2;

public:

    typedef typename first_or_second_pt<has_coordtype<exp1>::value,exp1,exp2>::coord_type coord_type;

    typedef orig orig_type;

    typedef int has_init;

    typedef int is_expression;

    typedef orig return_type;

    static const unsigned int nvals = r_type_dim<exp1::nvals,exp2::nvals,POINT_MUL>::value;

    __device__ __host__  inline point_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    __device__ __host__  inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename best_conv<typename std::remove_reference<decltype(o1.value(0))>::type,
                                       typename std::remove_reference<decltype(o2.value(0))>::type>::type >
    __device__ __host__  inline r_type value(size_t k) const
    {
        return o1.value(k) * o2.value(k);
    }

    template<typename T, typename test=typename boost::disable_if_c< std::is_same<T,orig>::value || exp1::nvals != 1 || exp2::nvals != 1 >::type  >
    __device__ __host__  operator T() const
    {
        init();
        return o1.value(0) * o2.value(0);
    }

};

template <typename orig, typename exp1, typename exp2>
class point_expression_op<orig,exp1,exp2,POINT_DIV>
{
    const exp1 o1;
    const exp2 o2;

public:


    typedef typename first_or_second_pt<has_coordtype<exp1>::value,exp1,exp2>::coord_type coord_type;

    typedef orig orig_type;

    typedef int is_expression;

    typedef int has_init;

    typedef orig return_type;

    static const unsigned int nvals = r_type_dim<exp1::nvals,exp2::nvals,POINT_DIV>::value;

    __device__ __host__  inline point_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    __device__ __host__  inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename best_conv<typename std::remove_reference<decltype(o1.value(0))>::type,
                                       typename std::remove_reference<decltype(o2.value(0))>::type>::type >
    __device__ __host__  inline r_type value(size_t k) const
    {
        return o1.value(k) / o2.value(k);
    }


    template<typename T, typename test=typename boost::disable_if_c< std::is_same<T,orig>::value || exp1::nvals != 1 || exp2::nvals != 1 >::type  > __device__ __host__  operator T() const
    {
        init();
        return o1.value(0) / o2.value(0);
    }
};

template<unsigned int dim, typename T> __device__ __host__  point_expression<T[dim]> getExprL(T (& a)[dim])
{
    return point_expression<T[dim]>(a);
}

template<unsigned int dim, typename T> __device__ __host__  point_expression<const T[dim]> getExprR(T (& a)[dim])
{
    return point_expression<const T[dim]>(a);
}

template<unsigned int dim, typename T>
struct ger
{
    template<typename vmpl> __device__ __host__
    __device__ __host__ static point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>
    getExprR(const openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl> & a)
    {
        return point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>(a);
    }

    __device__ __host__ static point_expression<const T[dim]>
    getExprR(T (& a)[dim])
    {
        return point_expression<const T[dim]>(a);
    }

    template<typename vmpl> __device__ __host__
    __device__ __host__ static point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>
    getExprL(const openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl> & a)
    {
        return point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>(a);
    }

    __device__ __host__ static point_expression<T[dim]>
    getExprL(T (& a)[dim])
    {
        return point_expression<T[dim]>(a);
    }
};


#define CREATE_POINT_OPERATOR(operator_name,OP_ID) \
\
\
template<unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<const T[dim]>,OP_ID>\
operator_name(const Point<dim,T> & va, const point_expression<const T[(unsigned int)dim]> & vb)\
{\
    point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<const T[dim]>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T, typename vmpl>\
__device__ __host__  inline point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,OP_ID>\
operator_name(const Point<dim,T> & va, const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & vb)\
{\
    point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<Point<dim,T>,point_expression<const T[dim]>,Point<dim,T>,OP_ID>\
operator_name(const point_expression<const T[(unsigned int)dim]> & va, const Point<dim,T> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<const T[dim]>,Point<dim,T>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T, typename vmpl>\
__device__ __host__  inline point_expression_op<Point<dim,T>,point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,Point<dim,T>,OP_ID>\
operator_name(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va,\
              const Point<dim,T> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,Point<dim,T>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<Point<dim,T>,point_expression<const T[dim]>,point_expression<double>,OP_ID>\
operator_name(const point_expression<const T[dim]> & va, double d)\
{\
    point_expression_op<Point<dim,T>,point_expression<const T[dim]>,point_expression<double>,OP_ID> exp_sum(va,point_expression<double>(d));\
\
    return exp_sum;\
}\
\
template<typename T, typename vmpl>\
__device__ __host__  inline point_expression_op<Point<subar_dim<vmpl>::type::value,T>,\
                                                point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                                                point_expression<double>,\
                                                OP_ID>\
operator_name(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va, double d)\
{\
    point_expression_op<Point<subar_dim<vmpl>::type::value,T>,\
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                        point_expression<double>,\
                        OP_ID>\
                        exp_sum(va,point_expression<double>(d));\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<Point<dim,T>,point_expression<double>,point_expression<const T[dim]>,OP_ID>\
operator_name(double d, const point_expression<const T[dim]> & va)\
{\
    point_expression_op<Point<dim,T>,point_expression<double>,point_expression<const T[dim]>,OP_ID> exp_sum(point_expression<double>(d),va);\
\
    return exp_sum;\
}\
\
template<typename T, typename vmpl>\
__device__ __host__  inline point_expression_op<Point<subar_dim<vmpl>::type::value,T>,\
                                                point_expression<double>,\
                                                point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                                                OP_ID>\
operator_name(double d, const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va)\
{\
    point_expression_op<Point<subar_dim<vmpl>::type::value,T>,\
                        point_expression<double>,\
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                        OP_ID>\
                        exp_sum(point_expression<double>(d),va);\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<Point<dim,T>,point_expression<const T[dim]>,point_expression<const T[dim]>,OP_ID>\
operator_name(const point_expression<const T[dim]> & va, const point_expression<const T[dim]> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<const T[dim]>,point_expression<const T[dim]>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename T, typename vmpl>\
__device__ __host__  inline point_expression_op<Point<subar_dim<vmpl>::type::value,T>,\
                                                point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                                                point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                                                OP_ID>\
operator_name(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va,\
              const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & vb)\
{\
    point_expression_op<Point<subar_dim<vmpl>::type::value,T>,\
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                        OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<Point<dim,T>,point_expression<const T[dim]>,point_expression<T[dim]>,OP_ID>\
operator_name(const point_expression<const T[dim]> & va, const point_expression<T[dim]> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<const T[dim]>,point_expression<T[dim]>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T[dim]>,OP_ID>\
operator_name(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression<T[dim]> & vb)\
{\
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T[dim]>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, typename T, typename vmpl>\
__device__ __host__  inline point_expression_op<orig,\
                                                point_expression_op<orig,exp1,exp2,op1>,\
                                                point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                                                OP_ID>\
operator_name(const point_expression_op<orig,exp1,exp2,op1> & va,\
              const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & vb)\
{\
    point_expression_op<orig,\
                        point_expression_op<orig,exp1,exp2,op1>,\
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                        OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, unsigned int dim, typename T>\
__device__ __host__ inline point_expression_op<orig,point_expression<T[dim]>,point_expression_op<orig,exp1,exp2,op1>,OP_ID>\
operator_name(const point_expression<T[dim]> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)\
{\
    point_expression_op<orig,point_expression<T[dim]>,point_expression_op<orig,exp1,exp2,op1>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, typename T, typename vmpl>\
__device__ __host__ inline point_expression_op<orig,\
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                                               point_expression_op<orig,exp1,exp2,op1>,\
                                               OP_ID>\
operator_name(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va,\
              const point_expression_op<orig,exp1,exp2,op1> & vb)\
{\
    point_expression_op<orig,\
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,\
                        point_expression_op<orig,exp1,exp2,op1>,\
                        OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
\
template<unsigned int dim, typename T>\
__device__ __host__  inline point_expression_op<Point<dim,T>,Point<dim,T>,Point<dim,T>,OP_ID>\
operator_name(const Point<dim,T> & va, const Point<dim,T> & vb)\
{\
    point_expression_op<Point<dim,T>,Point<dim,T>,Point<dim,T>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, unsigned int dim, typename T>\
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,OP_ID>\
operator_name(const point_expression_op<orig,exp1,exp2,op1> & va, const Point<dim,T> & vb)\
{\
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename orig,typename exp1 , typename exp2, unsigned int op1, unsigned int dim, typename T>\
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,OP_ID>\
operator_name(const Point<dim,T> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)\
{\
    point_expression_op<Point<dim,T>,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, typename T>\
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T>,OP_ID>\
operator_name(const point_expression_op<orig,exp1,exp2,op1> & va, T d)\
{\
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T>,OP_ID> exp_sum(va,d);\
\
    return exp_sum;\
}\
\
template<typename orig,typename exp1 , typename exp2, unsigned int op1, typename T>\
__device__ __host__ inline point_expression_op<orig,point_expression<T>,point_expression_op<orig,exp1,exp2,op1>,OP_ID>\
operator_name(T d, const point_expression_op<orig,exp1,exp2,op1> & vb)\
{\
    point_expression_op<orig,point_expression<T>,point_expression_op<orig,exp1,exp2,op1>,OP_ID> exp_sum(d,vb);\
\
    return exp_sum;\
}\
\
template<typename orig,typename exp1 , typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>\
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,OP_ID>\
operator_name(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression_op<orig,exp3,exp4,op2> & vb)\
{\
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,OP_ID> exp_sum(va,vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim , typename T, typename check=typename std::enable_if< !std::is_same<T,double>::value >::type >\
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<T>,OP_ID>\
operator_name(const Point<dim,T> & va, T d)\
{\
    point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<T>,OP_ID> exp_sum(va,point_expression<T>(d));\
\
    return exp_sum;\
}\
\
template<unsigned int dim , typename T>\
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<double>,OP_ID>\
operator_name(const Point<dim,T> & va, double d)\
{\
    point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<double>,OP_ID> exp_sum(va,point_expression<double>(d));\
\
    return exp_sum;\
}\
\
template<unsigned int dim , typename T, typename check=typename std::enable_if< !std::is_same<T,double>::value >::type >\
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<T>,Point<dim,T>,OP_ID>\
operator_name(T d, const Point<dim,T> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<T>,Point<dim,T>,OP_ID> exp_sum(point_expression<T>(d),vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim, typename T, typename check=typename std::enable_if< !std::is_same<T,double>::value >::type>\
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<T[dim]>,Point<dim,T>,OP_ID>\
operator_name(T (& d)[dim], const Point<dim,T> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<T[dim]>,Point<dim,T>,OP_ID> exp_sum(point_expression<T[dim]>(d),vb);\
\
    return exp_sum;\
}\
\
template<unsigned int dim , typename T>\
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<double>,Point<dim,T>,OP_ID>\
operator_name(double d, const Point<dim,T> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<double>,Point<dim,T>,OP_ID> exp_sum(point_expression<double>(d),vb);\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, typename T, typename check=typename std::enable_if< !std::is_same<T,double>::value >::type >\
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T>,OP_ID>\
operator_name(const point_expression_op<orig,exp1,exp2,op1> & va, T d)\
{\
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T>,OP_ID> exp_sum(va,point_expression<T>(d));\
\
    return exp_sum;\
}\
\
template<typename orig, typename exp1 , typename exp2, unsigned int op1, typename T>\
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<double>,OP_ID>\
operator_name(const point_expression_op<orig,exp1,exp2,op1> & va, double d)\
{\
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<double>,OP_ID> exp_sum(va,point_expression<double>(d));\
\
    return exp_sum;\
}


CREATE_POINT_OPERATOR(operator+,POINT_SUM)
CREATE_POINT_OPERATOR(operator-,POINT_SUB)
CREATE_POINT_OPERATOR(operator/,POINT_DIV)

/* \brief sum two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig, typename exp1 , typename exp2, unsigned int op1, typename T, typename check=typename std::enable_if< ! std::is_same<T,orig>::value >::type >
__device__ __host__ inline T &
operator+=(T & d, const point_expression_op<orig,exp1,exp2,op1> & va)
{
    va.init();
    d += va.value(0);

    return d;
}



/* \brief sum two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig, typename exp1 , typename exp2, unsigned int op1, typename T>
__device__ __host__ inline T &
operator+=(orig & d, const point_expression_op<orig,exp1,exp2,op1> & va)
{
    va.init();

    d = d + va;

    return d;
}

/* \brief minus points expression
 *
 * \param va point expression one
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig, typename exp1, typename exp2, unsigned int op1>
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,void, POINT_SUB_UNI >
operator-(const point_expression_op<orig,exp1,exp2,op1> & va)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,void,POINT_SUB_UNI> exp_sum(va);

    return exp_sum;
}

/* \brief minus points expression
 *
 * \param va point expression one
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,void, POINT_SUB_UNI >
operator-(const Point<dim,T> & va)
{
    point_expression_op<Point<dim,T>,Point<dim,T>,void,POINT_SUB_UNI> exp_sum(va);

    return exp_sum;
}




template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<T[dim]>,point_expression<double>,POINT_MUL>
operator*(const point_expression<T[dim]> & va, double d)
{
    point_expression_op<Point<dim,T>,point_expression<T[dim]>,point_expression<double>,POINT_MUL> exp_sum(va,point_expression<double>(d));

    return exp_sum;
}

template<typename T, typename vmpl>
__device__ __host__ inline point_expression_op<Point<subar_dim<vmpl>::type::value,T>,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               point_expression<double>,
                                               POINT_MUL>
operator*(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va,
          double d)
{
    point_expression_op<Point<subar_dim<vmpl>::type::value,T>,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        point_expression<double>,
                        POINT_MUL> exp_sum(va,point_expression<double>(d));

    return exp_sum;
}

template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<double>,point_expression<T[dim]>,POINT_MUL>
operator*(double d,
          const point_expression<T[dim]> & va)
{
    point_expression_op<Point<dim,T>,point_expression<double>,point_expression<T[dim]>,POINT_MUL> exp_sum(point_expression<double>(d),va);

    return exp_sum;
}

template<typename T, typename vmpl>
__device__ __host__ inline point_expression_op<Point<subar_dim<vmpl>::type::value,T>,
                                               point_expression<double>,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               POINT_MUL>
operator*(double d,
          const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va)
{
    point_expression_op<Point<subar_dim<vmpl>::type::value,T>,
                        point_expression<double>,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        POINT_MUL> exp_sum(point_expression<double>(d),va);

    return exp_sum;
}

template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<T[dim]>,point_expression<T[dim]>,POINT_MUL_POINT>
operator*(const point_expression<T[dim]> & va, const point_expression<T[dim]> & vb)
{
    point_expression_op<Point<dim,T>,point_expression<T[dim]>,point_expression<T[dim]>,POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

template<typename T, typename vmpl>
__device__ __host__ inline point_expression_op<Point<subar_dim<vmpl>::type::value,T>,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               POINT_MUL_POINT>
operator*(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va,
          const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & vb)
{
    point_expression_op<Point<subar_dim<vmpl>::type::value,T>,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

template<unsigned int dim, typename T, typename check=typename std::enable_if< !std::is_same<T,double>::value >::type>
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<T>,Point<dim,T>,POINT_MUL>
operator*(T d, const Point<dim,T> & vb)
{
    point_expression_op<Point<dim,T>,point_expression<T>,Point<dim,T>,POINT_MUL> exp_sum(point_expression<T>(d),vb);

    return exp_sum;
}


template<typename orig,
         typename exp1 ,
         typename exp2,
         unsigned int op1,
         unsigned int dim,
         typename T,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == dim >::type >
__device__ __host__ inline point_expression_op<orig,
                                               point_expression_op<orig,exp1,exp2,op1>,
                                               point_expression<T[dim]>,
                                               POINT_MUL_POINT>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression<T[dim]> & vb)
{
    point_expression_op<orig,
                        point_expression_op<orig,exp1,exp2,op1>,
                        point_expression<T[dim]>,
                        POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

template<typename orig,
         typename exp1 ,
         typename exp2,
         unsigned int op1,
         typename T,
         typename vmpl,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == subar_dim<vmpl>::type::value >::type >
__device__ __host__ inline point_expression_op<orig,
                                               point_expression_op<orig,exp1,exp2,op1>,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               POINT_MUL_POINT>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & vb)
{
    point_expression_op<orig,
                        point_expression_op<orig,exp1,exp2,op1>,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

template<typename orig,
         typename exp1 ,
         typename exp2,
         unsigned int op1,
         unsigned int dim,
         typename T,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == 1 >::type >
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T[dim]>,POINT_MUL>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression<T[dim]> & vb)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T[dim]>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}

template<typename orig,
         typename exp1 ,
         typename exp2,
         unsigned int op1,
         typename T,
         typename vmpl,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == 1 >::type >
__device__ __host__ inline point_expression_op<orig,
                                               point_expression_op<orig,exp1,exp2,op1>,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               POINT_MUL>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & vb)
{
    point_expression_op<orig,
                        point_expression_op<orig,exp1,exp2,op1>,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}


template<typename orig,
         typename exp1,
         typename exp2,
         unsigned int op1,
         unsigned int dim,
         typename T,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == dim >::type >
__device__ __host__ inline point_expression_op<orig,point_expression<T[dim]>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL_POINT>
operator*(const point_expression<T[dim]> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,point_expression<T[dim]>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

template<typename orig,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename T,
         typename vmpl,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == subar_dim<vmpl>::type::value >::type >
__device__ __host__ inline point_expression_op<orig,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               point_expression_op<orig,exp1,exp2,op1>,
                                               POINT_MUL_POINT>
operator*(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        point_expression_op<orig,exp1,exp2,op1>,
                        POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

template<typename orig,
         typename exp1,
         typename exp2,
         unsigned int op1,
         unsigned int dim,
         typename T,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == 1 >::type >
__device__ __host__ inline point_expression_op<orig,point_expression<T[dim]>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL>
operator*(const point_expression<T[dim]> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,point_expression<T[dim]>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}

template<typename orig,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename T,
         typename vmpl,
         typename sfinae = typename std::enable_if< point_expression_op<orig,exp1,exp2,op1>::nvals == 1 >::type >
__device__ __host__ inline point_expression_op<orig,
                                               point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                                               point_expression_op<orig,exp1,exp2,op1>,
                                               POINT_MUL>
operator*(const point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,
                        point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>,
                        point_expression_op<orig,exp1,exp2,op1>,
                        POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}


/* \brief Multiply two points expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,point_expression<double>,Point<dim,T>,POINT_MUL>
operator*(double d, const Point<dim,T> & vb)
{
    point_expression_op<Point<dim,T>,point_expression<double>,Point<dim,T>,POINT_MUL> exp_sum(point_expression<double>(d),vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int dim, typename T, typename check=typename std::enable_if< !std::is_same<T,double>::value >::type >
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<T>,POINT_MUL>
operator*(const Point<dim,T> & va, T d)
{
    point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<T>,POINT_MUL> exp_sum(va,point_expression<T>(d));

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<double>,POINT_MUL>
operator*(const Point<dim,T> & va, double d)
{
    point_expression_op<Point<dim,T>,Point<dim,T>,point_expression<double>,POINT_MUL> exp_sum(va,point_expression<double>(d));

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,Point<dim,T>,POINT_MUL_POINT>
operator*(const Point<dim,T> & va, const Point<dim,T> & vb)
{
    point_expression_op<Point<dim,T>,Point<dim,T>,Point<dim,T>,POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig,
         unsigned int dim,
         typename T,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename sfinae = typename std::enable_if<point_expression_op<orig,exp1,exp2,op1>::nvals != 1>::type >
__device__ __host__  inline point_expression_op<orig,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL_POINT>
operator*(const Point<dim,T> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig,
         unsigned int dim,
         typename T,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename sfinae = typename std::enable_if<point_expression_op<orig,exp1,exp2,op1>::nvals == 1>::type >
__device__ __host__ inline point_expression_op<orig,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL>
operator*(const Point<dim,T> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig,
         unsigned int dim,
         typename T,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename sfinae = typename std::enable_if<point_expression_op<orig,exp1,exp2,op1>::nvals != 1>::type >
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,POINT_MUL_POINT>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const Point<dim,T> & vb)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig,
         unsigned int dim,
         typename T,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename check = typename std::enable_if<point_expression_op<orig,exp1,exp2,op1>::nvals == 1>::type >
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,POINT_MUL>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const Point<dim,T> & vb)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply a point expression with a number
 *
 * \param d number
 * \param vb point expression
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig, typename T, typename exp1, typename exp2, unsigned int op1>
__device__ __host__ inline point_expression_op<orig,point_expression<T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL>
operator*(T d, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,point_expression<T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL> exp_sum(point_expression<T>(d),vb);

    return exp_sum;
}


/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param d constant value
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig, typename exp1, typename exp2, unsigned int op1, typename T>
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T>,POINT_MUL>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, T d)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression<T>,POINT_MUL> exp_sum(va,point_expression<T>(d));

    return exp_sum;
}


/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename exp3 ,
         typename exp4,
         unsigned int op2,
         typename check = typename std::enable_if<point_expression_op<orig,exp1,exp2,op1>::nvals != 1 && point_expression_op<orig,exp3,exp4,op2>::nvals != 1>::type >
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,POINT_MUL_POINT>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression_op<orig,exp3,exp4,op2> & vb)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,POINT_MUL_POINT> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig,
         typename exp1,
         typename exp2,
         unsigned int op1,
         typename exp3 ,
         typename exp4,
         unsigned int op2,
         typename check = typename std::enable_if<point_expression_op<orig,exp1,exp2,op1>::nvals == 1 || point_expression_op<orig,exp3,exp4,op2>::nvals == 1 >::type >
inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,POINT_MUL>
operator*(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression_op<orig,exp3,exp4,op2> & vb)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}


/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int dim, typename T>
__device__ __host__ inline point_expression_op<Point<dim,T>,Point<dim,T>,Point<dim,T>,POINT_MUL>
pmul(const Point<dim,T> & va, const Point<dim,T> & vb)
{
    point_expression_op<Point<dim,T>,Point<dim,T>,Point<dim,T>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig, unsigned int dim, typename T, typename exp1, typename exp2, unsigned int op1>
__device__ __host__ inline point_expression_op<orig,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL>
pmul(const Point<dim,T> & va, const point_expression_op<orig,exp1,exp2,op1> & vb)
{
    point_expression_op<orig,Point<dim,T>,point_expression_op<orig,exp1,exp2,op1>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig,unsigned int dim, typename T, typename exp1, typename exp2, unsigned int op1>
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,POINT_MUL>
pmul(const point_expression_op<orig,exp1,exp2,op1> & va, const Point<dim,T> & vb)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,Point<dim,T>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two points expression
 *
 * \param va point expression one
 * \param vb point expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename orig, typename exp1, typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
__device__ __host__ inline point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,POINT_MUL>
pmul(const point_expression_op<orig,exp1,exp2,op1> & va, const point_expression_op<orig,exp3,exp4,op2> & vb)
{
    point_expression_op<orig,point_expression_op<orig,exp1,exp2,op1>,point_expression_op<orig,exp3,exp4,op2>,POINT_MUL> exp_sum(va,vb);

    return exp_sum;
}


template<typename T, unsigned int dim>
class point_expression<T[dim]>
{
    T (& d)[dim];

public:

    typedef int has_init;

    typedef int is_expression;

    static const unsigned int nvals = dim;

    typedef T coord_type;

    inline point_expression(T (& d)[dim])
    :d(d)
    {
    }

    template<typename orig, typename exp1, typename exp2, unsigned int op>
    __device__ __host__  point_expression<T[dim]> & operator=(const point_expression_op<orig,exp1,exp2,op> & p_exp)
    {
        p_exp.init();

        for (size_t i = 0; i < dim ; i++)
        {d[i] = p_exp.value(i);}

        return *this;
    }

    template<typename T_>
    __device__ __host__  point_expression<T[dim]> & operator=(const point_expression<T_> & p_exp)
    {
        p_exp.init();

        for (size_t i = 0; i < dim ; i++)
        {d[i] = p_exp.value(i);}

        return *this;
    }

    __device__ __host__ T & operator[](int n)
    {
        return d[n];
    }

    __device__ __host__ const T & operator[](int n) const
    {
        return d[n];
    }

    __device__ __host__ inline void init() const
    {
    }

    __device__ __host__ inline T value(const size_t k) const
    {
        return d[k];
    }
};


template<typename T, unsigned int dim>
class point_expression<const T[dim]>
{
    const T (& d)[dim];

public:

    typedef int has_init;

    typedef int is_expression;

    static const unsigned int nvals = dim;

    typedef T coord_type;

    __device__ __host__ inline point_expression(const T (& d)[dim])
    :d(d)
    {
    }

    __device__ __host__ inline void init() const
    {
    }

    T operator[](int n) const
    {
        return d[n];
    }

    __device__ __host__ inline T value(const size_t k) const
    {
        return d[k];
    }
};

template<typename T, typename vmpl>
class point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>>
{
    mutable openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl> d;

public:

    typedef int has_init;

    typedef int is_expression;

    static const unsigned int nvals = subar_dim<vmpl>::type::value;

    typedef T coord_type;

    template<typename orig, typename exp1, typename exp2, unsigned int op>
    __device__ __host__  point_expression<openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl>> & 
    operator=(const point_expression_op<orig,exp1,exp2,op> & p_exp)
    {
        p_exp.init();

        for (size_t i = 0; i < nvals ; i++)
        {d[i] = p_exp.value(i);}

        return *this;
    }

    __device__ __host__ inline point_expression(const openfpm::detail::multi_array::sub_array_openfpm<T,1,vmpl> & d)
    :d(d.origin_mutable(),d.strides())
    {
    }

    __device__ __host__ inline void init() const
    {
    }

    __device__ __host__ inline T & operator[](int k)
    {
        return value(k);
    }

    __device__ __host__ inline T & operator[](int k) const
    {
        return value(k);
    }

    __device__ __host__ inline T & value(const int k)
    {
        return d[k];
    }

    __device__ __host__ inline T & value(const int k) const
    {
        return d[k];
    }
};

#include "Point_operators_functions.hpp"

#endif /* OPENFPM_DATA_SRC_SPACE_SHAPE_POINT_OPERATORS_HPP_ */