simd.h

Go to the documentation of this file.

// Copyright 2008-present Contributors to the OpenImageIO project.
// SPDX-License-Identifier: BSD-3-Clause
// https://github.com/OpenImageIO/oiio/blob/master/LICENSE.md

/// @file  simd.h
///
/// @brief Classes for SIMD processing.
///
/// Nice references for all the Intel intrinsics (SSE*, AVX*, etc.):
///   https://software.intel.com/sites/landingpage/IntrinsicsGuide/
///
/// Similar guide for ARM intrinsics:
///   https://developer.arm.com/architectures/instruction-sets/simd-isas/neon/intrinsics
///
/// It helped me a lot to peruse the source of these packages:
///   Syrah:     https://github.com/boulos/syrah
///   Embree:    https://github.com/embree
///   Vectorial: https://github.com/scoopr/vectorial
///
/// To find out which CPU features you have:
///   Linux: cat /proc/cpuinfo
///   OSX:   sysctl machdep.cpu.features
///
/// Additional web resources:
///   http://www.codersnotes.com/notes/maths-lib-2016/

// clang-format off

#pragma once

#include <algorithm>
#include <cstring>

#include <OpenImageIO/Imath.h>
#include <OpenImageIO/dassert.h>
#include <OpenImageIO/platform.h>


//////////////////////////////////////////////////////////////////////////
// Sort out which SIMD capabilities we have and set definitions
// appropriately. This is mostly for internal (within this file) use,
// but client applications using this header may find a few of the macros
// we define to be useful:
//
// OIIO_SIMD : Will be 0 if no hardware SIMD support is specified. If SIMD
//             hardware is available, this will hold the width in number of
//             float SIMD "lanes" of widest SIMD registers available. For
//             example, OIIO_SIMD will be 4 if vfloat4/vint4/vbool4 are
//             hardware accelerated, 8 if vfloat8/vint8/vbool8 are accelerated,
//             etc. Using SIMD classes wider than this should work (will be
//             emulated with narrower SIMD or scalar operations), but is not
//             expected to have high performance.
// OIIO_SIMD_SSE : if Intel SSE is supported, this will be nonzero,
//             specifically 2 for SSE2, 3 for SSSE3, 4 for SSE4.1 or
//             higher (including AVX).
// OIIO_SIMD_AVX : If Intel AVX is supported, this will be nonzero, and
//             specifically 1 for AVX (1.0), 2 for AVX2, 512 for AVX512f.
// OIIO_SIMD_NEON : If ARM NEON is supported, this will be nonzero.
// OIIO_SIMD_MAX_SIZE : holds the width in bytes of the widest SIMD
//             available (generally will be OIIO_SIMD*4).
// OIIO_SIMD4_ALIGN : macro for best alignment of 4-wide SIMD values in mem.
// OIIO_SIMD8_ALIGN : macro for best alignment of 8-wide SIMD values in mem.
// OIIO_SIMD16_ALIGN : macro for best alignment of 16-wide SIMD values in mem.
// OIIO_SIMD_HAS_MATRIX4 : nonzero if matrix44 is defined
// OIIO_SIMD_HAS_SIMD8 : nonzero if vfloat8, vint8, vbool8 are defined
// OIIO_SIMD_HAS_SIMD16 : nonzero if vfloat16, vint16, vbool16 are defined

#if defined(_WIN32)
#  include <intrin.h>
#elif defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__)) || defined(__e2k__)
#  include <x86intrin.h>
#elif defined(__GNUC__) && defined(__ARM_NEON__)
#  include <arm_neon.h>
#endif

// Disable SSE for 32 bit Windows patforms, it's unreliable and hard for us
// to test thoroughly. We presume that anybody needing high performance
// badly enough to want SIMD also is on a 64 bit CPU.
#if defined(_WIN32) && defined(__i386__) && !defined(__x86_64__) && !defined(OIIO_NO_SSE)
#define OIIO_NO_SSE 1
#endif

#if (defined(__SSE2__) || (_MSC_VER >= 1300 && !_M_CEE_PURE)) && !defined(OIIO_NO_SSE)
#  if (defined(__SSE4_1__) || defined(__SSE4_2__))
#    define OIIO_SIMD_SSE 4
      /* N.B. We consider both SSE4.1 and SSE4.2 to be "4". There are a few
       * instructions specific to 4.2, but they are all related to string
       * comparisons and CRCs, which don't currently seem relevant to OIIO,
       * so for simplicity, we sweep this difference under the rug.
       */
#  elif defined(__SSSE3__)
#    define OIIO_SIMD_SSE 3
     /* N.B. We only use OIIO_SIMD_SSE = 3 when fully at SSSE3. In theory,
      * there are a few older architectures that are SSE3 but not SSSE3,
      * and this simplification means that these particular old platforms
      * will only get SSE2 goodness out of our code. So be it. Anybody who
      * cares about performance is probably using a 64 bit machine that's
      * SSE 4.x or AVX by now.
      */
#  else
#    define OIIO_SIMD_SSE 2
#  endif
#  define OIIO_SIMD 4
#  define OIIO_SIMD_MAX_SIZE_BYTES 16
#  define OIIO_SIMD4_ALIGN OIIO_ALIGN(16)
#  define OIIO_SSE_ALIGN OIIO_ALIGN(16)
#else
#  define OIIO_SIMD_SSE 0
#endif

#if defined(__AVX__) && !defined(OIIO_NO_AVX)
   // N.B. Any machine with AVX will also have SSE
#  if defined(__AVX2__) && !defined(OIIO_NO_AVX2)
#    define OIIO_SIMD_AVX 2
#  else
#    define OIIO_SIMD_AVX 1
#  endif
#  undef OIIO_SIMD
#  define OIIO_SIMD 8
#  undef OIIO_SIMD_MAX_SIZE_BYTES
#  define OIIO_SIMD_MAX_SIZE_BYTES 32
#  define OIIO_SIMD8_ALIGN OIIO_ALIGN(32)
#  define OIIO_AVX_ALIGN OIIO_ALIGN(32)
#  if defined(__AVX512F__)
#    undef OIIO_SIMD_AVX
#    define OIIO_SIMD_AVX 512
#    undef OIIO_SIMD_MAX_SIZE_BYTES
#    define OIIO_SIMD_MAX_SIZE_BYTES 64
#    undef OIIO_SIMD
#    define OIIO_SIMD 16
#    define OIIO_SIMD16_ALIGN OIIO_ALIGN(64)
#    define OIIO_AVX512_ALIGN OIIO_ALIGN(64)
#    define OIIO_AVX512F_ENABLED 1
#  endif
#  if defined(__AVX512DQ__)
#    define OIIO_AVX512DQ_ENABLED 1   /* Doubleword and quadword */
#  else
#    define OIIO_AVX512DQ_ENABLED 0
#  endif
#  if defined(__AVX512PF__)
#    define OIIO_AVX512PF_ENABLED 1   /* Prefetch */
#  else
#    define OIIO_AVX512PF_ENABLED 0
#  endif
#  if defined(__AVX512ER__)
#    define OIIO_AVX512ER_ENABLED 1   /* Exponential & reciprocal */
#  else
#    define OIIO_AVX512ER_ENABLED 0
#  endif
#  if defined(__AVX512CD__)
#    define OIIO_AVX512CD_ENABLED 1   /* Conflict detection */
#  else
#    define OIIO_AVX512CD_ENABLED 0
#  endif
#  if defined(__AVX512BW__)
#    define OIIO_AVX512BW_ENABLED 1   /* Byte and word */
#  else
#    define OIIO_AVX512BW_ENABLED 0
#  endif
#  if defined(__AVX512VL__)
#    define OIIO_AVX512VL_ENABLED 1   /* Vector length extensions */
#  else
#    define OIIO_AVX512VL_ENABLED 0
#  endif
#else
#  define OIIO_SIMD_AVX 0
#  define OIIO_AVX512VL_ENABLED 0
#  define OIIO_AVX512DQ_ENABLED 0
#  define OIIO_AVX512PF_ENABLED 0
#  define OIIO_AVX512ER_ENABLED 0
#  define OIIO_AVX512CD_ENABLED 0
#  define OIIO_AVX512BW_ENABLED 0
#endif

#if defined(__FMA__)
#  define OIIO_FMA_ENABLED 1
#else
#  define OIIO_FMA_ENABLED 0
#endif
#if defined(__AVX512IFMA__)
#  define OIIO_AVX512IFMA_ENABLED 1
#else
#  define OIIO_AVX512IFMA_ENABLED 0
#endif

#if defined(__F16C__)
#  define OIIO_F16C_ENABLED 1
#else
#  define OIIO_F16C_ENABLED 0
#endif

// FIXME Future: support ARM Neon
// Uncomment this when somebody with Neon can verify it works
#if defined(__ARM_NEON__) && !defined(OIIO_NO_NEON)
#  define OIIO_SIMD 4
#  define OIIO_SIMD_NEON 1
#  define OIIO_SIMD_MAX_SIZE_BYTES 16
#  define OIIO_SIMD4_ALIGN OIIO_ALIGN(16)
#  define OIIO_SSE_ALIGN OIIO_ALIGN(16)
#else
#  define OIIO_SIMD_NEON 0
#endif

#ifndef OIIO_SIMD
   // No SIMD available
#  define OIIO_SIMD 0
#  define OIIO_SIMD4_ALIGN
#  define OIIO_SIMD_MAX_SIZE_BYTES 16
#endif

#ifndef OIIO_SIMD8_ALIGN
#  define OIIO_SIMD8_ALIGN OIIO_SIMD4_ALIGN
#endif
#ifndef OIIO_SIMD16_ALIGN
#  define OIIO_SIMD16_ALIGN OIIO_SIMD8_ALIGN
#endif


// General features that client apps may want to test for, for conditional
// compilation. Will add to this over time as needed. Note that just
// because a feature is present doesn't mean it's fast -- HAS_SIMD8 means
// the vfloat8 class (and friends) are in this version of simd.h, but that's
// different from OIIO_SIMD >= 8, which means it's supported in hardware.
#define OIIO_SIMD_HAS_MATRIX4 1  /* matrix44 defined */
#define OIIO_SIMD_HAS_FLOAT8 1   /* DEPRECATED(1.8) */
#define OIIO_SIMD_HAS_SIMD8 1    /* vfloat8, vint8, vbool8 defined */
#define OIIO_SIMD_HAS_SIMD16 1   /* vfloat16, vint16, vbool16 defined */


// Embarrassing hack: Xlib.h #define's True and False!
#ifdef True
#    undef True
#endif
#ifdef False
#    undef False
#endif



OIIO_NAMESPACE_BEGIN

namespace simd {

//////////////////////////////////////////////////////////////////////////
// Forward declarations of our main SIMD classes

class vbool4;
class vint4;
class vfloat4;
class vfloat3;
class matrix44;
class vbool8;
class vint8;
class vfloat8;
class vbool16;
class vint16;
class vfloat16;

// Deprecated names -- remove these in 1.9
typedef vbool4 mask4;    // old name
typedef vbool4 bool4;
typedef vbool8 bool8;
typedef vint4 int4;
typedef vint8 int8;
typedef vfloat3 float3;
typedef vfloat4 float4;
typedef vfloat8 float8;



//////////////////////////////////////////////////////////////////////////
// Template magic to determine the raw SIMD types involved, and other
// things helpful for metaprogramming.

template <typename T, int N> struct simd_raw_t { struct type { T val[N]; }; };
template <int N> struct simd_bool_t { struct type { int val[N]; }; };

#if OIIO_SIMD_SSE
template<> struct simd_raw_t<int,4> { typedef __m128i type; };
template<> struct simd_raw_t<float,4> { typedef __m128 type; };
template<> struct simd_bool_t<4> { typedef __m128 type; };
#endif

#if OIIO_SIMD_AVX
template<> struct simd_raw_t<int,8> { typedef __m256i type; };
template<> struct simd_raw_t<float,8> { typedef __m256 type; };
template<> struct simd_bool_t<8> { typedef __m256 type; };
#endif

#if OIIO_SIMD_AVX >= 512
template<> struct simd_raw_t<int,16> { typedef __m512i type; };
template<> struct simd_raw_t<float,16> { typedef __m512 type; };
template<> struct simd_bool_t<16> { typedef __mmask16 type; };
#else
// Note: change in strategy for 16-wide SIMD: instead of int[16] for
// vbool16, it's just a plain old bitmask, and __mask16 for actual HW.
template<> struct simd_bool_t<16> { typedef uint16_t type; };
#endif

#if OIIO_SIMD_NEON
template<> struct simd_raw_t<int,4> { typedef int32x4_t type; };
template<> struct simd_raw_t<float,4> { typedef float32x4_t type; };
template<> struct simd_bool_t<4> { typedef uint32x4_t type; };
#endif


/// Template to retrieve the vector type from the scalar. For example,
/// simd::VecType<int,4> will be vfloat4.
template<typename T,int elements> struct VecType {};
template<> struct VecType<int,1>   { typedef int type; };
template<> struct VecType<float,1> { typedef float type; };
template<> struct VecType<int,4>   { typedef vint4 type; };
template<> struct VecType<float,4>   { typedef vfloat4 type; };
template<> struct VecType<float,3> { typedef vfloat3 type; };
template<> struct VecType<bool,4>  { typedef vbool4 type; };
template<> struct VecType<int,8>   { typedef vint8 type; };
template<> struct VecType<float,8> { typedef vfloat8 type; };
template<> struct VecType<bool,8>  { typedef vbool8 type; };
template<> struct VecType<int,16>   { typedef vint16 type; };
template<> struct VecType<float,16> { typedef vfloat16 type; };
template<> struct VecType<bool,16>  { typedef vbool16 type; };

/// Template to retrieve the SIMD size of a SIMD type. Rigged to be 1 for
/// anything but our SIMD types.
template<typename T> struct SimdSize { static const int size = 1; };
template<> struct SimdSize<vint4>     { static const int size = 4; };
template<> struct SimdSize<vfloat4>   { static const int size = 4; };
template<> struct SimdSize<vfloat3>   { static const int size = 4; };
template<> struct SimdSize<vbool4>    { static const int size = 4; };
template<> struct SimdSize<vint8>     { static const int size = 8; };
template<> struct SimdSize<vfloat8>   { static const int size = 8; };
template<> struct SimdSize<vbool8>    { static const int size = 8; };
template<> struct SimdSize<vint16>    { static const int size = 16; };
template<> struct SimdSize<vfloat16>  { static const int size = 16; };
template<> struct SimdSize<vbool16>   { static const int size = 16; };

/// Template to retrieve the number of elements size of a SIMD type. Rigged
/// to be 1 for anything but our SIMD types.
template<typename T> struct SimdElements { static const int size = SimdSize<T>::size; };
template<> struct SimdElements<vfloat3>   { static const int size = 3; };

/// Template giving a printable name for each type
template<typename T> struct SimdTypeName { static const char *name() { return "unknown"; } };
template<> struct SimdTypeName<vfloat4>   { static const char *name() { return "vfloat4"; } };
template<> struct SimdTypeName<vint4>     { static const char *name() { return "vint4"; } };
template<> struct SimdTypeName<vbool4>    { static const char *name() { return "vbool4"; } };
template<> struct SimdTypeName<vfloat8>   { static const char *name() { return "vfloat8"; } };
template<> struct SimdTypeName<vint8>     { static const char *name() { return "vint8"; } };
template<> struct SimdTypeName<vbool8>    { static const char *name() { return "vbool8"; } };
template<> struct SimdTypeName<vfloat16>  { static const char *name() { return "vfloat16"; } };
template<> struct SimdTypeName<vint16>    { static const char *name() { return "vint16"; } };
template<> struct SimdTypeName<vbool16>   { static const char *name() { return "vbool16"; } };


//////////////////////////////////////////////////////////////////////////
// Macros helpful for making static constants in code.

# define OIIO_SIMD_FLOAT4_CONST(name,val) \
    static const OIIO_SIMD4_ALIGN float name[4] = { (val), (val), (val), (val) }
# define OIIO_SIMD_FLOAT4_CONST4(name,v0,v1,v2,v3) \
    static const OIIO_SIMD4_ALIGN float name[4] = { (v0), (v1), (v2), (v3) }
# define OIIO_SIMD_INT4_CONST(name,val) \
    static const OIIO_SIMD4_ALIGN int name[4] = { (val), (val), (val), (val) }
# define OIIO_SIMD_INT4_CONST4(name,v0,v1,v2,v3) \
    static const OIIO_SIMD4_ALIGN int name[4] = { (v0), (v1), (v2), (v3) }
# define OIIO_SIMD_UINT4_CONST(name,val) \
    static const OIIO_SIMD4_ALIGN uint32_t name[4] = { (val), (val), (val), (val) }
# define OIIO_SIMD_UINT4_CONST4(name,v0,v1,v2,v3) \
    static const OIIO_SIMD4_ALIGN uint32_t name[4] = { (v0), (v1), (v2), (v3) }

# define OIIO_SIMD_FLOAT8_CONST(name,val) \
    static const OIIO_SIMD8_ALIGN float name[8] = { (val), (val), (val), (val), \
                                                    (val), (val), (val), (val) }
# define OIIO_SIMD_FLOAT8_CONST8(name,v0,v1,v2,v3,v4,v5,v6,v7) \
    static const OIIO_SIMD8_ALIGN float name[8] = { (v0), (v1), (v2), (v3), \
                                                    (v4), (v5), (v6), (v7) }
# define OIIO_SIMD_INT8_CONST(name,val) \
    static const OIIO_SIMD8_ALIGN int name[8] = { (val), (val), (val), (val), \
                                                  (val), (val), (val), (val) }
# define OIIO_SIMD_INT8_CONST8(name,v0,v1,v2,v3,v4,v5,v6,v7) \
    static const OIIO_SIMD8_ALIGN int name[8] = { (v0), (v1), (v2), (v3), \
                                                  (v4), (v5), (v6), (v7) }
# define OIIO_SIMD_UINT8_CONST(name,val) \
    static const OIIO_SIMD8_ALIGN uint32_t name[8] = { (val), (val), (val), (val), \
                                                       (val), (val), (val), (val) }
# define OIIO_SIMD_UINT8_CONST8(name,v0,v1,v2,v3,v4,v5,v6,v7) \
    static const OIIO_SIMD8_ALIGN uint32_t name[8] = { (v0), (v1), (v2), (v3), \
                                                       (v4), (v5), (v6), (v7) }

# define OIIO_SIMD_VFLOAT16_CONST(name,val) \
    static const OIIO_SIMD16_ALIGN float name[16] = {           \
        (val), (val), (val), (val), (val), (val), (val), (val), \
        (val), (val), (val), (val), (val), (val), (val), (val) }
# define OIIO_SIMD_VFLOAT16_CONST16(name,v0,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15) \
    static const OIIO_SIMD16_ALIGN float name[16] = {           \
        (v0), (v1), (v2), (v3), (v4), (v5), (v6), (v7),         \
        (v8), (v9), (v10), (v11), (v12), (v13), (v14), (v15) }
# define OIIO_SIMD_INT16_CONST(name,val) \
    static const OIIO_SIMD16_ALIGN int name[16] = {             \
        (val), (val), (val), (val), (val), (val), (val), (val), \
        (val), (val), (val), (val), (val), (val), (val), (val) }
# define OIIO_SIMD_INT16_CONST16(name,v0,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15) \
    static const OIIO_SIMD16_ALIGN int name[16] = { \
        (v0), (v1), (v2), (v3), (v4), (v5), (v6), (v7), \
        (v8), (v9), (v10), (v11), (v12), (v13), (v14), (v15) }
# define OIIO_SIMD_UINT16_CONST(name,val) \
    static const OIIO_SIMD16_ALIGN uint32_t name[16] = {        \
        (val), (val), (val), (val), (val), (val), (val), (val), \
        (val), (val), (val), (val), (val), (val), (val), (val) }
# define OIIO_SIMD_UINT16_CONST16(name,v0,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15) \
    static const OIIO_SIMD16_ALIGN uint32_t name[16] = {        \
        (val), (val), (val), (val), (val), (val), (val), (val), \
        (val), (val), (val), (val), (val), (val), (val), (val) }


//////////////////////////////////////////////////////////////////////////
// Some macros just for use in this file (#undef-ed at the end) making
// it more succinct to express per-element operations.

#define SIMD_DO(x) for (int i = 0; i < elements; ++i) x
#define SIMD_CONSTRUCT(x) for (int i = 0; i < elements; ++i) m_val[i] = (x)
#define SIMD_CONSTRUCT_PAD(x) for (int i = 0; i < elements; ++i) m_val[i] = (x); \
                              for (int i = elements; i < paddedelements; ++i) m_val[i] = 0
#define SIMD_RETURN(T,x) T r; for (int i = 0; i < r.elements; ++i) r[i] = (x); return r
#define SIMD_RETURN_REDUCE(T,init,op) T r = init; for (int i = 0; i < v.elements; ++i) op; return r



//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
// The public declarations of the main SIMD classes follow: boolN, intN,
// floatN, matrix44.
//
// These class declarations are intended to be brief and self-documenting,
// and give all the information that users or client applications need to
// know to use these classes.
//
// No implementations are given inline except for the briefest, completely
// generic methods that don't have any architecture-specific overloads.
// After the class defintions, there will be an immense pile of full
// implementation definitions, which casual users are not expected to
// understand.
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////


/// vbool4: An 4-vector whose elements act mostly like bools, accelerated by
/// SIMD instructions when available. This is what is naturally produced by
/// SIMD comparison operators on the vfloat4 and vint4 types.
class vbool4 {
public:
    static const char* type_name() { return "vbool4"; }
    typedef bool value_t;        ///< Underlying equivalent scalar value type
    enum { elements = 4 };       ///< Number of scalar elements
    enum { paddedelements = 4 }; ///< Number of scalar elements for full pad
    enum { bits = elements*32 }; ///< Total number of bits
    typedef simd_bool_t<4>::type simd_t;  ///< the native SIMD type used

    /// Default constructor (contents undefined)
    vbool4 () { }

    /// Construct from a single value (store it in all slots)
    vbool4 (bool a) { load(a); }

    explicit vbool4 (const bool *a);

    /// Construct from 4 bool values
    vbool4 (bool a, bool b, bool c, bool d) { load (a, b, c, d); }

    /// Copy construct from another vbool4
    vbool4 (const vbool4 &other) { m_simd = other.m_simd; }

    /// Construct from 4 int values
    vbool4 (int a, int b, int c, int d) {
        load (bool(a), bool(b), bool(c), bool(d));
    }

    /// Construct from a SIMD int (is each element nonzero?)
    vbool4 (const vint4 &i);

    /// Construct from the underlying SIMD type
    vbool4 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Extract the bitmask
    int bitmask () const;

    /// Convert from integer bitmask to a true vbool4
    static vbool4 from_bitmask (int bitmask);

    /// Set all components to false
    void clear ();

    /// Return a vbool4 the is 'false' for all values
    static const vbool4 False ();

    /// Return a vbool4 the is 'true' for all values
    static const vbool4 True ();

    /// Assign one value to all components
    const vbool4 & operator= (bool a) { load(a); return *this; }

    /// Assignment of another vbool4
    const vbool4 & operator= (const vbool4 & other);

    /// Component access (get)
    int operator[] (int i) const;

    /// Component access (set).
    void setcomp (int i, bool value);

    /// Component access (set).
    /// NOTE: avoid this unsafe construct. It will go away some day.
    int& operator[] (int i);

    /// Helper: load a single value into all components.
    void load (bool a);

    /// Helper: load separate values into each component.
    void load (bool a, bool b, bool c, bool d);

    /// Helper: store the values into memory as bools.
    void store (bool *values) const;

    /// Store the first n values into memory.
    void store (bool *values, int n) const;

    /// Logical/bitwise operators, component-by-component
    friend vbool4 operator! (const vbool4& a);
    friend vbool4 operator& (const vbool4& a, const vbool4& b);
    friend vbool4 operator| (const vbool4& a, const vbool4& b);
    friend vbool4 operator^ (const vbool4& a, const vbool4& b);
    friend vbool4 operator~ (const vbool4& a);
    friend const vbool4& operator&= (vbool4& a, const vbool4& b);
    friend const vbool4& operator|= (vbool4& a, const vbool4& b);
    friend const vbool4& operator^= (vbool4& a, const vbool4& b);

    /// Comparison operators, component by component
    friend vbool4 operator== (const vbool4& a, const vbool4& b);
    friend vbool4 operator!= (const vbool4& a, const vbool4& b);

    /// Stream output
    friend std::ostream& operator<< (std::ostream& cout, const vbool4 & a);

private:
    // The actual data representation
    union {
        simd_t m_simd;
        int m_val[paddedelements];
    };
};



/// Helper: shuffle/swizzle with constant (templated) indices.
/// Example: shuffle<1,1,2,2>(vbool4(a,b,c,d)) returns (b,b,c,c)
template<int i0, int i1, int i2, int i3> vbool4 shuffle (const vbool4& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> vbool4 shuffle (const vbool4& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> bool extract (const vbool4& a);

/// Helper: substitute val for a[i]
template<int i> vbool4 insert (const vbool4& a, bool val);

/// Logical reduction across all components.
bool reduce_and (const vbool4& v);
bool reduce_or (const vbool4& v);

// Are all/any/no components true?
bool all (const vbool4& v);
bool any (const vbool4& v);
bool none (const vbool4& v);

// It's handy to have this defined for regular bool as well
inline bool all (bool v) { return v; }



/// vbool8: An 8-vector whose elements act mostly like bools, accelerated by
/// SIMD instructions when available. This is what is naturally produced by
/// SIMD comparison operators on the vfloat8 and vint8 types.
class vbool8 {
public:
    static const char* type_name() { return "vbool8"; }
    typedef bool value_t;        ///< Underlying equivalent scalar value type
    enum { elements = 8 };       ///< Number of scalar elements
    enum { paddedelements = 8 }; ///< Number of scalar elements for full pad
    enum { bits = elements*32 }; ///< Total number of bits
    typedef simd_bool_t<8>::type simd_t;  ///< the native SIMD type used

    /// Default constructor (contents undefined)
    vbool8 () { }

    /// Construct from a single value (store it in all slots)
    vbool8 (bool a) { load (a); }

    explicit vbool8 (const bool *values);

    /// Construct from 8 bool values
    vbool8 (bool a, bool b, bool c, bool d, bool e, bool f, bool g, bool h);

    /// Copy construct from another vbool8
    vbool8 (const vbool8 &other) { m_simd = other.m_simd; }

    /// Construct from 8 int values
    vbool8 (int a, int b, int c, int d, int e, int f, int g, int h);

    /// Construct from a SIMD int (is each element nonzero?)
    vbool8 (const vint8 &i);

    /// Construct from two vbool4's
    vbool8 (const vbool4 &lo, const vbool4 &hi);

    /// Construct from the underlying SIMD type
    vbool8 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Extract the bitmask
    int bitmask () const;

    /// Convert from integer bitmask to a true vbool8
    static vbool8 from_bitmask (int bitmask);

    /// Set all components to false
    void clear ();

    /// Return a vbool8 the is 'false' for all values
    static const vbool8 False ();

    /// Return a vbool8 the is 'true' for all values
    static const vbool8 True ();

    /// Assign one value to all components
    const vbool8 & operator= (bool a);

    /// Assignment of another vbool8
    const vbool8 & operator= (const vbool8 & other);

    /// Component access (get)
    int operator[] (int i) const;

    /// Component access (set).
    void setcomp (int i, bool value);

    /// Component access (set).
    /// NOTE: avoid this unsafe construct. It will go away some day.
    int& operator[] (int i);

    /// Extract the lower precision vbool4
    vbool4 lo () const;

    /// Extract the higher precision vbool4
    vbool4 hi () const;

    /// Helper: load a single value into all components.
    void load (bool a);

    /// Helper: load separate values into each component.
    void load (bool a, bool b, bool c, bool d,
               bool e, bool f, bool g, bool h);

    /// Helper: store the values into memory as bools.
    void store (bool *values) const;

    /// Store the first n values into memory.
    void store (bool *values, int n) const;

    /// Logical/bitwise operators, component-by-component
    friend vbool8 operator! (const vbool8& a);
    friend vbool8 operator& (const vbool8& a, const vbool8& b);
    friend vbool8 operator| (const vbool8& a, const vbool8& b);
    friend vbool8 operator^ (const vbool8& a, const vbool8& b);
    friend vbool8 operator~ (const vbool8& a);
    friend const vbool8& operator&= (vbool8& a, const vbool8& b);
    friend const vbool8& operator|= (vbool8& a, const vbool8& b);
    friend const vbool8& operator^= (vbool8& a, const vbool8& b);

    /// Comparison operators, component by component
    friend vbool8 operator== (const vbool8& a, const vbool8& b);
    friend vbool8 operator!= (const vbool8& a, const vbool8& b);

    /// Stream output
    friend std::ostream& operator<< (std::ostream& cout, const vbool8 & a);

private:
    // The actual data representation
    union {
        simd_t m_simd;
        int m_val[paddedelements];
        vbool4 m_4[2];
    };
};



/// Helper: shuffle/swizzle with constant (templated) indices.
/// Example: shuffle<1,1,2,2>(vbool4(a,b,c,d)) returns (b,b,c,c)
template<int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7>
vbool8 shuffle (const vbool8& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> vbool8 shuffle (const vbool8& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> bool extract (const vbool8& a);

/// Helper: substitute val for a[i]
template<int i> vbool8 insert (const vbool8& a, bool val);

/// Logical reduction across all components.
bool reduce_and (const vbool8& v);
bool reduce_or (const vbool8& v);

// Are all/any/no components true?
bool all (const vbool8& v);
bool any (const vbool8& v);
bool none (const vbool8& v);




/// vbool16: An 16-vector whose elements act mostly like bools, accelerated
/// by SIMD instructions when available. This is what is naturally produced
/// by SIMD comparison operators on the vfloat16 and vint16 types.
class vbool16 {
public:
    static const char* type_name() { return "vbool16"; }
    typedef bool value_t;        ///< Underlying equivalent scalar value type
    enum { elements = 16 };       ///< Number of scalar elements
    enum { paddedelements = 16 }; ///< Number of scalar elements for full pad
    enum { bits = 16 };           ///< Total number of bits
    typedef simd_bool_t<16>::type simd_t;  ///< the native SIMD type used

    /// Default constructor (contents undefined)
    vbool16 () { }

    /// Construct from a single value (store it in all slots)
    vbool16 (bool a) { load (a); }

    explicit vbool16 (int bitmask) { load_bitmask (bitmask); }

    explicit vbool16 (const bool *values);

    /// Construct from 16 bool values
    vbool16 (bool v0, bool v1, bool v2, bool v3, bool v4, bool v5, bool v6, bool v7,
            bool v8, bool v9, bool v10, bool v11, bool v12, bool v13, bool v14, bool v15);

    /// Copy construct from another vbool16
    vbool16 (const vbool16 &other) { m_simd = other.m_simd; }

    /// Construct from 16 int values
    vbool16 (int v0, int v1, int v2, int v3, int v4, int v5, int v6, int v7,
             int v8, int v9, int v10, int v11, int v12, int v13, int v14, int v15);

    /// Construct from a SIMD int (is each element nonzero?)
    vbool16 (const vint16 &i);

    /// Construct from two vbool8's
    vbool16 (const vbool8 &lo, const vbool8 &hi);

    /// Construct from four vbool4's
    vbool16 (const vbool4 &b4a, const vbool4 &b4b, const vbool4 &b4c, const vbool4 &b4d);

    /// Construct from the underlying SIMD type
    vbool16 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    int bitmask () const;

    /// Convert from integer bitmask to a true vbool16
    static vbool16 from_bitmask (int bitmask) { return vbool16(bitmask); }

    /// Set all components to false
    void clear ();

    /// Return a vbool16 the is 'false' for all values
    static const vbool16 False ();

    /// Return a vbool16 the is 'true' for all values
    static const vbool16 True ();

    /// Assign one value to all components
    const vbool16 & operator= (bool a);

    /// Assignment of another vbool16
    const vbool16 & operator= (const vbool16 & other);

    /// Component access (get)
    int operator[] (int i) const;

    /// Component access (set).
    void setcomp (int i, bool value);

    /// Extract the lower precision vbool8
    vbool8 lo () const;

    /// Extract the higher precision vbool8
    vbool8 hi () const;

    /// Helper: load a single value into all components.
    void load (bool a);

    /// Helper: load separate values into each component.
    void load (bool v0, bool v1, bool v2, bool v3, bool v4, bool v5, bool v6, bool v7,
               bool v8, bool v9, bool v10, bool v11, bool v12, bool v13, bool v14, bool v15);

    /// Helper: load all components from a bitmask in an int.
    void load_bitmask (int a);

    /// Helper: store the values into memory as bools.
    void store (bool *values) const;

    /// Store the first n values into memory.
    void store (bool *values, int n) const;

    /// Logical/bitwise operators, component-by-component
    friend vbool4 operator! (const vbool4& a);
    friend vbool16 operator! (const vbool16& a);
    friend vbool16 operator& (const vbool16& a, const vbool16& b);
    friend vbool16 operator| (const vbool16& a, const vbool16& b);
    friend vbool16 operator^ (const vbool16& a, const vbool16& b);
    friend vbool16 operator~ (const vbool16& a);
    friend const vbool16& operator&= (vbool16& a, const vbool16& b);
    friend const vbool16& operator|= (vbool16& a, const vbool16& b);
    friend const vbool16& operator^= (vbool16& a, const vbool16& b);

    /// Comparison operators, component by component
    friend vbool16 operator== (const vbool16& a, const vbool16& b);
    friend vbool16 operator!= (const vbool16& a, const vbool16& b);

    /// Stream output
    friend std::ostream& operator<< (std::ostream& cout, const vbool16 & a);

private:
    // The actual data representation
    union {
        simd_t   m_simd;
        uint16_t m_bits;
    };
};



/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> bool extract (const vbool16& a);

/// Helper: substitute val for a[i]
template<int i> vbool16 insert (const vbool16& a, bool val);

/// Logical reduction across all components.
bool reduce_and (const vbool16& v);
bool reduce_or (const vbool16& v);

// Are all/any/no components true?
bool all (const vbool16& v);
bool any (const vbool16& v);
bool none (const vbool16& v);





/// Integer 4-vector, accelerated by SIMD instructions when available.
class vint4 {
public:
    static const char* type_name() { return "vint4"; }
    typedef int value_t;      ///< Underlying equivalent scalar value type
    enum { elements = 4 };    ///< Number of scalar elements
    enum { paddedelements =4 }; ///< Number of scalar elements for full pad
    enum { bits = 128 };      ///< Total number of bits
    typedef simd_raw_t<int,elements>::type simd_t;  ///< the native SIMD type used
    typedef vbool4 vbool_t;   ///< bool type of the same length
    typedef vfloat4 vfloat_t; ///< float type of the same length
    typedef vint4 vint_t;     ///< int type of the same length
    typedef vbool4 bool_t;   // old name (deprecated 1.8)
    typedef vfloat4 float_t; // old name (deprecated 1.8)

    /// Default constructor (contents undefined)
    vint4 () { }

    /// Construct from a single value (store it in all slots)
    vint4 (int a);

    /// Construct from 2 values -- (a,a,b,b)
    vint4 (int a, int b);

    /// Construct from 4 values
    vint4 (int a, int b, int c, int d);

    /// Construct from a pointer to values
    vint4 (const int *vals);

    /// Construct from a pointer to unsigned short values
    explicit vint4 (const unsigned short *vals);

    /// Construct from a pointer to signed short values
    explicit vint4 (const short *vals);

    /// Construct from a pointer to unsigned char values (0 - 255)
    explicit vint4 (const unsigned char *vals);

    /// Construct from a pointer to signed char values (-128 - 127)
    explicit vint4 (const char *vals);

    /// Copy construct from another vint4
    vint4 (const vint4 & other) { m_simd = other.m_simd; }

    /// Convert a vfloat to an vint. Equivalent to i = (int)f;
    explicit vint4 (const vfloat4& f); // implementation below

    /// Construct from the underlying SIMD type
    vint4 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Return a pointer to the underlying scalar type
    const value_t* data () const { return (const value_t*)this; }
    value_t* data () { return (value_t*)this; }

    /// Sset all components to 0
    void clear () ;

    /// Return an vint4 with all components set to 0
    static const vint4 Zero ();

    /// Return an vint4 with all components set to 1
    static const vint4 One ();

    /// Return an vint4 with all components set to -1 (aka 0xffffffff)
    static const vint4 NegOne ();

    /// Return an vint4 with incremented components (e.g., 0,1,2,3).
    /// Optional arguments can give a non-zero starting point and step size.
    static const vint4 Iota (int start=0, int step=1);

    /// Return an vint4 with "geometric" iota: (1, 2, 4, 8).
    static const vint4 Giota ();

    /// Assign one value to all components.
    const vint4 & operator= (int a);

    /// Assignment from another vint4
    const vint4 & operator= (const vint4& other) ;

    /// Component access (get)
    int operator[] (int i) const;

    /// Component access (set)
    int& operator[] (int i);

    /// Component access (set).
    void setcomp (int i, int value);

    value_t x () const;
    value_t y () const;
    value_t z () const;
    value_t w () const;
    void set_x (value_t val);
    void set_y (value_t val);
    void set_z (value_t val);
    void set_w (value_t val);

    /// Helper: load a single int into all components
    void load (int a);

    /// Helper: load separate values into each component.
    void load (int a, int b, int c, int d);

    /// Load from an array of 4 values
    void load (const int *values);

    void load (const int *values, int n) ;

    /// Load from an array of 4 unsigned short values, convert to vint4
    void load (const unsigned short *values) ;

    /// Load from an array of 4 unsigned short values, convert to vint4
    void load (const short *values);

    /// Load from an array of 4 unsigned char values, convert to vint4
    void load (const unsigned char *values);

    /// Load from an array of 4 unsigned char values, convert to vint4
    void load (const char *values);

    /// Store the values into memory
    void store (int *values) const;

    /// Store the first n values into memory
    void store (int *values, int n) const;

    /// Store the least significant 16 bits of each element into adjacent
    /// unsigned shorts.
    void store (unsigned short *values) const;

    /// Store the least significant 8 bits of each element into adjacent
    /// unsigned chars.
    void store (unsigned char *values) const;

    /// Masked load -- read from values[] where mask is 1, load zero where
    /// mask is 0.
    void load_mask (int mask, const value_t *values);
    void load_mask (const vbool_t& mask, const value_t *values);

    /// Masked store -- write to values[] where mask is enabled, don't
    /// touch values[] where it's not.
    void store_mask (int mask, value_t *values) const;
    void store_mask (const vbool_t& mask, value_t *values) const;

    /// Load values from addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void gather (const value_t *baseptr, const vint_t& vindex);
    /// Gather elements defined by the mask, leave others unchanged.
    template<int scale=4>
    void gather_mask (const bool_t& mask, const value_t *baseptr, const vint_t& vindex);
    template<int scale=4>
    void gather_mask (int mask, const value_t *baseptr, const vint_t& vindex);

    /// Store values at addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void scatter (value_t *baseptr, const vint_t& vindex) const;
    /// Scatter elements defined by the mask
    template<int scale=4>
    void scatter_mask (const bool_t& mask, value_t *baseptr, const vint_t& vindex) const;
    template<int scale=4>
    void scatter_mask (int mask, value_t *baseptr, const vint_t& vindex) const;

    // Arithmetic operators (component-by-component)
    friend vint4 operator+ (const vint4& a, const vint4& b);
    friend vint4 operator- (const vint4& a);
    friend vint4 operator- (const vint4& a, const vint4& b);
    friend vint4 operator* (const vint4& a, const vint4& b);
    friend vint4 operator/ (const vint4& a, const vint4& b);
    friend vint4 operator% (const vint4& a, const vint4& b);
    friend const vint4 & operator+= (vint4& a, const vint4& b);
    friend const vint4 & operator-= (vint4& a, const vint4& b);
    friend const vint4 & operator*= (vint4& a, const vint4& b);
    friend const vint4 & operator/= (vint4& a, const vint4& b);
    friend const vint4 & operator%= (vint4& a, const vint4& b);
    // Bitwise operators (component-by-component)
    friend vint4 operator& (const vint4& a, const vint4& b);
    friend vint4 operator| (const vint4& a, const vint4& b);
    friend vint4 operator^ (const vint4& a, const vint4& b);
    friend const vint4& operator&= (vint4& a, const vint4& b);
    friend const vint4& operator|= (vint4& a, const vint4& b);
    friend const vint4& operator^= (vint4& a, const vint4& b);
    friend vint4 operator~ (const vint4& a);
    friend vint4 operator<< (const vint4& a, unsigned int bits);
    friend vint4 operator>> (const vint4& a, unsigned int bits);
    friend const vint4& operator<<= (vint4& a, unsigned int bits);
    friend const vint4& operator>>= (vint4& a, unsigned int bits);
    // Comparison operators (component-by-component)
    friend vbool4 operator== (const vint4& a, const vint4& b);
    friend vbool4 operator!= (const vint4& a, const vint4& b);
    friend vbool4 operator<  (const vint4& a, const vint4& b);
    friend vbool4 operator>  (const vint4& a, const vint4& b);
    friend vbool4 operator>= (const vint4& a, const vint4& b);
    friend vbool4 operator<= (const vint4& a, const vint4& b);

    /// Stream output
    friend std::ostream& operator<< (std::ostream& cout, const vint4 & a);

private:
    // The actual data representation
    union {
        simd_t  m_simd;
        value_t m_val[elements];
    };
};



// Shift right logical -- unsigned shift. This differs from operator>>
// in how it handles the sign bit.  (1<<31) >> 1 == (1<<31), but
// srl((1<<31),1) == 1<<30.
vint4 srl (const vint4& val, const unsigned int bits);

/// Helper: shuffle/swizzle with constant (templated) indices.
/// Example: shuffle<1,1,2,2>(vbool4(a,b,c,d)) returns (b,b,c,c)
template<int i0, int i1, int i2, int i3> vint4 shuffle (const vint4& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> vint4 shuffle (const vint4& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> int extract (const vint4& v);

/// The sum of all components, returned in all components.
vint4 vreduce_add (const vint4& v);

// Reduction across all components
int reduce_add (const vint4& v);
int reduce_and (const vint4& v);
int reduce_or (const vint4& v);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// and b (if mask[i] is true), i.e., mask[i] ? b[i] : a[i].
vint4 blend (const vint4& a, const vint4& b, const vbool4& mask);

/// Use a bool mask to select between `a` (if mask[i] is true) or 0 if
/// mask[i] is false), i.e., mask[i] ? a[i] : 0. Equivalent to
/// blend(0,a,mask).
vint4 blend0 (const vint4& a, const vbool4& mask);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// or 0 (if mask[i] is true), i.e., mask[i] ? 0 : a[i]. Equivalent to
/// blend(0,a,!mask), or blend(a,0,mask).
vint4 blend0not (const vint4& a, const vbool4& mask);

/// Select 'a' where mask is true, 'b' where mask is false. Sure, it's a
/// synonym for blend with arguments rearranged, but this is more clear
/// because the arguments are symmetric to scalar (cond ? a : b).
vint4 select (const vbool4& mask, const vint4& a, const vint4& b);

// Per-element math
vint4 abs (const vint4& a);
vint4 min (const vint4& a, const vint4& b);
vint4 max (const vint4& a, const vint4& b);

/// Circular bit rotate by s bits, for N values at once.
vint4 rotl (const vint4& x, const int s);
// DEPRECATED(2.1)
vint4 rotl32 (const vint4& x, const unsigned int k);

/// andnot(a,b) returns ((~a) & b)
vint4 andnot (const vint4& a, const vint4& b);

/// Bitcast back and forth to intN (not a convert -- move the bits!)
vint4 bitcast_to_int (const vbool4& x);
vint4 bitcast_to_int (const vfloat4& x);
vfloat4 bitcast_to_float (const vint4& x);

void transpose (vint4 &a, vint4 &b, vint4 &c, vint4 &d);
void transpose (const vint4& a, const vint4& b, const vint4& c, const vint4& d,
                vint4 &r0, vint4 &r1, vint4 &r2, vint4 &r3);

vint4 AxBxCxDx (const vint4& a, const vint4& b, const vint4& c, const vint4& d);

// safe_mod(a,b) is like a%b, but safely returns 0 when b==0.
vint4 safe_mod (const vint4& a, const vint4& b);
vint4 safe_mod (const vint4& a, int b);




/// Integer 8-vector, accelerated by SIMD instructions when available.
class vint8 {
public:
    static const char* type_name() { return "vint8"; }
    typedef int value_t;      ///< Underlying equivalent scalar value type
    enum { elements = 8 };    ///< Number of scalar elements
    enum { paddedelements =8 }; ///< Number of scalar elements for full pad
    enum { bits = elements*32 }; ///< Total number of bits
    typedef simd_raw_t<int,elements>::type simd_t;  ///< the native SIMD type used
    typedef vbool8 vbool_t;   ///< bool type of the same length
    typedef vfloat8 vfloat_t; ///< float type of the same length
    typedef vint8 vint_t;     ///< int type of the same length
    typedef vbool8 bool_t;   // old name (deprecated 1.8)
    typedef vfloat8 float_t; // old name (deprecated 1.8)

    /// Default constructor (contents undefined)
    vint8 () { }

    /// Construct from a single value (store it in all slots)
    vint8 (int a);

    /// Construct from 2 values -- (a,a,b,b)
    vint8 (int a, int b);

    /// Construct from 8 values (won't work for vint8)
    vint8 (int a, int b, int c, int d, int e, int f, int g, int h);

    /// Construct from a pointer to values
    vint8 (const int *vals);

    /// Construct from a pointer to unsigned short values
    explicit vint8 (const unsigned short *vals);

    /// Construct from a pointer to signed short values
    explicit vint8 (const short *vals);

    /// Construct from a pointer to unsigned char values (0 - 255)
    explicit vint8 (const unsigned char *vals);

    /// Construct from a pointer to signed char values (-128 - 127)
    explicit vint8 (const char *vals);

    /// Copy construct from another vint8
    vint8 (const vint8 & other) { m_simd = other.m_simd; }

    /// Convert a vfloat8 to an vint8. Equivalent to i = (int)f;
    explicit vint8 (const vfloat8& f); // implementation below

    /// Construct from two vint4's
    vint8 (const vint4 &lo, const vint4 &hi);

    /// Construct from the underlying SIMD type
    vint8 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Return a pointer to the underlying scalar type
    const value_t* data () const { return (const value_t*)this; }
    value_t* data () { return (value_t*)this; }

    /// Sset all components to 0
    void clear () ;

    /// Return an vint8 with all components set to 0
    static const vint8 Zero ();

    /// Return an vint8 with all components set to 1
    static const vint8 One ();

    /// Return an vint8 with all components set to -1 (aka 0xffffffff)
    static const vint8 NegOne ();

    /// Return an vint8 with incremented components (e.g., 0,1,2,3).
    /// Optional arguments can give a non-zero starting point and step size.
    static const vint8 Iota (int start=0, int step=1);

    /// Return an vint8 with "geometric" iota: (1, 2, 4, 8, ...).
    static const vint8 Giota ();

    /// Assign one value to all components.
    const vint8 & operator= (int a);

    /// Assignment from another vint8
    const vint8 & operator= (const vint8& other) ;

    /// Component access (get)
    int operator[] (int i) const;

    /// Component access (set)
    int& operator[] (int i);

    /// Component access (set).
    void setcomp (int i, int value);

    value_t x () const;
    value_t y () const;
    value_t z () const;
    value_t w () const;
    void set_x (value_t val);
    void set_y (value_t val);
    void set_z (value_t val);
    void set_w (value_t val);

    /// Extract the lower precision vint4
    vint4 lo () const;

    /// Extract the higher precision vint4
    vint4 hi () const;

    /// Helper: load a single int into all components
    void load (int a);

    /// Load separate values into each component.
    void load (int a, int b, int c, int d, int e, int f, int g, int h);

    /// Load from an array of 8 values
    void load (const int *values);

    void load (const int *values, int n) ;

    /// Load from an array of 8 unsigned short values, convert to vint8
    void load (const unsigned short *values) ;

    /// Load from an array of 8 unsigned short values, convert to vint8
    void load (const short *values);

    /// Load from an array of 8 unsigned char values, convert to vint8
    void load (const unsigned char *values);

    /// Load from an array of 8 unsigned char values, convert to vint8
    void load (const char *values);

    /// Store the values into memory
    void store (int *values) const;

    /// Store the first n values into memory
    void store (int *values, int n) const;

    /// Store the least significant 16 bits of each element into adjacent
    /// unsigned shorts.
    void store (unsigned short *values) const;

    /// Store the least significant 8 bits of each element into adjacent
    /// unsigned chars.
    void store (unsigned char *values) const;

    /// Masked load -- read from values[] where mask is 1, load zero where
    /// mask is 0.
    void load_mask (int mask, const value_t *values);
    void load_mask (const vbool_t& mask, const value_t *values);

    /// Masked store -- write to values[] where mask is enabled, don't
    /// touch values[] where it's not.
    void store_mask (int mask, value_t *values) const;
    void store_mask (const vbool_t& mask, value_t *values) const;

    /// Load values from addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void gather (const value_t *baseptr, const vint_t& vindex);
    /// Gather elements defined by the mask, leave others unchanged.
    template<int scale=4>
    void gather_mask (const bool_t& mask, const value_t *baseptr, const vint_t& vindex);
    template<int scale=4>
    void gather_mask (int mask, const value_t *baseptr, const vint_t& vindex);

    /// Store values at addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void scatter (value_t *baseptr, const vint_t& vindex) const;
    /// Scatter elements defined by the mask
    template<int scale=4>
    void scatter_mask (const bool_t& mask, value_t *baseptr, const vint_t& vindex) const;
    template<int scale=4>
    void scatter_mask (int mask, value_t *baseptr, const vint_t& vindex) const;

    // Arithmetic operators (component-by-component)
    friend vint8 operator+ (const vint8& a, const vint8& b);
    friend vint8 operator- (const vint8& a);
    friend vint8 operator- (const vint8& a, const vint8& b);
    friend vint8 operator* (const vint8& a, const vint8& b);
    friend vint8 operator/ (const vint8& a, const vint8& b);
    friend vint8 operator% (const vint8& a, const vint8& b);
    friend const vint8 & operator+= (vint8& a, const vint8& b);
    friend const vint8 & operator-= (vint8& a, const vint8& b);
    friend const vint8 & operator*= (vint8& a, const vint8& b);
    friend const vint8 & operator/= (vint8& a, const vint8& b);
    friend const vint8 & operator%= (vint8& a, const vint8& b);
    // Bitwise operators (component-by-component)
    friend vint8 operator& (const vint8& a, const vint8& b);
    friend vint8 operator| (const vint8& a, const vint8& b);
    friend vint8 operator^ (const vint8& a, const vint8& b);
    friend const vint8& operator&= (vint8& a, const vint8& b);
    friend const vint8& operator|= (vint8& a, const vint8& b);
    friend const vint8& operator^= (vint8& a, const vint8& b);
    friend vint8 operator~ (const vint8& a);
    friend vint8 operator<< (const vint8& a, unsigned int bits);
    friend vint8 operator>> (const vint8& a, unsigned int bits);
    friend const vint8& operator<<= (vint8& a, unsigned int bits);
    friend const vint8& operator>>= (vint8& a, unsigned int bits);
    // Comparison operators (component-by-component)
    friend vbool8 operator== (const vint8& a, const vint8& b);
    friend vbool8 operator!= (const vint8& a, const vint8& b);
    friend vbool8 operator<  (const vint8& a, const vint8& b);
    friend vbool8 operator>  (const vint8& a, const vint8& b);
    friend vbool8 operator>= (const vint8& a, const vint8& b);
    friend vbool8 operator<= (const vint8& a, const vint8& b);

    /// Stream output
    friend std::ostream& operator<< (std::ostream& cout, const vint8& a);

private:
    // The actual data representation
    union {
        simd_t  m_simd;
        value_t m_val[elements];
        vint4 m_4[2];
    };
};



// Shift right logical -- unsigned shift. This differs from operator>>
// in how it handles the sign bit.  (1<<31) >> 1 == (1<<31), but
// srl((1<<31),1) == 1<<30.
vint8 srl (const vint8& val, const unsigned int bits);

/// Helper: shuffle/swizzle with constant (templated) indices.
/// Example: shuffle<1,1,2,2>(vbool4(a,b,c,d)) returns (b,b,c,c)
template<int i0, int i1, int i2, int i3,
         int i4, int i5, int i6, int i7> vint8 shuffle (const vint8& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> vint8 shuffle (const vint8& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> int extract (const vint8& v);

/// Helper: substitute val for a[i]
template<int i> vint8 insert (const vint8& a, int val);

/// The sum of all components, returned in all components.
vint8 vreduce_add (const vint8& v);

// Reduction across all components
int reduce_add (const vint8& v);
int reduce_and (const vint8& v);
int reduce_or (const vint8& v);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// and b (if mask[i] is true), i.e., mask[i] ? b[i] : a[i].
vint8 blend (const vint8& a, const vint8& b, const vbool8& mask);

/// Use a bool mask to select between `a` (if mask[i] is true) or 0 if
/// mask[i] is false), i.e., mask[i] ? a[i] : 0. Equivalent to
/// blend(0,a,mask).
vint8 blend0 (const vint8& a, const vbool8& mask);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// or 0 (if mask[i] is true), i.e., mask[i] ? 0 : a[i]. Equivalent to
/// blend(0,a,!mask), or blend(a,0,mask).
vint8 blend0not (const vint8& a, const vbool8& mask);

/// Select 'a' where mask is true, 'b' where mask is false. Sure, it's a
/// synonym for blend with arguments rearranged, but this is more clear
/// because the arguments are symmetric to scalar (cond ? a : b).
vint8 select (const vbool8& mask, const vint8& a, const vint8& b);

// Per-element math
vint8 abs (const vint8& a);
vint8 min (const vint8& a, const vint8& b);
vint8 max (const vint8& a, const vint8& b);

/// Circular bit rotate by s bits, for N values at once.
vint8 rotl (const vint8& x, const int s);
// DEPRECATED(2.1)
vint8 rotl32 (const vint8& x, const unsigned int k);

/// andnot(a,b) returns ((~a) & b)
vint8 andnot (const vint8& a, const vint8& b);

/// Bitcast back and forth to intN (not a convert -- move the bits!)
vint8 bitcast_to_int (const vbool8& x);
vint8 bitcast_to_int (const vfloat8& x);
vfloat8 bitcast_to_float (const vint8& x);

// safe_mod(a,b) is like a%b, but safely returns 0 when b==0.
vint8 safe_mod (const vint8& a, const vint8& b);
vint8 safe_mod (const vint8& a, int b);





/// Integer 16-vector, accelerated by SIMD instructions when available.
class vint16 {
public:
    static const char* type_name() { return "vint16"; }
    typedef int value_t;      ///< Underlying equivalent scalar value type
    enum { elements = 16 };    ///< Number of scalar elements
    enum { paddedelements =16 }; ///< Number of scalar elements for full pad
    enum { bits = 128 };      ///< Total number of bits
    typedef simd_raw_t<int,elements>::type simd_t;  ///< the native SIMD type used
    typedef vbool16 vbool_t;   ///< bool type of the same length
    typedef vfloat16 vfloat_t; ///< float type of the same length
    typedef vint16 vint_t;     ///< int type of the same length
    typedef vbool16 bool_t;   // old name (deprecated 1.8)
    typedef vfloat16 float_t; // old name (deprecated 1.8)

    /// Default constructor (contents undefined)
    vint16 () { }

    /// Construct from a single value (store it in all slots)
    vint16 (int a);

    /// Construct from 16 values (won't work for vint16)
    vint16 (int v0, int v1, int v2, int v3, int v4, int v5, int v6, int v7,
           int v8, int v9, int v10, int v11, int v12, int v13, int v14, int v15);

    /// Construct from a pointer to values
    vint16 (const int *vals);

    /// Construct from a pointer to unsigned short values
    explicit vint16 (const unsigned short *vals);

    /// Construct from a pointer to signed short values
    explicit vint16 (const short *vals);

    /// Construct from a pointer to unsigned char values (0 - 255)
    explicit vint16 (const unsigned char *vals);

    /// Construct from a pointer to signed char values (-128 - 127)
    explicit vint16 (const char *vals);

    /// Copy construct from another vint16
    vint16 (const vint16 & other) { m_simd = other.m_simd; }

    /// Convert a vfloat16 to an vint16. Equivalent to i = (int)f;
    explicit vint16 (const vfloat16& f); // implementation below

    /// Construct from two vint8's
    vint16 (const vint8 &lo, const vint8 &hi);

    /// Construct from four vint4's
    vint16 (const vint4 &a, const vint4 &b, const vint4 &c, const vint4 &d);

    /// Construct from the underlying SIMD type
    vint16 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Return a pointer to the underlying scalar type
    const value_t* data () const { return (const value_t*)this; }
    value_t* data () { return (value_t*)this; }

    /// Sset all components to 0
    void clear () ;

    /// Return an vint16 with all components set to 0
    static const vint16 Zero ();

    /// Return an vint16 with all components set to 1
    static const vint16 One ();

    /// Return an vint16 with all components set to -1 (aka 0xffffffff)
    static const vint16 NegOne ();

    /// Return an vint16 with incremented components (e.g., 0,1,2,3).
    /// Optional arguments can give a non-zero starting point and step size.
    static const vint16 Iota (int start=0, int step=1);

    /// Return an vint16 with "geometric" iota: (1, 2, 4, 8, ...).
    static const vint16 Giota ();

    /// Assign one value to all components.
    const vint16 & operator= (int a);

    /// Assignment from another vint16
    const vint16 & operator= (const vint16& other) ;

    /// Component access (get)
    int operator[] (int i) const;

    /// Component access (set)
    int& operator[] (int i);

    /// Component access (set).
    void setcomp (int i, int value);

    value_t x () const;
    value_t y () const;
    value_t z () const;
    value_t w () const;
    void set_x (value_t val);
    void set_y (value_t val);
    void set_z (value_t val);
    void set_w (value_t val);

    /// Extract the lower precision vint8
    vint8 lo () const;

    /// Extract the higher precision vint8
    vint8 hi () const;

    /// Helper: load a single int into all components
    void load (int a);

    /// Load separate values into each component.
    void load (int v0, int v1, int v2, int v3, int v4, int v5, int v6, int v7,
               int v8, int v9, int v10, int v11, int v12, int v13, int v14, int v15);

    /// Load from an array of 16 values
    void load (const int *values);

    void load (const int *values, int n) ;

    /// Load from an array of 16 unsigned short values, convert to vint16
    void load (const unsigned short *values) ;

    /// Load from an array of 16 unsigned short values, convert to vint16
    void load (const short *values);

    /// Load from an array of 16 unsigned char values, convert to vint16
    void load (const unsigned char *values);

    /// Load from an array of 16 unsigned char values, convert to vint16
    void load (const char *values);

    /// Store the values into memory
    void store (int *values) const;

    /// Store the first n values into memory
    void store (int *values, int n) const;

    /// Store the least significant 16 bits of each element into adjacent
    /// unsigned shorts.
    void store (unsigned short *values) const;

    /// Store the least significant 8 bits of each element into adjacent
    /// unsigned chars.
    void store (unsigned char *values) const;

    /// Masked load -- read from values[] where mask is 1, load zero where
    /// mask is 0.
    void load_mask (const vbool_t &mask, const value_t *values);
    void load_mask (int mask, const value_t *values) { load_mask(vbool_t(mask), values); }

    /// Masked store -- write to values[] where mask is enabled, don't
    /// touch values[] where it's not.
    void store_mask (const vbool_t &mask, value_t *values) const;
    void store_mask (int mask, value_t *values) const { store_mask(vbool_t(mask), values); }

    /// Load values from addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void gather (const value_t *baseptr, const vint_t& vindex);
    /// Gather elements defined by the mask, leave others unchanged.
    template<int scale=4>
    void gather_mask (const bool_t& mask, const value_t *baseptr, const vint_t& vindex);
    template<int scale=4>
    void gather_mask (int mask, const value_t *baseptr, const vint_t& vindex) {
        gather_mask<scale> (vbool_t(mask), baseptr, vindex);
    }

    /// Store values at addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void scatter (value_t *baseptr, const vint_t& vindex) const;
    /// Scatter elements defined by the mask
    template<int scale=4>
    void scatter_mask (const bool_t& mask, value_t *baseptr, const vint_t& vindex) const;
    template<int scale=4>
    void scatter_mask (int mask, value_t *baseptr, const vint_t& vindex) const {
        scatter_mask<scale> (vbool_t(mask), baseptr, vindex);
    }

    // Arithmetic operators (component-by-component)
    friend vint16 operator+ (const vint16& a, const vint16& b);
    friend vint16 operator- (const vint16& a);
    friend vint16 operator- (const vint16& a, const vint16& b);
    friend vint16 operator* (const vint16& a, const vint16& b);
    friend vint16 operator/ (const vint16& a, const vint16& b);
    friend vint16 operator% (const vint16& a, const vint16& b);
    friend const vint16 & operator+= (vint16& a, const vint16& b);
    friend const vint16 & operator-= (vint16& a, const vint16& b);
    friend const vint16 & operator*= (vint16& a, const vint16& b);
    friend const vint16 & operator/= (vint16& a, const vint16& b);
    friend const vint16 & operator%= (vint16& a, const vint16& b);
    // Bitwise operators (component-by-component)
    friend vint16 operator& (const vint16& a, const vint16& b);
    friend vint16 operator| (const vint16& a, const vint16& b);
    friend vint16 operator^ (const vint16& a, const vint16& b);
    friend const vint16& operator&= (vint16& a, const vint16& b);
    friend const vint16& operator|= (vint16& a, const vint16& b);
    friend const vint16& operator^= (vint16& a, const vint16& b);
    friend vint16 operator~ (const vint16& a);
    friend vint16 operator<< (const vint16& a, unsigned int bits);
    friend vint16 operator>> (const vint16& a, unsigned int bits);
    friend const vint16& operator<<= (vint16& a, unsigned int bits);
    friend const vint16& operator>>= (vint16& a, unsigned int bits);
    // Comparison operators (component-by-component)
    friend vbool16 operator== (const vint16& a, const vint16& b);
    friend vbool16 operator!= (const vint16& a, const vint16& b);
    friend vbool16 operator<  (const vint16& a, const vint16& b);
    friend vbool16 operator>  (const vint16& a, const vint16& b);
    friend vbool16 operator>= (const vint16& a, const vint16& b);
    friend vbool16 operator<= (const vint16& a, const vint16& b);

    /// Stream output
    friend std::ostream& operator<< (std::ostream& cout, const vint16& a);

private:
    // The actual data representation
    union {
        simd_t  m_simd;
        value_t m_val[elements];
        vint8    m_8[2];
    };
};



/// Shift right logical -- unsigned shift. This differs from operator>>
/// in how it handles the sign bit.  (1<<31) >> 1 == (1<<31), but
/// srl((1<<31),1) == 1<<30.
vint16 srl (const vint16& val, const unsigned int bits);

/// Shuffle groups of 4
template<int i0, int i1, int i2, int i3>
vint16 shuffle4 (const vint16& a);

/// shuffle4<i>(a) is the same as shuffle4<i,i,i,i>(a)
template<int i> vint16 shuffle4 (const vint16& a);

/// Shuffle within each group of 4
template<int i0, int i1, int i2, int i3>
vint16 shuffle (const vint16& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> vint16 shuffle (const vint16& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> int extract (const vint16& v);

/// Helper: substitute val for a[i]
template<int i> vint16 insert (const vint16& a, int val);

/// The sum of all components, returned in all components.
vint16 vreduce_add (const vint16& v);

// Reduction across all components
int reduce_add (const vint16& v);
int reduce_and (const vint16& v);
int reduce_or  (const vint16& v);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// and b (if mask[i] is true), i.e., mask[i] ? b[i] : a[i].
vint16 blend (const vint16& a, const vint16& b, const vbool16& mask);

/// Use a bool mask to select between `a` (if mask[i] is true) or 0 if
/// mask[i] is false), i.e., mask[i] ? a[i] : 0. Equivalent to
/// blend(0,a,mask).
vint16 blend0 (const vint16& a, const vbool16& mask);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// or 0 (if mask[i] is true), i.e., mask[i] ? 0 : a[i]. Equivalent to
/// blend(0,a,!mask), or blend(a,0,mask).
vint16 blend0not (const vint16& a, const vbool16& mask);

/// Select 'a' where mask is true, 'b' where mask is false. Sure, it's a
/// synonym for blend with arguments rearranged, but this is more clear
/// because the arguments are symmetric to scalar (cond ? a : b).
vint16 select (const vbool16& mask, const vint16& a, const vint16& b);

// Per-element math
vint16 abs (const vint16& a);
vint16 min (const vint16& a, const vint16& b);
vint16 max (const vint16& a, const vint16& b);

/// Circular bit rotate by s bits, for N values at once.
vint16 rotl (const vint16& x, const int s);
// DEPRECATED(2.1)
vint16 rotl32 (const vint16& x, const unsigned int k);

/// andnot(a,b) returns ((~a) & b)
vint16 andnot (const vint16& a, const vint16& b);

/// Bitcast back and forth to intN (not a convert -- move the bits!)
vint16 bitcast_to_int (const vbool16& x);
vint16 bitcast_to_int (const vfloat16& x);
vfloat16 bitcast_to_float (const vint16& x);

// safe_mod(a,b) is like a%b, but safely returns 0 when b==0.
vint16 safe_mod (const vint16& a, const vint16& b);
vint16 safe_mod (const vint16& a, int b);





/// Floating point 4-vector, accelerated by SIMD instructions when
/// available.
class vfloat4 {
public:
    static const char* type_name() { return "vfloat4"; }
    typedef float value_t;    ///< Underlying equivalent scalar value type
    enum { elements = 4 };    ///< Number of scalar elements
    enum { paddedelements = 4 }; ///< Number of scalar elements for full pad
    enum { bits = elements*32 }; ///< Total number of bits
    typedef simd_raw_t<float,4>::type simd_t;  ///< the native SIMD type used
    typedef vfloat4 vfloat_t; ///< SIMD int type
    typedef vint4 vint_t;     ///< SIMD int type
    typedef vbool4 vbool_t;   ///< SIMD bool type
    typedef vint4 int_t;      // old name (deprecated 1.8)
    typedef vbool4 bool_t;    // old name (deprecated 1.8)

    /// Default constructor (contents undefined)
    vfloat4 () { }

    /// Construct from a single value (store it in all slots)
    vfloat4 (float a) { load(a); }

    /// Construct from 3 or 4 values
    vfloat4 (float a, float b, float c, float d=0.0f) { load(a,b,c,d); }

    /// Construct from a pointer to 4 values
    vfloat4 (const float *f) { load (f); }

    /// Copy construct from another vfloat4
    vfloat4 (const vfloat4 &other) { m_simd = other.m_simd; }

    /// Construct from an vint4 (promoting all components to float)
    explicit vfloat4 (const vint4& ival);

    /// Construct from the underlying SIMD type
    vfloat4 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Return a pointer to the underlying scalar type
    const value_t* data () const { return (const value_t*)this; }
    value_t* data () { return (value_t*)this; }

    /// Construct from a Imath::V3f
    explicit vfloat4 (const Imath::V3f &v) { load (v[0], v[1], v[2]); }

    /// Cast to a Imath::V3f
    const Imath::V3f& V3f () const { return *(const Imath::V3f*)this; }

    /// Construct from a Imath::V4f
    explicit vfloat4 (const Imath::V4f &v) { load ((const float *)&v); }

    /// Cast to a Imath::V4f
    const Imath::V4f& V4f () const { return *(const Imath::V4f*)this; }

    /// Construct from a pointer to 4 unsigned short values
    explicit vfloat4 (const unsigned short *vals) { load(vals); }

    /// Construct from a pointer to 4 short values
    explicit vfloat4 (const short *vals) { load(vals); }

    /// Construct from a pointer to 4 unsigned char values
    explicit vfloat4 (const unsigned char *vals) { load(vals); }

    /// Construct from a pointer to 4 char values
    explicit vfloat4 (const char *vals) { load(vals); }

#ifdef _HALF_H_
    /// Construct from a pointer to 4 half (16 bit float) values
    explicit vfloat4 (const half *vals) { load(vals); }
#endif

    /// Assign a single value to all components
    const vfloat4 & operator= (float a) { load(a); return *this; }

    /// Assign a vfloat4
    const vfloat4 & operator= (vfloat4 other) {
        m_simd = other.m_simd;
        return *this;
    }

    /// Return a vfloat4 with all components set to 0.0
    static const vfloat4 Zero ();

    /// Return a vfloat4 with all components set to 1.0
    static const vfloat4 One ();

    /// Return a vfloat4 with incremented components (e.g., 0.0,1.0,2.0,3.0).
    /// Optional argument can give a non-zero starting point and non-1 step.
    static const vfloat4 Iota (float start=0.0f, float step=1.0f);

    /// Set all components to 0.0
    void clear ();

    /// Assign from a Imath::V4f
    const vfloat4 & operator= (const Imath::V4f &v);

    /// Assign from a Imath::V3f
    const vfloat4 & operator= (const Imath::V3f &v);

    /// Component access (get)
    float operator[] (int i) const;
    /// Component access (set)
    float& operator[] (int i);

    /// Component access (set).
    void setcomp (int i, float value);

    value_t x () const;
    value_t y () const;
    value_t z () const;
    value_t w () const;
    void set_x (value_t val);
    void set_y (value_t val);
    void set_z (value_t val);
    void set_w (value_t val);

    /// Helper: load a single value into all components
    void load (float val);

    /// Helper: load 3 or 4 values. (If 3 are supplied, the 4th will be 0.)
    void load (float a, float b, float c, float d=0.0f);

    /// Load from an array of 4 values
    void load (const float *values);

    /// Load from a partial array of <=4 values. Unassigned values are
    /// undefined.
    void load (const float *values, int n);

    /// Load from an array of 4 unsigned short values, convert to float
    void load (const unsigned short *values);

    /// Load from an array of 4 short values, convert to float
    void load (const short *values);

    /// Load from an array of 4 unsigned char values, convert to float
    void load (const unsigned char *values);

    /// Load from an array of 4 char values, convert to float
    void load (const char *values);

#ifdef _HALF_H_
    /// Load from an array of 4 half values, convert to float
    void load (const half *values);
#endif /* _HALF_H_ */

    void store (float *values) const;

    /// Store the first n values into memory
    void store (float *values, int n) const;

#ifdef _HALF_H_
    void store (half *values) const;
#endif

    /// Masked load -- read from values[] where mask is 1, load zero where
    /// mask is 0.
    void load_mask (int mask, const value_t *values);
    void load_mask (const vbool_t& mask, const value_t *values);

    /// Masked store -- write to values[] where mask is enabled, don't
    /// touch values[] where it's not.
    void store_mask (int mask, value_t *values) const;
    void store_mask (const vbool_t& mask, value_t *values) const;

    /// Load values from addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void gather (const value_t *baseptr, const vint_t& vindex);
    /// Gather elements defined by the mask, leave others unchanged.
    template<int scale=4>
    void gather_mask (const bool_t& mask, const value_t *baseptr, const vint_t& vindex);
    template<int scale=4>
    void gather_mask (int mask, const value_t *baseptr, const vint_t& vindex);

    /// Store values at addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void scatter (value_t *baseptr, const vint_t& vindex) const;
    /// Scatter elements defined by the mask
    template<int scale=4>
    void scatter_mask (const bool_t& mask, value_t *baseptr, const vint_t& vindex) const;
    template<int scale=4>
    void scatter_mask (int mask, value_t *baseptr, const vint_t& vindex) const;

    // Arithmetic operators
    friend vfloat4 operator+ (const vfloat4& a, const vfloat4& b);
    const vfloat4 & operator+= (const vfloat4& a);
    vfloat4 operator- () const;
    friend vfloat4 operator- (const vfloat4& a, const vfloat4& b);
    const vfloat4 & operator-= (const vfloat4& a);
    friend vfloat4 operator* (const vfloat4& a, const vfloat4& b);
    friend vfloat4 operator* (const vfloat4& a, float b);
    friend vfloat4 operator* (float a, const vfloat4& b);
    const vfloat4 & operator*= (const vfloat4& a);
    const vfloat4 & operator*= (float val);
    friend vfloat4 operator/ (const vfloat4& a, const vfloat4& b);
    const vfloat4 & operator/= (const vfloat4& a);
    const vfloat4 & operator/= (float val);

    // Comparison operations
    friend vbool4 operator== (const vfloat4& a, const vfloat4& b);
    friend vbool4 operator!= (const vfloat4& a, const vfloat4& b);
    friend vbool4 operator<  (const vfloat4& a, const vfloat4& b);
    friend vbool4 operator>  (const vfloat4& a, const vfloat4& b);
    friend vbool4 operator>= (const vfloat4& a, const vfloat4& b);
    friend vbool4 operator<= (const vfloat4& a, const vfloat4& b);

    // Some oddball items that are handy

    /// Combine the first two components of A with the first two components
    /// of B.
    friend vfloat4 AxyBxy (const vfloat4& a, const vfloat4& b);

    /// Combine the first two components of A with the first two components
    /// of B, but interleaved.
    friend vfloat4 AxBxAyBy (const vfloat4& a, const vfloat4& b);

    /// Return xyz components, plus 0 for w
    vfloat4 xyz0 () const;

    /// Return xyz components, plus 1 for w
    vfloat4 xyz1 () const;

    /// Stream output
    friend inline std::ostream& operator<< (std::ostream& cout, const vfloat4& val);

protected:
    // The actual data representation
    union {
        simd_t  m_simd;
        value_t m_val[paddedelements];
    };
};


/// Helper: shuffle/swizzle with constant (templated) indices.
/// Example: shuffle<1,1,2,2>(vbool4(a,b,c,d)) returns (b,b,c,c)
template<int i0, int i1, int i2, int i3> vfloat4 shuffle (const vfloat4& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> vfloat4 shuffle (const vfloat4& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> float extract (const vfloat4& a);

/// Helper: substitute val for a[i]
template<int i> vfloat4 insert (const vfloat4& a, float val);

/// The sum of all components, returned in all components.
vfloat4 vreduce_add (const vfloat4& v);

/// The sum of all components, returned as a scalar.
float reduce_add (const vfloat4& v);

/// Return the float dot (inner) product of a and b in every component.
vfloat4 vdot (const vfloat4 &a, const vfloat4 &b);

/// Return the float dot (inner) product of a and b.
float dot (const vfloat4 &a, const vfloat4 &b);

/// Return the float 3-component dot (inner) product of a and b in
/// all components.
vfloat4 vdot3 (const vfloat4 &a, const vfloat4 &b);

/// Return the float 3-component dot (inner) product of a and b.
float dot3 (const vfloat4 &a, const vfloat4 &b);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// and b (if mask[i] is true), i.e., mask[i] ? b[i] : a[i].
vfloat4 blend (const vfloat4& a, const vfloat4& b, const vbool4& mask);

/// Use a bool mask to select between `a` (if mask[i] is true) or 0 if
/// mask[i] is false), i.e., mask[i] ? a[i] : 0. Equivalent to
/// blend(0,a,mask).
vfloat4 blend0 (const vfloat4& a, const vbool4& mask);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// or 0 (if mask[i] is true), i.e., mask[i] ? 0 : a[i]. Equivalent to
/// blend(0,a,!mask), or blend(a,0,mask).
vfloat4 blend0not (const vfloat4& a, const vbool4& mask);

/// "Safe" divide of vfloat4/vfloat4 -- for any component of the divisor
/// that is 0, return 0 rather than Inf.
vfloat4 safe_div (const vfloat4 &a, const vfloat4 &b);

/// Homogeneous divide to turn a vfloat4 into a vfloat3.
vfloat3 hdiv (const vfloat4 &a);

/// Select 'a' where mask is true, 'b' where mask is false. Sure, it's a
/// synonym for blend with arguments rearranged, but this is more clear
/// because the arguments are symmetric to scalar (cond ? a : b).
vfloat4 select (const vbool4& mask, const vfloat4& a, const vfloat4& b);

// Per-element math
vfloat4 abs (const vfloat4& a);    ///< absolute value (float)
vfloat4 sign (const vfloat4& a);   ///< 1.0 when value >= 0, -1 when negative
vfloat4 ceil (const vfloat4& a);
vfloat4 floor (const vfloat4& a);
vint4 ifloor (const vfloat4& a);    ///< (int)floor
inline vint4 floori (const vfloat4& a) { return ifloor(a); }  // DEPRECATED(1.8) alias

/// Per-element round to nearest integer.
/// CAVEAT: the rounding when mid-way between integers may differ depending
/// on hardware. Intel SSE/AVX does "banker's founding" (to nearest even
/// integer) but std::round() says to round away from 0 regardless of
/// current rounding mode (but that is multiple instructions on x64).
/// USE WITH CAUTION, and maybe avoid this if it is critical to exactly
/// match std::round().
vfloat4 round (const vfloat4& a);

/// Per-element round to nearest integer (equivalent to vint(round(a))).
/// CAVEAT: On SSE/AVX this uses banker's rounding, which may differ from
/// C++ std::rint() which says to use the current rounding mode.
/// USE WITH CAUTION, and maybe avoid this if it is critical to exactly
/// match std::rint().
vint4 rint (const vfloat4& a);

vfloat4 rcp_fast (const vfloat4 &a);  ///< Fast, approximate 1/a
vfloat4 sqrt (const vfloat4 &a);
vfloat4 rsqrt (const vfloat4 &a);   ///< Fully accurate 1/sqrt
vfloat4 rsqrt_fast (const vfloat4 &a);  ///< Fast, approximate 1/sqrt
vfloat4 min (const vfloat4& a, const vfloat4& b); ///< Per-element min
vfloat4 max (const vfloat4& a, const vfloat4& b); ///< Per-element max
template <typename T> T exp (const T& v);  // template for all SIMD variants
template <typename T> T log (const T& v);

/// andnot(a,b) returns ((~a) & b)
vfloat4 andnot (const vfloat4& a, const vfloat4& b);

// Fused multiply and add (or subtract):
vfloat4 madd (const vfloat4& a, const vfloat4& b, const vfloat4& c); // a*b + c
vfloat4 msub (const vfloat4& a, const vfloat4& b, const vfloat4& c); // a*b - c
vfloat4 nmadd (const vfloat4& a, const vfloat4& b, const vfloat4& c); // -a*b + c
vfloat4 nmsub (const vfloat4& a, const vfloat4& b, const vfloat4& c); // -a*b - c

/// Transpose the rows and columns of the 4x4 matrix [a b c d].
/// In the end, a will have the original (a[0], b[0], c[0], d[0]),
/// b will have the original (a[1], b[1], c[1], d[1]), and so on.
void transpose (vfloat4 &a, vfloat4 &b, vfloat4 &c, vfloat4 &d);
void transpose (const vfloat4& a, const vfloat4& b, const vfloat4& c, const vfloat4& d,
                vfloat4 &r0, vfloat4 &r1, vfloat4 &r2, vfloat4 &r3);

/// Make a vfloat4 consisting of the first element of each of 4 vfloat4's.
vfloat4 AxBxCxDx (const vfloat4& a, const vfloat4& b,
                 const vfloat4& c, const vfloat4& d);



/// Floating point 3-vector, aligned to be internally identical to a vfloat4.
/// The way it differs from vfloat4 is that all of he load functions only
/// load three values, and all the stores only store 3 values. The vast
/// majority of ops just fall back to the vfloat4 version, and so will
/// operate on the 4th component, but we won't care about that results.
class vfloat3 : public vfloat4 {
public:
    static const char* type_name() { return "vfloat3"; }
    enum { elements = 3 };    ///< Number of scalar elements
    enum { paddedelements = 4 }; ///< Number of scalar elements for full pad

    /// Default constructor (contents undefined)
    vfloat3 () { }

    /// Construct from a single value (store it in all slots)
    vfloat3 (float a) { load(a); }

    /// Construct from 3 or 4 values
    vfloat3 (float a, float b, float c) { vfloat4::load(a,b,c); }

    /// Construct from a pointer to 4 values
    vfloat3 (const float *f) { load (f); }

    /// Copy construct from another vfloat3
    vfloat3 (const vfloat3 &other);

    /// Construct from a vfloat4. Note: it will not zero out the internal
    /// 4th component, but rather accept on faith that the vfloat4 you are
    /// giving it is a valid vfloat3. Be careful!
    explicit vfloat3 (const vfloat4 &other);

#if OIIO_SIMD
    /// Construct from the underlying SIMD type. Note: it will not zero out
    /// the internal 4th component, but rather accept on faith that the
    /// vfloat4 you are giving it is a valid vfloat3. Be careful!
    explicit vfloat3 (const simd_t& m) : vfloat4(m) { }
#endif

    /// Construct from a Imath::V3f
    vfloat3 (const Imath::V3f &v) : vfloat4(v) { }

    /// Cast to a Imath::V3f
    const Imath::V3f& V3f () const { return *(const Imath::V3f*)this; }

    /// Construct from a pointer to 4 unsigned short values
    explicit vfloat3 (const unsigned short *vals) { load(vals); }

    /// Construct from a pointer to 4 short values
    explicit vfloat3 (const short *vals) { load(vals); }

    /// Construct from a pointer to 4 unsigned char values
    explicit vfloat3 (const unsigned char *vals) { load(vals); }

    /// Construct from a pointer to 4 char values
    explicit vfloat3 (const char *vals) { load(vals); }

#ifdef _HALF_H_
    /// Construct from a pointer to 4 half (16 bit float) values
    explicit vfloat3 (const half *vals) { load(vals); }
#endif

    /// Assign a single value to all components
    const vfloat3 & operator= (float a) { load(a); return *this; }

    /// Return a vfloat3 with all components set to 0.0
    static const vfloat3 Zero ();

    /// Return a vfloat3 with all components set to 1.0
    static const vfloat3 One ();

    /// Return a vfloat3 with incremented components (e.g., 0.0,1.0,2.0).
    /// Optional argument can give a non-zero starting point and non-1 step.
    static const vfloat3 Iota (float start=0.0f, float step=1.0f);

    /// Helper: load a single value into all components
    void load (float val);

    /// Load from an array of 4 values
    void load (const float *values);

    /// Load from an array of 4 values
    void load (const float *values, int n);

    /// Load from an array of 4 unsigned short values, convert to float
    void load (const unsigned short *values);

    /// Load from an array of 4 short values, convert to float
    void load (const short *values);

    /// Load from an array of 4 unsigned char values, convert to float
    void load (const unsigned char *values);

    /// Load from an array of 4 char values, convert to float
    void load (const char *values);

#ifdef _HALF_H_
    /// Load from an array of 4 half values, convert to float
    void load (const half *values);
#endif /* _HALF_H_ */

    void store (float *values) const;

    void store (float *values, int n) const;

#ifdef _HALF_H_
    void store (half *values) const;
#endif

    /// Store into an Imath::V3f reference.
    void store (Imath::V3f &vec) const;

    // Math operators -- define in terms of vfloat3.
    friend vfloat3 operator+ (const vfloat3& a, const vfloat3& b);
    const vfloat3 & operator+= (const vfloat3& a);
    vfloat3 operator- () const;
    friend vfloat3 operator- (const vfloat3& a, const vfloat3& b);
    const vfloat3 & operator-= (const vfloat3& a);
    friend vfloat3 operator* (const vfloat3& a, const vfloat3& b);
    friend vfloat3 operator* (const vfloat3& a, float b);
    friend vfloat3 operator* (float a, const vfloat3& b);
    const vfloat3 & operator*= (const vfloat3& a);
    const vfloat3 & operator*= (float a);
    friend vfloat3 operator/ (const vfloat3& a, const vfloat3& b);
    const vfloat3 & operator/= (const vfloat3& a);
    const vfloat3 & operator/= (float a);

    /// Square of the length of the vector
    float length2() const;
    /// Length of the vector
    float length() const;

    /// Return a normalized version of the vector.
    vfloat3 normalized () const;
    /// Return a fast, approximate normalized version of the vector.
    vfloat3 normalized_fast () const;
    /// Normalize in place.
    void normalize() { *this = normalized(); }

    /// Stream output
    friend inline std::ostream& operator<< (std::ostream& cout, const vfloat3& val);
};



// Per-element math on float3
vfloat3 abs (const vfloat3& a);
vfloat3 sign (const vfloat3& a);
vfloat3 ceil (const vfloat3& a);
vfloat3 floor (const vfloat3& a);
vfloat3 round (const vfloat3& a);



/// SIMD-based 4x4 matrix. This is guaranteed to have memory layout (when
/// not in registers) isomorphic to Imath::M44f.
class matrix44 {
public:
    // Uninitialized
    OIIO_FORCEINLINE matrix44 ()
#ifndef OIIO_SIMD_SSE
        : m_mat(Imath::UNINITIALIZED)
#endif
    { }

    /// Construct from a reference to an Imath::M44f
    OIIO_FORCEINLINE explicit matrix44 (const Imath::M44f &M) {
#if OIIO_SIMD_SSE
        m_row[0].load (M[0]);
        m_row[1].load (M[1]);
        m_row[2].load (M[2]);
        m_row[3].load (M[3]);
#else
        m_mat = M;
#endif
    }

    /// Construct from a float array
    OIIO_FORCEINLINE explicit matrix44 (const float *f) {
#if OIIO_SIMD_SSE
        m_row[0].load (f+0);
        m_row[1].load (f+4);
        m_row[2].load (f+8);
        m_row[3].load (f+12);
#else
        m_mat = *(const Imath::M44f*)f;
#endif
    }

    /// Construct from 4 vfloat4 rows
    OIIO_FORCEINLINE explicit matrix44 (const vfloat4& a, const vfloat4& b,
                                        const vfloat4& c, const vfloat4& d) {
#if OIIO_SIMD_SSE
        m_row[0] = a; m_row[1] = b; m_row[2] = c; m_row[3] = d;
#else
        a.store (m_mat[0]);
        b.store (m_mat[1]);
        c.store (m_mat[2]);
        d.store (m_mat[3]);
#endif
    }
    /// Construct from 4 float[4] rows
    OIIO_FORCEINLINE explicit matrix44 (const float *a, const float *b,
                                        const float *c, const float *d) {
#if OIIO_SIMD_SSE
        m_row[0].load(a); m_row[1].load(b); m_row[2].load(c); m_row[3].load(d);
#else
        memcpy (m_mat[0], a, 4*sizeof(float));
        memcpy (m_mat[1], b, 4*sizeof(float));
        memcpy (m_mat[2], c, 4*sizeof(float));
        memcpy (m_mat[3], d, 4*sizeof(float));
#endif
    }

    /// Construct from 16 floats
    OIIO_FORCEINLINE matrix44 (float f00, float f01, float f02, float f03,
                               float f10, float f11, float f12, float f13,
                               float f20, float f21, float f22, float f23,
                               float f30, float f31, float f32, float f33)
    {
#if OIIO_SIMD_SSE
        m_row[0].load (f00, f01, f02, f03);
        m_row[1].load (f10, f11, f12, f13);
        m_row[2].load (f20, f21, f22, f23);
        m_row[3].load (f30, f31, f32, f33);
#else
        m_mat[0][0] = f00; m_mat[0][1] = f01; m_mat[0][2] = f02; m_mat[0][3] = f03;
        m_mat[1][0] = f10; m_mat[1][1] = f11; m_mat[1][2] = f12; m_mat[1][3] = f13;
        m_mat[2][0] = f20; m_mat[2][1] = f21; m_mat[2][2] = f22; m_mat[2][3] = f23;
        m_mat[3][0] = f30; m_mat[3][1] = f31; m_mat[3][2] = f32; m_mat[3][3] = f33;
#endif
    }

    /// Present as an Imath::M44f
    const Imath::M44f& M44f() const;

    /// Return one row
    vfloat4 operator[] (int i) const;

    /// Return the transposed matrix
    matrix44 transposed () const;

    /// Transform 3-point V by 4x4 matrix M.
    vfloat3 transformp (const vfloat3 &V) const;

    /// Transform 3-vector V by 4x4 matrix M.
    vfloat3 transformv (const vfloat3 &V) const;

    /// Transform 3-vector V by the transpose of 4x4 matrix M.
    vfloat3 transformvT (const vfloat3 &V) const;

    friend vfloat4 operator* (const vfloat4 &V, const matrix44& M);
    friend vfloat4 operator* (const matrix44& M, const vfloat4 &V);

    bool operator== (const matrix44& m) const;

    bool operator== (const Imath::M44f& m) const ;
    friend bool operator== (const Imath::M44f& a, const matrix44 &b);

    bool operator!= (const matrix44& m) const;

    bool operator!= (const Imath::M44f& m) const;
    friend bool operator!= (const Imath::M44f& a, const matrix44 &b);

    /// Return the inverse of the matrix.
    matrix44 inverse() const;

    /// Stream output
    friend inline std::ostream& operator<< (std::ostream& cout, const matrix44 &M);

private:
#if OIIO_SIMD_SSE
    vfloat4 m_row[4];
#else
    Imath::M44f m_mat;
#endif
};

/// Transform 3-point V by 4x4 matrix M.
vfloat3 transformp (const matrix44 &M, const vfloat3 &V);
vfloat3 transformp (const Imath::M44f &M, const vfloat3 &V);

/// Transform 3-vector V by 4x4 matrix M.
vfloat3 transformv (const matrix44 &M, const vfloat3 &V);
vfloat3 transformv (const Imath::M44f &M, const vfloat3 &V);

// Transform 3-vector by the transpose of 4x4 matrix M.
vfloat3 transformvT (const matrix44 &M, const vfloat3 &V);
vfloat3 transformvT (const Imath::M44f &M, const vfloat3 &V);




/// Floating point 8-vector, accelerated by SIMD instructions when
/// available.
class vfloat8 {
public:
    static const char* type_name() { return "vfloat8"; }
    typedef float value_t;    ///< Underlying equivalent scalar value type
    enum { elements = 8 };    ///< Number of scalar elements
    enum { paddedelements = 8 }; ///< Number of scalar elements for full pad
    enum { bits = elements*32 }; ///< Total number of bits
    typedef simd_raw_t<float,8>::type simd_t;  ///< the native SIMD type used
    typedef vfloat8 vfloat_t; ///< SIMD int type
    typedef vint8 vint_t;     ///< SIMD int type
    typedef vbool8 vbool_t;   ///< SIMD bool type
    typedef vint8 int_t;      // old name (deprecated 1.8)
    typedef vbool8 bool_t;    // old name (deprecated 1.8)

    /// Default constructor (contents undefined)
    vfloat8 () { }

    /// Construct from a single value (store it in all slots)
    vfloat8 (float a) { load(a); }

    /// Construct from 8 values
    vfloat8 (float a, float b, float c, float d,
            float e, float f, float g, float h) { load(a,b,c,d,e,f,g,h); }

    /// Construct from a pointer to 8 values
    vfloat8 (const float *f) { load (f); }

    /// Copy construct from another vfloat8
    vfloat8 (const vfloat8 &other) { m_simd = other.m_simd; }

    /// Construct from an int vector (promoting all components to float)
    explicit vfloat8 (const vint8& ival);

    /// Construct from two vfloat4's
    vfloat8 (const vfloat4 &lo, const vfloat4 &hi);

    /// Construct from the underlying SIMD type
    vfloat8 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Return a pointer to the underlying scalar type
    const value_t* data () const { return (const value_t*)this; }
    value_t* data () { return (value_t*)this; }

    /// Construct from a pointer to unsigned short values
    explicit vfloat8 (const unsigned short *vals) { load(vals); }

    /// Construct from a pointer to short values
    explicit vfloat8 (const short *vals) { load(vals); }

    /// Construct from a pointer to unsigned char values
    explicit vfloat8 (const unsigned char *vals) { load(vals); }

    /// Construct from a pointer to char values
    explicit vfloat8 (const char *vals) { load(vals); }

#ifdef _HALF_H_
    /// Construct from a pointer to half (16 bit float) values
    explicit vfloat8 (const half *vals) { load(vals); }
#endif

    /// Assign a single value to all components
    const vfloat8& operator= (float a) { load(a); return *this; }

    /// Assign a vfloat8
    const vfloat8& operator= (vfloat8 other) {
        m_simd = other.m_simd;
        return *this;
    }

    /// Return a vfloat8 with all components set to 0.0
    static const vfloat8 Zero ();

    /// Return a vfloat8 with all components set to 1.0
    static const vfloat8 One ();

    /// Return a vfloat8 with incremented components (e.g., 0,1,2,3,...)
    /// Optional argument can give a non-zero starting point and non-1 step.
    static const vfloat8 Iota (float start=0.0f, float step=1.0f);

    /// Set all components to 0.0
    void clear ();

    /// Component access (get)
    float operator[] (int i) const;
    /// Component access (set)
    float& operator[] (int i);

    /// Component access (set).
    void setcomp (int i, float value);

    value_t x () const;
    value_t y () const;
    value_t z () const;
    value_t w () const;
    void set_x (value_t val);
    void set_y (value_t val);
    void set_z (value_t val);
    void set_w (value_t val);

    /// Extract the lower precision vfloat4
    vfloat4 lo () const;

    /// Extract the higher precision vfloat4
    vfloat4 hi () const;

    /// Helper: load a single value into all components
    void load (float val);

    /// Helper: load 8 values
    void load (float a, float b, float c, float d,
               float e, float f, float g, float h);

    /// Load from an array of values
    void load (const float *values);

    /// Load from a partial array of <=8 values. Unassigned values are
    /// undefined.
    void load (const float *values, int n);

    /// Load from an array of 8 unsigned short values, convert to float
    void load (const unsigned short *values);

    /// Load from an array of 8 short values, convert to float
    void load (const short *values);

    /// Load from an array of 8 unsigned char values, convert to float
    void load (const unsigned char *values);

    /// Load from an array of 8 char values, convert to float
    void load (const char *values);

#ifdef _HALF_H_
    /// Load from an array of 8 half values, convert to float
    void load (const half *values);
#endif /* _HALF_H_ */

    void store (float *values) const;

    /// Store the first n values into memory
    void store (float *values, int n) const;

#ifdef _HALF_H_
    void store (half *values) const;
#endif

    /// Masked load -- read from values[] where mask is 1, load zero where
    /// mask is 0.
    void load_mask (int mask, const value_t *values);
    void load_mask (const vbool_t& mask, const value_t *values);

    /// Masked store -- write to values[] where mask is enabled, don't
    /// touch values[] where it's not.
    void store_mask (int mask, value_t *values) const;
    void store_mask (const vbool_t& mask, value_t *values) const;

    /// Load values from addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void gather (const value_t *baseptr, const vint_t& vindex);
    template<int scale=4>
    // Fastest way to fill with all 1 bits is to cmp any value to itself.
    void gather_mask (const bool_t& mask, const value_t *baseptr, const vint_t& vindex);
    template<int scale=4>
    void gather_mask (int mask, const value_t *baseptr, const vint_t& vindex);

    /// Store values at addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void scatter (value_t *baseptr, const vint_t& vindex) const;
    /// Scatter elements defined by the mask
    template<int scale=4>
    void scatter_mask (const bool_t& mask, value_t *baseptr, const vint_t& vindex) const;
    template<int scale=4>
    void scatter_mask (int mask, value_t *baseptr, const vint_t& vindex) const;

    // Arithmetic operators (component-by-component)
    friend vfloat8 operator+ (const vfloat8& a, const vfloat8& b);
    friend vfloat8 operator- (const vfloat8& a);
    friend vfloat8 operator- (const vfloat8& a, const vfloat8& b);
    friend vfloat8 operator* (const vfloat8& a, const vfloat8& b);
    friend vfloat8 operator* (const vfloat8& a, float b);
    friend vfloat8 operator* (float a, const vfloat8& b);
    friend vfloat8 operator/ (const vfloat8& a, const vfloat8& b);
    friend vfloat8 operator% (const vfloat8& a, const vfloat8& b);
    friend const vfloat8 & operator+= (vfloat8& a, const vfloat8& b);
    friend const vfloat8 & operator-= (vfloat8& a, const vfloat8& b);
    friend const vfloat8 & operator*= (vfloat8& a, const vfloat8& b);
    friend const vfloat8 & operator/= (vfloat8& a, const vfloat8& b);

    // Comparison operations
    friend vbool8 operator== (const vfloat8& a, const vfloat8& b);
    friend vbool8 operator!= (const vfloat8& a, const vfloat8& b);
    friend vbool8 operator<  (const vfloat8& a, const vfloat8& b);
    friend vbool8 operator>  (const vfloat8& a, const vfloat8& b);
    friend vbool8 operator>= (const vfloat8& a, const vfloat8& b);
    friend vbool8 operator<= (const vfloat8& a, const vfloat8& b);

    // Some oddball items that are handy

    /// Stream output
    friend inline std::ostream& operator<< (std::ostream& cout, const vfloat8& val);

protected:
    // The actual data representation
    union {
        simd_t  m_simd;
        value_t m_val[paddedelements];
        vfloat4 m_4[2];
    };
};


/// Helper: shuffle/swizzle with constant (templated) indices.
/// Example: shuffle<1,1,2,2>(vbool4(a,b,c,d)) returns (b,b,c,c)
template<int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7>
vfloat8 shuffle (const vfloat8& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i,...>(a)
template<int i> vfloat8 shuffle (const vfloat8& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> float extract (const vfloat8& a);

/// Helper: substitute val for a[i]
template<int i> vfloat8 insert (const vfloat8& a, float val);

/// The sum of all components, returned in all components.
vfloat8 vreduce_add (const vfloat8& v);

/// The sum of all components, returned as a scalar.
float reduce_add (const vfloat8& v);

/// Return the float dot (inner) product of a and b in every component.
vfloat8 vdot (const vfloat8 &a, const vfloat8 &b);

/// Return the float dot (inner) product of a and b.
float dot (const vfloat8 &a, const vfloat8 &b);

/// Return the float 3-component dot (inner) product of a and b in
/// all components.
vfloat8 vdot3 (const vfloat8 &a, const vfloat8 &b);

/// Return the float 3-component dot (inner) product of a and b.
float dot3 (const vfloat8 &a, const vfloat8 &b);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// and b (if mask[i] is true), i.e., mask[i] ? b[i] : a[i].
vfloat8 blend (const vfloat8& a, const vfloat8& b, const vbool8& mask);

/// Use a bool mask to select between `a` (if mask[i] is true) or 0 if
/// mask[i] is false), i.e., mask[i] ? a[i] : 0. Equivalent to
/// blend(0,a,mask).
vfloat8 blend0 (const vfloat8& a, const vbool8& mask);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// or 0 (if mask[i] is true), i.e., mask[i] ? 0 : a[i]. Equivalent to
/// blend(0,a,!mask), or blend(a,0,mask).
vfloat8 blend0not (const vfloat8& a, const vbool8& mask);

/// "Safe" divide of vfloat8/vfloat8 -- for any component of the divisor
/// that is 0, return 0 rather than Inf.
vfloat8 safe_div (const vfloat8 &a, const vfloat8 &b);

/// Select 'a' where mask is true, 'b' where mask is false. Sure, it's a
/// synonym for blend with arguments rearranged, but this is more clear
/// because the arguments are symmetric to scalar (cond ? a : b).
vfloat8 select (const vbool8& mask, const vfloat8& a, const vfloat8& b);

// Per-element math
vfloat8 abs (const vfloat8& a);    ///< absolute value (float)
vfloat8 sign (const vfloat8& a);   ///< 1.0 when value >= 0, -1 when negative
vfloat8 ceil (const vfloat8& a);
vfloat8 floor (const vfloat8& a);
vint8 ifloor (const vfloat8& a);    ///< (int)floor
inline vint8 floori (const vfloat8& a) { return ifloor(a); }  // DEPRECATED(1.8) alias

/// Per-element round to nearest integer.
/// CAVEAT: the rounding when mid-way between integers may differ depending
/// on hardware. Intel SSE/AVX does "banker's founding" (to nearest even
/// integer) but std::round() says to round away from 0 regardless of
/// current rounding mode (but that is multiple instructions on x64).
/// USE WITH CAUTION, and maybe avoid this if it is critical to exactly
/// match std::round().
vfloat8 round (const vfloat8& a);

/// Per-element round to nearest integer (equivalent to vint(round(a))).
/// CAVEAT: On SSE/AVX this uses banker's rounding, which may differ from
/// C++ std::rint() which says to use the current rounding mode.
/// USE WITH CAUTION, and maybe avoid this if it is critical to exactly
/// match std::rint().
vint8 rint (const vfloat8& a);

vfloat8 rcp_fast (const vfloat8 &a);  ///< Fast, approximate 1/a
vfloat8 sqrt (const vfloat8 &a);
vfloat8 rsqrt (const vfloat8 &a);   ///< Fully accurate 1/sqrt
vfloat8 rsqrt_fast (const vfloat8 &a);  ///< Fast, approximate 1/sqrt
vfloat8 min (const vfloat8& a, const vfloat8& b); ///< Per-element min
vfloat8 max (const vfloat8& a, const vfloat8& b); ///< Per-element max
// vfloat8 exp (const vfloat8& v);  // See template with vfloat4
// vfloat8 log (const vfloat8& v);  // See template with vfloat4

/// andnot(a,b) returns ((~a) & b)
vfloat8 andnot (const vfloat8& a, const vfloat8& b);

// Fused multiply and add (or subtract):
vfloat8 madd (const vfloat8& a, const vfloat8& b, const vfloat8& c); // a*b + c
vfloat8 msub (const vfloat8& a, const vfloat8& b, const vfloat8& c); // a*b - c
vfloat8 nmadd (const vfloat8& a, const vfloat8& b, const vfloat8& c); // -a*b + c
vfloat8 nmsub (const vfloat8& a, const vfloat8& b, const vfloat8& c); // -a*b - c



/// Floating point 16-vector, accelerated by SIMD instructions when
/// available.
class vfloat16 {
public:
    static const char* type_name() { return "vfloat16"; }
    typedef float value_t;    ///< Underlying equivalent scalar value type
    enum { elements = 16 };    ///< Number of scalar elements
    enum { paddedelements = 16 }; ///< Number of scalar elements for full pad
    enum { bits = elements*32 }; ///< Total number of bits
    typedef simd_raw_t<float,16>::type simd_t;  ///< the native SIMD type used
    typedef vfloat16 vfloat_t; ///< SIMD int type
    typedef vint16 vint_t;     ///< SIMD int type
    typedef vbool16 vbool_t;   ///< SIMD bool type
    typedef vint16 int_t;      // old name (deprecated 1.8)
    typedef vbool16 bool_t;    // old name (deprecated 1.8)

    /// Default constructor (contents undefined)
    vfloat16 () { }

    /// Construct from a single value (store it in all slots)
    vfloat16 (float a) { load(a); }

    /// Construct from 16 values
    vfloat16 (float v0, float v1, float v2, float v3,
             float v4, float v5, float v6, float v7,
             float v8, float v9, float v10, float v11,
             float v12, float v13, float v14, float v15);

    /// Construct from a pointer to 16 values
    vfloat16 (const float *f) { load (f); }

    /// Copy construct from another vfloat16
    vfloat16 (const vfloat16 &other) { m_simd = other.m_simd; }

    /// Construct from an int vector (promoting all components to float)
    explicit vfloat16 (const vint16& ival);

    /// Construct from two vfloat8's
    vfloat16 (const vfloat8 &lo, const vfloat8 &hi);

    /// Construct from four vfloat4's
    vfloat16 (const vfloat4 &a, const vfloat4 &b, const vfloat4 &c, const vfloat4 &d);

    /// Construct from the underlying SIMD type
    vfloat16 (const simd_t& m) : m_simd(m) { }

    /// Return the raw SIMD type
    operator simd_t () const { return m_simd; }
    simd_t simd () const { return m_simd; }
    simd_t& simd () { return m_simd; }

    /// Return a pointer to the underlying scalar type
    const value_t* data () const { return (const value_t*)this; }
    value_t* data () { return (value_t*)this; }

    /// Construct from a pointer to unsigned short values
    explicit vfloat16 (const unsigned short *vals) { load(vals); }

    /// Construct from a pointer to short values
    explicit vfloat16 (const short *vals) { load(vals); }

    /// Construct from a pointer to unsigned char values
    explicit vfloat16 (const unsigned char *vals) { load(vals); }

    /// Construct from a pointer to char values
    explicit vfloat16 (const char *vals) { load(vals); }

#ifdef _HALF_H_
    /// Construct from a pointer to half (16 bit float) values
    explicit vfloat16 (const half *vals) { load(vals); }
#endif

    /// Assign a single value to all components
    const vfloat16& operator= (float a) { load(a); return *this; }

    /// Assign a vfloat16
    const vfloat16& operator= (vfloat16 other) {
        m_simd = other.m_simd;
        return *this;
    }

    /// Return a vfloat16 with all components set to 0.0
    static const vfloat16 Zero ();

    /// Return a vfloat16 with all components set to 1.0
    static const vfloat16 One ();

    /// Return a vfloat16 with incremented components (e.g., 0,1,2,3,...)
    /// Optional argument can give a non-zero starting point and non-1 step.
    static const vfloat16 Iota (float start=0.0f, float step=1.0f);

    /// Set all components to 0.0
    void clear ();

    /// Component access (get)
    float operator[] (int i) const;
    /// Component access (set)
    float& operator[] (int i);

    /// Component access (set).
    void setcomp (int i, float value);

    value_t x () const;
    value_t y () const;
    value_t z () const;
    value_t w () const;
    void set_x (value_t val);
    void set_y (value_t val);
    void set_z (value_t val);
    void set_w (value_t val);

    /// Extract the lower precision vfloat8
    vfloat8 lo () const;

    /// Extract the higher precision vfloat8
    vfloat8 hi () const;

    /// Helper: load a single value into all components
    void load (float val);

    /// Load separate values into each component.
    void load (float v0, float v1, float v2, float v3,
               float v4, float v5, float v6, float v7,
               float v8, float v9, float v10, float v11,
               float v12, float v13, float v14, float v15);

    /// Load from an array of values
    void load (const float *values);

    /// Load from a partial array of <=16 values. Unassigned values are
    /// undefined.
    void load (const float *values, int n);

    /// Load from an array of 16 unsigned short values, convert to float
    void load (const unsigned short *values);

    /// Load from an array of 16 short values, convert to float
    void load (const short *values);

    /// Load from an array of 16 unsigned char values, convert to float
    void load (const unsigned char *values);

    /// Load from an array of 16 char values, convert to float
    void load (const char *values);

#ifdef _HALF_H_
    /// Load from an array of 16 half values, convert to float
    void load (const half *values);
#endif /* _HALF_H_ */

    void store (float *values) const;

    /// Store the first n values into memory
    void store (float *values, int n) const;

#ifdef _HALF_H_
    void store (half *values) const;
#endif

    /// Masked load -- read from values[] where mask is 1, load zero where
    /// mask is 0.
    void load_mask (const vbool_t &mask, const value_t *values);
    void load_mask (int mask, const value_t *values) { load_mask(vbool_t(mask), values); }

    /// Masked store -- write to values[] where mask is enabled, don't
    /// touch values[] where it's not.
    void store_mask (const vbool_t &mask, value_t *values) const;
    void store_mask (int mask, value_t *values) const { store_mask(vbool_t(mask), values); }

    /// Load values from addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void gather (const value_t *baseptr, const vint_t& vindex);
    /// Gather elements defined by the mask, leave others unchanged.
    template<int scale=4>
    void gather_mask (const bool_t& mask, const value_t *baseptr, const vint_t& vindex);
    template<int scale=4>
    void gather_mask (int mask, const value_t *baseptr, const vint_t& vindex) {
        gather_mask<scale> (vbool_t(mask), baseptr, vindex);
    }

    /// Store values at addresses  (char*)basepatr + vindex[i]*scale
    template<int scale=4>
    void scatter (value_t *baseptr, const vint_t& vindex) const;
    /// Scatter elements defined by the mask
    template<int scale=4>
    void scatter_mask (const bool_t& mask, value_t *baseptr, const vint_t& vindex) const;
    template<int scale=4>
    void scatter_mask (int mask, value_t *baseptr, const vint_t& vindex) const {
        scatter_mask<scale> (vbool_t(mask), baseptr, vindex);
    }

    // Arithmetic operators (component-by-component)
    friend vfloat16 operator+ (const vfloat16& a, const vfloat16& b);
    friend vfloat16 operator- (const vfloat16& a);
    friend vfloat16 operator- (const vfloat16& a, const vfloat16& b);
    friend vfloat16 operator* (const vfloat16& a, const vfloat16& b);
    friend vfloat16 operator* (const vfloat16& a, float b);
    friend vfloat16 operator* (float a, const vfloat16& b);
    friend vfloat16 operator/ (const vfloat16& a, const vfloat16& b);
    friend vfloat16 operator% (const vfloat16& a, const vfloat16& b);
    friend const vfloat16 & operator+= (vfloat16& a, const vfloat16& b);
    friend const vfloat16 & operator-= (vfloat16& a, const vfloat16& b);
    friend const vfloat16 & operator*= (vfloat16& a, const vfloat16& b);
    friend const vfloat16 & operator/= (vfloat16& a, const vfloat16& b);

    // Comparison operations
    friend vbool16 operator== (const vfloat16& a, const vfloat16& b);
    friend vbool16 operator!= (const vfloat16& a, const vfloat16& b);
    friend vbool16 operator<  (const vfloat16& a, const vfloat16& b);
    friend vbool16 operator>  (const vfloat16& a, const vfloat16& b);
    friend vbool16 operator>= (const vfloat16& a, const vfloat16& b);
    friend vbool16 operator<= (const vfloat16& a, const vfloat16& b);

    // Some oddball items that are handy

    /// Stream output
    friend inline std::ostream& operator<< (std::ostream& cout, const vfloat16& val);

protected:
    // The actual data representation
    union {
        simd_t  m_simd;
        value_t m_val[paddedelements];
        vfloat8  m_8[2];
    };
};


/// Shuffle groups of 4
template<int i0, int i1, int i2, int i3>
vfloat16 shuffle4 (const vfloat16& a);

/// shuffle4<i>(a) is the same as shuffle4<i,i,i,i>(a)
template<int i> vfloat16 shuffle4 (const vfloat16& a);

/// Shuffle within each group of 4
template<int i0, int i1, int i2, int i3>
vfloat16 shuffle (const vfloat16& a);

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> vfloat16 shuffle (const vfloat16& a);

/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i> float extract (const vfloat16& a);

/// Helper: substitute val for a[i]
template<int i> vfloat16 insert (const vfloat16& a, float val);

/// The sum of all components, returned in all components.
vfloat16 vreduce_add (const vfloat16& v);

/// The sum of all components, returned as a scalar.
float reduce_add (const vfloat16& v);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// and b (if mask[i] is true), i.e., mask[i] ? b[i] : a[i].
vfloat16 blend (const vfloat16& a, const vfloat16& b, const vbool4& mask);

/// Use a bool mask to select between `a` (if mask[i] is true) or 0 if
/// mask[i] is false), i.e., mask[i] ? a[i] : 0. Equivalent to
/// blend(0,a,mask).
vfloat16 blend0 (const vfloat16& a, const vbool4& mask);

/// Use a bool mask to select between components of a (if mask[i] is false)
/// or 0 (if mask[i] is true), i.e., mask[i] ? 0 : a[i]. Equivalent to
/// blend(0,a,!mask), or blend(a,0,mask).
vfloat16 blend0not (const vfloat16& a, const vbool4& mask);

/// "Safe" divide of vfloat16/vfloat16 -- for any component of the divisor
/// that is 0, return 0 rather than Inf.
vfloat16 safe_div (const vfloat16 &a, const vfloat16 &b);

/// Select 'a' where mask is true, 'b' where mask is false. Sure, it's a
/// synonym for blend with arguments rearranged, but this is more clear
/// because the arguments are symmetric to scalar (cond ? a : b).
vfloat16 select (const vbool16& mask, const vfloat16& a, const vfloat16& b);

// Per-element math
vfloat16 abs (const vfloat16& a);    ///< absolute value (float)
vfloat16 sign (const vfloat16& a);   ///< 1.0 when value >= 0, -1 when negative
vfloat16 ceil (const vfloat16& a);
vfloat16 floor (const vfloat16& a);
vint16 ifloor (const vfloat16& a);    ///< (int)floor
inline vint16 floori (const vfloat16& a) { return ifloor(a); }  // DEPRECATED(1.8) alias

/// Per-element round to nearest integer.
/// CAVEAT: the rounding when mid-way between integers may differ depending
/// on hardware. Intel SSE/AVX does "banker's founding" (to nearest even
/// integer) but std::round() says to round away from 0 regardless of
/// current rounding mode (but that is multiple instructions on x64).
/// USE WITH CAUTION, and maybe avoid this if it is critical to exactly
/// match std::round().
vfloat16 round (const vfloat16& a);

/// Per-element round to nearest integer (equivalent to vint(round(a))).
/// CAVEAT: On SSE/AVX this uses banker's rounding, which may differ from
/// C++ std::rint() which says to use the current rounding mode.
/// USE WITH CAUTION, and maybe avoid this if it is critical to exactly
/// match std::rint().
vint16 rint (const vfloat16& a);

vfloat16 rcp_fast (const vfloat16 &a);  ///< Fast, approximate 1/a
vfloat16 sqrt (const vfloat16 &a);
vfloat16 rsqrt (const vfloat16 &a);   ///< Fully accurate 1/sqrt
vfloat16 rsqrt_fast (const vfloat16 &a);  ///< Fast, approximate 1/sqrt
vfloat16 min (const vfloat16& a, const vfloat16& b); ///< Per-element min
vfloat16 max (const vfloat16& a, const vfloat16& b); ///< Per-element max
// vfloat16 exp (const vfloat16& v);  // See template with vfloat4
// vfloat16 log (const vfloat16& v);  // See template with vfloat4

/// andnot(a,b) returns ((~a) & b)
vfloat16 andnot (const vfloat16& a, const vfloat16& b);

// Fused multiply and add (or subtract):
vfloat16 madd (const vfloat16& a, const vfloat16& b, const vfloat16& c); // a*b + c
vfloat16 msub (const vfloat16& a, const vfloat16& b, const vfloat16& c); // a*b - c
vfloat16 nmadd (const vfloat16& a, const vfloat16& b, const vfloat16& c); // -a*b + c
vfloat16 nmsub (const vfloat16& a, const vfloat16& b, const vfloat16& c); // -a*b - c



// Odds and ends, other CPU hardware tricks

// Try to set the flush_zero_mode CPU flag on x86. Return true if we are
// able, otherwise false (because it's not available on that platform,
// or because it's gcc 4.8 which has a bug that lacks this intrinsic).
inline bool set_flush_zero_mode (bool on) {
#if (defined(__x86_64__) || defined(__i386__)) && (OIIO_GNUC_VERSION == 0 || OIIO_GNUC_VERSION > 40900)
    _MM_SET_FLUSH_ZERO_MODE (on ? _MM_FLUSH_ZERO_ON : _MM_FLUSH_ZERO_OFF);
    return true;
#endif
    return false;
}

// Try to set the denorms_zero_mode CPU flag on x86. Return true if we are
// able, otherwise false (because it's not available on that platform,
// or because it's gcc 4.8 which has a bug that lacks this intrinsic).
inline bool set_denorms_zero_mode (bool on) {
#if (defined(__x86_64__) || defined(__i386__)) && (OIIO_GNUC_VERSION == 0 || OIIO_GNUC_VERSION > 40900)
    _MM_SET_DENORMALS_ZERO_MODE (on ? _MM_DENORMALS_ZERO_ON : _MM_DENORMALS_ZERO_OFF);
    return true;
#endif
    return false;
}

// Get the flush_zero_mode CPU flag on x86.
inline bool get_flush_zero_mode () {
#if (defined(__x86_64__) || defined(__i386__)) && (OIIO_GNUC_VERSION == 0 || OIIO_GNUC_VERSION > 40900)
    return _MM_GET_FLUSH_ZERO_MODE() == _MM_FLUSH_ZERO_ON;
#endif
    return false;
}

// Get the denorms_zero_mode CPU flag on x86.
inline bool get_denorms_zero_mode () {
#if (defined(__x86_64__) || defined(__i386__)) && (OIIO_GNUC_VERSION == 0 || OIIO_GNUC_VERSION > 40900)
    return _MM_GET_DENORMALS_ZERO_MODE() == _MM_DENORMALS_ZERO_ON;
#endif
    return false;
}






//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//
// Gory implementation details follow.
//
// ^^^ All declarations and documention is above ^^^
//
// vvv Below is the implementation, often considerably cluttered with
//     #if's for each architeture, and unapologitic use of intrinsics and
//     every manner of dirty trick we can think of to make things fast.
//     Some of this isn't pretty. We won't recapitulate comments or
//     documentation of what the functions are supposed to do, please
//     consult the declarations above for that.
//
//     Here be dragons.
//
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////



//////////////////////////////////////////////////////////////////////
// vbool4 implementation


OIIO_FORCEINLINE int vbool4::operator[] (int i) const {
    OIIO_DASSERT(i >= 0 && i < elements);
#if OIIO_SIMD_SSE
    return ((_mm_movemask_ps(m_simd) >> i) & 1) ? -1 : 0;
#else
    return m_val[i];
#endif
}

OIIO_FORCEINLINE int& vbool4::operator[] (int i) {
    OIIO_DASSERT(i >= 0 && i < elements);
    return m_val[i];
}


OIIO_FORCEINLINE void vbool4::setcomp (int i, bool value) {
    OIIO_DASSERT(i >= 0 && i < elements);
    m_val[i] = value ? -1 : 0;
}


OIIO_FORCEINLINE std::ostream& operator<< (std::ostream& cout, const vbool4& a) {
    cout << a[0];
    for (int i = 1; i < a.elements; ++i)
        cout << ' ' << a[i];
    return cout;
}


OIIO_FORCEINLINE void vbool4::load (bool a) {
#if OIIO_SIMD_SSE
    m_simd = _mm_castsi128_ps(_mm_set1_epi32(-int(a)));
#elif OIIO_SIMD_NEON
    m_simd = vdupq_n_u32(a ? 0xffffffff : 0);
#else
    int val = -int(a);
    SIMD_CONSTRUCT (val);
#endif
}


OIIO_FORCEINLINE void vbool4::load (bool a, bool b, bool c, bool d) {
#if OIIO_SIMD_SSE
    // N.B. -- we need to reverse the order because of our convention
    // of storing a,b,c,d in the same order in memory.
    m_simd = _mm_castsi128_ps(_mm_set_epi32(-int(d), -int(c), -int(b), -int(a)));
// #elif OIIO_SIMD_NEON
// FIXME
#else
    m_val[0] = -int(a);
    m_val[1] = -int(b);
    m_val[2] = -int(c);
    m_val[3] = -int(d);
#endif
}

OIIO_FORCEINLINE vbool4::vbool4 (const bool *a) {
    load (a[0], a[1], a[2], a[3]);
}

OIIO_FORCEINLINE const vbool4& vbool4::operator= (const vbool4 & other) {
    m_simd = other.m_simd;
    return *this;
}


OIIO_FORCEINLINE int vbool4::bitmask () const {
#if OIIO_SIMD_SSE
    return _mm_movemask_ps(m_simd);
#else
    int r = 0;
    for (int i = 0; i < elements; ++i)
        if (m_val[i])
            r |= 1<<i;
    return r;
#endif
}


OIIO_FORCEINLINE vbool4
vbool4::from_bitmask (int bitmask) {
    // I think this is a fast conversion from int bitmask to vbool4
    return (vint4::Giota() & vint4(bitmask)) != vint4::Zero();
}


OIIO_FORCEINLINE void vbool4::clear () {
#if OIIO_SIMD_SSE
    m_simd = _mm_setzero_ps();
#else
    *this = false;
#endif
}


OIIO_FORCEINLINE const vbool4 vbool4::False () {
#if OIIO_SIMD_SSE
    return _mm_setzero_ps();
#else
    return false;
#endif
}

OIIO_FORCEINLINE const vbool4 vbool4::True () {
    // Fastest way to fill with all 1 bits is to cmp any value to itself.
#if OIIO_SIMD_SSE
# if OIIO_SIMD_AVX && (OIIO_GNUC_VERSION > 50000)
    __m128i anyval = _mm_undefined_si128();
# else
    __m128i anyval = _mm_setzero_si128();
# endif
    return _mm_castsi128_ps (_mm_cmpeq_epi8 (anyval, anyval));
#else
    return true;
#endif
}

OIIO_FORCEINLINE void vbool4::store (bool *values) const {
    SIMD_DO (values[i] = m_val[i] ? true : false);
}

OIIO_FORCEINLINE void vbool4::store (bool *values, int n) const {
    OIIO_DASSERT (n >= 0 && n <= elements);
    for (int i = 0; i < n; ++i)
        values[i] = m_val[i] ? true : false;
}



OIIO_FORCEINLINE vbool4 operator! (const vbool4 & a) {
#if OIIO_SIMD_SSE
    return _mm_xor_ps (a.simd(), vbool4::True());
#elif OIIO_SIMD_NEON
    return vmvnq_u32(a.simd());
#else
    SIMD_RETURN (vbool4, a[i] ^ (-1));
#endif
}

OIIO_FORCEINLINE vbool4 operator& (const vbool4 & a, const vbool4 & b) {
#if OIIO_SIMD_SSE
    return _mm_and_ps (a.simd(), b.simd());
#elif OIIO_SIMD_NEON
    return vandq_u32(a.simd(), b.simd());
#else
    SIMD_RETURN (vbool4, a[i] & b[i]);
#endif
}

OIIO_FORCEINLINE vbool4 operator| (const vbool4 & a, const vbool4 & b) {
#if OIIO_SIMD_SSE
    return _mm_or_ps (a.simd(), b.simd());
#elif OIIO_SIMD_NEON
    return vorrq_u32(a.simd(), b.simd());
#else
    SIMD_RETURN (vbool4, a[i] | b[i]);
#endif
}

OIIO_FORCEINLINE vbool4 operator^ (const vbool4& a, const vbool4& b) {
#if OIIO_SIMD_SSE
    return _mm_xor_ps (a.simd(), b.simd());
#elif OIIO_SIMD_NEON
    return veorq_u32(a.simd(), b.simd());
#else
    SIMD_RETURN (vbool4, a[i] ^ b[i]);
#endif
}


OIIO_FORCEINLINE const vbool4& operator&= (vbool4& a, const vbool4 &b) {
    return a = a & b;
}

OIIO_FORCEINLINE const vbool4& operator|= (vbool4& a, const vbool4& b) {
    return a = a | b;
}

OIIO_FORCEINLINE const vbool4& operator^= (vbool4& a, const vbool4& b) {
    return a = a ^ b;
}

OIIO_FORCEINLINE vbool4 operator~ (const vbool4& a) {
#if OIIO_SIMD_SSE
    // Fastest way to bit-complement in SSE is to xor with 0xffffffff.
    return _mm_xor_ps (a.simd(), vbool4::True());
#elif OIIO_SIMD_NEON
    return vmvnq_u32(a.m_simd);
#else
    SIMD_RETURN (vbool4, ~a[i]);
#endif
}

OIIO_FORCEINLINE vbool4 operator== (const vbool4 & a, const vbool4 & b) {
#if OIIO_SIMD_SSE
    return _mm_castsi128_ps (_mm_cmpeq_epi32 (_mm_castps_si128 (a), _mm_castps_si128(b)));
#elif OIIO_SIMD_NEON
    return vceqq_u32 (a.m_simd, b.m_simd);
#else
    SIMD_RETURN (vbool4, a[i] == b[i] ? -1 : 0);
#endif
}

OIIO_FORCEINLINE vbool4 operator!= (const vbool4 & a, const vbool4 & b) {
#if OIIO_SIMD_SSE
    return _mm_xor_ps (a, b);
#elif OIIO_SIMD_NEON
    return !(a == b);
#else
    SIMD_RETURN (vbool4, a[i] != b[i] ? -1 : 0);
#endif
}




#if OIIO_SIMD_SSE
// Shuffling. Use like this:  x = shuffle<3,2,1,0>(b)
template<int i0, int i1, int i2, int i3>
OIIO_FORCEINLINE __m128i shuffle_sse (__m128i v) {
    return _mm_shuffle_epi32(v, _MM_SHUFFLE(i3, i2, i1, i0));
}
#endif

#if OIIO_SIMD_SSE >= 3
// SSE3 has intrinsics for a few special cases
template<> OIIO_FORCEINLINE __m128i shuffle_sse<0, 0, 2, 2> (__m128i a) {
    return _mm_castps_si128(_mm_moveldup_ps(_mm_castsi128_ps(a)));
}
template<> OIIO_FORCEINLINE __m128i shuffle_sse<1, 1, 3, 3> (__m128i a) {
    return _mm_castps_si128(_mm_movehdup_ps(_mm_castsi128_ps(a)));
}
template<> OIIO_FORCEINLINE __m128i shuffle_sse<0, 1, 0, 1> (__m128i a) {
    return _mm_castpd_si128(_mm_movedup_pd(_mm_castsi128_pd(a)));
}
#endif

#if OIIO_SIMD_SSE
template<int i0, int i1, int i2, int i3>
OIIO_FORCEINLINE __m128 shuffle_sse (__m128 a) {
    return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(a), _MM_SHUFFLE(i3, i2, i1, i0)));
}
#endif

#if OIIO_SIMD_SSE >= 3
// SSE3 has intrinsics for a few special cases
template<> OIIO_FORCEINLINE __m128 shuffle_sse<0, 0, 2, 2> (__m128 a) {
    return _mm_moveldup_ps(a);
}
template<> OIIO_FORCEINLINE __m128 shuffle_sse<1, 1, 3, 3> (__m128 a) {
    return _mm_movehdup_ps(a);
}
template<> OIIO_FORCEINLINE __m128 shuffle_sse<0, 1, 0, 1> (__m128 a) {
    return _mm_castpd_ps(_mm_movedup_pd(_mm_castps_pd(a)));
}
#endif


/// Helper: shuffle/swizzle with constant (templated) indices.
/// Example: shuffle<1,1,2,2>(vbool4(a,b,c,d)) returns (b,b,c,c)
template<int i0, int i1, int i2, int i3>
OIIO_FORCEINLINE vbool4 shuffle (const vbool4& a) {
#if OIIO_SIMD_SSE
    return shuffle_sse<i0,i1,i2,i3> (a.simd());
#else
    return vbool4 (a[i0], a[i1], a[i2], a[i3]);
#endif
}

/// shuffle<i>(a) is the same as shuffle<i,i,i,i>(a)
template<int i> OIIO_FORCEINLINE vbool4 shuffle (const vbool4& a) {
    return shuffle<i,i,i,i>(a);
}


/// Helper: as rapid as possible extraction of one component, when the
/// index is fixed.
template<int i>
OIIO_FORCEINLINE bool extract (const vbool4& a) {
#if OIIO_SIMD_SSE >= 4
    return _mm_extract_epi32(_mm_castps_si128(a.simd()), i);  // SSE4.1 only
#else
    return a[i];
#endif
}

/// Helper: substitute val for a[i]
template<int i>
OIIO_FORCEINLINE vbool4 insert (const vbool4& a, bool val) {
#if OIIO_SIMD_SSE >= 4
    int ival = -int(val);
    return _mm_castsi128_ps (_mm_insert_epi32 (_mm_castps_si128(a), ival, i));
#else
    vbool4 tmp = a;
    tmp[i] = -int(val);
    return tmp;
#endif
}

OIIO_FORCEINLINE bool reduce_and (const vbool4& v) {
#if OIIO_SIMD_AVX
    return _mm_testc_ps (v, vbool4(true)) != 0;
#elif OIIO_SIMD_SSE
    return _mm_movemask_ps(v.simd()) == 0xf;
#else
    SIMD_RETURN_REDUCE (bool, true, r &= (v[i] != 0));
#endif
}

OIIO_FORCEINLINE bool reduce_or (const vbool4& v) {
#if OIIO_SIMD_AVX
    return ! _mm_testz_ps (v, v);
#elif OIIO_SIMD_SSE
    return _mm_movemask_ps(v) != 0;
#else
    SIMD_RETURN_REDUCE (bool, false, r |= (v[i] != 0));
#endif
}

OIIO_FORCEINLINE bool all (const vbool4& v) { return reduce_and(v) == true; }
OIIO_FORCEINLINE bool any (const vbool4& v) { return reduce_or(v) == true; }
OIIO_FORCEINLINE bool none (const vbool4& v) { return reduce_or(v) == false; }



//////////////////////////////////////////////////////////////////////
// vbool8 implementation


OIIO_FORCEINLINE int vbool8::operator[] (int i) const {
    OIIO_DASSERT(i >= 0 && i < elements);
#if OIIO_SIMD_AVX
    return ((_mm256_movemask_ps(m_simd) >> i) & 1) ? -1 : 0;
#else
    return m_val[i];
#endif
}

OIIO_FORCEINLINE void vbool8::setcomp (int i, bool value) {
    OIIO_DASSERT(i >= 0 && i < elements);
    m_val[i] = value ? -1 : 0;
}

OIIO_FORCEINLINE int& vbool8::operator[] (int i) {
    OIIO_DASSERT(i >= 0 && i < elements);
    return m_val[i];
}


OIIO_FORCEINLINE std::ostream& operator<< (std::ostream& cout, const vbool8& a) {
    cout << a[0];
    for (int i = 1; i < a.elements; ++i)
        cout << ' ' << a[i];
    return cout;
}


OIIO_FORCEINLINE void vbool8::load (bool a) {
#if OIIO_SIMD_AVX
    m_simd = _mm256_castsi256_ps(_mm256_set1_epi32(-int(a)));
#elif OIIO_SIMD_SSE || OIIO_SIMD_NEON
    m_4[0].load(a);
    m_4[1].load(a);
#else
    int val = -int(a);
    SIMD_CONSTRUCT (val);
#endif
}


OIIO_FORCEINLINE void vbool8::load (bool a, bool b, bool c, bool d,
                                   bool e, bool f, bool g, bool h) {
#if OIIO_SIMD_AVX
    // N.B. -- we need to reverse the order because of our convention
    // of storing a,b,c,d in the same order in memory.
    m_simd = _mm256_castsi256_ps(_mm256_set_epi32(-int(h), -int(g), -int(f), -int(e),
                                                  -int(d), -int(c), -int(b), -int(a)));
#elif OIIO_SIMD_SSE || OIIO_SIMD_NEON
    m_4[0].load(a, b, c, d);
    m_4[1].load(e, f, g, h);
#else
    m_val[0] = -int(a);
    m_val[1] = -int(b);
    m_val[2] = -int(c);
    m_val[3] = -int(d);
    m_val[4] = -int(e);
    m_val[5] = -int(f);
    m_val[6] = -int(g);
    m_val[7] = -int(h);
#endif
}

OIIO_FORCEINLINE vbool8::vbool8 (bool a, bool b, bool c, bool d,
                               bool e, bool f, bool g, bool h) {
    load (a, b, c, d, e, f, g, h);
}

OIIO_FORCEINLINE vbool8::vbool8 (int a, int b, int c, int d,
                                 int e, int f, int g, int h) {
    load (bool(a), bool(b), bool(c), bool(d),
          bool(e), bool(f), bool(g), bool(h));
}

OIIO_FORCEINLINE vbool8::vbool8 (const bool *a) {
    load (a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7]);
}


OIIO_FORCEINLINE const vbool8& vbool8::operator= (bool a) {
    load(a);
    return *this;
}

OIIO_FORCEINLINE const vbool8& vbool8::operator= (const vbool8 & other) {
    m_simd = other.m_simd;
    return *this;
}

OIIO_FORCEINLINE int vbool8::bitmask () const {
#if OIIO_SIMD_AVX
    return _mm256_movemask_ps(m_simd);
#else
    return lo().bitmask() | (hi().bitmask() << 4);
#endif
}


OIIO_FORCEINLINE vbool8
vbool8::from_bitmask (int bitmask) {
    // I think this is a fast conversion from int bitmask to vbool8
    return (vint8::Giota() & vint8(bitmask)) != vint8::Zero();
}


OIIO_FORCEINLINE void vbool8::clear () {
#if OIIO_SIMD_AVX
    m_simd = _mm256_setzero_ps();
#else
    *this = false;
#endif
}

OIIO_FORCEINLINE const vbool8 vbool8::False () {
#if OIIO_SIMD_AVX
    return _mm256_setzero_ps();
#else
    return false;
#endif
}


OIIO_FORCEINLINE const vbool8 vbool8::True () {
#if OIIO_SIMD_AVX
# if OIIO_SIMD_AVX >= 2 && (OIIO_GNUC_VERSION > 50000)
    // Fastest way to fill with all 1 bits is to cmp any value to itself.
    __m256i anyval = _mm256_undefined_si256();
    return _mm256_castsi256_ps (_mm256_cmpeq_epi8 (anyval, anyval));
# else
    return _mm256_castsi256_ps (_mm256_set1_epi32 (-1));
# endif
#else
    return true;
#endif
}


OIIO_FORCEINLINE void vbool8::store (bool *values) const {
    SIMD_DO (values[i] = m_val[i] ? true : false);
}

OIIO_FORCEINLINE void vbool8::store (bool *values, int n) const {
    OIIO_DASSERT (n >= 0 && n <= elements);
    for (int i = 0; i < n; ++i)
        values[i] = m_val[i] ? true : false;
}


OIIO_FORCEINLINE vbool4 vbool8::lo () const {
#if OIIO_SIMD_AVX
    return _mm256_castps256_ps128 (simd());
#else
    return m_4[0];
#endif
}

OIIO_FORCEINLINE vbool4 vbool8::hi () const {
#if OIIO_SIMD_AVX
    return _mm256_extractf128_ps (simd(), 1);
#else
    return m_4[1];
#endif
}


OIIO_FORCEINLINE vbool8::vbool8 (const vbool4& lo, const vbool4 &hi) {
#if OIIO_SIMD_AVX
    __m256 r = _mm256_castps128_ps256 (lo);
    m_simd = _mm256_insertf128_ps (r, hi, 1);
    // N.B. equivalent, if available: m_simd = _mm256_set_m128 (hi, lo);
#else
    m_4[0] = lo;
    m_4[1] = hi;
#endif
}


OIIO_FORCEINLINE vbool8 operator! (const vbool8 & a) {
#if OIIO_SIMD_AVX
    return _mm256_xor_ps (a.simd(), vbool8::True());
#else
    SIMD_RETURN (vbool8, a[i] ^ (-1));
#endif
}

OIIO_FORCEINLINE vbool8 operator& (const vbool8 & a, const vbool8 & b) {
#if OIIO_SIMD_AVX
    return _mm256_and_ps (a.simd(), b.simd());
#else
    SIMD_RETURN (vbool8, a[i] & b[i]);
#endif
}

OIIO_FORCEINLINE vbool8 operator| (const vbool8 & a, const vbool8 & b) {
#if OIIO_SIMD_AVX
    return _mm256_or_ps (a.simd(), b.simd());
#else
    SIMD_RETURN (vbool8, a[i] | b[i]);
#endif
}

OIIO_FORCEINLINE vbool8 operator^ (const vbool8& a, const vbool8& b) {
#if OIIO_SIMD_AVX
    return _mm256_xor_ps (a.simd(), b.simd());
#else
    SIMD_RETURN (vbool8, a[i] ^ b[i]);
#endif
}


OIIO_FORCEINLINE const vbool8& operator&= (vbool8& a, const vbool8 &b) {
    return a = a & b;
}

OIIO_FORCEINLINE const vbool8& operator|= (vbool8& a, const vbool8& b) {
    return a = a | b;
}

OIIO_FORCEINLINE const vbool8& operator^= (vbool8& a, const vbool8& b) {
    return a = a ^ b;
}


OIIO_FORCEINLINE vbool8 operator~ (const vbool8& a) {
#if OIIO_SIMD_AVX
    // Fastest way to bit-complement in SSE is to xor with 0xffffffff.
    return _mm256_xor_ps (a.simd(), vbool8::True());
#else
    SIMD_RETURN (vbool8, ~a[i]);
#endif
}


OIIO_FORCEINLINE vbool8 operator== (const vbool8 & a, const vbool8 & b) {
#if OIIO_SIMD_AVX >= 2
    return _mm256_castsi256_ps (_mm256_cmpeq_epi32 (_mm256_castps_si256 (a), _mm256_castps_si256(b)));
#elif OIIO_SIMD_AVX
    return _mm256_cmp_ps (a, b, _CMP_EQ_UQ);
#else
    SIMD_RETURN (vbool8, a[i] == b[i] ? -1 : 0);
#endif
}

OIIO_FORCEINLINE vbool8 operator!= (const vbool8 & a, const vbool8 & b) {
#if OIIO_SIMD_AVX
    return _mm256_xor_ps (a, b);
#else
    SIMD_RETURN (vbool8, a[i] != b[i] ? -1 : 0);
#endif
}


template<int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7>
OIIO_FORCEINLINE vbool8 shuffle (const vbool8& a) {
#if OIIO_SIMD_AVX >= 2
    vint8 index (i0, i1, i2, i3, i4, i5, i6, i7);
    return _mm256_permutevar8x32_ps (a.simd(), index.simd());
#else
    return vbool8 (a[i0], a[i1], a[i2], a[i3], a[i4], a[i5], a[i6], a[i7]);
#endif
}

template<int i> OIIO_FORCEINLINE vbool8 shuffle (const vbool8& a) {
    return shuffle<i,i,i,i,i,i,i,i>(a);
}


template<int i>
OIIO_FORCEINLINE bool extract (const vbool8& a) {
#if OIIO_SIMD_AVX && !_WIN32
    return _mm256_extract_epi32(_mm256_castps_si256(a.simd()), i);  // SSE4.1 only
#else
    return a[i];
#endif
}

template<int i>
OIIO_FORCEINLINE vbool8 insert (const vbool8& a, bool val) {
#if OIIO_SIMD_AVX && !_WIN32
    int ival = -int(val);
    return _mm256_castsi256_ps (_mm256_insert_epi32 (_mm256_castps_si256(a.simd()), ival, i));
#else
    vbool8 tmp = a;
    tmp[i] = -int(val);
    return tmp;
#endif
}


OIIO_FORCEINLINE bool reduce_and (const vbool8& v) {
#if OIIO_SIMD_AVX
    return _mm256_testc_ps (v, vbool8(true)) != 0;
    // return _mm256_movemask_ps(v.simd()) == 0xff;
#else
    SIMD_RETURN_REDUCE (bool, true, r &= bool(v[i]));
#endif
}

OIIO_FORCEINLINE bool reduce_or (const vbool8& v) {
#if OIIO_SIMD_AVX
    return ! _mm256_testz_ps (v, v);   // FIXME? Not in all immintrin.h !
    // return _mm256_movemask_ps(v) != 0;
#else
    SIMD_RETURN_REDUCE (bool, false, r |= bool(v[i]));
#endif
}


OIIO_FORCEINLINE bool all (const vbool8& v) { return reduce_and(v) == true; }
OIIO_FORCEINLINE bool any (const vbool8& v) { return reduce_or(v) == true; }
OIIO_FORCEINLINE bool none (const vbool8& v) { return reduce_or(v) == false; }



//////////////////////////////////////////////////////////////////////
// vbool16 implementation


OIIO_FORCEINLINE int vbool16::operator[] (int i) const {
    OIIO_DASSERT(i >= 0 && i < elements);
#if OIIO_SIMD_AVX >= 512
    return (int(m_simd) >> i) & 1;
#else
    return (m_bits >> i) & 1;
#endif
}

OIIO_FORCEINLINE void vbool16::setcomp (int i, bool value) {
    OIIO_DASSERT(i >= 0 && i < elements);
    int bits = m_bits;
    bits &= (0xffff ^ (1<<i));
    bits |= (int(value)<<i);
    m_bits = bits;
}


OIIO_FORCEINLINE std::ostream& operator<< (std::ostream& cout, const vbool16& a) {
    cout << a[0];
    for (int i = 1; i < a.elements; ++i)
        cout << ' ' << a[i];
    return cout;
}


OIIO_FORCEINLINE void vbool16::load (bool a) {
    m_simd = a ? 0xffff : 0;
}


OIIO_FORCEINLINE void vbool16::load_bitmask (int a) {
    m_simd = simd_t(a);
}


OIIO_FORCEINLINE void vbool16::load (bool v0, bool v1, bool v2, bool v3,
                                    bool v4, bool v5, bool v6, bool v7,
                                    bool v8, bool v9, bool v10, bool v11,
                                    bool v12, bool v13, bool v14, bool v15) {
    m_simd = simd_t((int(v0) << 0) |
                    (int(v1) << 1) |
                    (int(v2) << 2) |
                    (int(v3) << 3) |
                    (int(v4) << 4) |
                    (int(v5) << 5) |
                    (int(v6) << 6) |
                    (int(v7) << 7) |
                    (int(v8) << 8) |
                    (int(v9) << 9) |
                    (int(v10) << 10) |
                    (int(v11) << 11) |
                    (int(v12) << 12) |
                    (int(v13) << 13) |
                    (int(v14) << 14) |
                    (int(v15) << 15));
}

OIIO_FORCEINLINE vbool16::vbool16 (bool v0, bool v1, bool v2, bool v3,
                                 bool v4, bool v5, bool v6, bool v7,
                                 bool v8, bool v9, bool v10, bool v11,
                                 bool v12, bool v13, bool v14, bool v15) {
    load (v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13, v14, v15);
}

OIIO_FORCEINLINE vbool16::vbool16 (int v0, int v1, int v2, int v3,
                                   int v4, int v5, int v6, int v7,
                                   int v8, int v9, int v10, int v11,
                                   int v12, int v13, int v14, int v15) {
    load (bool(v0), bool(v1), bool(v2), bool(v3),
          bool(v4), bool(v5), bool(v6), bool(v7),
          bool(v8), bool(v9), bool(v10), bool(v11),
          bool(v12), bool(v13), bool(v14), bool(v15));
}

OIIO_FORCEINLINE vbool16::vbool16 (const vbool8& a, const vbool8& b) {
    load_bitmask (a.bitmask() | (b.bitmask() << 8));
}

OIIO_FORCEINLINE vbool16::vbool16 (const bool *a) {
    load (a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7],
          a[8], a[9], a[10], a[11], a[12], a[13], a[14], a[15]);
}


OIIO_FORCEINLINE const vbool16& vbool16::operator= (bool a) {
    load(a);
    return *this;
}

OIIO_FORCEINLINE const vbool16& vbool16::operator= (const vbool16 & other) {
    m_simd = other.m_simd;
    return *this;
}


OIIO_FORCEINLINE int vbool16::bitmask () const {
#if OIIO_SIMD_AVX >= 512
    return int(m_simd);
#else
    return int(m_bits);
#endif
}


OIIO_FORCEINLINE void vbool16::clear () {
    m_simd = simd_t(0);
}

OIIO_FORCEINLINE const vbool16 vbool16::False () {
    return simd_t(0);
}


OIIO_FORCEINLINE const vbool16 vbool16::True () {
    return simd_t(0xffff);
}


OIIO_FORCEINLINE void vbool16::store (bool *values) const {
    SIMD_DO (values[i] = m_bits & (1<<i));
}

OIIO_FORCEINLINE void vbool16::store (bool *values, int n) const {
    OIIO_DASSERT (n >= 0 && n <= elements);
    for (int i = 0; i < n; ++i)
        values[i] = m_bits & (1<<i);
}



OIIO_FORCEINLINE vbool8 vbool16::lo () const {
#if OIIO_SIMD_AVX >= 512
    return _mm256_castsi256_ps (_mm256_maskz_set1_epi32 (bitmask()&0xff, -1));
#else
    SIMD_RETURN (vbool8, (*this)[i] ? -1 : 0);
#endif
}

OIIO_FORCEINLINE vbool8 vbool16::hi () const {
#if OIIO_SIMD_AVX >= 512
    return _mm256_castsi256_ps (_mm256_maskz_set1_epi32 (bitmask()>>8, -1));
#else
    SIMD_RETURN (vbool8, (*this)[i+8] ? -1 : 0);
#endif
}


OIIO_FORCEINLINE vbool16 operator! (const vbool16 & a) {
#if OIIO_SIMD_AVX >= 512
    return _mm512_knot (a.simd());
#else
    return vbool16 (a.m_bits ^ 0xffff);
#endif
}

OIIO_FORCEINLINE vbool16 operator& (const vbool16 & a, const vbool16 & b) {
#if OIIO_SIMD_AVX >= 512
    return _mm512_kand (a.simd(), b.simd());
#else
    return vbool16 (a.m_bits & b.m_bits);
#endif
}

OIIO_FORCEINLINE vbool16 operator| (const vbool16 & a, const vbool16 & b) {
#if OIIO_SIMD_AVX >= 512
    return _mm512_kor (a.simd(), b.simd());
#else
    return vbool16 (a.m_bits | b.m_bits);
#endif
}

OIIO_FORCEINLINE vbool16 operator^ (const vbool16& a, const vbool16& b) {
#if OIIO_SIMD_AVX >= 512
    return _mm512_kxor (a.simd(), b.simd());
#else
    return vbool16 (a.m_bits ^ b.m_bits);
#endif
}


OIIO_FORCEINLINE const vbool16& operator&= (vbool16& a, const vbool16 &b) {
    return a = a & b;
}

OIIO_FORCEINLINE const vbool16& operator|= (vbool16& a, const vbool16& b) {
    return a = a | b;
}

OIIO_FORCEINLINE const vbool16& operator^= (vbool16& a, const vbool16& b) {
    return a = a ^ b;
}


OIIO_FORCEINLINE vbool16 operator~ (const vbool16& a) {
    return a ^ vbool16::True();
}


OIIO_FORCEINLINE vbool16 operator== (const vbool16 & a, const vbool16 & b) {
#if OIIO_SIMD_AVX >= 512
    return _mm512_kxnor (a.simd(), b.simd());
#else
    return vbool16 (!(a.m_bits ^ b.m_bits));
#endif
}

OIIO_FORCEINLINE vbool16 operator!= (const vbool16 & a, const vbool16 & b) {
#if OIIO_SIMD_AVX >= 512
    return _mm512_kxor (a.simd(), b.simd());
#else
    return vbool16 (a.m_bits ^ b.m_bits);
#endif
}


template<int i>
OIIO_FORCEINLINE bool extract (const vbool16& a) {
    return a[i];
}

template<int i>
OIIO_FORCEINLINE vbool16 insert (const vbool16& a, bool val) {