vector.h

Go to the documentation of this file.

/***************************************************************************************************
 * Copyright 2022 NVIDIA Corporation. All rights reserved.
 **************************************************************************************************/

#ifndef MI_MATH_VECTOR_H
#define MI_MATH_VECTOR_H

#include <mi/base/types.h>
#include <mi/math/assert.h>
#include <mi/math/function.h>

namespace mi {

namespace math {

enum From_iterator_tag {
    FROM_ITERATOR 
};


// Color and Vector can be converted into each other. To avoid cyclic dependencies among the
// headers, the Color_struct class is already defined here.

//------- POD struct that provides storage for color elements --------

struct Color_struct
{
    Float32 r;
    Float32 g;
    Float32 b;
    Float32 a;
};


//------- POD struct that provides storage for vector elements --------

template <typename T, Size DIM>
struct Vector_struct
{
    T elements[DIM]; 
};

template <typename T> struct Vector_struct<T, 1>
{
    T x;  
};

template <typename T> struct Vector_struct<T, 2>
{
    T x;  
    T y;  
};

template <typename T> struct Vector_struct<T, 3>
{
    T x;  
    T y;  
    T z;  
};

template <typename T> struct Vector_struct<T, 4>
{
    T x;  
    T y;  
    T z;  
    T w;  
};


//------ Indirect access to vector storage base ptr to keep Vector_struct a POD --

template <typename T, Size DIM>
inline T* vector_base_ptr( Vector_struct<T,DIM>& vec)
{
    return vec.elements;
}

template <typename T, Size DIM>
inline const T* vector_base_ptr( const Vector_struct<T,DIM>& vec)
{
    return vec.elements;
}

template <typename T>
inline T* vector_base_ptr( Vector_struct<T,1>& vec)
{
    return &vec.x;
}

template <typename T>
inline const T* vector_base_ptr( const Vector_struct<T,1>& vec)
{
    return &vec.x;
}

template <typename T>
inline T* vector_base_ptr( Vector_struct<T,2>& vec)
{
    return &vec.x;
}

template <typename T>
inline const T* vector_base_ptr( const Vector_struct<T,2>& vec)
{
    return &vec.x;
}

template <typename T>
inline T* vector_base_ptr( Vector_struct<T,3>& vec)
{
    return &vec.x;
}

template <typename T>
inline const T* vector_base_ptr( const Vector_struct<T,3>& vec)
{
    return &vec.x;
}

template <typename T>
inline T* vector_base_ptr( Vector_struct<T,4>& vec)
{
    return &vec.x;
}

template <typename T>
inline const T* vector_base_ptr( const Vector_struct<T,4>& vec)
{
    return &vec.x;
}

 // end group mi_math_vector_struct


//------ Generic Vector Class -------------------------------------------------

template < class T, Size DIM>
class Vector : public Vector_struct<T, DIM> //-V690 PVS
{
public:
    typedef Vector_struct<T,DIM> Pod_type;      
    typedef Vector_struct<T,DIM> storage_type;  

    typedef T           value_type;             
    typedef Size        size_type;              
    typedef Difference  difference_type;        
    typedef T *         pointer;                
    typedef const T *   const_pointer;          
    typedef T &         reference;              
    typedef const T &   const_reference;        

    static const Size DIMENSION = DIM;          
    static const Size SIZE      = DIM;          

    static inline Size size()     { return SIZE; }

    static inline Size max_size() { return SIZE; }

    inline T* begin() { return mi::math::vector_base_ptr( *this); }

    inline const T* begin() const { return mi::math::vector_base_ptr( *this); }

    inline T* end() { return begin() + DIM; }

    inline const T* end() const { return begin() + DIM; }

    inline Vector()
    {
#if defined(DEBUG) || (defined(_MSC_VER) && _MSC_VER <= 1310)
        // In debug mode, default-constructed vectors are initialized with signaling NaNs or, if not
        // applicable, with a maximum value to increase the chances of diagnosing incorrect use of
        // an uninitialized vector.
        //
        // When compiling with Visual C++ 7.1 or earlier, this code is enabled in all variants to
        // work around a very obscure compiler bug that causes the compiler to crash.
        typedef mi::base::numeric_traits<T> Traits;
        T v = (Traits::has_signaling_NaN) ? Traits::signaling_NaN()
                                          : Traits::max MI_PREVENT_MACRO_EXPAND ();
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = v;
#endif
    }

#if (__cplusplus >= 201103L)
    Vector( const Vector<T,DIM>& vec ) = default;
#endif

    inline Vector( const Vector_struct<T,DIM>& vec )
    {
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = mi::math::vector_base_ptr(vec)[i];
    }

    inline explicit Vector(T v)
    {
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = v;
    }

    template <typename Iterator>
    inline Vector(From_iterator_tag, Iterator p)
    {
        for( Size i(0u); i < DIM; ++i, ++p)
            (*this)[i] = *p;
    }

    template <typename T2>
    inline explicit Vector( T2 const (& array)[DIM])
    {
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = array[i];
    }

    template <typename T2>
    inline explicit Vector( const Vector<T2,DIM>& other)
    {
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = T(other[i]);
    }

    template <typename T2>
    inline explicit Vector( const Vector_struct<T2,DIM>& other)
    {
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = T(mi::math::vector_base_ptr(other)[i]);
    }

    inline Vector(T v1, T v2)
    {
        mi_static_assert(DIM == 2);
        begin()[0] = v1;
        begin()[1] = v2;
    }

    inline Vector(T v1, T v2, T v3)
    {
        mi_static_assert(DIM == 3);
        begin()[0] = v1;
        begin()[1] = v2;
        begin()[2] = v3;
    }

    inline Vector(T v1, const Vector<T,2>& v2)
    {
        mi_static_assert(DIM == 3);
        begin()[0] = v1;
        begin()[1] = v2.x;
        begin()[2] = v2.y;
    }

    inline Vector(const Vector<T,2>& v1, T v2)
    {
        mi_static_assert(DIM == 3);
        begin()[0] = v1.x;
        begin()[1] = v1.y;
        begin()[2] = v2;
    }

    inline Vector(T v1, T v2, T v3, T v4)
    {
        mi_static_assert(DIM == 4);
        begin()[0] = v1;
        begin()[1] = v2;
        begin()[2] = v3;
        begin()[3] = v4;
    }

    inline Vector(T v1, T v2, const Vector<T,2>& v3)
    {
        mi_static_assert(DIM == 4);
        begin()[0] = v1;
        begin()[1] = v2;
        begin()[2] = v3.x;
        begin()[3] = v3.y;
    }


    inline Vector(T v1, const Vector<T,2>& v2, T v3)
    {
        mi_static_assert(DIM == 4);
        begin()[0] = v1;
        begin()[1] = v2.x;
        begin()[2] = v2.y;
        begin()[3] = v3;
    }

    inline Vector(const Vector<T,2>& v1, T v2, T v3)
    {
        mi_static_assert(DIM == 4);
        begin()[0] = v1.x;
        begin()[1] = v1.y;
        begin()[2] = v2;
        begin()[3] = v3;
    }

    inline Vector(const Vector<T,2>& v1, const Vector<T,2>& v2)
    {
        mi_static_assert(DIM == 4);
        begin()[0] = v1.x;
        begin()[1] = v1.y;
        begin()[2] = v2.x;
        begin()[3] = v2.y;
    }

    inline Vector(T v1, const Vector<T,3>& v2)
    {
        mi_static_assert(DIM == 4);
        begin()[0] = v1;
        begin()[1] = v2.x;
        begin()[2] = v2.y;
        begin()[3] = v2.z;
    }

    inline Vector(const Vector<T,3>& v1, T v2)
    {
        mi_static_assert(DIM == 4);
        begin()[0] = v1.x;
        begin()[1] = v1.y;
        begin()[2] = v1.z;
        begin()[3] = v2;
    }

    inline explicit Vector( const Color_struct& color)
    {
        mi_static_assert(DIM == 4);
        this->x = color.r;
        this->y = color.g;
        this->z = color.b;
        this->w = color.a;
    }

    inline Vector& operator= ( const Vector& other)
    {
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = other[i];
        return *this;
    }

    inline Vector& operator= ( T s)
    {
        for( Size i(0u); i < DIM; ++i)
            (*this)[i] = s;
        return *this;
    }
    inline Vector& operator= ( const Color_struct& color)
    {
        mi_static_assert(DIM == 4);
        this->x = color.r;
        this->y = color.g;
        this->z = color.b;
        this->w = color.a;
        return *this;
    }

    inline T& operator[] (Size i)
    {
        mi_math_assert_msg(i < DIM, "precondition");
        return begin()[i];
    }

    inline const T& operator[] (Size i) const
    {
        mi_math_assert_msg(i < DIM, "precondition");
        return begin()[i];
    }

    inline const T& get(Size i) const
    {
        mi_math_assert_msg(i < DIM, "precondition");
        return begin()[i];
    }

    inline void set(Size i, T value)
    {
        mi_math_assert_msg(i < DIM, "precondition");
        begin()[i] = value;
    }



    inline bool normalize()
    {
        const T rec_length = T(1) / length( *this);
        const bool result = isfinite( rec_length); //-V601 PVS
        if( result)
            (*this) *= rec_length;
        return result;
    }


    //------ Free comparison operators ==, !=, <, <=, >, >= for vectors --------

    inline bool operator==( Vector<T,DIM> rhs) const
    {
        return is_equal( *this, rhs);
    }

    inline bool operator!=( Vector<T,DIM> rhs) const
    {
        return is_not_equal( *this, rhs);
    }

    inline bool operator<( Vector<T,DIM> rhs) const
    {
        return lexicographically_less( *this, rhs);
    }

    inline bool operator<=( Vector<T,DIM> rhs) const
    {
        return lexicographically_less_or_equal( *this, rhs);
    }

    inline bool operator>( Vector<T,DIM> rhs) const
    {
        return lexicographically_greater( *this, rhs);
    }

    inline bool operator>=( Vector<T,DIM> rhs) const
    {
        return lexicographically_greater_or_equal( *this, rhs);
    }
};


//------ Free operators +=, -=, *=, /=, +, -, *, and / for vectors -------------

template <typename T, Size DIM>
inline Vector<T,DIM>& operator+=(
    Vector<T,DIM>&       lhs,
    const Vector<T,DIM>& rhs)
{
    for( Size i(0u); i < DIM; ++i)
        lhs[i] += rhs[i];
    return lhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM>& operator-=(
    Vector<T,DIM>&       lhs,
    const Vector<T,DIM>& rhs)
{
    for( Size i(0u); i < DIM; ++i)
        lhs[i] -= rhs[i];
    return lhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM>& operator*=(
    Vector<T,DIM>&       lhs,
    const Vector<T,DIM>& rhs)
{
    for( Size i(0u); i < DIM; ++i)
        lhs[i] *= rhs[i];
    return lhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM>& operator%=(
    Vector<T,DIM>&       lhs,
    const Vector<T,DIM>& rhs)
{
    for( Size i(0u); i < DIM; ++i)
        lhs[i] %= rhs[i];
    return lhs;
}

template <typename T, typename U, Size DIM>
inline Vector<T,DIM>& operator/=(
    Vector<T,DIM>&       lhs,
    const Vector<U,DIM>& rhs)
{
    for( Size i(0u); i < DIM; ++i)
        lhs[i] = T(lhs[i] / rhs[i]);
    return lhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM> operator+(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<T,DIM> tmp( lhs);
    return tmp += rhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM> operator-(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<T,DIM> tmp( lhs);
    return tmp -= rhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM> operator*(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<T,DIM> tmp( lhs);
    return tmp *= rhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM> operator%(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<T,DIM> tmp( lhs);
    return tmp %= rhs;
}

template <typename T, typename U, Size DIM>
inline Vector<T,DIM> operator/(
    const Vector<T,DIM>& lhs,
    const Vector<U,DIM>& rhs)
{
    Vector<T,DIM> tmp(lhs);
    return tmp /= rhs;
}

template <typename T, Size DIM>
inline Vector<T,DIM> operator-( const Vector<T,DIM>& v)
{
    Vector<T,DIM> tmp;
    for( Size i(0u); i < DIM; ++i)
        tmp[i] = -v[i];
    return tmp;
}


//------ Free operator *=, /=, *, and / definitions for scalars ---------------

template <typename T, typename TT, Size DIM>
inline Vector<T,DIM>& operator*=(
    Vector<T,DIM>& v,
    TT             s)
{
    for( Size i(0u); i < DIM; ++i)
        v[i] = T(v[i] * s);
    return v;
}

template <typename T, typename TT, Size DIM>
inline Vector<T,DIM>& operator%=(
    Vector<T,DIM>& v,
    TT             s)
{
    for( Size i(0u); i < DIM; ++i)
        v[i] = T(v[i] % s);
    return v;
}

template <typename T, typename TT, Size DIM>
inline Vector<T,DIM>& operator/=(
    Vector<T,DIM>&  v,
    TT              s)
{
    for( Size i(0u); i < DIM; ++i)
        v[i] = T(v[i] / s);
    return v;
}

template <typename T, typename TT, Size DIM>
inline Vector<T,DIM> operator*(
    const Vector<T,DIM>& v,
    TT                   s)
{
    Vector<T,DIM> tmp( v);
    return tmp *= s;
}

template <typename T, typename TT, Size DIM>
inline Vector<T,DIM> operator*(
    TT                   s,
    const Vector<T,DIM>& v)
{
    Vector<T,DIM> tmp(v);
    return tmp *= s;
}

template <typename T, typename TT, Size DIM>
inline Vector<T,DIM> operator%(
    const Vector<T,DIM>& v,
    TT                   s)
{
    Vector<T,DIM> tmp(v);
    return tmp %= s;
}

template <typename T, typename TT, Size DIM>
inline Vector<T,DIM> operator/(
    const Vector<T,DIM>& v,
    TT                   s)
{
    Vector<T,DIM> tmp( v);
    return tmp /= s;
}


//------ Free operator ++, -- for vectors -------------------------------------

template <typename T, Size DIM>
inline Vector<T,DIM>& operator++( Vector<T,DIM>& vec)
{
    general::for_each( vec, functor::Operator_pre_incr());
    return vec;
}

template <typename T, Size DIM>
inline Vector<T,DIM>& operator--( Vector<T,DIM>& vec)
{
    general::for_each( vec, functor::Operator_pre_decr());
    return vec;
}


//------ Free operators !, &&, ||, ^ for bool vectors and bool scalars ---------

template <Size DIM>
inline Vector<bool,DIM> operator&&(
    const Vector<bool,DIM>& lhs,
    const Vector<bool,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result, functor::Operator_and_and());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator&&(
    bool                    lhs,
    const Vector<bool,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform_left_scalar( lhs, rhs, result, functor::Operator_and_and());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator&&(
    const Vector<bool,DIM>& lhs,
    bool                    rhs)
{
    Vector<bool,DIM> result;
    general::transform_right_scalar( lhs, rhs, result, functor::Operator_and_and());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator||(
    const Vector<bool,DIM>& lhs,
    const Vector<bool,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform(lhs, rhs, result, functor::Operator_or_or());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator||(
    bool                    lhs,
    const Vector<bool,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform_left_scalar( lhs, rhs, result, functor::Operator_or_or());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator||(
    const Vector<bool,DIM>& lhs,
    bool                     rhs)
{
    Vector<bool,DIM> result;
    general::transform_right_scalar( lhs, rhs, result, functor::Operator_or_or());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator^(
    const Vector<bool,DIM>& lhs,
    const Vector<bool,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result, functor::Operator_xor());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator^(
    bool                    lhs,
    const Vector<bool,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform_left_scalar( lhs, rhs, result, functor::Operator_xor());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator^(
    const Vector<bool,DIM>& lhs,
    bool                    rhs)
{
    Vector<bool,DIM> result;
    general::transform_right_scalar( lhs, rhs, result, functor::Operator_xor());
    return result;
}

template <Size DIM>
inline Vector<bool,DIM> operator!(
    const Vector<bool,DIM>& vec)
{
    Vector<bool,DIM> result;
    general::transform( vec, result, functor::Operator_not());
    return result;
}


//------ Elementwise comparison operators returning a bool vector. ------------

template <typename T, Size DIM>
inline Vector<bool,DIM> elementwise_is_equal(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result,functor::Operator_equal_equal());
    return result;
}

template <typename T, Size DIM>
inline Vector<bool,DIM> elementwise_is_not_equal(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result,functor::Operator_not_equal());
    return result;
}

template <typename T, Size DIM>
inline Vector<bool,DIM> elementwise_is_less_than(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result,functor::Operator_less());
    return result;
}

template <typename T, Size DIM>
inline Vector<bool,DIM> elementwise_is_less_than_or_equal(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result,functor::Operator_less_equal());
    return result;
}

template <typename T, Size DIM>
inline Vector<bool,DIM> elementwise_is_greater_than(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result,functor::Operator_greater());
    return result;
}

template <typename T, Size DIM>
inline Vector<bool,DIM> elementwise_is_greater_than_or_equal(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<bool,DIM> result;
    general::transform( lhs, rhs, result,functor::Operator_greater_equal());
    return result;
}


//------ Function Overloads for Vector Algorithms -----------------------------

template <typename T, Size DIM>
inline Vector<T,DIM> abs( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = abs( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> acos( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = acos( v[i]);
    return result;
}

template <typename T, Size DIM>
inline bool all( const Vector<T,DIM>& v)
{
    for( Size i = 0; i != DIM; ++i)
        if( !all(v[i]))
            return false;
    return true;
}

template <typename T, Size DIM>
inline bool any( const Vector<T,DIM>& v)
{
    for( Size i = 0; i != DIM; ++i)
        if( any(v[i]))
           return true;
    return false;
}

template <typename T, Size DIM>
inline Vector<T,DIM> asin( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = asin( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> atan( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = atan( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> atan2( const Vector<T,DIM>& v,  const Vector<T,DIM>& w)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = atan2( v[i], w[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> ceil( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = ceil( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> clamp(
    const Vector<T,DIM>&  v,
    const Vector<T,DIM>&  low,
    const Vector<T,DIM>&  high)
{
    Vector<T,DIM> result;
    for( Size i = 0u; i < DIM; ++i)
        result[i] = clamp( v[i], low[i], high[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> clamp(
    const Vector<T,DIM>& v,
    const Vector<T,DIM>& low,
    T                    high)
{
    Vector<T,DIM> result;
    for( Size i = 0u; i < DIM; ++i)
        result[i] = clamp( v[i], low[i], high);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> clamp(
    const Vector<T,DIM>& v,
    T                    low,
    const Vector<T,DIM>& high)
{
    Vector<T,DIM> result;
    for( Size i = 0u; i < DIM; ++i)
        result[i] = clamp( v[i], low, high[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> clamp(
    const Vector<T,DIM>& v,
    T                    low,
    T                    high)
{
    Vector<T,DIM> result;
    for( Size i = 0u; i < DIM; ++i)
        result[i] = clamp( v[i], low, high);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> cos( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = cos( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> degrees( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = degrees( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> elementwise_max(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<T,DIM> r;
    for( Size i(0u); i < Vector<T,DIM>::DIMENSION; ++i)
        r[i] = base::max MI_PREVENT_MACRO_EXPAND ( lhs[i], rhs[i] );
    return r;
}

template <typename T, Size DIM>
inline Vector<T,DIM> elementwise_min(
    const Vector<T,DIM>& lhs,
    const Vector<T,DIM>& rhs)
{
    Vector<T,DIM> r;
    for( Size i(0u); i < Vector<T,DIM>::DIMENSION; ++i)
        r[i] = base::min MI_PREVENT_MACRO_EXPAND ( lhs[i], rhs[i] );
    return r;
}

template <typename T, Size DIM>
inline Vector<T,DIM> exp( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = exp( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> exp2( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = exp2( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> floor( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = floor( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> fmod( const Vector<T,DIM>& a, const Vector<T,DIM>& b)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = fmod( a[i], b[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> fmod( const Vector<T,DIM>& a, T b)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = fmod( a[i], b);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> frac( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = frac( v[i]);
    return result;
}

template <typename T, Size DIM>
inline bool is_approx_equal(
    const Vector<T,DIM>& left,
    const Vector<T,DIM>& right,
    T                    e)
{
    for( Size i = 0u; i < DIM; ++i)
        if( !is_approx_equal( left[i], right[i], e))
            return false;
    return true;
}

template <typename T, Size DIM>
inline Vector<T,DIM> lerp(
    const Vector<T,DIM>& v1,  
    const Vector<T,DIM>& v2,  
    const Vector<T,DIM>& t)   
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = v1[i] * (T(1)-t[i]) + v2[i] * t[i];
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> lerp(
    const Vector<T,DIM>& v1,  
    const Vector<T,DIM>& v2,  
    T          t)             
{
    // equivalent to: return v1 * (T(1)-t) + v2 * t;
    Vector<T,DIM> result;
    T t2 = T(1) - t;
    for( Size i = 0; i != DIM; ++i)
        result[i] = v1[i] * t2 + v2[i] * t;
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> log( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = log( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> log2 MI_PREVENT_MACRO_EXPAND ( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = log2 MI_PREVENT_MACRO_EXPAND ( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> log10( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = log10( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> modf( const Vector<T,DIM>& v, Vector<T,DIM>& i)
{
    Vector<T,DIM> result;
    for( Size j = 0; j != DIM; ++j)
        result[j] = modf( v[j], i[j]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> pow( const Vector<T,DIM>& a, const Vector<T,DIM>& b)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = pow( a[i], b[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> pow( const Vector<T,DIM>& a, T b)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = pow( a[i], b);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> radians( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = radians( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> round( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = round( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> rsqrt( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = rsqrt( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> saturate( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = saturate( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> sign( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = sign( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> sin( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = sin( v[i]);
    return result;
}

template <typename T, Size DIM>
inline void sincos( const Vector<T,DIM>& a, Vector<T,DIM>& s, Vector<T,DIM>& c)
{
    for( Size i = 0; i != DIM; ++i)
        sincos( a[i], s[i], c[i]);
}

template <typename T, Size DIM>
inline Vector<T,DIM> smoothstep(
    const Vector<T,DIM>& a,
    const Vector<T,DIM>& b,
    const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = smoothstep( a[i], b[i], v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> smoothstep(
    const Vector<T,DIM>& a,
    const Vector<T,DIM>& b,
    T x)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = smoothstep( a[i], b[i], x);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> sqrt( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = sqrt( v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> step( const Vector<T,DIM>& a, const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = step( a[i], v[i]);
    return result;
}

template <typename T, Size DIM>
inline Vector<T,DIM> tan( const Vector<T,DIM>& v)
{
    Vector<T,DIM> result;
    for( Size i = 0; i != DIM; ++i)
        result[i] = tan( v[i]);
    return result;
}


//------ Geometric Vector Algorithms ------------------------------------------

template <typename T>
inline T cross(
    const Vector<T,2>& lhs,
    const Vector<T,2>& rhs)
{
    return lhs.x * rhs.y - lhs.y * rhs.x;
}

template <typename T>
inline Vector<T,3> cross(
    const Vector<T,3>& lhs,
    const Vector<T,3>& rhs)
{
    return Vector<T,3>( lhs.y * rhs.z - lhs.z * rhs.y,
                        lhs.z * rhs.x - lhs.x * rhs.z,
                        lhs.x * rhs.y - lhs.y * rhs.x);
}

template <typename T>
inline void make_basis(
    const Vector<T,3>& n,  
    Vector<T,3>*       u,  
    Vector<T,3>*       v)  
{
#ifdef mi_base_assert_enabled
    const T eps    = 1e-6f;       // smallest resolvable factor
#endif

    mi_math_assert_msg( u != 0, "precondition");
    mi_math_assert_msg( v != 0, "precondition");
    // Sanity check: the normal vector must be unit length.
    mi_math_assert_msg( abs( length(n) - 1.0f) < eps, "precondition");

    // Compute u.
    if( abs(n.x) < abs(n.y)) {
        // u = cross(x, n), x = (1, 0, 0)
        u->x = T(0);
        u->y = -n.z;
        u->z =  n.y;
    } else {
        // u = cross(y, n), y = (0, 1, 0)
        u->x =  n.z;
        u->y = T(0);
        u->z = -n.x;
    }
    u->normalize();

    // Compute v. Since *u and n are orthogonal and unit-length,
    // there is no need to normalize *v.
    *v = cross( *u, n);

    // Sanity check: make sure (u, n, v) is an orthogonal basis.
    mi_math_assert_msg( abs( dot( *u,  n)) < eps, "postcondition");
    mi_math_assert_msg( abs( dot( *u, *v)) < eps, "postcondition");
    mi_math_assert_msg( abs( dot(  n, *v)) < eps, "postcondition");
    // Sanity check: make sure u and v are unit length.
    mi_math_assert_msg( abs( length( *u) - T(1)) < eps, "postcondition");
    mi_math_assert_msg( abs( length( *v) - T(1)) < eps, "postcondition");
}

template <typename T>
inline void make_basis(
    const Vector<T,3>& n,     
    const Vector<T,3>& u,     
    const Vector<T,3>& v,     
    Vector<T,3>*       t,     
    Vector<T,3>*       b)     
{
    const T eps    = 1e-6f;       // smallest resolvable factor
    (void)eps;

    mi_math_assert_msg( t != 0, "precondition");
    mi_math_assert_msg( b != 0, "precondition");
    // Sanity check: the normal vector must be unit length.
    mi_math_assert_msg( abs( length( n) - 1.0f) < eps, "precondition");
    // Sanity check: the other vector lengths should be finite and non-zero
    mi_math_assert_msg( length( u) > 0., "precondition");
    mi_math_assert_msg( length( v) > 0., "precondition");
    mi_math_assert_msg( isfinite( length( u)), "precondition");
    mi_math_assert_msg( isfinite( length( v)), "precondition");

    // Compute b
    *b = cross(u,n);
    b->normalize();

    // Compute t. Since *b and n are orthogonal and unit-length,
    // there is no need to normalize *t.
    *t = cross(n,*b);

    // Check that b has the same orientation of v
    if( dot( *b,v) < T(0))
        *b = -*b;

    // Sanity check: make sure *u and t have the same orientation.
    mi_math_assert_msg( dot( u, *t) > T(0), "postcondition");
    // Sanity check: make sure (t, n, b) is an orthogonal basis.
    // We use a scaled epsilon in order to avoid false positives.
    mi_math_assert_msg( abs( dot( *t,  n)) < 20*eps, "postcondition");
    mi_math_assert_msg( abs( dot( *t, *b)) < 20*eps, "postcondition");
    mi_math_assert_msg( abs( dot(  n, *b)) < 20*eps, "postcondition");
    // Sanity check: make sure t and b are unit length.
    mi_math_assert_msg( abs( length( *t) - T(1)) < eps, "postcondition");
    mi_math_assert_msg( abs( length( *b) - T(1)) < eps, "postcondition");
}

template <typename T2, Size DIM2, typename T1, Size DIM1>
inline Vector<T2, DIM2> convert_vector(
    const Vector<T1, DIM1>& v,
    const T2& fill = T2(0))
{
    const Size dim_min = base::min MI_PREVENT_MACRO_EXPAND ( DIM1, DIM2 );
    Vector<T2, DIM2> result;
    for( Size i = 0; i < dim_min; ++i)
        result[i] = T2(v[i]);
    for( Size i = dim_min; i < DIM2; ++i)
        result[i] = fill;
    return result;
}
 // end group mi_math_vector

} // namespace math

} // namespace mi

#endif // MI_MATH_VECTOR_H