/* -*- c++ -*- vim:et:sw=4:ft=cpp
*/

#ifndef POINT_H
#define POINT_H 1

/*
  Template of a pure class interface definition. It contains no method
  definitions, these are supposed to be provided by specialisation.
*/
template <typename R>
struct Point3I: public R
{
public:

    Point3I();
    Point3I(const float x,const float y,const float z);

    Point3I(const Point3I<R> &that);
    Point3I<R>     	&operator =(const Point3I<R> &that);

    float   	    	xGet() const;
    float       	yGet() const;
    float       	zGet() const;

    float       	xPut(const float x);
    float      		yPut(const float y);
    float       	zPut(const float z);

    Point3I<R>     	&sub(const Point3I<R> &that);

    float       	dot(const Point3I<R> &that) const;
    float       	norm2() const;
};

/*
  Representation class for an unaligned point without SSE acceleration.
*/
struct Point3C
{
protected:
    float x,y,z;

    Point3C
    (
        register const float x,
        register const float y,
        register const float z
    )
        : x(x), y(y), z(z)
    {}
};

/*
  Definition of 'Point3I<Point3C>' methods.
*/

template<>
inline Point3I<Point3C>::Point3I()
    : Point3C(0.0f,0.0f,0.0f)
{
}

template<>
inline Point3I<Point3C>::Point3I
(
    register const float x,register const float y,register const float z
)
    : Point3C(x,y,z)
{
}

template<>
inline Point3I<Point3C>::Point3I(register const Point3I<Point3C> &that)
    : Point3C(that.x,that.y,that.z)
{
}

template<> inline float Point3I<Point3C>::xGet() const { return this->x; }
template<> inline float Point3I<Point3C>::yGet() const { return this->y; }
template<> inline float Point3I<Point3C>::zGet() const { return this->z; }

template<> inline float Point3I<Point3C>::xPut(register const float x)
{
    return this->x = x;
}

template<> inline float Point3I<Point3C>::yPut(register const float y)
{
    return this->y = y;
}

template<> inline float Point3I<Point3C>::zPut(register const float z)
{
    return this->z = z;
}

template<>
inline Point3I<Point3C> &Point3I<Point3C>::operator =
(
    register const Point3I<Point3C> &that
)
{
    this->x = that.x, this->y = that.y, this->z = that.z;

    return *this;
}

template<>
inline Point3I<Point3C> &Point3I<Point3C>::sub
(
    register const Point3I<Point3C> &that
)
{
    this->x -= that.x, this->y -= that.y, this->z -= that.z;

    return *this;
}

template<>
inline float Point3I<Point3C>::dot(register const Point3I<Point3C> &that) const
{
    return this->x*that.x + this->y*that.y + this->z*that.z;
}

template<>
inline float Point3I<Point3C>::norm2() const
{
    return this->dot(*this);
}

#if defined __SSE__
#include <xmmintrin.h>

/*
  Representation class for an aligned point with SSE acceleration.
*/
struct Point3S
{
    union
    {
        __m128        	q;
        float          	v[4];
        struct	        { float w,x,y,z; };
    };

    Point3S(register const __m128 q) : q(q) {}
};

/*
  Definition of 'Point3I<Point3S>' methods, using SSE.
*/

template<>
inline Point3I<Point3S>::Point3I()
    : Point3S(_mm_xor_ps(this->q,this->q))
{
}

template<>
inline Point3I<Point3S>::Point3I
(
    register const float x,register const float y,register const float z
)
    : Point3S(_mm_set_ps(z,y,x,0.0f))
{
}

template<>
inline Point3I<Point3S>::Point3I(register const Point3I<Point3S> &that)
    : Point3S(that.q)
{
}

template<> inline float Point3I<Point3S>::xGet() const { return this->x; }
template<> inline float Point3I<Point3S>::yGet() const { return this->y; }
template<> inline float Point3I<Point3S>::zGet() const { return this->z; }

template<>
inline float Point3I<Point3S>::xPut(register const float x)
{
    return this->x = x;
}

template<>
inline float Point3I<Point3S>::yPut(register const float y)
{
    return this->y = y;
}

template<>
inline float Point3I<Point3S>::zPut(register const float z)
{
    return this->z = z;
}

template<>
inline Point3I<Point3S> &Point3I<Point3S>::operator =
(
    register const Point3I<Point3S> &that
)
{
    this->q = that.q;

    return *this;
}

template<>
inline Point3I<Point3S> &Point3I<Point3S>::sub
(
    register const Point3I<Point3S> &that
)
{
    this->q = _mm_sub_ps(this->q,that.q);

    return *this;
}

/**
 * Add together the first 3 elements of the parameter,
 * and return the value as element 0 of the result.
 *
 * Horizontal adds are part of SSE2 but not SSE.
 */
static inline __m128 _mm_local_hadd_ps(register const __m128 q)
{
    return
        _mm_add_ss(
            _mm_shuffle_ps(q,q,_MM_SHUFFLE(1,1,1,1)),
            _mm_add_ss(
                _mm_shuffle_ps(q,q,_MM_SHUFFLE(2,2,2,2)),
                _mm_shuffle_ps(q,q,_MM_SHUFFLE(3,3,3,3))));
}

template<>
inline float Point3I<Point3S>::dot(register const Point3I<Point3S> &that) const
{
    register float d;

    _mm_store_ss(&d,_mm_local_hadd_ps(_mm_mul_ps(this->q,that.q)));

    return d;
}

template<>
inline float Point3I<Point3S>::norm2() const
{
    return this->dot(*this);
}
#endif

#endif