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_

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

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;

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

    inline void init() const
    {
    }

    inline T value(const size_t k) const
    {
        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 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;

    inline point_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(o1.value(0))>::type > 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  > inline 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>
{
    const exp1 o1;
    const exp2 o2;

public:

    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;

    inline point_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(o1.value(0))>::type > 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  > 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;

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

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

    template<typename r_type=typename std::remove_reference<decltype(o1.value(0))>::type > 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  >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 orig orig_type;

    typedef int has_init;

    typedef int is_expression;

    typedef typename orig::coord_type return_type;

    static const unsigned int nvals = 1;

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

    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 std::remove_reference<decltype(o1.value(0))>::type > 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 > 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 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;

    inline point_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(o1.value(0))>::type > 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  > 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 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;

    inline point_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(o1.value(0))>::type > 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  > operator T() const
    {
        init();
        return o1.value(0) / o2.value(0);
    }
};

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

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

#define CREATE_POINT_OPERATOR(operator_name,OP_ID) \
\
\
template<unsigned int dim, typename T>\
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>\
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>\
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<unsigned int dim, typename T>\
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<unsigned int dim, typename T>\
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<unsigned int dim, typename T>\
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<unsigned int dim, typename T>\
inline point_expression_op<Point<dim,T>,point_expression<T[dim]>,point_expression<const T[dim]>,OP_ID>\
operator_name(const point_expression<T[dim]> & va, const point_expression<const T[dim]> & vb)\
{\
    point_expression_op<Point<dim,T>,point_expression<T[dim]>,point_expression<const 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>\
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, unsigned int dim, typename T>\
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<unsigned int dim, typename T>\
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>\
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>\
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>\
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>\
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>\
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 >\
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>\
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 >\
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>\
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>\
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 >\
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>\
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 >
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>
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>
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>
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>
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<unsigned int dim, typename T>
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<unsigned int dim, typename T>
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<unsigned int dim, typename T, typename check=typename std::enable_if< !std::is_same<T,double>::value >::type>
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>
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, unsigned int dim, typename T>
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;
}

/* \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>
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 >
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>
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>
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 >
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 >
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 >
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 >
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>
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>
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 >
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>
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>
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>
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>
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;

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

    template<typename orig, typename exp1, typename exp2, unsigned int op> 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;
    }

    inline void init() const
    {
    }

    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;

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

    inline void init() const
    {
    }

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

#include "Point_operators_functions.hpp"

#endif /* OPENFPM_DATA_SRC_SPACE_SHAPE_POINT_OPERATORS_HPP_ */