// ======================================================================== //
// Copyright 2009-2017 Intel Corporation                                    //
//                                                                          //
// Licensed under the Apache License, Version 2.0 (the "License");          //
// you may not use this file except in compliance with the License.         //
// You may obtain a copy of the License at                                  //
//                                                                          //
//     http://www.apache.org/licenses/LICENSE-2.0                           //
//                                                                          //
// Unless required by applicable law or agreed to in writing, software      //
// distributed under the License is distributed on an "AS IS" BASIS,        //
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
// See the License for the specific language governing permissions and      //
// limitations under the License.                                           //
// ======================================================================== //

#pragma once

namespace embree
{ 
  /* 16-wide AVX-512 unsigned integer type */
  template<>
    struct vuint<16>
  {
    typedef vboolf16 Bool;
    typedef vuint16  UInt;
    typedef vfloat16 Float;

    enum  { size = 16 }; // number of SIMD elements
    union {              // data
      __m512i v; 
      unsigned int i[16]; 
    };
    
    ////////////////////////////////////////////////////////////////////////////////
    /// Constructors, Assignment & Cast Operators
    ////////////////////////////////////////////////////////////////////////////////
       
    __forceinline vuint() {}
    __forceinline vuint(const vuint16& t) { v = t.v; }
    __forceinline vuint16& operator=(const vuint16& f) { v = f.v; return *this; }

    __forceinline vuint(const __m512i& t) { v = t; }
    __forceinline operator __m512i () const { return v; }
    __forceinline operator __m256i () const { return _mm512_castsi512_si256(v); }

    __forceinline vuint(const unsigned int i) {
      v = _mm512_set1_epi32(i);
    }
    
    __forceinline vuint(const unsigned int a, const unsigned int b, const unsigned int c, const unsigned int d) {
      v = _mm512_set4_epi32(d,c,b,a);      
    }

    __forceinline vuint(const unsigned int a0 , const unsigned int a1 , const unsigned int a2 , const unsigned int a3, 
                        const unsigned int a4 , const unsigned int a5 , const unsigned int a6 , const unsigned int a7, 
                        const unsigned int a8 , const unsigned int a9 , const unsigned int a10, const unsigned int a11, 
                        const unsigned int a12, const unsigned int a13, const unsigned int a14, const unsigned int a15)
    {
      v = _mm512_set_epi32(a15,a14,a13,a12,a11,a10,a9,a8,a7,a6,a5,a4,a3,a2,a1,a0);
    }
   
    __forceinline explicit vuint(const __m512& f) {
      v = _mm512_cvtps_epu32(f);
    }
    
    ////////////////////////////////////////////////////////////////////////////////
    /// Constants
    ////////////////////////////////////////////////////////////////////////////////
    
    __forceinline vuint( ZeroTy   ) : v(_mm512_setzero_epi32()) {}
    __forceinline vuint( OneTy    ) : v(_mm512_set1_epi32(1)) {}
    __forceinline vuint( StepTy   ) : v(_mm512_set_epi32(15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0)) {}
    __forceinline vuint( ReverseStepTy )   : v(_mm512_setr_epi32(15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0)) {}

    __forceinline static vuint16 zero() { return _mm512_setzero_epi32(); }
    __forceinline static vuint16 one () { return _mm512_set1_epi32(1); }
    __forceinline static vuint16 neg_one () { return _mm512_set1_epi32(-1); }

    ////////////////////////////////////////////////////////////////////////////////
    /// Loads and Stores
    ////////////////////////////////////////////////////////////////////////////////

    static __forceinline void store_nt(void *__restrict__ ptr, const vuint16& a) {
      _mm512_stream_si512((__m512i*)ptr,a);
    }

    static __forceinline vuint16 loadu(const void* addr)
    {
      return _mm512_loadu_si512(addr);
    }

    static __forceinline vuint16 load(const vuint16* addr) {
      return _mm512_load_si512(addr);
    }

    static __forceinline vuint16 load(const unsigned int* addr) {
      return _mm512_load_si512(addr);
    }

    static __forceinline void store(void* ptr, const vuint16& v) {
      _mm512_store_si512(ptr,v);
    }

    static __forceinline void storeu(void* ptr, const vuint16& v ) {
      _mm512_storeu_si512(ptr,v);
    }

    static __forceinline void storeu(const vboolf16& mask, int* ptr, const vuint16& f ) {
      _mm512_mask_storeu_epi32(ptr,mask,f);
    }

    static __forceinline void store(const vboolf16& mask, void* addr, const vuint16& v2) {
      _mm512_mask_store_epi32(addr,mask,v2);
    }

  /* pass by value to avoid compiler generating inefficient code */
    static __forceinline void storeu_compact(const vboolf16 mask,void * addr, const vuint16 reg) {
      _mm512_mask_compressstoreu_epi32(addr,mask,reg);
    }

    static __forceinline vuint16 compact64bit(const vboolf16& mask, vuint16 &v) {
      return _mm512_mask_compress_epi64(v,mask,v);
    }

    static __forceinline vuint16 compact(const vboolf16& mask, vuint16 &v) {
      return _mm512_mask_compress_epi32(v,mask,v);
    }

    static __forceinline vuint16 compact(const vboolf16& mask, const vuint16 &a, vuint16 &b) {
      return _mm512_mask_compress_epi32(a,mask,b);
    }

    static __forceinline vuint16 broadcast64bit(size_t v) {
      return _mm512_set1_epi64(v);
    }

    static __forceinline size_t extract64bit(const vuint16& v)
    {
      return _mm_cvtsi128_si64(_mm512_castsi512_si128(v));
    }


    ////////////////////////////////////////////////////////////////////////////////
    /// Array Access
    ////////////////////////////////////////////////////////////////////////////////
    
    __forceinline       unsigned int& operator[](const size_t index)       { assert(index < 16); return i[index]; }
    __forceinline const unsigned int& operator[](const size_t index) const { assert(index < 16); return i[index]; }

    __forceinline unsigned int uint    (const size_t index) const { assert(index < 16); return ((unsigned int*)i)[index]; }
    __forceinline size_t&      uint64_t(const size_t index) const { assert(index < 8);  return ((size_t*)i)[index]; }
  };
  
  ////////////////////////////////////////////////////////////////////////////////
  /// Unary Operators
  ////////////////////////////////////////////////////////////////////////////////

  __forceinline const vuint16 asUInt    ( const __m512& a ) { return _mm512_castps_si512(a); }
  __forceinline const vuint16 operator +( const vuint16& a ) { return a; }
  __forceinline const vuint16 operator -( const vuint16& a ) { return _mm512_sub_epi32(_mm512_setzero_epi32(), a); }

  ////////////////////////////////////////////////////////////////////////////////
  /// Binary Operators
  ////////////////////////////////////////////////////////////////////////////////

  __forceinline const vuint16 operator +( const vuint16& a, const vuint16& b )     { return _mm512_add_epi32(a, b); }
  __forceinline const vuint16 operator +( const vuint16& a, const unsigned int b ) { return a + vuint16(b); }
  __forceinline const vuint16 operator +( const unsigned int a, const vuint16& b ) { return vuint16(a) + b; }

  __forceinline const vuint16 operator -( const vuint16& a, const vuint16& b     ) { return _mm512_sub_epi32(a, b); }
  __forceinline const vuint16 operator -( const vuint16& a, const unsigned int b ) { return a - vuint16(b); }
  __forceinline const vuint16 operator -( const unsigned int a, const vuint16& b ) { return vuint16(a) - b; }

  __forceinline const vuint16 operator *( const vuint16& a, const vuint16& b )     { return _mm512_mul_epu32(a, b); }
  __forceinline const vuint16 operator *( const vuint16& a, const unsigned int b ) { return a * vuint16(b); }
  __forceinline const vuint16 operator *( const unsigned int a, const vuint16& b ) { return vuint16(a) * b; }

  __forceinline const vuint16 operator &( const vuint16& a, const vuint16& b )     { return _mm512_and_epi32(a, b); }
  __forceinline const vuint16 operator &( const vuint16& a, const unsigned int b ) { return a & vuint16(b); }
  __forceinline const vuint16 operator &( const unsigned int a, const vuint16& b ) { return vuint16(a) & b; }

  __forceinline const vuint16 operator |( const vuint16& a, const vuint16& b )     { return _mm512_or_epi32(a, b); }
  __forceinline const vuint16 operator |( const vuint16& a, const unsigned int b ) { return a | vuint16(b); }
  __forceinline const vuint16 operator |( const unsigned int a, const vuint16& b ) { return vuint16(a) | b; }

  __forceinline const vuint16 operator ^( const vuint16& a, const vuint16& b )     { return _mm512_xor_epi32(a, b); }
  __forceinline const vuint16 operator ^( const vuint16& a, const unsigned int b ) { return a ^ vuint16(b); }
  __forceinline const vuint16 operator ^( const unsigned int a, const vuint16& b ) { return vuint16(a) ^ b; }

  __forceinline const vuint16 operator <<( const vuint16& a, const unsigned int n ) { return _mm512_slli_epi32(a, n); }
  __forceinline const vuint16 operator >>( const vuint16& a, const unsigned int n ) { return _mm512_srli_epi32(a, n); }

  __forceinline const vuint16 operator <<( const vuint16& a, const vuint16& n ) { return _mm512_sllv_epi32(a, n); }
  __forceinline const vuint16 operator >>( const vuint16& a, const vuint16& n ) { return _mm512_srlv_epi32(a, n); }

  __forceinline const vuint16 sll ( const vuint16& a, const unsigned int b ) { return _mm512_slli_epi32(a, b); }
  __forceinline const vuint16 sra ( const vuint16& a, const unsigned int b ) { return _mm512_srai_epi32(a, b); }
  __forceinline const vuint16 srl ( const vuint16& a, const unsigned int b ) { return _mm512_srli_epi32(a, b); }
  
  __forceinline const vuint16 min( const vuint16& a, const vuint16& b )     { return _mm512_min_epu32(a, b); }
  __forceinline const vuint16 min( const vuint16& a, const unsigned int b ) { return min(a,vuint16(b)); }
  __forceinline const vuint16 min( const unsigned int a, const vuint16& b ) { return min(vuint16(a),b); }

  __forceinline const vuint16 max( const vuint16& a, const vuint16& b )     { return _mm512_max_epu32(a, b); }
  __forceinline const vuint16 max( const vuint16& a, const unsigned int b ) { return max(a,vuint16(b)); }
  __forceinline const vuint16 max( const unsigned int a, const vuint16& b ) { return max(vuint16(a),b); }
  
  __forceinline const vuint16 mask_add(const vboolf16& mask, vuint16& c, const vuint16& a, const vuint16& b) { return _mm512_mask_add_epi32(c,mask,a,b); }
  __forceinline const vuint16 mask_sub(const vboolf16& mask, vuint16& c, const vuint16& a, const vuint16& b) { return _mm512_mask_sub_epi32(c,mask,a,b); }

  __forceinline const vuint16 mask_and(const vboolf16& m,vuint16& c, const vuint16& a, const vuint16& b) { return _mm512_mask_and_epi32(c,m,a,b); }
  __forceinline const vuint16 mask_or (const vboolf16& m,vuint16& c, const vuint16& a, const vuint16& b) { return _mm512_mask_or_epi32(c,m,a,b); }

  ////////////////////////////////////////////////////////////////////////////////
  /// Assignment Operators
  ////////////////////////////////////////////////////////////////////////////////

  __forceinline vuint16& operator +=( vuint16& a, const vuint16& b ) { return a = a + b; }
  __forceinline vuint16& operator +=( vuint16& a, const unsigned int     b ) { return a = a + b; }
  
  __forceinline vuint16& operator -=( vuint16& a, const vuint16& b ) { return a = a - b; }
  __forceinline vuint16& operator -=( vuint16& a, const unsigned int     b ) { return a = a - b; }

  __forceinline vuint16& operator *=( vuint16& a, const vuint16& b ) { return a = a * b; }
  __forceinline vuint16& operator *=( vuint16& a, const unsigned int     b ) { return a = a * b; }
  
  __forceinline vuint16& operator &=( vuint16& a, const vuint16& b ) { return a = a & b; }
  __forceinline vuint16& operator &=( vuint16& a, const unsigned int     b ) { return a = a & b; }
  
  __forceinline vuint16& operator |=( vuint16& a, const vuint16& b ) { return a = a | b; }
  __forceinline vuint16& operator |=( vuint16& a, const unsigned int     b ) { return a = a | b; }
  
  __forceinline vuint16& operator <<=( vuint16& a, const unsigned int b ) { return a = a << b; }
  __forceinline vuint16& operator >>=( vuint16& a, const unsigned int b ) { return a = a >> b; }


  ////////////////////////////////////////////////////////////////////////////////
  /// Comparison Operators + Select
  ////////////////////////////////////////////////////////////////////////////////

  __forceinline const vboolf16 operator ==( const vuint16& a, const vuint16& b ) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_EQ); }
  __forceinline const vboolf16 operator ==( const vuint16& a, const unsigned int b ) { return a == vuint16(b); }
  __forceinline const vboolf16 operator ==( const unsigned int a, const vuint16& b ) { return vuint16(a) == b; }
  
  __forceinline const vboolf16 operator !=( const vuint16& a, const vuint16& b ) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_NE); }
  __forceinline const vboolf16 operator !=( const vuint16& a, const unsigned int b ) { return a != vuint16(b); }
  __forceinline const vboolf16 operator !=( const unsigned int a, const vuint16& b ) { return vuint16(a) != b; }
  
  __forceinline const vboolf16 operator < ( const vuint16& a, const vuint16& b ) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_LT); }
  __forceinline const vboolf16 operator < ( const vuint16& a, const unsigned int b ) { return a <  vuint16(b); }
  __forceinline const vboolf16 operator < ( const unsigned int a, const vuint16& b ) { return vuint16(a) <  b; }
  
  __forceinline const vboolf16 operator >=( const vuint16& a, const vuint16& b ) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_GE); }
  __forceinline const vboolf16 operator >=( const vuint16& a, const unsigned int b ) { return a >= vuint16(b); }
  __forceinline const vboolf16 operator >=( const unsigned int a, const vuint16& b ) { return vuint16(a) >= b; }

  __forceinline const vboolf16 operator > ( const vuint16& a, const vuint16& b ) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_GT); }
  __forceinline const vboolf16 operator > ( const vuint16& a, const unsigned int b ) { return a >  vuint16(b); }
  __forceinline const vboolf16 operator > ( const unsigned int a, const vuint16& b ) { return vuint16(a) >  b; }

  __forceinline const vboolf16 operator <=( const vuint16& a, const vuint16& b ) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_LE); }
  __forceinline const vboolf16 operator <=( const vuint16& a, const unsigned int b ) { return a <= vuint16(b); }
  __forceinline const vboolf16 operator <=( const unsigned int a, const vuint16& b ) { return vuint16(a) <= b; }

  __forceinline vboolf16 eq(const vuint16& a, const vuint16& b) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_EQ); }
  __forceinline vboolf16 ne(const vuint16& a, const vuint16& b) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_NE); }
  __forceinline vboolf16 lt(const vuint16& a, const vuint16& b) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_LT); }
  __forceinline vboolf16 ge(const vuint16& a, const vuint16& b) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_GE); }
  __forceinline vboolf16 gt(const vuint16& a, const vuint16& b) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_GT); }
  __forceinline vboolf16 le(const vuint16& a, const vuint16& b) { return _mm512_cmp_epu32_mask(a,b,_MM_CMPINT_LE); }

  __forceinline vboolf16 eq(const vboolf16 mask, const vuint16& a, const vuint16& b) { return _mm512_mask_cmp_epu32_mask(mask,a,b,_MM_CMPINT_EQ); }
  __forceinline vboolf16 ne(const vboolf16 mask, const vuint16& a, const vuint16& b) { return _mm512_mask_cmp_epu32_mask(mask,a,b,_MM_CMPINT_NE); }
  __forceinline vboolf16 lt(const vboolf16 mask, const vuint16& a, const vuint16& b) { return _mm512_mask_cmp_epu32_mask(mask,a,b,_MM_CMPINT_LT); }
  __forceinline vboolf16 ge(const vboolf16 mask, const vuint16& a, const vuint16& b) { return _mm512_mask_cmp_epu32_mask(mask,a,b,_MM_CMPINT_GE); }
  __forceinline vboolf16 gt(const vboolf16 mask, const vuint16& a, const vuint16& b) { return _mm512_mask_cmp_epu32_mask(mask,a,b,_MM_CMPINT_GT); }
  __forceinline vboolf16 le(const vboolf16 mask, const vuint16& a, const vuint16& b) { return _mm512_mask_cmp_epu32_mask(mask,a,b,_MM_CMPINT_LE); }
    
 
  __forceinline const vuint16 select( const vboolf16& m, const vuint16& t, const vuint16& f ) {
    return _mm512_mask_or_epi32(f,m,t,t); 
  }

  __forceinline void xchg(const vboolf16& m, vuint16& a, vuint16& b) {
    const vuint16 c = a; a = select(m,b,a); b = select(m,c,b);
  }

  __forceinline vboolf16 test(const vboolf16& m, const vuint16& a, const vuint16& b) {
    return _mm512_mask_test_epi32_mask(m,a,b);
  }

  __forceinline vboolf16 test(const vuint16& a, const vuint16& b) {
    return _mm512_test_epi32_mask(a,b);
  }

  ////////////////////////////////////////////////////////////////////////////////
  // Movement/Shifting/Shuffling Functions
  ////////////////////////////////////////////////////////////////////////////////

  template<size_t i> 
    __forceinline const vuint16 shuffle( const vuint16& a ) {
    return _mm512_castps_si512(_mm512_permute_ps(_mm512_castsi512_ps(a), _MM_SHUFFLE(i, i, i, i)));
  }

  template<int A, int B, int C, int D> 
    __forceinline vuint16 shuffle (const vuint16& v) {
    return _mm512_castps_si512(_mm512_permute_ps(_mm512_castsi512_ps(v),_MM_SHUFFLE(D,C,B,A))); 
  }

  template<int i>
    __forceinline vuint16 shuffle4(const vuint16& x) { 
    return _mm512_castps_si512(_mm512_shuffle_f32x4(_mm512_castsi512_ps(x),_mm512_castsi512_ps(x),_MM_SHUFFLE(i,i,i,i)));
  }

  template<int A, int B, int C, int D>
    __forceinline vuint16 shuffle4(const vuint16& x) { 
    return _mm512_castps_si512(_mm512_shuffle_f32x4(_mm512_castsi512_ps(x),_mm512_castsi512_ps(x),_MM_SHUFFLE(D,C,B,A)));
  }

  template<int i>
    __forceinline vuint16 align_shift_right(const vuint16& a, const vuint16& b)
  {
    return _mm512_alignr_epi32(a,b,i); 
  };

  __forceinline unsigned int toScalar(const vuint16& a)
  {
    return _mm_cvtsi128_si32(_mm512_castsi512_si128(a));
  }

  ////////////////////////////////////////////////////////////////////////////////
  /// Reductions
  ////////////////////////////////////////////////////////////////////////////////

  __forceinline unsigned int reduce_add(vuint16 a) { return _mm512_reduce_add_epi32(a); }
  __forceinline unsigned int reduce_min(vuint16 a) { return _mm512_reduce_min_epu32(a); }
  __forceinline unsigned int reduce_max(vuint16 a) { return _mm512_reduce_max_epu32(a); }
  __forceinline unsigned int reduce_and(vuint16 a) { return _mm512_reduce_and_epi32(a); }

  __forceinline vuint16 vreduce_min2(vuint16 x) {                      return min(x,shuffle<1,0,3,2>(x)); }
  __forceinline vuint16 vreduce_min4(vuint16 x) { x = vreduce_min2(x); return min(x,shuffle<2,3,0,1>(x)); }
  __forceinline vuint16 vreduce_min8(vuint16 x) { x = vreduce_min4(x); return min(x,shuffle4<1,0,3,2>(x)); }
  __forceinline vuint16 vreduce_min (vuint16 x) { x = vreduce_min8(x); return min(x,shuffle4<2,3,0,1>(x)); }

  __forceinline vuint16 vreduce_max2(vuint16 x) {                      return max(x,shuffle<1,0,3,2>(x)); }
  __forceinline vuint16 vreduce_max4(vuint16 x) { x = vreduce_max2(x); return max(x,shuffle<2,3,0,1>(x)); }
  __forceinline vuint16 vreduce_max8(vuint16 x) { x = vreduce_max4(x); return max(x,shuffle4<1,0,3,2>(x)); }
  __forceinline vuint16 vreduce_max (vuint16 x) { x = vreduce_max8(x); return max(x,shuffle4<2,3,0,1>(x)); }

  __forceinline vuint16 vreduce_and2(vuint16 x) {                      return x & shuffle<1,0,3,2>(x); }
  __forceinline vuint16 vreduce_and4(vuint16 x) { x = vreduce_and2(x); return x & shuffle<2,3,0,1>(x); }
  __forceinline vuint16 vreduce_and8(vuint16 x) { x = vreduce_and4(x); return x & shuffle4<1,0,3,2>(x); }
  __forceinline vuint16 vreduce_and (vuint16 x) { x = vreduce_and8(x); return x & shuffle4<2,3,0,1>(x); }

  __forceinline vuint16 vreduce_or2(vuint16 x) {                     return x | shuffle<1,0,3,2>(x); }
  __forceinline vuint16 vreduce_or4(vuint16 x) { x = vreduce_or2(x); return x | shuffle<2,3,0,1>(x); }
  __forceinline vuint16 vreduce_or8(vuint16 x) { x = vreduce_or4(x); return x | shuffle4<1,0,3,2>(x); }
  __forceinline vuint16 vreduce_or (vuint16 x) { x = vreduce_or8(x); return x | shuffle4<2,3,0,1>(x); }

  __forceinline vuint16 vreduce_add2(vuint16 x) {                      return x + shuffle<1,0,3,2>(x); }
  __forceinline vuint16 vreduce_add4(vuint16 x) { x = vreduce_add2(x); return x + shuffle<2,3,0,1>(x); }
  __forceinline vuint16 vreduce_add8(vuint16 x) { x = vreduce_add4(x); return x + shuffle4<1,0,3,2>(x); }
  __forceinline vuint16 vreduce_add (vuint16 x) { x = vreduce_add8(x); return x + shuffle4<2,3,0,1>(x); }

  ////////////////////////////////////////////////////////////////////////////////
  /// Memory load and store operations
  ////////////////////////////////////////////////////////////////////////////////

  template<int scale = 4>
  __forceinline vuint16 gather16i(const vboolf16& mask, const unsigned int *const ptr, const vuint16& index) {
    return _mm512_mask_i32gather_epi32(_mm512_undefined_epi32(),mask,index,ptr,scale);
  }
  
  template<int scale = 4>
    __forceinline vuint16 gather16i(const vboolf16& mask, vuint16& dest, const unsigned int *const ptr, const vuint16& index) {
    return _mm512_mask_i32gather_epi32(dest,mask,index,ptr,scale);
  }
  
  template<int scale = 4>
    __forceinline void scatter16i(const vboolf16& mask,unsigned int *const ptr, const vuint16& index,const vuint16& v) {
    _mm512_mask_i32scatter_epi32((int*)ptr,mask,index,v,scale);
  }

  __forceinline void compactustore16i_low(const vboolf16 mask, void *addr, const vuint16& reg) {
    _mm512_mask_compressstoreu_epi32(addr,mask,reg);
  }
  
  __forceinline vuint16 permute(vuint16 v,vuint16 index) {
    return _mm512_permutexvar_epi32(index,v);  
  }

  __forceinline vuint16 reverse(const vuint16 &a) {
    return permute(a,vuint16(reverse_step));
  }

  __forceinline vuint16 prefix_sum(const vuint16& a) 
  {
    const vuint16 z(zero);
    vuint16 v = a;
    v = v + align_shift_right<16-1>(v,z);
    v = v + align_shift_right<16-2>(v,z);
    v = v + align_shift_right<16-4>(v,z);
    v = v + align_shift_right<16-8>(v,z);
    return v;  
  }

  __forceinline vuint16 reverse_prefix_sum(const vuint16& a) 
  {
    const vuint16 z(zero);
    vuint16 v = a;
    v = v + align_shift_right<1>(z,v);
    v = v + align_shift_right<2>(z,v);
    v = v + align_shift_right<4>(z,v);
    v = v + align_shift_right<8>(z,v);
    return v;  
  }

  ////////////////////////////////////////////////////////////////////////////////
  /// Output Operators
  ////////////////////////////////////////////////////////////////////////////////
  
  __forceinline std::ostream& operator<<(std::ostream& cout, const vuint16& v)
  {
    cout << "<" << v[0];
    for (int i=1; i<16; i++) cout << ", " << v[i];
    cout << ">";
    return cout;
  }
}