function.h

Go to the documentation of this file.

/***************************************************************************************************
 * Copyright 2024 NVIDIA Corporation. All rights reserved.
 **************************************************************************************************/
 
#ifndef MI_MATH_FUNCTION_H
#define MI_MATH_FUNCTION_H
 
#include <mi/base/assert.h>
#include <mi/base/types.h>
 
#ifdef MI_PLATFORM_WINDOWS
#include <intrin.h>
#pragma intrinsic(_BitScanReverse)
#ifdef MI_ARCH_64BIT
#pragma intrinsic(_BitScanReverse64)
#endif
#endif
 
namespace mi {
 
namespace math {
 
namespace functor {
struct Operator_equal_equal {
    template <typename T1, typename T2>
    inline bool operator()( const T1& t1, const T2& t2) const { return t1 == t2; }
};
 
struct Operator_not_equal {
    template <typename T1, typename T2>
    inline bool operator()( const T1& t1, const T2& t2) const { return t1 != t2; }
};
 
struct Operator_less {
    template <typename T1, typename T2>
    inline bool operator()( const T1& t1, const T2& t2) const { return t1 < t2; }
};
 
struct Operator_less_equal {
    template <typename T1, typename T2>
    inline bool operator()( const T1& t1, const T2& t2) const { return t1 <= t2; }
};
 
struct Operator_greater {
    template <typename T1, typename T2>
    inline bool operator()( const T1& t1, const T2& t2) const { return t1 > t2; }
};
 
struct Operator_greater_equal {
    template <typename T1, typename T2>
    inline bool operator()( const T1& t1, const T2& t2) const { return t1 >= t2; }
};
 
struct Operator_plus {
    template <typename T>
    inline T operator()( const T& t1, const T& t2) const { return t1 + t2; }
};
 
struct Operator_minus {
    template <typename T>
    inline T operator()( const T& t)               const { return - t; }
    template <typename T>
    inline T operator()( const T& t1, const T& t2) const { return t1 - t2; }
};
 
struct Operator_multiply {
    template <typename T>
    inline T operator()( const T& t1, const T& t2) const { return t1 * t2; }
};
 
struct Operator_divide {
    template <typename T>
    inline T operator()( const T& t1, const T& t2) const { return t1 / t2; }
};
 
struct Operator_and_and {
    template <typename T>
    inline bool operator()( const T& t1, const T& t2) const { return t1 && t2; }
};
 
struct Operator_or_or {
    template <typename T>
    inline bool operator()( const T& t1, const T& t2) const { return t1 || t2; }
};
 
struct Operator_xor {
    template <typename T>
    inline bool operator()( const T& t1, const T& t2) const { return t1 ^ t2; }
};
 
struct Operator_not {
    template <typename T>
    inline bool operator()( const T& t) const { return ! t; }
};
 
struct Operator_pre_incr {
    template <typename T>
    inline T operator()( T& t) const { return ++ t; }
};
 
struct Operator_post_incr {
    template <typename T>
    inline T operator()( T& t) const { return t ++; }
};
 
struct Operator_pre_decr {
    template <typename T>
    inline T operator()( T& t) const { return -- t; }
};
 
struct Operator_post_decr {
    template <typename T>
    inline T operator()( T& t) const { return t --; }
};
 // end group mi_math_functor
 
} // namespace functor
 
namespace general {
template <class Vector, class ResultVector, class UnaryFunctor>
inline void transform( const Vector& vec, ResultVector& result, UnaryFunctor f)
{
    mi_static_assert( Vector::SIZE == ResultVector::SIZE);
    for( Size i = 0; i != Vector::SIZE; ++i)
        result.set( i, f( vec.get(i)));
}
 
template <class Vector1, class Vector2, class ResultVector, class BinaryFunctor>
inline void transform(
    const Vector1& vec1, const Vector2& vec2, ResultVector& result, BinaryFunctor f)
{
    mi_static_assert( Vector1::SIZE == Vector2::SIZE);
    mi_static_assert( Vector1::SIZE == ResultVector::SIZE);
    for( Size i = 0; i != Vector1::SIZE; ++i)
        result.set( i, f( vec1.get(i), vec2.get(i)));
}
 
template <class Scalar, class Vector, class ResultVector, class BinaryFunctor>
inline void transform_left_scalar(
    const Scalar& s, const Vector& vec, ResultVector& result, BinaryFunctor f)
{
    mi_static_assert( Vector::SIZE == ResultVector::SIZE);
    for( Size i = 0; i != Vector::SIZE; ++i)
        result.set( i, f( s, vec.get(i)));
}
 
template <class Scalar, class Vector, class ResultVector, class BinaryFunctor>
inline void transform_right_scalar(
    const Vector& vec, const Scalar& s, ResultVector& result, BinaryFunctor f)
{
    mi_static_assert( Vector::SIZE == ResultVector::SIZE);
    for( Size i = 0; i != Vector::SIZE; ++i)
        result.set( i, f( vec.get(i), s));
}
 
template <class Vector, class UnaryFunctor>
inline void for_each( Vector& vec, UnaryFunctor f)
{
    for( Size i = 0; i != Vector::SIZE; ++i)
        f( vec.begin()[i]);
}
 
template <class Vector1, class Vector2, class BinaryFunctor>
inline void for_each( Vector1& vec1, const Vector2& vec2, BinaryFunctor f)
{
    mi_static_assert( Vector1::SIZE == Vector2::SIZE);
    for( Size i = 0; i != Vector1::SIZE; ++i)
        f( vec1.begin()[i], vec2.begin()[i]);
}
 // end group mi_math_functor
 
} // namespace general
 
inline Float32 exp( Float32 s) { return std::exp(s); }
inline Float64 exp( Float64 s) { return std::exp(s); }
 
inline Float32 log( Float32 s) { return std::log(s); }
inline Float64 log( Float64 s) { return std::log(s); }
 
 
#ifndef __CUDACC__
using mi::base::min;
using mi::base::max;
#endif
using mi::base::abs;
 
 
inline Float32 fast_sqrt( Float32 i)
{
    int tmp = base::binary_cast<int>( i);
    tmp -= 1 << 23;     // Remove last bit to not let it go to mantissa
    // tmp is now an approximation to logbase2(i)
    tmp = tmp >> 1;     // Divide by 2
    tmp += 1 << 29;     // Add 64 to the exponent: (e+127)/2 = (e/2)+63
                        // that represents (e/2)-64 but we want e/2
    return base::binary_cast<Float32>( tmp);
}
 
inline Float32 fast_exp( Float32 x)
{
    const Float32 EXP_C   = 8388608.0f; // 2^23
    const Float32 LOG_2_E = 1.4426950408889634073599246810019f; // 1 / log(2)
 
    x *= LOG_2_E;
    Float32 y = x - std::floor(x);
    y = (y - y*y) * 0.33971f;
    const Float32 z = max MI_PREVENT_MACRO_EXPAND ((x + 127.0f - y) * EXP_C, 0.f);
    return base::binary_cast<Float32>( static_cast<int>( z));
}
 
inline Float32 fast_pow2( Float32 x)
{
    const Float32 EXP_C = 8388608.0f; // 2^23
 
    Float32 y = x - std::floor(x);
    y = (y - y*y) * 0.33971f;
    const Float32 z = max MI_PREVENT_MACRO_EXPAND ((x + 127.0f - y) * EXP_C, 0.f);
    return base::binary_cast<Float32>( static_cast<int>( z));
}
 
inline Float32 fast_log2( Float32 i)
{
    const Float32 LOG_C = 0.00000011920928955078125f; // 2^(-23)
 
    const Float32 x = static_cast<Float32>( base::binary_cast<int>( i)) * LOG_C - 127.f;
    const Float32 y = x - std::floor(x);
    return x + (y - y*y) * 0.346607f;
}
 
inline Float32 fast_pow(
    Float32 b,  
    Float32 e)  
{
    if( b == 0.0f)
        return 0.0f;
 
    const Float32 LOG_C = 0.00000011920928955078125f; // 2^(-23)
    const Float32 EXP_C = 8388608.0f; // 2^23
 
    const Float32 x  = static_cast<Float32>( base::binary_cast<int>( b)) * LOG_C - 127.f;
    const Float32 y  = x - std::floor(x);
    const Float32 fl = e * (x + (y - y*y) * 0.346607f);
    const Float32 y2 = fl - std::floor(fl);
    const Float32 z  = max(
        fl*EXP_C + static_cast<Float32>( 127.0*EXP_C)
                 - (y2 - y2*y2) * static_cast<Float32>( 0.33971*EXP_C), 0.f);
 
    return base::binary_cast<Float32>( static_cast<int>( z));
}
 // end group mi_math_approx_function
 
 
inline Float32 acos( Float32 s) { return std::acos(s); }
inline Float64 acos( Float64 s) { return std::acos(s); }
 
inline bool all( Uint8   v) { return Uint8  (0) != v; }
inline bool all( Uint16  v) { return Uint16 (0) != v; }
inline bool all( Uint32  v) { return Uint32 (0) != v; }
inline bool all( Uint64  v) { return Uint64 (0) != v; }
inline bool all( Sint8   v) { return Sint8  (0) != v; }
inline bool all( Sint16  v) { return Sint16 (0) != v; }
inline bool all( Sint32  v) { return Sint32 (0) != v; }
inline bool all( Sint64  v) { return Sint64 (0) != v; }
inline bool all( Float32 v) { return Float32(0) != v; }
inline bool all( Float64 v) { return Float64(0) != v; }
 
inline bool any( Uint8   v) { return Uint8  (0) != v; } //-V524 PVS
inline bool any( Uint16  v) { return Uint16 (0) != v; } //-V524 PVS
inline bool any( Uint32  v) { return Uint32 (0) != v; } //-V524 PVS
inline bool any( Uint64  v) { return Uint64 (0) != v; } //-V524 PVS
inline bool any( Sint8   v) { return Sint8  (0) != v; } //-V524 PVS
inline bool any( Sint16  v) { return Sint16 (0) != v; } //-V524 PVS
inline bool any( Sint32  v) { return Sint32 (0) != v; } //-V524 PVS
inline bool any( Sint64  v) { return Sint64 (0) != v; } //-V524 PVS
inline bool any( Float32 v) { return Float32(0) != v; } //-V524 PVS
inline bool any( Float64 v) { return Float64(0) != v; } //-V524 PVS
 
inline Float32 asin( Float32 s) { return std::asin(s); }
inline Float64 asin( Float64 s) { return std::asin(s); }
 
inline Float32 atan( Float32 s) { return std::atan(s); }
inline Float64 atan( Float64 s) { return std::atan(s); }
 
inline Float32 atan2( Float32 s, Float32 t) { return std::atan2(s,t); }
inline Float64 atan2( Float64 s, Float64 t) { return std::atan2(s,t); }
 
inline Float32 ceil( Float32 s) { return std::ceil(s); }
inline Float64 ceil( Float64 s) { return std::ceil(s); }
 
inline  Uint8 clamp(  Uint8 s,  Uint8 low,  Uint8 high)
{
    return min(high, max(low, s));
}
inline Uint16 clamp( Uint16 s, Uint16 low, Uint16 high)
{
    return min(high, max(low, s));
}
inline Uint32 clamp( Uint32 s, Uint32 low, Uint32 high)
{
    return min(high, max(low, s));
}
inline Uint64 clamp( Uint64 s, Uint64 low, Uint64 high)
{
    return min(high, max(low, s));
}
inline  Sint8 clamp(  Sint8 s,  Sint8 low,  Sint8 high)
{
    return min(high, max(low, s));
}
inline Sint16 clamp( Sint16 s, Sint16 low, Sint16 high)
{
    return min(high, max(low, s));
}
inline Sint32 clamp( Sint32 s, Sint32 low, Sint32 high)
{
    return min(high, max(low, s));
}
inline Sint64 clamp( Sint64 s, Sint64 low, Sint64 high)
{
    return min(high, max(low, s));
}
inline Float32 clamp( Float32 s, Float32 low, Float32 high)
{
    return min(high, max(low, s));
}
inline Float64 clamp( Float64 s, Float64 low, Float64 high)
{
    return min(high, max(low, s));
}
 
inline Float32 cos( Float32 a) { return std::cos(a); }
inline Float64 cos( Float64 a) { return std::cos(a); }
 
inline Float32 degrees( Float32 r) { return r * Float32(180.0/MI_PI); }
inline Float64 degrees( Float64 r) { return r * Float64(180.0/MI_PI); }
 
inline Float32 exp2( Float32 s) { return fast_pow2(s); }
inline Float64 exp2( Float64 s)
{
    return std::exp(s * 0.69314718055994530941723212145818 /* log(2) */ );
}
 
inline Float32 floor( Float32 s) { return std::floor(s); }
inline Float64 floor( Float64 s) { return std::floor(s); }
 
inline Float32 fmod( Float32 a, Float32 b) { return std::fmod(a,b); }
inline Float64 fmod( Float64 a, Float64 b) { return std::fmod(a,b); }
 
inline Float32 frac( Float32 s)
{
    Float32 dummy;
    if( s >= 0.0f)
        return std::modf( s, &dummy);
    else
        return 1.0f + std::modf( s, &dummy);
}
inline Float64 frac( Float64 s)
{
    Float64 dummy;
    if( s >= 0.0f)
        return std::modf( s, &dummy);
    else
        return 1.0f + std::modf( s, &dummy);
}
 
inline bool is_approx_equal(
    Float32 left,
    Float32 right,
    Float32 e)
{
    return abs( left - right ) <= e;
}
 
inline bool is_approx_equal(
    Float64 left,
    Float64 right,
    Float64 e)
{
    return abs( left - right ) <= e;
}
 
inline Uint32 leading_zeros(Uint32 v) {
    // This implementation tries to use built-in functions if available. For the fallback
    // method, see Henry Warren: "Hacker's Delight" for reference.
#if defined(MI_COMPILER_MSC)
    unsigned long index;
    const unsigned char valid = _BitScanReverse(&index, v);
    return (valid != 0) ? 31 - index : 32;
#elif defined(MI_COMPILER_ICC) || defined(MI_COMPILER_GCC)
    return (v != 0) ? __builtin_clz(v) : 32;
#else
    // use fallback
    if (v == 0) return 32;
    Uint32 n = 1;
    if ((v >> 16) == 0) { n += 16; v <<= 16; };
    if ((v >> 24) == 0) { n += 8;  v <<= 8; };
    if ((v >> 28) == 0) { n += 4;  v <<= 4; };
    if ((v >> 30) == 0) { n += 2;  v <<= 2; };
    n -= Uint32(v >> 31);
    return n;
#endif
}
 
inline Uint32 leading_zeros(Uint64 v) {
    // This implementation tries to use built-in functions if available. For the fallback
    // method, see Henry Warren: "Hacker's Delight" for reference.
#if defined(MI_COMPILER_MSC)
#if defined(MI_ARCH_64BIT)
    unsigned long index;
    const unsigned char valid = _BitScanReverse64(&index, v);
    return (valid != 0) ? 63 - index : 64;
#else
    unsigned long index_h, index_l;
    const unsigned char valid_h = _BitScanReverse(&index_h,(Uint32)(v >> 32));
    const unsigned char valid_l = _BitScanReverse(&index_l,(Uint32)(v & 0xFFFFFFFF));
    if (valid_h == 0)
    return (valid_l != 0) ? 63 - index_l : 64;
    return 63 - index_h + 32;
#endif
#elif defined(MI_COMPILER_ICC) || defined(MI_COMPILER_GCC)
    return (v != 0) ? __builtin_clzll(v) : 64;
#else
    // use fallback
    if (v == 0) return 64;
    Uint32 n = 1;
    if ((v >> 32) == 0) { n += 32; v <<= 32; };
    if ((v >> 48) == 0) { n += 16; v <<= 16; };
    if ((v >> 56) == 0) { n += 8;  v <<= 8; };
    if ((v >> 60) == 0) { n += 4;  v <<= 4; };
    if ((v >> 62) == 0) { n += 2;  v <<= 2; };
    n -= Uint32(v >> 63);
    return n;
#endif
}
 
inline Float32 lerp(
    Float32 s1,  
    Float32 s2,  
    Float32 t)   
{
    return s1 * (Float32(1)-t) + s2 * t;
}
 
inline Float64 lerp(
    Float64 s1,  
    Float64 s2,  
    Float64 t)   
{
    return s1 * (Float64(1)-t) + s2 * t;
}
 
inline Float32 log2 MI_PREVENT_MACRO_EXPAND ( Float32 s)
{ return std::log(s) * 1.4426950408889634073599246810019f /* log(2) */; }
inline Float64 log2 MI_PREVENT_MACRO_EXPAND ( Float64 s)
{ return std::log(s) * 1.4426950408889634073599246810019 /* log(2) */; }
 
inline Sint32 log2_int( const Uint32 v) { return (v != 0) ? 31 - leading_zeros(v) : 0; }
 
inline Sint32 log2_int( const Uint64 v) { return (v != 0) ? 63 - leading_zeros(v) : 0; }
 
inline Sint32 log2_int( const Float32 v) { 
    return (mi::base::binary_cast<Uint32>(v) >> 23) - 127; 
}
 
inline Sint32 log2_int( const Float64 v) { 
    return static_cast<Sint32>(mi::base::binary_cast<Uint64>(v) >> 52) - 1023;
}
 
template<typename Integer>
inline Sint32 log2_int_ceil( const Integer v) {
    // See Henry Warren: "Hacker's Delight" for reference.
    return (v > 1) ? log2_int(v - 1) + 1 : 0;
}
 
inline Float32 log10( Float32 s) { return std::log10(s); }
inline Float64 log10( Float64 s) { return std::log10(s); }
 
inline Float32 modf( Float32 s, Float32& i) { return std::modf( s, &i); }
inline Float64 modf( Float64 s, Float64& i) { return std::modf( s, &i); }
 
inline Uint32  pow( Uint32  a, Uint32  b) { return Uint32(std::pow(double(a), int(b))); }
inline Uint64  pow( Uint64  a, Uint64  b) { return Uint64(std::pow(double(a), int(b))); }
inline Sint32  pow( Sint32  a, Sint32  b) { return Sint32(std::pow(double(a), int(b))); }
inline Sint64  pow( Sint64  a, Sint64  b) { return Sint64(std::pow(double(a), int(b))); }
inline Float32 pow( Float32 a, Float32 b) { return std::pow( a, b); }
inline Float64 pow( Float64 a, Float64 b) { return std::pow( a, b); }
 
inline Float32 radians( Float32 d) { return d * Float32(MI_PI/180.0); }
inline Float64 radians( Float64 d) { return d * Float64(MI_PI/180.0); }
 
inline Float32 round( Float32 s) { return std::floor(s + 0.5f); }
inline Float64 round( Float64 s) { return std::floor(s + 0.5); }
 
inline Float32 rsqrt( Float32 s) { return 1.0f / std::sqrt(s); }
inline Float64 rsqrt( Float64 s) { return 1.0  / std::sqrt(s); }
 
inline Float32 saturate( Float32 s) { return min(1.f, max(0.f, s)); }
inline Float64 saturate( Float64 s) { return min(1. , max(0. , s)); }
 
inline Sint8   sign( Sint8   s)
{ int r = (s < 0 ) ? -1 : (s > 0) ? 1 : 0; return static_cast<Sint8>(  r); }
inline Sint16  sign( Sint16  s)
{ int r = (s < 0 ) ? -1 : (s > 0) ? 1 : 0; return static_cast<Sint16>( r); }
inline Sint32  sign( Sint32  s) { return (s < 0  )  ? -1 : (s > 0)   ? 1   : 0;   }
inline Sint64  sign( Sint64  s) { return (s < 0  )  ? -1 : (s > 0)   ? 1   : 0;   }
inline Float32 sign( Float32 s) { return (s < 0.f) ? -1.f: (s > 0.f) ? 1.f : 0.f; }
inline Float64 sign( Float64 s) { return (s < 0. ) ? -1. : (s > 0.)  ? 1.  : 0.;  }
 
inline bool   sign_bit( Sint8    s) { return s < 0; }
inline bool   sign_bit( Sint16   s) { return s < 0; }
inline bool   sign_bit( Sint32   s) { return s < 0; }
inline bool   sign_bit( Sint64   s) { return s < 0; }
 
inline bool sign_bit( Float32 s)
{
    return (base::binary_cast<Uint32>(s) & (1U << 31)) != 0U;
}
 
inline bool sign_bit( Float64 s)
{
    return (base::binary_cast<Uint64>(s) & (1ULL << 63)) != 0ULL;
}
 
#if (__cplusplus < 201103L)
inline bool isnan MI_PREVENT_MACRO_EXPAND (const Float32 x)
{
    // interpret as Uint32 value
    const Uint32 f = base::binary_cast<Uint32>(x);
 
    // check bit pattern
    return (f << 1) > 0xFF000000U; // shift sign bit, 8bit exp == 2^8-1, fraction != 0
}
 
inline bool isnan MI_PREVENT_MACRO_EXPAND (const Float64 x)
{
    // interpret as Uint64 value
    const Uint64 f = base::binary_cast<Uint64>(x);
 
    return (f << 1) > 0xFFE0000000000000ULL; // shift sign bit, 11bit exp == 2^11-1, fraction != 0
}
#else
using std::isnan;
#endif
 
 
#ifndef __CUDA_ARCH__
inline float __uint_as_float(const unsigned v)
{ return base::binary_cast<float>(v);}
 
inline unsigned __float_as_uint(const float v)
{ return base::binary_cast<unsigned>(v);}
#endif
 
 
inline bool isinfinite MI_PREVENT_MACRO_EXPAND (const Float32 x)
{
    const Uint32 exponent_mask = 0x7F800000; // 8 bit exponent
    const Uint32 fraction_mask = 0x7FFFFF;   // 23 bit fraction
 
    // interpret as Uint32 value
    const Uint32 f = base::binary_cast<Uint32>(x);
 
    // check bit pattern
    return ((f & exponent_mask) == exponent_mask) && // exp == 2^8 - 1
           ((f & fraction_mask) == 0);               // fraction == 0
}
 
inline bool isinfinite MI_PREVENT_MACRO_EXPAND (const Float64 x)
{
    const Uint64 exponent_mask = 0x7FF0000000000000ULL; // 11 bit exponent
    const Uint64 fraction_mask = 0xFFFFFFFFFFFFFULL;    // 52 bit fraction
 
    // interpret as Uint64 value
    const Uint64 f = base::binary_cast<Uint64>(x);
 
    // check bit pattern
    return ((f & exponent_mask) == exponent_mask) && // exp == 2^11 - 1
           ((f & fraction_mask) == 0);               // fraction == 0
}
 
 
#if (__cplusplus < 201103L)
inline bool isfinite MI_PREVENT_MACRO_EXPAND (const Float32 x)
{
    const Uint32 exponent_mask = 0x7F800000; // 8 bit exponent
 
    // interpret as Uint32 value
    const Uint32 f = base::binary_cast<Uint32>(x);
 
    // check exponent bits
    return ((f & exponent_mask) != exponent_mask); // exp != 2^8 - 1
}
 
inline bool isfinite MI_PREVENT_MACRO_EXPAND (const Float64 x)
{
    const Uint64 exponent_mask = 0x7FF0000000000000ULL; // 11 bit exponent
 
    // interpret as Uint64 value
    const Uint64 f = base::binary_cast<Uint64>(x);
 
    // check exponent bits
    return ((f & exponent_mask) != exponent_mask); // exp != 2^11 - 1
}
#else
using std::isfinite;
#endif
 
 
inline Float32 sin( Float32 a) { return std::sin(a); }
inline Float64 sin( Float64 a) { return std::sin(a); }
 
inline void sincos( Float32 a, Float32& s, Float32& c)
{
    s = std::sin(a);
    c = std::cos(a);
}
inline void sincos( Float64 a, Float64& s, Float64& c)
{
    s = std::sin(a);
    c = std::cos(a);
}
 
inline Float32 smoothstep( Float32 a, Float32 b, Float32 x)
{
    if(x < a)
        return 0.0f;
    if(b < x)
        return 1.0f;
    Float32 t = (x - a) / (b - a);
    return t * t * (3.0f - 2.0f * t);
}
inline Float64 smoothstep( Float64 a, Float64 b, Float64 x)
{
    if(x < a)
        return 0.0;
    if(b < x)
        return 1.0;
    Float64 t = (x - a) / (b - a);
    return t * t * (3.0 - 2.0 * t);
}
 
inline Float32 sqrt( Float32 s) { return std::sqrt(s); }
inline Float64 sqrt( Float64 s) { return std::sqrt(s); }
 
inline Float32 step( Float32 a, Float32 x) { return (x < a) ? 0.0f : 1.0f; }
inline Float64 step( Float64 a, Float64 x) { return (x < a) ? 0.0 : 1.0; }
 
inline Float32 tan( Float32 a) { return std::tan(a); }
inline Float64 tan( Float64 a) { return std::tan(a); }
 
 
MI_HOST_DEVICE_INLINE void to_rgbe( const Float32 color[3], Uint32& rgbe)
{
    Float32 c[3];
    c[0] = max( color[0], 0.0f);
    c[1] = max( color[1], 0.0f);
    c[2] = max( color[2], 0.0f);
 
    const Float32 m = max( max( c[0], c[1]), c[2]);
 
    // should actually be -126 or even -128, but avoid precision problems / denormalized numbers
    if( m <= 7.5231631727e-37f) // ~2^(-120)
        rgbe = 0;
    else if( m >= 1.7014118346046923173168730371588e+38f) // 2^127
        rgbe = 0xFFFFFFFFu;
    else {
        const Uint32  e = __float_as_uint( m) & 0x7F800000u;
        const Float32 v = __uint_as_float( 0x82800000u - e);
 
        rgbe =  Uint32( c[0] * v)
             | (Uint32( c[1] * v) <<  8)
             | (Uint32( c[2] * v) << 16)
             | (e * 2 + (2 << 24));
    }
}
 
MI_HOST_DEVICE_INLINE void to_rgbe( const Float32 color[3], Uint8 rgbe[4])
{
    Float32 c[3];
    c[0] = max( color[0], 0.0f);
    c[1] = max( color[1], 0.0f);
    c[2] = max( color[2], 0.0f);
 
    const Float32 m = max( max( c[0], c[1]), c[2]);
 
    // should actually be -126 or even -128, but avoid precision problems / denormalized numbers
    if( m <= 7.5231631727e-37f) // ~2^(-120)
        rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 0;
    else if( m >= 1.7014118346046923173168730371588e+38f) // 2^127
        rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 255;
    else {
        const Uint32  e = __float_as_uint( m) & 0x7F800000u;
        const Float32 v = __uint_as_float( 0x82800000u - e);
 
        rgbe[0] = Uint8( c[0] * v);
        rgbe[1] = Uint8( c[1] * v);
        rgbe[2] = Uint8( c[2] * v);
        rgbe[3] = Uint8( (e >> 23) + 2);
    }
}
 
MI_HOST_DEVICE_INLINE void from_rgbe( const Uint8 rgbe[4], Float32 color[3])
{
    if( rgbe[3] == 0) {
        color[0] = color[1] = color[2] = 0.0f;
        return;
    }
 
    const Uint32  e = (static_cast<Uint32>( rgbe[3]) << 23) - 0x800000u;
    const Float32 v = __uint_as_float( e);
    const Float32 c = static_cast<Float32>( 1.0 - 0.5/256.0) * v;
 
    color[0] = __uint_as_float( e | (static_cast<Uint32>( rgbe[0]) << 15)) - c;
    color[1] = __uint_as_float( e | (static_cast<Uint32>( rgbe[1]) << 15)) - c;
    color[2] = __uint_as_float( e | (static_cast<Uint32>( rgbe[2]) << 15)) - c;
}
 
MI_HOST_DEVICE_INLINE void from_rgbe( const Uint32 rgbe, Float32 color[3])
{
    const Uint32 rgbe3 = rgbe & 0xFF000000u;
    if( rgbe3 == 0) {
        color[0] = color[1] = color[2] = 0.0f;
        return;
    }
 
    const Uint32  e = (rgbe3 >> 1) - 0x800000u;
    const Float32 v = __uint_as_float( e);
    const Float32 c = static_cast<Float32>( 1.0 - 0.5/256.0) * v;
 
    color[0] = __uint_as_float( e | ((rgbe << 15) & 0x7F8000u)) - c;
    color[1] = __uint_as_float( e | ((rgbe <<  7) & 0x7F8000u)) - c;
    color[2] = __uint_as_float( e | ((rgbe >>  1) & 0x7F8000u)) - c;
}
 
 
//------ Generic Vector Algorithms --------------------------------------------
 
// overloads for 1D vectors (scalars)
 
inline Sint32 dot( Sint32 a, Sint32 b) { return a * b; }
inline Float32 dot( Float32 a, Float32 b) { return a * b; }
inline Float64 dot( Float64 a, Float64 b) { return a * b; }
 
template <class V>
inline typename V::value_type dot( const V& lhs, const V& rhs)
{
    typename V::value_type v(0);
    for( Size i(0u); i < V::SIZE; ++i)
        v += lhs.get(i) * rhs.get(i);
    return v;
}
 
template <class V>
inline typename V::value_type square_length( const V& v)
{
    return dot( v, v);
}
 
// base case for scalars
 
inline Float32 length( Float32 a) { return abs(a); }
inline Float64 length( Float64 a) { return abs(a); }
 
 
template <class V>
inline typename V::value_type length( const V& v)
{
    return sqrt( square_length(v));
}
 
template <class V>
inline typename V::value_type square_euclidean_distance( const V& lhs, const V& rhs)
{
    return square_length( lhs - rhs);
}
 
template <class V>
inline typename V::value_type euclidean_distance(
    const V& lhs, const V& rhs)
{
    return length( lhs - rhs);
}
 
template <class V>
inline void set_bounds( V& v, const V& low, const V& high)
{
    for( Size i(0u); i < V::SIZE; ++i)
        v[i] = min MI_PREVENT_MACRO_EXPAND (
                   max MI_PREVENT_MACRO_EXPAND (v.get(i), low.get(i)), high.get(i));
}
 
template <class V>
inline bool is_equal( const V& lhs, const V& rhs)
{
    for( Size i(0u); i < V::SIZE; ++i)
        if( ! (lhs.get(i) == rhs.get(i)))
            return false;
    return true;
}
 
template <class V>
inline bool is_not_equal( const V& lhs, const V& rhs)
{
    for( Size i(0u); i < V::SIZE; ++i)
        if( lhs.get(i) != rhs.get(i))
            return true;
    return false;
}
 
template <class V>
inline bool lexicographically_less( const V& lhs, const V& rhs)
{
    for( Size i(0u); i < V::SIZE-1; ++i) {
        if( lhs.get(i) < rhs.get(i))
            return true;
        if( lhs.get(i) > rhs.get(i))
            return false;
    }
    return lhs.get(V::SIZE-1) < rhs.get(V::SIZE-1);
}
 
template <class V>
inline bool lexicographically_less_or_equal( const V& lhs, const V& rhs)
{
    for( Size i(0u); i < V::SIZE-1; ++i) {
        if( lhs.get(i) < rhs.get(i))
            return true;
        if( lhs.get(i) > rhs.get(i))
            return false;
    }
    return lhs.get(V::SIZE-1) <= rhs.get(V::SIZE-1);
}
 
template <class V>
inline bool lexicographically_greater( const V& lhs, const V& rhs)
{
    for( Size i(0u); i < V::SIZE-1; ++i) {
        if( lhs.get(i) > rhs.get(i))
            return true;
        if( lhs.get(i) < rhs.get(i))
            return false;
    }
    return lhs.get(V::SIZE-1) > rhs.get(V::SIZE-1);
}
 
template <class V>
inline bool lexicographically_greater_or_equal( const V& lhs, const V& rhs)
{
    for( Size i(0u); i < V::SIZE-1; ++i) {
        if( lhs.get(i) > rhs.get(i))
            return true;
        if( lhs.get(i) < rhs.get(i))
            return false;
    }
    return lhs.get(V::SIZE-1) >= rhs.get(V::SIZE-1);
}
 
template <class V>
inline Comparison_result lexicographically_compare( const V& lhs, const V& rhs)
{
    for( Size i(0u); i < V::SIZE; ++i) {
        Comparison_result result = three_valued_compare( lhs.get(i), rhs.get(i));
        if( result != EQUAL)
            return result;
    }
    return EQUAL;
}
 // end group mi_math_function
 
} // namespace math
 
} // namespace mi
 
#endif // MI_MATH_FUNCTION_H