Vector.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

#ifndef __math_Vector_h__
#define __math_Vector_h__

namespace fe
{

/**************************************************************************//**
    @brief Dense vector - size fixed by template

    @ingroup math
*//***************************************************************************/
template<int N, class T>
class Vector
{
    public:
                    Vector(void);
                    Vector(const Vector<N,T> &other);
                    template<typename U>
                    Vector(const Vector<N,U> &other);
                    template<typename U>
                    Vector(const Vector<N+1,U> &other);
                    template<typename U>
                    Vector(const Vector<N+2,U> &other);
                    template<typename U>
                    Vector(const Vector<N-1,U> &other);
                    template<typename U>
                    Vector(const Vector<N-2,U> &other);
                    Vector(const T *array);

//                  Vector(T x,T y,T z=(T)0,T w=(T)0,T v=(T)0,T u=(T)0);
                    Vector(T x,T y,T z,T w,T v,T u);
                    Vector(T x,T y,T z,T w,T v);
                    Vector(T x,T y,T z,T w);
                    Vector(T x,T y,T z);
                    Vector(T x,T y);
                    Vector(T x);

//                  Vector(const T &all);

const   T&          operator[](unsigned int index) const;
        T&          operator[](unsigned int index);

        T           *raw(void);
const   T           *raw(void) const;

                    // like raw, but only valid if used on the call
const   T           *temp(void) const;

        Vector<N,T> &operator=(const T *array);
        Vector<N,T> &operator=(const Vector<N,T> &other);
        template<typename U>
        Vector<N,T> &operator=(const Vector<N,U> &other);
        Vector<N,T> &operator=(const Vector<N+1,T> &other);
        template<typename U>
        Vector<N,T> &operator=(const Vector<N+1,U> &other);
        template<typename U>
        Vector<N,T> &operator=(const Vector<N+2,U> &other);
        template<typename U>
        Vector<N,T> &operator=(const Vector<N-1,U> &other);
        template<typename U>
        Vector<N,T> &operator=(const Vector<N-2,U> &other);

        bool        operator==(const Vector<N,T> &other) const;
        bool        operator!=(const Vector<N,T> &other) const;

    public:
        T   m_data[N];
};

/// @internal
template<typename T>
inline
BWORD checkValid(const T& a_value)
{
    return TRUE;
}

/// @internal
template<typename T>
inline
void confirmValid(const T& a_value)
{
}

/// @internal
inline
BWORD checkValid(const F32& a_value)
{
    return !FE_INVALID_SCALAR(a_value);
}

/// @internal
inline
void confirmValid(const F32& a_value)
{
    if(!checkValid(a_value))
        feX(e_corrupt,"confirmValid(F32)","scalar invalid");
}

/// @internal
inline
BWORD checkValid(const F64& a_value)
{
    return !FE_INVALID_SCALAR(a_value);
}

/// @internal
inline
void confirmValid(const F64& a_value)
{
    if(!checkValid(a_value))
        feX(e_corrupt,"confirmValid(F64)","scalar invalid");
}

/// @internal
template<int N,typename T>
inline
BWORD checkValid(const Vector<N,T>& a_vector)
{
    return TRUE;
}

/// @internal
template<int N,typename T>
inline
void confirmValid(const Vector<N,T>& a_vector)
{
}

/// @internal
template<int N>
inline
BWORD checkValid(const Vector<N,F32>& a_vector)
{
    const F32* data=a_vector.raw();
    if(N>0 && FE_INVALID_SCALAR(data[0]))
        return FALSE;
    if(N>1 && FE_INVALID_SCALAR(data[1]))
        return FALSE;
    if(N>2 && FE_INVALID_SCALAR(data[2]))
        return FALSE;
    if(N>3 && FE_INVALID_SCALAR(data[3]))
        return FALSE;
    if(N>4 && FE_INVALID_SCALAR(data[4]))
        return FALSE;
    if(N>5 && FE_INVALID_SCALAR(data[5]))
        return FALSE;
    return TRUE;
}

/// @internal
template<int N>
inline
void confirmValid(const Vector<N,F32>& a_vector)
{
    const F32* data=a_vector.raw();
    if(N>0 && FE_INVALID_SCALAR(data[0]))
        feX(e_corrupt,"confirmValid(Vector<N,F32>)","vector element 0 invalid");
    if(N>1 && FE_INVALID_SCALAR(data[1]))
        feX(e_corrupt,"confirmValid(Vector<N,F32>)","vector element 1 invalid");
    if(N>2 && FE_INVALID_SCALAR(data[2]))
        feX(e_corrupt,"confirmValid(Vector<N,F32>)","vector element 2 invalid");
    if(N>3 && FE_INVALID_SCALAR(data[3]))
        feX(e_corrupt,"confirmValid(Vector<N,F32>)","vector element 3 invalid");
    if(N>4 && FE_INVALID_SCALAR(data[4]))
        feX(e_corrupt,"confirmValid(Vector<N,F32>)","vector element 4 invalid");
    if(N>5 && FE_INVALID_SCALAR(data[5]))
        feX(e_corrupt,"confirmValid(Vector<N,F32>)","vector element 5 invalid");
}

/// @internal
template<int N>
inline
BWORD checkValid(const Vector<N,F64>& a_vector)
{
    const F64* data=a_vector.raw();
    if(N>0 && FE_INVALID_SCALAR(data[0]))
        return FALSE;
    if(N>1 && FE_INVALID_SCALAR(data[1]))
        return FALSE;
    if(N>2 && FE_INVALID_SCALAR(data[2]))
        return FALSE;
    if(N>3 && FE_INVALID_SCALAR(data[3]))
        return FALSE;
    if(N>4 && FE_INVALID_SCALAR(data[4]))
        return FALSE;
    if(N>5 && FE_INVALID_SCALAR(data[5]))
        return FALSE;
    return TRUE;
}

/// @internal
template<int N>
inline
void confirmValid(const Vector<N,F64>& a_vector)
{
    const F64* data=a_vector.raw();
    if(N>0 && FE_INVALID_SCALAR(data[0]))
        feX(e_corrupt,"confirmValid(Vector<N,F64>)","vector element 0 invalid");
    if(N>1 && FE_INVALID_SCALAR(data[1]))
        feX(e_corrupt,"confirmValid(Vector<N,F64>)","vector element 1 invalid");
    if(N>2 && FE_INVALID_SCALAR(data[2]))
        feX(e_corrupt,"confirmValid(Vector<N,F64>)","vector element 2 invalid");
    if(N>3 && FE_INVALID_SCALAR(data[3]))
        feX(e_corrupt,"confirmValid(Vector<N,F64>)","vector element 3 invalid");
    if(N>4 && FE_INVALID_SCALAR(data[4]))
        feX(e_corrupt,"confirmValid(Vector<N,F64>)","vector element 4 invalid");
    if(N>5 && FE_INVALID_SCALAR(data[5]))
        feX(e_corrupt,"confirmValid(Vector<N,F64>)","vector element 5 invalid");
}

/// @internal
template<typename T>
inline
void maybeConfirmValid(const T& a_value)
{
#if FE_VEC_CHECK_VALID
    confirmValid(a_value);
#endif
}

template<int N, class T>
inline Vector<N,T>::Vector(void)
{
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = T(0);
    }
    maybeConfirmValid(*this);
}

template<int N, class T>
inline Vector<N,T>::Vector(const Vector<N,T> &other)
{
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = other[i];
    }
    maybeConfirmValid(*this);
}

template<int N, class T>
template<typename U>
inline Vector<N,T>::Vector(const Vector<N,U> &other)
{
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = (T)other[i];
    }
    maybeConfirmValid(*this);
}

template<int N, class T>
template<typename U>
inline Vector<N,T>::Vector(const Vector<N+1,U> &other)
{
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = other[i];
    }
    maybeConfirmValid(*this);
}

template<int N, class T>
template<typename U>
inline Vector<N,T>::Vector(const Vector<N+2,U> &other)
{
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = other[i];
    }
    maybeConfirmValid(*this);
}

template<int N, class T>
template<typename U>
inline Vector<N,T>::Vector(const Vector<N-1,U> &other)
{
    for(unsigned int i = 0; i < N-1; i++)
    {
        m_data[i] = other[i];
    }
    if(N>1)
    {
        m_data[N-1]=T(0);
    }
    maybeConfirmValid(*this);
}

template<int N, class T>
template<typename U>
inline Vector<N,T>::Vector(const Vector<N-2,U> &other)
{
    for(unsigned int i = 0; i < N-2; i++)
    {
        m_data[i] = other[i];
    }
    if(N>1)
    {
        m_data[N-1]=T(0);
        if(N>2)
        {
            m_data[N-2]=T(0);
        }
    }
    maybeConfirmValid(*this);
}

template<int N, class T>
inline Vector<N,T>::Vector(const T *array)
{
    operator=(array);
}

#if FALSE
template<int N,class T>
inline Vector<N,T>::Vector(T x, T y, T z, T w, T v, T u)
{
    m_data[0] = x;
    if(N >=2) m_data[1] = y;
    if(N >=3) m_data[2] = z;
    if(N >=4) m_data[3] = w;
    if(N >=5) m_data[4] = v;
    if(N >=6) m_data[5] = u;

    maybeConfirmValid(*this);
}
#else
template<int N,class T>
inline Vector<N,T>::Vector(T x, T y, T z, T w, T v, T u)
{
    m_data[0] = x;
    m_data[1] = y;
    m_data[2] = z;
    m_data[3] = w;
    m_data[4] = v;
    m_data[5] = u;
    for(I32 n=6;n<N;n++) m_data[n]=T(0);

    maybeConfirmValid(*this);
}
template<int N,class T>
inline Vector<N,T>::Vector(T x, T y, T z, T w, T v)
{
    m_data[0] = x;
    m_data[1] = y;
    m_data[2] = z;
    m_data[3] = w;
    m_data[4] = v;
    for(I32 n=5;n<N;n++) m_data[n]=T(0);

    maybeConfirmValid(*this);
}
template<int N,class T>
inline Vector<N,T>::Vector(T x, T y, T z, T w)
{
    m_data[0] = x;
    m_data[1] = y;
    m_data[2] = z;
    m_data[3] = w;
    for(I32 n=4;n<N;n++) m_data[n]=T(0);

    maybeConfirmValid(*this);
}
template<int N,class T>
inline Vector<N,T>::Vector(T x, T y, T z)
{
    m_data[0] = x;
    m_data[1] = y;
    m_data[2] = z;
    for(I32 n=3;n<N;n++) m_data[n]=T(0);

    maybeConfirmValid(*this);
}
template<int N,class T>
inline Vector<N,T>::Vector(T x, T y)
{
    m_data[0] = x;
    m_data[1] = y;
    for(I32 n=2;n<N;n++) m_data[n]=T(0);

    maybeConfirmValid(*this);
}
template<int N,class T>
inline Vector<N,T>::Vector(T x)
{
    m_data[0] = x;
    for(I32 n=1;n<N;n++) m_data[n]=T(0);

    maybeConfirmValid(*this);
}
#endif

#if FALSE
template<int N,class T>
inline Vector<N,T>::Vector(const T &all)
{
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = all;
    }
}
#endif

template<int N, class T>
inline const T& Vector<N,T>::operator[](unsigned int index) const
{
#if FE_BOUNDSCHECK
    if(index >= N) { feX(e_invalidRange); }
#endif
    return m_data[index];
}

template<int N, class T>
inline T &Vector<N,T>::operator[](unsigned int index)
{
#if FE_BOUNDSCHECK
    if(index >= N) { feX(e_invalidRange); }
#endif
    return m_data[index];
}

template<int N, class T>
inline T *Vector<N,T>::raw(void)
{
    return m_data;
}

template<int N, class T>
inline const T *Vector<N,T>::raw(void) const
{
    return m_data;
}

template<int N, class T>
inline const T *Vector<N,T>::temp(void) const
{
    maybeConfirmValid(*this);
    return m_data;
}

template<int N, class T>
inline bool Vector<N,T>::operator==(const Vector<N,T> &other) const
{
    maybeConfirmValid(*this);
    maybeConfirmValid(other);
    if((void*)this != (void*)&other)
    {
        for(unsigned int i = 0; i < N; i++)
        {
            if(m_data[i] != other[i]) { return false; }
        }
    }
    return true;
}

template<int N, class T>
inline bool Vector<N,T>::operator!=(const Vector<N,T> &other) const
{
    return !operator==(other);
}

template<int N, class T>
inline Vector<N,T> &Vector<N,T>::operator=(const Vector<N,T> &other)
{
    maybeConfirmValid(other);
    if((void*)this != (void*)&other)
    {
        for(unsigned int i = 0; i < N; i++)
        {
            m_data[i] = other[i];
        }
    }
    maybeConfirmValid(*this);
    return *this;
}

template<int N, class T>
template<typename U>
inline Vector<N,T> &Vector<N,T>::operator=(const Vector<N,U> &other)
{
    maybeConfirmValid(other);
    if((void*)this != (void*)&other)
    {
        for(unsigned int i = 0; i < N; i++)
        {
            m_data[i] = (T)other[i];
        }
    }
    maybeConfirmValid(*this);
    return *this;
}

template<int N, class T>
inline Vector<N,T> &Vector<N,T>::operator=(const Vector<N+1,T> &other)
{
    maybeConfirmValid(other);
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = other[i];
    }
    maybeConfirmValid(*this);
    return *this;
}

template<int N, class T>
template<typename U>
inline Vector<N,T> &Vector<N,T>::operator=(const Vector<N+1,U> &other)
{
    maybeConfirmValid(other);
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = (T)other[i];
    }
    maybeConfirmValid(*this);
    return *this;
}

template<int N, class T>
template<typename U>
inline Vector<N,T> &Vector<N,T>::operator=(const Vector<N+2,U> &other)
{
    maybeConfirmValid(other);
    for(unsigned int i = 0; i < N; i++)
    {
        m_data[i] = (T)other[i];
    }
    maybeConfirmValid(*this);
    return *this;
}

template<int N, class T>
template<typename U>
inline Vector<N,T> &Vector<N,T>::operator=(const Vector<N-1,U> &other)
{
    maybeConfirmValid(other);
    for(unsigned int i = 0; i < N-1; i++)
    {
        m_data[i] = (T)other[i];
    }
    if(N>1)
    {
        m_data[N-1]=T(0);
    }
    maybeConfirmValid(*this);
    return *this;
}

template<int N, class T>
template<typename U>
inline Vector<N,T>& Vector<N,T>::operator=(const Vector<N-2,U> &other)
{
    maybeConfirmValid(other);
    for(unsigned int i = 0; i < N-2; i++)
    {
        m_data[i] = (T)other[i];
    }
    if(N>1)
    {
        m_data[N-1]=T(0);
        if(N>2)
        {
            m_data[N-2]=T(0);
        }
    }
    maybeConfirmValid(*this);
    return *this;
}

template<int N, class T>
inline Vector<N,T> &Vector<N,T>::operator=(const T *array)
{
    for(unsigned int i = 0; i < N; i++)
    {
        maybeConfirmValid(array[i]);
        m_data[i] = array[i];
    }
    maybeConfirmValid(*this);
    return *this;
}

/* non member operations */

/** Set all the elements to zero
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T>& set(Vector<N,T> &lhs)
{
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i] = T(0);
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** Set all the elements to the given value
    @relates Vector
    */
template<int N, class T, class U>
inline Vector<N,T>& setAll(Vector<N,T> &lhs,const U value)
{
    maybeConfirmValid(value);
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i] = value;
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** Set the value at the index
    @relates Vector
    */
template<int N, class T,class U>
inline Vector<N,T> &setAt(Vector<N,T> &lhs, U32 index,const U value)
{
    maybeConfirmValid(value);
    lhs[index]=value;
    maybeConfirmValid(lhs);
    return lhs;
}

/** Return the number of elements
    @relates Vector
    */
template<int N, class T>
inline U32 size(const Vector<N,T> &lhs)
{
    maybeConfirmValid(lhs);
    return N;
}

/** In place add operator
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> &operator+=(Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(rhs);
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i] += rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** In place add operator
    @relates Vector
    */
template<int M, int N, class T>
inline Vector<N,T> &operator+=(Vector<N,T> &lhs, const Vector<M,T> &rhs)
{
    maybeConfirmValid(rhs);
    const unsigned int limit=(M<N)? M: N;
    for(unsigned int i = 0; i < limit; i++)
    {
        lhs[i] += rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** In place subtract operator
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> &operator-=(Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(rhs);
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i] -= rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** In place subtract operator
    @relates Vector
    */
template<int M, int N, class T>
inline Vector<N,T> &operator-=(Vector<N,T> &lhs, const Vector<M,T> &rhs)
{
    maybeConfirmValid(rhs);
    const int limit=(M<N)? M: N;
    for(unsigned int i = 0; i < limit; i++)
    {
        lhs[i] -= rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** Negate operation
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> operator-(const Vector<N,T> &rhs)
{
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        v[i] = -rhs[i];
    }
    return v;
}

/** In place piecewise multiply operator
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> &operator*=(Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(rhs);
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i] = lhs[i] * rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** In place piecewise multiply operator
    @relates Vector
    */
template<int M, int N, class T>
inline Vector<N,T> &operator*=(Vector<N,T> &lhs, const Vector<M,T> &rhs)
{
    maybeConfirmValid(rhs);
    const int limit=(M<N)? M: N;
    for(unsigned int i = 0; i < limit; i++)
    {
        lhs[i] *= rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** In place piecewise scale operator
    @relates Vector
    */
template<int N, class T,class U>
inline Vector<N,T> &operator*=(Vector<N,T> &lhs, U scale)
{
    maybeConfirmValid(scale);
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i] *= scale;
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** Dot (inner) product
    @relates Vector
    */
template<int N, class T>
inline T dot(const Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    T t = (T)0;
    for(unsigned int i = 0; i < N; i++)
    {
        t += lhs[i] * rhs[i];
    }
    maybeConfirmValid(t);
    return t;
}

/** Frobenius norm operation
    @relates Vector
    */
template<int N, class T>
inline T magnitude(const Vector<N,T> &rhs)
{
    return sqrt(magnitudeSquared(rhs));
}

/** Square of the length
    @relates Vector
    */
template<int N, class T>
inline T magnitudeSquared(const Vector<N,T> &rhs)
{
    maybeConfirmValid(rhs);
    T mag = 0.0;
    for(unsigned int i = 0; i < N; i++)
    {
        mag += (T)rhs[i] * (T)rhs[i];
    }
    maybeConfirmValid(mag);
    return mag;
}

/** Return vector scaled to unit length
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> unit(const Vector<N,T> &vec)
{
    maybeConfirmValid(vec);
    T mag=magnitude(vec);
    if(mag==T(0))
    {
        feX("unit(Vector<N,T>)","attempt to normalize zero magnitude vector");
    }
    T inv=T(1)/mag;
    Vector<N,T> v;
    for(U32 i=0;i<N;i++)
    {
        v[i]=inv*vec[i];
    }
    maybeConfirmValid(v);
    return v;
}

/** Return vector scaled to unit length with zero check
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> unitSafe(const Vector<N,T> &vec)
{
    maybeConfirmValid(vec);
    T mag=magnitude(vec);
    if(mag>T(0))
    {
        T inv=T(1)/mag;
        Vector<N,T> v;
        for(U32 i=0;i<N;i++)
        {
            v[i]=inv*vec[i];
        }
        maybeConfirmValid(v);
        return v;
    }
    return vec;
}

/** In place normalize operator
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> &normalize(Vector<N,T> &vec)
{
    maybeConfirmValid(vec);
    T mag=magnitude(vec);
    if(mag==T(0))
    {
        feX("Vector<N,T>::normal",
            "attempt to normalize zero magnitude vector");
    }
    T inv=T(1)/mag;
    for(U32 i=0;i<N;i++)
    {
        vec[i]*=inv;
    }
    maybeConfirmValid(vec);
    return vec;
}

/** In place normalize operator with zero length check
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> &normalizeSafe(Vector<N,T> &vec)
{
    maybeConfirmValid(vec);
    T mag=magnitude(vec);
    if(mag>T(0))
    {
        T inv=T(1)/mag;
        for(U32 i=0;i<N;i++)
        {
            vec[i]*=inv;
        }
        maybeConfirmValid(vec);
    }
    return vec;
}

/** Print to a string
    @relates Vector
    */
template<int N, class T>
inline String print(const Vector<N,T> &vec)
{
    maybeConfirmValid(vec);
    String s = "";
    for(unsigned int i = 0; i < N; i++)
    {
        s.cat(print(vec[i]));
        if(i < N-1)
        {
            s.cat(" ");
        }
    }
    return s;
}

/** add operation
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> operator+(const Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        v[i] = lhs[i] + rhs[i];
    }
    maybeConfirmValid(v);
    return v;
}

/** add operation
    @relates Vector
    */
template<int N, int M, class T>
inline Vector<N,T> operator+(const Vector<N,T> &lhs, const Vector<M,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        if(i >= M)
        {
            v[i] = lhs[i];
        }
        else
        {
            v[i] = lhs[i] + rhs[i];
        }
    }
    maybeConfirmValid(v);
    return v;
}


/** subtract operation
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> operator-(const Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        v[i] = lhs[i] - rhs[i];
    }
    maybeConfirmValid(v);
    return v;
}

/** subtract operation
    @relates Vector
    */
template<int N, int M, class T>
inline Vector<N,T> operator-(const Vector<N,T> &lhs, const Vector<M,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        if(i >= M)
        {
            v[i] = lhs[i];
        }
        else
        {
            v[i] = lhs[i] - rhs[i];
        }
    }
    maybeConfirmValid(v);
    return v;
}

/** equality test
    @relates Vector
    */
template<int N, class T>
inline bool operator==(const Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    for(unsigned int i = 0; i < N; i++)
    {
        if(rhs[i] != lhs[i])
        {
            return false;
        }
    }
    return true;
}

/** Equivalence test within the given tolerance @em margin
    @relates Vector
    */
template<int N, class T,class U>
inline bool equivalent(const Vector<N,T> &lhs, const Vector<N,T> &rhs, U margin)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    maybeConfirmValid(margin);
    for(unsigned int i = 0; i < N; i++)
    {
        if(FE_INVALID_SCALAR(lhs[i]) ||
                FE_INVALID_SCALAR(rhs[i]) ||
                fabs(rhs[i] - lhs[i]) > margin)
        {
            return false;
        }
    }
    return true;
}

/** Piecewise multiply operation
    @relates Vector
    */
template<int N, class T>
inline Vector<N,T> operator*(const Vector<N,T> &lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        v[i] = lhs[i] * rhs[i];
    }
    maybeConfirmValid(v);
    return v;
}

/** Scale operation
    @relates Vector
    */
template<int N, class T,class U>
inline Vector<N,T> operator*(const U lhs, const Vector<N,T> &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        v[i] = rhs[i] * lhs;
    }
    maybeConfirmValid(v);
    return v;
}

/** @brief Scale operation

    @relates Vector
    */
template<int N, class T,class U>
inline Vector<N,T> operator*(const Vector<N,T> &lhs, const U &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    Vector<N,T> v;
    for(unsigned int i = 0; i < N; i++)
    {
        v[i] = lhs[i] * rhs;
    }
    maybeConfirmValid(v);
    return v;
}

/** @brief Inverse Scale operation

    @relates Vector
    */
template<int N, class T,class U>
inline Vector<N,T> operator/(const Vector<N,T> &lhs, const U &rhs)
{
    maybeConfirmValid(lhs);
    maybeConfirmValid(rhs);
    if(rhs==0.0f)
        feX(e_unsolvable,"operator/(Vector<N,T>,T)","divide by zero");
    return lhs*(1.0f/rhs);
}

/** @brief Add with scaling

    lhs = lhs + rhs * scalar

    @relates Vector
    */
template<int N, class T,class U>
inline Vector<N,T> &addScaled(Vector<N,T> &lhs,U scalar,
        const Vector<N,T> &rhs)
{
    maybeConfirmValid(rhs);
    maybeConfirmValid(scalar);
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i]+=scalar*rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

/** @brief Scale then add

    lhs = lhs * scalar + rhs

    @relates Vector
    */
template<int N, class T,class U>
inline Vector<N,T> &scaleAndAdd(Vector<N,T> &lhs,U scalar,
        const Vector<N,T> &rhs)
{
    maybeConfirmValid(rhs);
    maybeConfirmValid(scalar);
    for(unsigned int i = 0; i < N; i++)
    {
        lhs[i]=lhs[i]*scalar+rhs[i];
    }
    maybeConfirmValid(lhs);
    return lhs;
}

} /* namespace */

#endif /* __math_Vector_h__ */