vectori512.h

/****************************  vectori512.h   *******************************
* Author:        Agner Fog
* Date created:  2014-07-23
* Last modified: 2019-08-01
* Version:       1.40.00
* Project:       vector class library
* Description:
* Header file defining 512-bit integer vector classes for 32 and 64 bit integers.
* For x86 microprocessors with AVX512F and later instruction sets.
*
* Instructions: see vcl_manual.pdf
*
* The following vector classes are defined here:
* Vec16i    Vector of  16  32-bit signed   integers
* Vec16ui   Vector of  16  32-bit unsigned integers
* Vec16ib   Vector of  16  Booleans for use with Vec16i and Vec16ui
* Vec8q     Vector of   8  64-bit signed   integers
* Vec8uq    Vector of   8  64-bit unsigned integers
* Vec8qb    Vector of   8  Booleans for use with Vec8q and Vec8uq
* Other 512-bit integer vectors are defined in Vectori512s.h
*
* Each vector object is represented internally in the CPU as a 512-bit register.
* This header file defines operators and functions for these vectors.
*
* (c) Copyright 2012-2019 Agner Fog.
* Apache License version 2.0 or later.
*****************************************************************************/

#ifndef VECTORCLASS_H
#include "vectorclass.h"
#endif

#if VECTORCLASS_H > 19999
#error Incompatible versions of vector class library mixed
#endif

#ifndef VECTORI512_H
#define VECTORI512_H  1

// check combination of header files
#ifdef VECTORI512E_H
#error Two different versions of vectori512.h included
#endif


#ifdef VCL_NAMESPACE
namespace VCL_NAMESPACE {
#endif

// Bug fix for missing intrinsics: !
// _mm512_cmpgt_epu32_mask, _mm512_cmpgt_epu64_mask
// all typecast intrinsics
// Fix expected in GCC version 4.9.3 or 5.0. https://gcc.gnu.org/bugzilla/show_bug.cgi?id=61878

// questionable
// _mm512_mask_mov_epi32 check select(). 
// Doc at https://software.intel.com/en-us/node/513888 is wrong. Bug report filed
// Doc at https://software.intel.com/en-us/cpp-compiler-developer-guide-and-reference-intrinsics-for-integer-move-operations is still wrong. Filed another bug report 2019


#if defined (GCC_VERSION) && GCC_VERSION < 50000 && !defined(__INTEL_COMPILER) && !defined(__clang__)

static inline  __m512i _mm512_castsi256_si512(__m256i x) {
    union {
        __m512i a;
        __m256i b;
    } u;
    u.b = x;
    return u.a;
}

static inline  __m256i _mm512_castsi512_si256(__m512i x) {
    union {
        __m512i a;
        __m256i b;
    } u;
    u.a = x;
    return u.b;
}

static inline  __m512i _mm512_castsi128_si512(__m128i x) {
    union {
        __m128i a;
        __m512i b;
    } u;
    u.a = x;
    return u.b;
}

static inline  __m128i _mm512_castsi512_si128(__m512i x) {
    union {
        __m512i a;
        __m128i b;
    } u;
    u.a = x;
    return u.b;
}

#endif


/*****************************************************************************
*
*          Generate compile-time constant vector
*
*****************************************************************************/
// Generate a constant vector of 8 integers stored in memory.
// Can be converted to any integer vector type
template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7, 
int i8, int i9, int i10, int i11, int i12, int i13, int i14, int i15>
static inline __m512i constant16i() {
    static const union {
        int32_t i[16];
        __m512i zmm;
    } u = {{i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15}};
    return u.zmm;
}
template <uint32_t i0, uint32_t i1, uint32_t i2, uint32_t i3, uint32_t i4, uint32_t i5, uint32_t i6, uint32_t i7, 
    uint32_t i8, uint32_t i9, uint32_t i10, uint32_t i11, uint32_t i12, uint32_t i13, uint32_t i14, uint32_t i15>
    static inline __m512i constant16ui() {
    static const union {
        uint32_t i[16];
        __m512i zmm;
    } u = {{i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15}};
    return u.zmm;
}


/*****************************************************************************
*
*          Boolean vector base classes for AVX512
*
*****************************************************************************/

class Vec16b {
protected:
    __mmask16  m16; // Boolean vector
public:
    // Default constructor:
    Vec16b () {
    }
    // Constructor to convert from type __mmask16 used in intrinsics:
    Vec16b (__mmask16 x) {
        m16 = x;
    }
    // Constructor to build from all elements:
    Vec16b(bool b0, bool b1, bool b2, bool b3, bool b4, bool b5, bool b6, bool b7, 
    bool b8, bool b9, bool b10, bool b11, bool b12, bool b13, bool b14, bool b15) {
        m16 = uint16_t(
			(uint16_t)b0 | (uint16_t)b1<<1 | (uint16_t)b2<<2 | (uint16_t)b3<<3 | 
			(uint16_t)b4<<4 | (uint16_t)b5<<5 | (uint16_t)b6<<6 | (uint16_t)b7<<7 |
			(uint16_t)b8<<8 | (uint16_t)b9<<9 | (uint16_t)b10<<10 | (uint16_t)b11<<11 | 
			(uint16_t)b12<<12 | (uint16_t)b13<<13 | (uint16_t)b14<<14 | (uint16_t)b15<<15);
    }
    // Constructor to broadcast single value:
    Vec16b(bool b) {
        m16 = __mmask16(-int16_t(b));
    }
    // Constructor to make from two halves
    Vec16b (Vec8ib const & x0, Vec8ib const & x1) {
        // = Vec16i(x0,x1) != 0;  (not defined yet)
        __m512i z = _mm512_inserti64x4(_mm512_castsi256_si512(x0), x1, 1);
        m16 = _mm512_cmpneq_epi32_mask(z, _mm512_setzero_epi32());
    }        
    // Assignment operator to convert from type __mmask16 used in intrinsics:
    Vec16b & operator = (__mmask16 x) {
        m16 = x;
        return *this;
    }
    // Assignment operator to broadcast scalar value:
    Vec16b & operator = (bool b) {
        m16 = Vec16b(b);
        return *this;
    }
    // Type cast operator to convert to __mmask16 used in intrinsics
    operator __mmask16() const {
        return m16;
    }
    // split into two halves
    Vec8ib get_low() const {
        return Vec8ib().load_bits(uint8_t(m16));
    }
    Vec8ib get_high() const {
        return Vec8ib().load_bits(uint8_t((uint16_t)m16 >> 8u));
    }
    // Member function to change a single element in vector
    Vec16b const & insert(int index, bool value) {
        m16 = __mmask16(((uint16_t)m16 & ~(1 << index)) | (int)value << index);
        return *this;
    }
    // Member function extract a single element from vector
    bool extract(int index) const {
        return ((uint32_t)m16 >> index) & 1;
    }
    // Extract a single element. Operator [] can only read an element, not write.
    bool operator [] (int index) const {
        return extract(index);
    }
    // Member function to change a bitfield to a boolean vector
    Vec16b & load_bits(uint16_t a) {
        m16 = __mmask16(a);
        return *this;
    }
    // Number of elements
    static int size() {
        return 16;
    }
    // Type of elements
    static int elementtype() {
        return 2;
    }
private: // Prevent constructing from int, etc.
    Vec16b(int b);
    Vec16b & operator = (int x);
};

// Define operators for this class

// vector operator & : and
static inline Vec16b operator & (Vec16b a, Vec16b b) {
    return _mm512_kand(__mmask16(a), __mmask16(b));
}
static inline Vec16b operator && (Vec16b a, Vec16b b) {
    return a & b;
}

// vector operator | : or
static inline Vec16b operator | (Vec16b a, Vec16b b) {
    return _mm512_kor(__mmask16(a), __mmask16(b));
}
static inline Vec16b operator || (Vec16b a, Vec16b b) {
    return a | b;
}

// vector operator ^ : xor
static inline Vec16b operator ^ (Vec16b a, Vec16b b) {
    return _mm512_kxor(__mmask16(a), __mmask16(b));
}

// vector operator ~ : not
static inline Vec16b operator ~ (Vec16b a) {
    return _mm512_knot(__mmask16(a));
}

// vector operator ! : element not
static inline Vec16b operator ! (Vec16b a) {
    return ~a;
}

// vector operator &= : and
static inline Vec16b & operator &= (Vec16b & a, Vec16b b) {
    a = a & b;
    return a;
}

// vector operator |= : or
static inline Vec16b & operator |= (Vec16b & a, Vec16b b) {
    a = a | b;
    return a;
}

// vector operator ^= : xor
static inline Vec16b & operator ^= (Vec16b & a, Vec16b b) {
    a = a ^ b;
    return a;
}


/*****************************************************************************
*
*          Functions for boolean vectors
*
*****************************************************************************/

// function andnot: a & ~ b
static inline Vec16b andnot (Vec16b a, Vec16b b) {
    return _mm512_kandn(b, a);
}

// horizontal_and. Returns true if all bits are 1
static inline bool horizontal_and (Vec16b const & a) {
    return (uint16_t)(__mmask16)a == 0xFFFF;
}

// horizontal_or. Returns true if at least one bit is 1
static inline bool horizontal_or (Vec16b const & a) {
    return (uint16_t)(__mmask16)a != 0;
}


/*****************************************************************************
*
*          Vec16ib: Vector of 16 Booleans for use with Vec16i and Vec16ui
*
*****************************************************************************/

class Vec16ib : public Vec16b {
public:
    // Default constructor:
    Vec16ib () {
    }
    Vec16ib (Vec16b x) {
        m16 = x;
    }
    // Constructor to build from all elements:
    Vec16ib(bool x0, bool x1, bool x2, bool x3, bool x4, bool x5, bool x6, bool x7,
        bool x8, bool x9, bool x10, bool x11, bool x12, bool x13, bool x14, bool x15) :
        Vec16b(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15) {
    }
    // Constructor to convert from type __mmask16 used in intrinsics:
    Vec16ib (__mmask16 x) {
        m16 = x;
    }
    // Constructor to broadcast single value:
    Vec16ib(bool b) : Vec16b(b) {}
    // Constructor to make from two halves
    Vec16ib (Vec8ib const & x0, Vec8ib const & x1) {
        m16 = Vec16b(x0, x1);
    }
    // Assignment operator to convert from type __mmask16 used in intrinsics:
    Vec16ib & operator = (__mmask16 x) {
        m16 = x;
        return *this;
    }
    // Assignment operator to broadcast scalar value:
    Vec16ib & operator = (bool b) {
        m16 = Vec16b(b);
        return *this;
    }
private: // Prevent constructing from int, etc.
    Vec16ib(int b);
    Vec16ib & operator = (int x);
};

// Define operators for Vec16ib

// vector operator & : and
static inline Vec16ib operator & (Vec16ib a, Vec16ib b) {
    return Vec16b(a) & Vec16b(b);
}
static inline Vec16ib operator && (Vec16ib a, Vec16ib b) {
    return a & b;
}

// vector operator | : or
static inline Vec16ib operator | (Vec16ib a, Vec16ib b) {
    return Vec16b(a) | Vec16b(b);
}
static inline Vec16ib operator || (Vec16ib a, Vec16ib b) {
    return a | b;
}

// vector operator ^ : xor
static inline Vec16ib operator ^ (Vec16ib a, Vec16ib b) {
    return Vec16b(a) ^ Vec16b(b);
}

// vector operator == : xnor
static inline Vec16ib operator == (Vec16ib const & a, Vec16ib const & b) {
    return _mm512_kxnor(__mmask16(a), __mmask16(b));
}

// vector operator != : xor
static inline Vec16ib operator != (Vec16ib const & a, Vec16ib const & b) {
    return Vec16ib(a ^ b);
}

// vector operator ~ : not
static inline Vec16ib operator ~ (Vec16ib a) {
    return ~Vec16b(a);
}

// vector operator ! : not
static inline Vec16ib operator ! (Vec16ib a) {
    return ~a;
}

// vector operator &= : and
static inline Vec16ib & operator &= (Vec16ib & a, Vec16ib b) {
    a = a & b;
    return a;
}

// vector operator |= : or
static inline Vec16ib & operator |= (Vec16ib & a, Vec16ib b) {
    a = a | b;
    return a;
}

// vector operator ^= : xor
static inline Vec16ib & operator ^= (Vec16ib & a, Vec16ib b) {
    a = a ^ b;
    return a;
}

// vector function andnot
static inline Vec16ib andnot (Vec16ib a, Vec16ib b) {
    return Vec16ib(andnot(Vec16b(a), Vec16b(b)));
}


/*****************************************************************************
*
*          Vec8b: Base class vector of 8 Booleans
*
*****************************************************************************/

class Vec8b : public Vec16b {
public:
    // Default constructor:
    Vec8b () {
    }
    // Constructor to convert from type __mmask16 used in intrinsics:
    Vec8b (__mmask16 x) {
        m16 = x;
    }
    // Constructor to build from all elements:
    Vec8b(bool b0, bool b1, bool b2, bool b3, bool b4, bool b5, bool b6, bool b7) {
        m16 = uint16_t((uint16_t)b0 | (uint16_t)b1<<1 | (uint16_t)b2<<2 | (uint16_t)b3<<3 | 
			(uint16_t)b4<<4 | (uint16_t)b5<<5 | (uint16_t)b6<<6 | (uint16_t)b7<<7);
    }
    Vec8b (Vec16b const & x) {
        m16 = __mmask16(x);
    }
    // Constructor to broadcast single value:
    Vec8b(bool b) {
        m16 = __mmask16(-int8_t(b));
    }
    // Assignment operator to convert from type __mmask16 used in intrinsics:
    Vec8b & operator = (__mmask16 x) {
        m16 = x;
        return *this;
    }
	// Type cast operator to convert to __mmask8 used in intrinsics
	operator __mmask8() const {
		return (__mmask8)m16;
	}
    // split into two halves
    Vec4qb get_low() const {
        return Vec4qb(Vec4q(_mm512_castsi512_si256(_mm512_maskz_set1_epi64(__mmask8(m16), -1LL))));
    }
    Vec4qb get_high() const {
        return Vec8b(__mmask16(m16 >> 4)).get_low();
    }
    // Member function to change a bitfield to a boolean vector
    Vec8b & load_bits(uint8_t a) {
        m16 = __mmask16(uint16_t(a));
        return *this;
    }
    static int size() {
        return 8;
    }
    static int elementtype() {
        return 2;
    } 
private: // Prevent constructing from int etc. because of ambiguity
    Vec8b(int b);
    Vec8b & operator = (int x);
};

/*****************************************************************************
*
*          Functions for boolean vectors
*
*****************************************************************************/

// function andnot: a & ~ b
static inline Vec8b andnot (Vec8b a, Vec8b b) {
    return _mm512_kandn(b, a);
}

// horizontal_and. Returns true if all bits are 1
static inline bool horizontal_and (Vec8b const & a) {
    return (uint8_t)(__mmask16)a == 0xFF;
}

// horizontal_or. Returns true if at least one bit is 1
static inline bool horizontal_or (Vec8b const & a) {
    return (uint8_t)(__mmask16)a != 0;
}


/*****************************************************************************
*
*          Vec8qb: Vector of 8 Booleans for use with Vec8q and Vec8qu
*
*****************************************************************************/

class Vec8qb : public Vec8b {
public:
    // Default constructor:
    Vec8qb () {
    }
    Vec8qb (Vec16b x) {
        m16 = x;
    }
    // Constructor to build from all elements:
    Vec8qb(bool x0, bool x1, bool x2, bool x3, bool x4, bool x5, bool x6, bool x7) :
        Vec8b(x0, x1, x2, x3, x4, x5, x6, x7) {
    }
    // Constructor to convert from type __mmask8 used in intrinsics:
    Vec8qb (__mmask8 x) {
        m16 = (__mmask16)x;
    }
    // Constructor to convert from type __mmask16 used in intrinsics:
    Vec8qb (__mmask16 x) {
        m16 = x;
    }
    // Assignment operator to convert from type __mmask16 used in intrinsics:
    Vec8qb & operator = (__mmask16 x) {
        m16 = x;
        return *this;
    }
    // Constructor to broadcast single value:
    Vec8qb(bool b) : Vec8b(b) {}
    // Assignment operator to broadcast scalar:
    Vec8qb & operator = (bool b) {
        m16 = Vec8b(b);
        return *this;
    }
private: // Prevent constructing from int, etc.
    Vec8qb(int b);
    Vec8qb & operator = (int x);
public:
    // Constructor to make from two halves
    Vec8qb (Vec4qb const & x0, Vec4qb const & x1) {
        // = Vec8q(x0,x1) != 0;  (not defined yet)
        __m512i z = _mm512_inserti64x4(_mm512_castsi256_si512(x0), x1, 1);
        m16 = _mm512_cmpneq_epi64_mask(z, _mm512_setzero_si512());
    }        
};

// Define operators for Vec8qb

// vector operator & : and
static inline Vec8qb operator & (Vec8qb a, Vec8qb b) {
    return Vec16b(a) & Vec16b(b);
}
static inline Vec8qb operator && (Vec8qb a, Vec8qb b) {
    return a & b;
}

// vector operator | : or
static inline Vec8qb operator | (Vec8qb a, Vec8qb b) {
    return Vec16b(a) | Vec16b(b);
}
static inline Vec8qb operator || (Vec8qb a, Vec8qb b) {
    return a | b;
}

// vector operator ^ : xor
static inline Vec8qb operator ^ (Vec8qb a, Vec8qb b) {
    return Vec16b(a) ^ Vec16b(b);
}

// vector operator == : xnor
static inline Vec8qb operator == (Vec8qb const & a, Vec8qb const & b) {
    return _mm512_kxnor(a, b);
}

// vector operator != : xor
static inline Vec8qb operator != (Vec8qb const & a, Vec8qb const & b) {
    return Vec8qb(a ^ b);
}

// vector operator ~ : not
static inline Vec8qb operator ~ (Vec8qb a) {
    return ~Vec16b(a);
}

// vector operator ! : not
static inline Vec8qb operator ! (Vec8qb a) {
    return ~a;
}

// vector operator &= : and
static inline Vec8qb & operator &= (Vec8qb & a, Vec8qb b) {
    a = a & b;
    return a;
}

// vector operator |= : or
static inline Vec8qb & operator |= (Vec8qb & a, Vec8qb b) {
    a = a | b;
    return a;
}

// vector operator ^= : xor
static inline Vec8qb & operator ^= (Vec8qb & a, Vec8qb b) {
    a = a ^ b;
    return a;
}

// to_bits: convert to integer bitfield
static inline uint32_t to_bits(Vec8qb a) {
    return (uint8_t)(__mmask16)a;
}

// vector function andnot
static inline Vec8qb andnot (Vec8qb a, Vec8qb b) {
    return Vec8qb(andnot(Vec16b(a), Vec16b(b)));
}


/*****************************************************************************
*
*          Vector of 512 1-bit unsigned integers (base class for Vec16i)
*
*****************************************************************************/
class Vec512b {
protected:
    __m512i zmm; // Integer vector
public:
    // Default constructor:
    Vec512b() {
    }
    // Constructor to build from two Vec256b:
    Vec512b(Vec256b const & a0, Vec256b const & a1) {
        zmm = _mm512_inserti64x4(_mm512_castsi256_si512(a0), a1, 1);
    }
    // Constructor to convert from type __m512i used in intrinsics:
    Vec512b(__m512i const & x) {
        zmm = x;
    }
    // Assignment operator to convert from type __m512i used in intrinsics:
    Vec512b & operator = (__m512i const & x) {
        zmm = x;
        return *this;
    }
    // Type cast operator to convert to __m512i used in intrinsics
    operator __m512i() const {
        return zmm;
    }
    // Member function to load from array (unaligned)
    Vec512b & load(void const * p) {
        zmm = _mm512_loadu_si512(p);
        return *this;
    }
    // Member function to load from array, aligned by 64
    // You may use load_a instead of load if you are certain that p points to an address
    // divisible by 64, but there is hardly any speed advantage of load_a on modern processors
    Vec512b & load_a(void const * p) {
        zmm = _mm512_load_si512(p);
        return *this;
    }
    // Member function to store into array (unaligned)
    void store(void * p) const {
        _mm512_storeu_si512(p, zmm);
    }
    // Member function to store into array, aligned by 64
    // You may use store_a instead of store if you are certain that p points to an address
    // divisible by 64, but there is hardly any speed advantage of store_a on modern processors
    void store_a(void * p) const {
        _mm512_store_si512(p, zmm);
    }
    // Member function to change a single bit, mainly for test purposes
    // DEPRECATED
    Vec512b const & set_bit(int index, int value) {
        static uint64_t m[16] = {0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0};
        int wi = (index >> 6) & 7;               // qword index
        int bi = index & 0x3F;                   // bit index within qword w

        __m512i mask = Vec512b().load(m+8-wi);   // 1 in qword number wi
        mask = _mm512_sll_epi64(mask,_mm_cvtsi32_si128(bi)); // mask with bit number b set
        if (value & 1) {
            zmm = _mm512_or_si512(mask,zmm);
        }
        else {
            zmm = _mm512_andnot_si512(mask,zmm);
        }
        return *this;
    }
    // Member function to get a single bit, mainly for test purposes
    // DEPRECATED
    int get_bit(int index) const {
        union {
            __m512i z;
            uint8_t i[64];
        } u;
        u.z = zmm; 
        int wi = (index >> 3) & 0x3F;            // byte index
        int bi = index & 7;                      // bit index within byte w
        return (u.i[wi] >> bi) & 1;
    }
    // Member functions to split into two Vec256b:
    Vec256b get_low() const {
        return _mm512_castsi512_si256(zmm);
    }
    Vec256b get_high() const {
        return _mm512_extracti64x4_epi64(zmm,1);
    }
    static int size() {
        return 512;
    }
    static int elementtype() {
        return 1;
    }
};


// Define operators for this class

// vector operator & : bitwise and
static inline Vec512b operator & (Vec512b const & a, Vec512b const & b) {
    return _mm512_and_epi32(a, b);
}
static inline Vec512b operator && (Vec512b const & a, Vec512b const & b) {
    return a & b;
}

// vector operator | : bitwise or
static inline Vec512b operator | (Vec512b const & a, Vec512b const & b) {
    return _mm512_or_epi32(a, b);
}
static inline Vec512b operator || (Vec512b const & a, Vec512b const & b) {
    return a | b;
}

// vector operator ^ : bitwise xor
static inline Vec512b operator ^ (Vec512b const & a, Vec512b const & b) {
    return _mm512_xor_epi32(a, b);
}

// vector operator ~ : bitwise not
static inline Vec512b operator ~ (Vec512b const & a) {
    return _mm512_xor_epi32(a, _mm512_set1_epi32(-1));
}

// vector operator &= : bitwise and
static inline Vec512b & operator &= (Vec512b & a, Vec512b const & b) {
    a = a & b;
    return a;
}

// vector operator |= : bitwise or
static inline Vec512b & operator |= (Vec512b & a, Vec512b const & b) {
    a = a | b;
    return a;
}

// vector operator ^= : bitwise xor
static inline Vec512b & operator ^= (Vec512b & a, Vec512b const & b) {
    a = a ^ b;
    return a;
}

// Define functions for this class

// function andnot: a & ~ b
static inline Vec512b andnot (Vec512b const & a, Vec512b const & b) {
    return _mm512_andnot_epi32(b, a);
}


/*****************************************************************************
*
*          Vector of 16 32-bit signed integers
*
*****************************************************************************/

class Vec16i: public Vec512b {
public:
    // Default constructor:
    Vec16i() {
    }
    // Constructor to broadcast the same value into all elements:
    Vec16i(int i) {
        zmm = _mm512_set1_epi32(i);
    }
    // Constructor to build from all elements:
    Vec16i(int32_t i0, int32_t i1, int32_t i2, int32_t i3, int32_t i4, int32_t i5, int32_t i6, int32_t i7,
        int32_t i8, int32_t i9, int32_t i10, int32_t i11, int32_t i12, int32_t i13, int32_t i14, int32_t i15) {
        zmm = _mm512_setr_epi32(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15);
    }
    // Constructor to build from two Vec8i:
    Vec16i(Vec8i const & a0, Vec8i const & a1) {
        zmm = _mm512_inserti64x4(_mm512_castsi256_si512(a0), a1, 1);
    }
    // Constructor to convert from type __m512i used in intrinsics:
    Vec16i(__m512i const & x) {
        zmm = x;
    }
    // Assignment operator to convert from type __m512i used in intrinsics:
    Vec16i & operator = (__m512i const & x) {
        zmm = x;
        return *this;
    }
    // Constructor to convert from type Vec512b used in emulation:
    Vec16i(Vec512b const & x) {
        zmm = x;
    }
    // Type cast operator to convert to __m512i used in intrinsics
    operator __m512i() const {
        return zmm;
    }
    // Member function to load from array (unaligned)
    Vec16i & load(void const * p) {
        zmm = _mm512_loadu_si512(p);
        return *this;
    }
    // Member function to load from array, aligned by 64
    Vec16i & load_a(void const * p) {
        zmm = _mm512_load_si512(p);
        return *this;
    }
    // Partial load. Load n elements and set the rest to 0
    Vec16i & load_partial(int n, void const * p) {
        zmm = _mm512_maskz_loadu_epi32(__mmask16((1u << n) - 1), p);
        return *this;
    }
    // Partial store. Store n elements
    void store_partial(int n, void * p) const {
        _mm512_mask_storeu_epi32(p, __mmask16((1u << n) - 1), zmm);
    }
    // cut off vector to n elements. The last 16-n elements are set to zero
    Vec16i & cutoff(int n) {
        zmm = _mm512_maskz_mov_epi32(__mmask16((1u << n) - 1), zmm);
        return *this;
    }
    // Member function to change a single element in vector
    Vec16i const & insert(int index, int32_t value) {
        zmm = _mm512_mask_set1_epi32(zmm, __mmask16(1u << index), value);
        return *this;
    }
    // Member function extract a single element from vector
    int32_t extract(int index) const {
#if INSTRSET >= 10
        __m512i x = _mm512_maskz_compress_epi32(__mmask16(1u << index), zmm);
        return _mm_cvtsi128_si32(_mm512_castsi512_si128(x));        
#else 
        int32_t a[16];
        store(a);
        return a[index & 15];
#endif
    }
    // Extract a single element. Use store function if extracting more than one element.
    // Operator [] can only read an element, not write.
    int32_t operator [] (int index) const {
        return extract(index);
    }
    // Member functions to split into two Vec8i:
    Vec8i get_low() const {
        return _mm512_castsi512_si256(zmm);
    }
    Vec8i get_high() const {
        return _mm512_extracti64x4_epi64(zmm,1);
    }
    static int size () {
        return 16;
    }
    static int elementtype() {
        return 8;
    }
};


// Define operators for Vec16i

// vector operator + : add element by element
static inline Vec16i operator + (Vec16i const & a, Vec16i const & b) {
    return _mm512_add_epi32(a, b);
}

// vector operator += : add
static inline Vec16i & operator += (Vec16i & a, Vec16i const & b) {
    a = a + b;
    return a;
}

// postfix operator ++
static inline Vec16i operator ++ (Vec16i & a, int) {
    Vec16i a0 = a;
    a = a + 1;
    return a0;
}

// prefix operator ++
static inline Vec16i & operator ++ (Vec16i & a) {
    a = a + 1;
    return a;
}

// vector operator - : subtract element by element
static inline Vec16i operator - (Vec16i const & a, Vec16i const & b) {
    return _mm512_sub_epi32(a, b);
}

// vector operator - : unary minus
static inline Vec16i operator - (Vec16i const & a) {
    return _mm512_sub_epi32(_mm512_setzero_epi32(), a);
}

// vector operator -= : subtract
static inline Vec16i & operator -= (Vec16i & a, Vec16i const & b) {
    a = a - b;
    return a;
}

// postfix operator --
static inline Vec16i operator -- (Vec16i & a, int) {
    Vec16i a0 = a;
    a = a - 1;
    return a0;
}

// prefix operator --
static inline Vec16i & operator -- (Vec16i & a) {
    a = a - 1;
    return a;
}

// vector operator * : multiply element by element
static inline Vec16i operator * (Vec16i const & a, Vec16i const & b) {
    return _mm512_mullo_epi32(a, b);
}

// vector operator *= : multiply
static inline Vec16i & operator *= (Vec16i & a, Vec16i const & b) {
    a = a * b;
    return a;
}

// vector operator / : divide all elements by same integer
// See bottom of file


// vector operator << : shift left
static inline Vec16i operator << (Vec16i const & a, int32_t b) {
    return _mm512_sll_epi32(a, _mm_cvtsi32_si128(b));
}

// vector operator <<= : shift left
static inline Vec16i & operator <<= (Vec16i & a, int32_t b) {
    a = a << b;
    return a;
}

// vector operator >> : shift right arithmetic
static inline Vec16i operator >> (Vec16i const & a, int32_t b) {
    return _mm512_sra_epi32(a, _mm_cvtsi32_si128(b));
}

// vector operator >>= : shift right arithmetic
static inline Vec16i & operator >>= (Vec16i & a, int32_t b) {
    a = a >> b;
    return a;
}

// vector operator == : returns true for elements for which a == b
static inline Vec16ib operator == (Vec16i const & a, Vec16i const & b) {
    return _mm512_cmpeq_epi32_mask(a, b);
}

// vector operator != : returns true for elements for which a != b
static inline Vec16ib operator != (Vec16i const & a, Vec16i const & b) {
    return _mm512_cmpneq_epi32_mask(a, b);
}
  
// vector operator > : returns true for elements for which a > b
static inline Vec16ib operator > (Vec16i const & a, Vec16i const & b) {
    return  _mm512_cmpgt_epi32_mask(a, b);
}

// vector operator < : returns true for elements for which a < b
static inline Vec16ib operator < (Vec16i const & a, Vec16i const & b) {
    return b > a;
}

// vector operator >= : returns true for elements for which a >= b (signed)
static inline Vec16ib operator >= (Vec16i const & a, Vec16i const & b) {
    return _mm512_cmpge_epi32_mask(a, b);
}

// vector operator <= : returns true for elements for which a <= b (signed)
static inline Vec16ib operator <= (Vec16i const & a, Vec16i const & b) {
    return b >= a;
}

// vector operator & : bitwise and
static inline Vec16i operator & (Vec16i const & a, Vec16i const & b) {
    return _mm512_and_epi32(a, b);
}

// vector operator &= : bitwise and
static inline Vec16i & operator &= (Vec16i & a, Vec16i const & b) {
    a = a & b;
    return a;
}

// vector operator | : bitwise or
static inline Vec16i operator | (Vec16i const & a, Vec16i const & b) {
    return _mm512_or_epi32(a, b);
}

// vector operator |= : bitwise or
static inline Vec16i & operator |= (Vec16i & a, Vec16i const & b) {
    a = a | b;
    return a;
}

// vector operator ^ : bitwise xor
static inline Vec16i operator ^ (Vec16i const & a, Vec16i const & b) {
    return _mm512_xor_epi32(a, b);
}

// vector operator ^= : bitwise xor
static inline Vec16i & operator ^= (Vec16i & a, Vec16i const & b) {
    a = a ^ b;
    return a;
}

// vector operator ~ : bitwise not
static inline Vec16i operator ~ (Vec16i const & a) {
    return a ^ Vec16i(-1);
}

// Functions for this class

// Select between two operands. Corresponds to this pseudocode:
// for (int i = 0; i < 16; i++) result[i] = s[i] ? a[i] : b[i];
static inline Vec16i select (Vec16ib const & s, Vec16i const & a, Vec16i const & b) {
    return _mm512_mask_mov_epi32(b, s, a);  // conditional move may be optimized better by the compiler than blend
    // return _mm512_mask_blend_epi32(s, b, a);
}

// Conditional add: For all vector elements i: result[i] = f[i] ? (a[i] + b[i]) : a[i]
static inline Vec16i if_add (Vec16ib const & f, Vec16i const & a, Vec16i const & b) {
    return _mm512_mask_add_epi32(a, f, a, b);
}

// Conditional sub: For all vector elements i: result[i] = f[i] ? (a[i] - b[i]) : a[i]
static inline Vec16i if_sub (Vec16ib const & f, Vec16i const & a, Vec16i const & b) {
    return _mm512_mask_sub_epi32(a, f, a, b);
}

// Conditional mul: For all vector elements i: result[i] = f[i] ? (a[i] * b[i]) : a[i]
static inline Vec16i if_mul (Vec16ib const & f, Vec16i const & a, Vec16i const & b) {
    return _mm512_mask_mullo_epi32(a, f, a, b);
}

// Horizontal add: Calculates the sum of all vector elements.
// Overflow will wrap around
static inline int32_t horizontal_add (Vec16i const & a) {
#if defined(__INTEL_COMPILER)
    return _mm512_reduce_add_epi32(a);
#else
    return horizontal_add(a.get_low() + a.get_high());
#endif
}

// function add_saturated: add element by element, signed with saturation
// (is it faster to up-convert to 64 bit integers, and then downconvert the sum with saturation?)
static inline Vec16i add_saturated(Vec16i const & a, Vec16i const & b) {
    __m512i sum    = _mm512_add_epi32(a, b);                  // a + b
    __m512i axb    = _mm512_xor_epi32(a, b);                  // check if a and b have different sign
    __m512i axs    = _mm512_xor_epi32(a, sum);                // check if a and sum have different sign
    __m512i ovf1   = _mm512_andnot_epi32(axb,axs);            // check if sum has wrong sign
    __m512i ovf2   = _mm512_srai_epi32(ovf1,31);              // -1 if overflow
    __mmask16 ovf3 = _mm512_cmpneq_epi32_mask(ovf2, _mm512_setzero_epi32()); // same, as mask
    __m512i asign  = _mm512_srli_epi32(a,31);                 // 1  if a < 0
    __m512i sat1   = _mm512_srli_epi32(ovf2,1);               // 7FFFFFFF if overflow
    __m512i sat2   = _mm512_add_epi32(sat1,asign);            // 7FFFFFFF if positive overflow 80000000 if negative overflow
    return _mm512_mask_blend_epi32(ovf3, sum, sat2);          // sum if not overflow, else sat2
}

// function sub_saturated: subtract element by element, signed with saturation
static inline Vec16i sub_saturated(Vec16i const & a, Vec16i const & b) {
    __m512i diff   = _mm512_sub_epi32(a, b);                  // a + b
    __m512i axb    = _mm512_xor_si512(a, b);                  // check if a and b have different sign
    __m512i axs    = _mm512_xor_si512(a, diff);               // check if a and sum have different sign
    __m512i ovf1   = _mm512_and_si512(axb,axs);               // check if sum has wrong sign
    __m512i ovf2   = _mm512_srai_epi32(ovf1,31);              // -1 if overflow
    __mmask16 ovf3 = _mm512_cmpneq_epi32_mask(ovf2, _mm512_setzero_epi32()); // same, as mask
    __m512i asign  = _mm512_srli_epi32(a,31);                 // 1  if a < 0
    __m512i sat1   = _mm512_srli_epi32(ovf2,1);               // 7FFFFFFF if overflow
    __m512i sat2   = _mm512_add_epi32(sat1,asign);            // 7FFFFFFF if positive overflow 80000000 if negative overflow
    return _mm512_mask_blend_epi32(ovf3, diff, sat2);         // sum if not overflow, else sat2
}

// function max: a > b ? a : b
static inline Vec16i max(Vec16i const & a, Vec16i const & b) {
    return _mm512_max_epi32(a,b);
}

// function min: a < b ? a : b
static inline Vec16i min(Vec16i const & a, Vec16i const & b) {
    return _mm512_min_epi32(a,b);
}

// function abs: a >= 0 ? a : -a
static inline Vec16i abs(Vec16i const & a) {
    return _mm512_abs_epi32(a);
}

// function abs_saturated: same as abs, saturate if overflow
static inline Vec16i abs_saturated(Vec16i const & a) {
    return _mm512_min_epu32(abs(a), Vec16i(0x7FFFFFFF));
}

// function rotate_left all elements
// Use negative count to rotate right
static inline Vec16i rotate_left(Vec16i const & a, int b) {
    return _mm512_rolv_epi32(a, Vec16i(b));
}


/*****************************************************************************
*
*          Vector of 16 32-bit unsigned integers
*
*****************************************************************************/


class Vec16ui : public Vec16i {
public:
    // Default constructor:
    Vec16ui() {
    }
    // Constructor to broadcast the same value into all elements:
    Vec16ui(uint32_t i) {
        zmm = _mm512_set1_epi32((int32_t)i);
    }
    // Constructor to build from all elements:
    Vec16ui(uint32_t i0, uint32_t i1, uint32_t i2, uint32_t i3, uint32_t i4, uint32_t i5, uint32_t i6, uint32_t i7, 
        uint32_t i8, uint32_t i9, uint32_t i10, uint32_t i11, uint32_t i12, uint32_t i13, uint32_t i14, uint32_t i15) {
        zmm = _mm512_setr_epi32((int32_t)i0, (int32_t)i1, (int32_t)i2, (int32_t)i3, (int32_t)i4, (int32_t)i5, (int32_t)i6, (int32_t)i7, 
            (int32_t)i8, (int32_t)i9, (int32_t)i10, (int32_t)i11, (int32_t)i12, (int32_t)i13, (int32_t)i14, (int32_t)i15);
    }
    // Constructor to build from two Vec8ui:
    Vec16ui(Vec8ui const & a0, Vec8ui const & a1) {
        zmm = Vec16i(Vec8i(a0), Vec8i(a1));
    }
    // Constructor to convert from type __m512i used in intrinsics:
    Vec16ui(__m512i const & x) {
        zmm = x;
    }
    // Assignment operator to convert from type __m512i used in intrinsics:
    Vec16ui & operator = (__m512i const & x) {
        zmm = x;
        return *this;
    }
    // Member function to load from array (unaligned)
    Vec16ui & load(void const * p) {
        Vec16i::load(p);
        return *this;
    }
    // Member function to load from array, aligned by 64
    Vec16ui & load_a(void const * p) {
        Vec16i::load_a(p);
        return *this;
    }
    // Member function to change a single element in vector
    Vec16ui const & insert(int index, uint32_t value) {
        Vec16i::insert(index, (int32_t)value);
        return *this;
    }
    // Member function extract a single element from vector
    uint32_t extract(int index) const {
        return (uint32_t)Vec16i::extract(index);
    }
    // Extract a single element. Use store function if extracting more than one element.
    // Operator [] can only read an element, not write.
    uint32_t operator [] (int index) const {
        return extract(index);
    }
    // Member functions to split into two Vec4ui:
    Vec8ui get_low() const {
        return Vec8ui(Vec16i::get_low());
    }
    Vec8ui get_high() const {
        return Vec8ui(Vec16i::get_high());
    }
    static int elementtype() {
        return 9;
    }
};

// Define operators for this class

// vector operator + : add
static inline Vec16ui operator + (Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui (Vec16i(a) + Vec16i(b));
}

// vector operator - : subtract
static inline Vec16ui operator - (Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui (Vec16i(a) - Vec16i(b));
}

// vector operator * : multiply
static inline Vec16ui operator * (Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui (Vec16i(a) * Vec16i(b));
}

// vector operator / : divide
// See bottom of file

// vector operator >> : shift right logical all elements
static inline Vec16ui operator >> (Vec16ui const & a, uint32_t b) {
    return _mm512_srl_epi32(a, _mm_cvtsi32_si128((int32_t)b)); 
}
static inline Vec16ui operator >> (Vec16ui const & a, int32_t b) {
    return a >> (uint32_t)b;
}

// vector operator >>= : shift right logical
static inline Vec16ui & operator >>= (Vec16ui & a, uint32_t b) {
    a = a >> b;
    return a;
} 

// vector operator >>= : shift right logical
static inline Vec16ui & operator >>= (Vec16ui & a, int32_t b) {
    a = a >> uint32_t(b);
    return a;
}

// vector operator << : shift left all elements
static inline Vec16ui operator << (Vec16ui const & a, uint32_t b) {
    return Vec16ui ((Vec16i)a << (int32_t)b);
}

// vector operator << : shift left all elements
static inline Vec16ui operator << (Vec16ui const & a, int32_t b) {
    return Vec16ui ((Vec16i)a << (int32_t)b);
}

// vector operator < : returns true for elements for which a < b (unsigned)
static inline Vec16ib operator < (Vec16ui const & a, Vec16ui const & b) {
    return _mm512_cmplt_epu32_mask(a, b);
}

// vector operator > : returns true for elements for which a > b (unsigned)
static inline Vec16ib operator > (Vec16ui const & a, Vec16ui const & b) {
    return b < a;
}


// vector operator >= : returns true for elements for which a >= b (unsigned)
static inline Vec16ib operator >= (Vec16ui const & a, Vec16ui const & b) {
    return  _mm512_cmpge_epu32_mask(a, b);
}            

// vector operator <= : returns true for elements for which a <= b (unsigned)
static inline Vec16ib operator <= (Vec16ui const & a, Vec16ui const & b) {
    return b >= a;
}

// vector operator & : bitwise and
static inline Vec16ui operator & (Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui(Vec16i(a) & Vec16i(b));
}

// vector operator | : bitwise or
static inline Vec16ui operator | (Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui(Vec16i(a) | Vec16i(b));
}

// vector operator ^ : bitwise xor
static inline Vec16ui operator ^ (Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui(Vec16i(a) ^ Vec16i(b));
}

// vector operator ~ : bitwise not
static inline Vec16ui operator ~ (Vec16ui const & a) {
    return Vec16ui( ~ Vec16i(a));
}

// Functions for this class

// Select between two operands. Corresponds to this pseudocode:
// for (int i = 0; i < 16; i++) result[i] = s[i] ? a[i] : b[i];
static inline Vec16ui select (Vec16ib const & s, Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui(select(s, Vec16i(a), Vec16i(b)));
}

// Conditional add: For all vector elements i: result[i] = f[i] ? (a[i] + b[i]) : a[i]
static inline Vec16ui if_add (Vec16ib const & f, Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui(if_add(f, Vec16i(a), Vec16i(b)));
}

// Conditional sub: For all vector elements i: result[i] = f[i] ? (a[i] - b[i]) : a[i]
static inline Vec16ui if_sub (Vec16ib const & f, Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui(if_sub(f, Vec16i(a), Vec16i(b)));
}

// Conditional mul: For all vector elements i: result[i] = f[i] ? (a[i] * b[i]) : a[i]
static inline Vec16ui if_mul (Vec16ib const & f, Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui(if_mul(f, Vec16i(a), Vec16i(b)));
}

// Horizontal add: Calculates the sum of all vector elements.
// Overflow will wrap around
static inline uint32_t horizontal_add (Vec16ui const & a) {
    return (uint32_t)horizontal_add((Vec16i)a);
}

// horizontal_add_x: Horizontal add extended: Calculates the sum of all vector elements. Defined later in this file

// function add_saturated: add element by element, unsigned with saturation
static inline Vec16ui add_saturated(Vec16ui const & a, Vec16ui const & b) {
    Vec16ui sum      = a + b;
    Vec16ib overflow = sum < (a | b);                  // overflow if (a + b) < (a | b)
    return _mm512_mask_set1_epi32(sum, overflow, -1);  // 0xFFFFFFFF if overflow
}

// function sub_saturated: subtract element by element, unsigned with saturation
static inline Vec16ui sub_saturated(Vec16ui const & a, Vec16ui const & b) {
    Vec16ui diff      = a - b;
    return _mm512_maskz_mov_epi32(diff <= a, diff);   // underflow if diff > a gives zero
}

// function max: a > b ? a : b
static inline Vec16ui max(Vec16ui const & a, Vec16ui const & b) {
    return _mm512_max_epu32(a,b);
}

// function min: a < b ? a : b
static inline Vec16ui min(Vec16ui const & a, Vec16ui const & b) {
    return _mm512_min_epu32(a,b);
}


/*****************************************************************************
*
*          Vector of 8 64-bit signed integers
*
*****************************************************************************/

class Vec8q : public Vec512b {
public:
    // Default constructor:
    Vec8q() {
    }
    // Constructor to broadcast the same value into all elements:
    Vec8q(int64_t i) {
        zmm = _mm512_set1_epi64(i);
    }
    // Constructor to build from all elements:
    Vec8q(int64_t i0, int64_t i1, int64_t i2, int64_t i3, int64_t i4, int64_t i5, int64_t i6, int64_t i7) {
        zmm = _mm512_setr_epi64(i0, i1, i2, i3, i4, i5, i6, i7);
    }
    // Constructor to build from two Vec4q:
    Vec8q(Vec4q const & a0, Vec4q const & a1) {
        zmm = _mm512_inserti64x4(_mm512_castsi256_si512(a0), a1, 1);
    }
    // Constructor to convert from type __m512i used in intrinsics:
    Vec8q(__m512i const & x) {
        zmm = x;
    }
    // Assignment operator to convert from type __m512i used in intrinsics:
    Vec8q & operator = (__m512i const & x) {
        zmm = x;
        return *this;
    }
    // Type cast operator to convert to __m512i used in intrinsics
    operator __m512i() const {
        return zmm;
    }
    // Member function to load from array (unaligned)
    Vec8q & load(void const * p) {
        zmm = _mm512_loadu_si512(p);
        return *this;
    }
    // Member function to load from array, aligned by 64
    Vec8q & load_a(void const * p) {
        zmm = _mm512_load_si512(p);
        return *this;
    }
    // Partial load. Load n elements and set the rest to 0
    Vec8q & load_partial(int n, void const * p) {
        zmm = _mm512_maskz_loadu_epi64(__mmask16((1 << n) - 1), p);
        return *this;
    }
    // Partial store. Store n elements
    void store_partial(int n, void * p) const {
        _mm512_mask_storeu_epi64(p, __mmask16((1 << n) - 1), zmm);
    }
    // cut off vector to n elements. The last 8-n elements are set to zero
    Vec8q & cutoff(int n) {
        zmm = _mm512_maskz_mov_epi64(__mmask16((1 << n) - 1), zmm);
        return *this;
    }
    // Member function to change a single element in vector
    Vec8q const & insert(int index, int64_t value) {
        zmm = _mm512_mask_set1_epi64(zmm, __mmask16(1 << index), value);
        return *this;
    }
    // Member function extract a single element from vector
    int64_t extract(int index) const {
#if INSTRSET >= 10
        __m512i x = _mm512_maskz_compress_epi64(__mmask8(1u << index), zmm);
        return _emulate_movq(_mm512_castsi512_si128(x));        
#else 
        int64_t a[8];
        store (a);
        return a[index & 7];
#endif
    }
    // Extract a single element. Use store function if extracting more than one element.
    // Operator [] can only read an element, not write.
    int64_t operator [] (int index) const {
        return extract(index);
    }
    // Member functions to split into two Vec2q:
    Vec4q get_low() const {
        return _mm512_castsi512_si256(zmm);
    }
    Vec4q get_high() const {
        return _mm512_extracti64x4_epi64(zmm,1);
    }
    static int size () {
        return 8;
    }
    static int elementtype() {
        return 10;
    }
};


// Define operators for Vec8q

// vector operator + : add element by element
static inline Vec8q operator + (Vec8q const & a, Vec8q const & b) {
    return _mm512_add_epi64(a, b);
}

// vector operator += : add
static inline Vec8q & operator += (Vec8q & a, Vec8q const & b) {
    a = a + b;
    return a;
}

// postfix operator ++
static inline Vec8q operator ++ (Vec8q & a, int) {
    Vec8q a0 = a;
    a = a + 1;
    return a0;
}

// prefix operator ++
static inline Vec8q & operator ++ (Vec8q & a) {
    a = a + 1;
    return a;
}

// vector operator - : subtract element by element
static inline Vec8q operator - (Vec8q const & a, Vec8q const & b) {
    return _mm512_sub_epi64(a, b);
}

// vector operator - : unary minus
static inline Vec8q operator - (Vec8q const & a) {
    return _mm512_sub_epi64(_mm512_setzero_epi32(), a);
}

// vector operator -= : subtract
static inline Vec8q & operator -= (Vec8q & a, Vec8q const & b) {
    a = a - b;
    return a;
}

// postfix operator --
static inline Vec8q operator -- (Vec8q & a, int) {
    Vec8q a0 = a;
    a = a - 1;
    return a0;
}

// prefix operator --
static inline Vec8q & operator -- (Vec8q & a) {
    a = a - 1;
    return a;
}

// vector operator * : multiply element by element
static inline Vec8q operator * (Vec8q const & a, Vec8q const & b) {
#if INSTRSET >= 10  // __AVX512DQ__
    return _mm512_mullo_epi64(a, b);
#elif defined (__INTEL_COMPILER)
    return _mm512_mullox_epi64(a, b);                      // _mm512_mullox_epi64 missing in gcc
#else
    // instruction does not exist. Split into 32-bit multiplies
    //__m512i ahigh = _mm512_shuffle_epi32(a, 0xB1);       // swap H<->L
    __m512i ahigh   = _mm512_srli_epi64(a, 32);            // high 32 bits of each a
    __m512i bhigh   = _mm512_srli_epi64(b, 32);            // high 32 bits of each b
    __m512i prodahb = _mm512_mul_epu32(ahigh, b);          // ahigh*b
    __m512i prodbha = _mm512_mul_epu32(bhigh, a);          // bhigh*a
    __m512i prodhl  = _mm512_add_epi64(prodahb, prodbha);  // sum of high*low products
    __m512i prodhi  = _mm512_slli_epi64(prodhl, 32);       // same, shifted high
    __m512i prodll  = _mm512_mul_epu32(a, b);              // alow*blow = 64 bit unsigned products
    __m512i prod    = _mm512_add_epi64(prodll, prodhi);    // low*low+(high*low)<<32
    return  prod;
#endif
}

// vector operator *= : multiply
static inline Vec8q & operator *= (Vec8q & a, Vec8q const & b) {
    a = a * b;
    return a;
}

// vector operator << : shift left
static inline Vec8q operator << (Vec8q const & a, int32_t b) {
    return _mm512_sll_epi64(a, _mm_cvtsi32_si128(b));
}

// vector operator <<= : shift left
static inline Vec8q & operator <<= (Vec8q & a, int32_t b) {
    a = a << b;
    return a;
}

// vector operator >> : shift right arithmetic
static inline Vec8q operator >> (Vec8q const & a, int32_t b) {
    return _mm512_sra_epi64(a, _mm_cvtsi32_si128(b));
}

// vector operator >>= : shift right arithmetic
static inline Vec8q & operator >>= (Vec8q & a, int32_t b) {
    a = a >> b;
    return a;
}

// vector operator == : returns true for elements for which a == b
static inline Vec8qb operator == (Vec8q const & a, Vec8q const & b) {
    return Vec8qb(_mm512_cmpeq_epi64_mask(a, b));
}

// vector operator != : returns true for elements for which a != b
static inline Vec8qb operator != (Vec8q const & a, Vec8q const & b) {
    return Vec8qb(_mm512_cmpneq_epi64_mask(a, b));
}
  
// vector operator < : returns true for elements for which a < b
static inline Vec8qb operator < (Vec8q const & a, Vec8q const & b) {
    return Vec8qb(_mm512_cmplt_epi64_mask(a, b));
}

// vector operator > : returns true for elements for which a > b
static inline Vec8qb operator > (Vec8q const & a, Vec8q const & b) {
    return b < a;
}

// vector operator >= : returns true for elements for which a >= b (signed)
static inline Vec8qb operator >= (Vec8q const & a, Vec8q const & b) {
    return Vec8qb(_mm512_cmpge_epi64_mask(a, b));
}

// vector operator <= : returns true for elements for which a <= b (signed)
static inline Vec8qb operator <= (Vec8q const & a, Vec8q const & b) {
    return b >= a;
}

// vector operator & : bitwise and
static inline Vec8q operator & (Vec8q const & a, Vec8q const & b) {
    return _mm512_and_epi32(a, b);
}

// vector operator &= : bitwise and
static inline Vec8q & operator &= (Vec8q & a, Vec8q const & b) {
    a = a & b;
    return a;
}

// vector operator | : bitwise or
static inline Vec8q operator | (Vec8q const & a, Vec8q const & b) {
    return _mm512_or_epi32(a, b);
}

// vector operator |= : bitwise or
static inline Vec8q & operator |= (Vec8q & a, Vec8q const & b) {
    a = a | b;
    return a;
}

// vector operator ^ : bitwise xor
static inline Vec8q operator ^ (Vec8q const & a, Vec8q const & b) {
    return _mm512_xor_epi32(a, b);
}
// vector operator ^= : bitwise xor
static inline Vec8q & operator ^= (Vec8q & a, Vec8q const & b) {
    a = a ^ b;
    return a;
}

// vector operator ~ : bitwise not
static inline Vec8q operator ~ (Vec8q const & a) {
    return Vec8q(~ Vec16i(a));
}

// Functions for this class

// Select between two operands. Corresponds to this pseudocode:
// for (int i = 0; i < 4; i++) result[i] = s[i] ? a[i] : b[i];
static inline Vec8q select (Vec8qb const & s, Vec8q const & a, Vec8q const & b) {
    return _mm512_mask_mov_epi64(b, s, a);
    //return _mm512_mask_blend_epi64(s, b, a);
}

// Conditional add: For all vector elements i: result[i] = f[i] ? (a[i] + b[i]) : a[i]
static inline Vec8q if_add (Vec8qb const & f, Vec8q const & a, Vec8q const & b) {
    return _mm512_mask_add_epi64(a, f, a, b);
}

// Conditional sub: For all vector elements i: result[i] = f[i] ? (a[i] - b[i]) : a[i]
static inline Vec8q if_sub (Vec8qb const & f, Vec8q const & a, Vec8q const & b) {
    return _mm512_mask_sub_epi64(a, f, a, b);
}

// Conditional mul: For all vector elements i: result[i] = f[i] ? (a[i] * b[i]) : a[i]
static inline Vec8q if_mul (Vec8qb const & f, Vec8q const & a, Vec8q const & b) {
#if INSTRSET >= 10
    return _mm512_mask_mullo_epi64(a, f, a, b);  // AVX512DQ
#else
    return select(f, a*b, a);
#endif
}

// Horizontal add: Calculates the sum of all vector elements.
// Overflow will wrap around
static inline int64_t horizontal_add (Vec8q const & a) {
#if defined(__INTEL_COMPILER)
    return _mm512_reduce_add_epi64(a);
#else
    return horizontal_add(a.get_low()+a.get_high());
#endif
}

// Horizontal add extended: Calculates the sum of all vector elements
// Elements are sign extended before adding to avoid overflow
static inline int64_t horizontal_add_x (Vec16i const & x) {
    Vec8q a = _mm512_cvtepi32_epi64(x.get_low());
    Vec8q b = _mm512_cvtepi32_epi64(x.get_high());
    return horizontal_add(a+b);
}

// Horizontal add extended: Calculates the sum of all vector elements
// Elements are zero extended before adding to avoid overflow
static inline uint64_t horizontal_add_x (Vec16ui const & x) {
    Vec8q a = _mm512_cvtepu32_epi64(x.get_low());
    Vec8q b = _mm512_cvtepu32_epi64(x.get_high());
    return (uint64_t)horizontal_add(a+b);
}

// function max: a > b ? a : b
static inline Vec8q max(Vec8q const & a, Vec8q const & b) {
    return _mm512_max_epi64(a, b);
}

// function min: a < b ? a : b
static inline Vec8q min(Vec8q const & a, Vec8q const & b) {
    return _mm512_min_epi64(a, b);
}

// function abs: a >= 0 ? a : -a
static inline Vec8q abs(Vec8q const & a) {
    return _mm512_abs_epi64(a);
}

// function abs_saturated: same as abs, saturate if overflow
static inline Vec8q abs_saturated(Vec8q const & a) {
    return _mm512_min_epu64(abs(a), Vec8q(0x7FFFFFFFFFFFFFFF));
}

// function rotate_left all elements
// Use negative count to rotate right
static inline Vec8q rotate_left(Vec8q const & a, int b) {
    return _mm512_rolv_epi64(a, Vec8q(b));
}


/*****************************************************************************
*
*          Vector of 8 64-bit unsigned integers
*
*****************************************************************************/

class Vec8uq : public Vec8q {
public:
    // Default constructor:
    Vec8uq() {
    }
    // Constructor to broadcast the same value into all elements:
    Vec8uq(uint64_t i) {
        zmm = Vec8q((int64_t)i);
    }
    // Constructor to convert from Vec8q:
    Vec8uq(Vec8q const & x) {
        zmm = x;
    }
    // Constructor to convert from type __m512i used in intrinsics:
    Vec8uq(__m512i const & x) {
        zmm = x;
    }
    // Constructor to build from all elements:
    Vec8uq(uint64_t i0, uint64_t i1, uint64_t i2, uint64_t i3, uint64_t i4, uint64_t i5, uint64_t i6, uint64_t i7) {
        zmm = Vec8q((int64_t)i0, (int64_t)i1, (int64_t)i2, (int64_t)i3, (int64_t)i4, (int64_t)i5, (int64_t)i6, (int64_t)i7);
    }
    // Constructor to build from two Vec4uq:
    Vec8uq(Vec4uq const & a0, Vec4uq const & a1) {
        zmm = Vec8q(Vec4q(a0), Vec4q(a1));
    }
    // Assignment operator to convert from Vec8q:
    Vec8uq  & operator = (Vec8q const & x) {
        zmm = x;
        return *this;
    }
    // Assignment operator to convert from type __m512i used in intrinsics:
    Vec8uq & operator = (__m512i const & x) {
        zmm = x;
        return *this;
    }
    // Member function to load from array (unaligned)
    Vec8uq & load(void const * p) {
        Vec8q::load(p);
        return *this;
    }
    // Member function to load from array, aligned by 32
    Vec8uq & load_a(void const * p) {
        Vec8q::load_a(p);
        return *this;
    }
    // Member function to change a single element in vector
    Vec8uq const & insert(int index, uint64_t value) {
        Vec8q::insert(index, (int64_t)value);
        return *this;
    }
    // Member function extract a single element from vector
    uint64_t extract(int index) const {
        return (uint64_t)Vec8q::extract(index);
    }
    // Extract a single element. Use store function if extracting more than one element.
    // Operator [] can only read an element, not write.
    uint64_t operator [] (int index) const {
        return extract(index);
    }
    // Member functions to split into two Vec2uq:
    Vec4uq get_low() const {
        return Vec4uq(Vec8q::get_low());
    }
    Vec4uq get_high() const {
        return Vec4uq(Vec8q::get_high());
    }
    static int elementtype() {
        return 10;
    }
};

// Define operators for this class

// vector operator + : add
static inline Vec8uq operator + (Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq (Vec8q(a) + Vec8q(b));
}

// vector operator - : subtract
static inline Vec8uq operator - (Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq (Vec8q(a) - Vec8q(b));
}

// vector operator * : multiply element by element
static inline Vec8uq operator * (Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq (Vec8q(a) * Vec8q(b));
}

// vector operator >> : shift right logical all elements
static inline Vec8uq operator >> (Vec8uq const & a, uint32_t b) {
    return _mm512_srl_epi64(a,_mm_cvtsi32_si128((int32_t)b)); 
}
static inline Vec8uq operator >> (Vec8uq const & a, int32_t b) {
    return a >> (uint32_t)b;
}

// vector operator >>= : shift right artihmetic
static inline Vec8uq & operator >>= (Vec8uq & a, uint32_t b) {
    a = a >> b;
    return a;
}

// vector operator >>= : shift right logical
static inline Vec8uq & operator >>= (Vec8uq & a, int32_t b) {
    a = a >> uint32_t(b);
    return a;
}

// vector operator << : shift left all elements
static inline Vec8uq operator << (Vec8uq const & a, uint32_t b) {
    return Vec8uq ((Vec8q)a << (int32_t)b);
}

// vector operator << : shift left all elements
static inline Vec8uq operator << (Vec8uq const & a, int32_t b) {
    return Vec8uq ((Vec8q)a << b);
}

// vector operator < : returns true for elements for which a < b (unsigned)
static inline Vec8qb operator < (Vec8uq const & a, Vec8uq const & b) {
    return _mm512_cmplt_epu64_mask(a, b);
}

// vector operator > : returns true for elements for which a > b (unsigned)
static inline Vec8qb operator > (Vec8uq const & a, Vec8uq const & b) {
    return b < a;
}

// vector operator >= : returns true for elements for which a >= b (unsigned)
static inline Vec8qb operator >= (Vec8uq const & a, Vec8uq const & b) {
    return _mm512_cmpge_epu64_mask(a, b);
}

// vector operator <= : returns true for elements for which a <= b (unsigned)
static inline Vec8qb operator <= (Vec8uq const & a, Vec8uq const & b) {
    return b >= a;
}

// vector operator & : bitwise and
static inline Vec8uq operator & (Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq(Vec8q(a) & Vec8q(b));
}

// vector operator | : bitwise or
static inline Vec8uq operator | (Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq(Vec8q(a) | Vec8q(b));
}

// vector operator ^ : bitwise xor
static inline Vec8uq operator ^ (Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq(Vec8q(a) ^ Vec8q(b));
}

// Functions for this class

// Select between two operands. Corresponds to this pseudocode:
// for (int i = 0; i < 4; i++) result[i] = s[i] ? a[i] : b[i];
static inline Vec8uq select (Vec8qb const & s, Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq(select(s, Vec8q(a), Vec8q(b)));
}

// Conditional add: For all vector elements i: result[i] = f[i] ? (a[i] + b[i]) : a[i]
static inline Vec8uq if_add (Vec8qb const & f, Vec8uq const & a, Vec8uq const & b) {
    return _mm512_mask_add_epi64(a, f, a, b);
}

// Conditional sub: For all vector elements i: result[i] = f[i] ? (a[i] - b[i]) : a[i]
static inline Vec8uq if_sub (Vec8qb const & f, Vec8uq const & a, Vec8uq const & b) {
    return _mm512_mask_sub_epi64(a, f, a, b);
}

// Conditional mul: For all vector elements i: result[i] = f[i] ? (a[i] * b[i]) : a[i]
static inline Vec8uq if_mul (Vec8qb const & f, Vec8uq const & a, Vec8uq const & b) {
#if INSTRSET >= 10
    return _mm512_mask_mullo_epi64(a, f, a, b);  // AVX512DQ
#else
    return select(f, a*b, a);
#endif
}

// Horizontal add: Calculates the sum of all vector elements.
// Overflow will wrap around
static inline uint64_t horizontal_add (Vec8uq const & a) {
    return (uint64_t)horizontal_add(Vec8q(a));
}

// function max: a > b ? a : b
static inline Vec8uq max(Vec8uq const & a, Vec8uq const & b) {
    return _mm512_max_epu64(a, b);
}

// function min: a < b ? a : b
static inline Vec8uq min(Vec8uq const & a, Vec8uq const & b) {
    return _mm512_min_epu64(a, b);
}


/*****************************************************************************
*
*          Vector permute functions
*
******************************************************************************
*
* These permute functions can reorder the elements of a vector and optionally
* set some elements to zero. 
*
* The indexes are inserted as template parameters in <>. These indexes must be
* constants. Each template parameter is an index to the element you want to select.
* An index of -1 will generate zero. An index of -256 means don't care.
*
* Example:
* Vec8q a(10,11,12,13,14,15,16,17);      // a is (10,11,12,13,14,15,16,17)
* Vec8q b;
* b = permute8<0,2,7,7,-1,-1,1,1>(a);   // b is (10,12,17,17, 0, 0,11,11)
*
* A lot of the code here is metaprogramming aiming to find the instructions
* that best fit the template parameters and instruction set. The metacode
* will be reduced out to leave only a few vector instructions in release
* mode with optimization on.
*****************************************************************************/

// Permute vector of 8 64-bit integers.
// Index -1 gives 0, index -256 means don't care.
template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7>
static inline Vec8q permute8(Vec8q const & a) {

    // Combine indexes into a single bitfield, with 4 bits for each
    const int m1 = (i0&7) | (i1&7)<<4 | (i2&7)<< 8 | (i3&7)<<12 | (i4&7)<<16 | (i5&7)<<20 | (i6&7)<<24 | (i7&7)<<28;

    // Mask to zero out negative indexes
    const int mz = (i0<0?0:0xF) | (i1<0?0:0xF0) | (i2<0?0:0xF00) | (i3<0?0:0xF000) | (i4<0?0:0xF0000) | (i5<0?0:0xF00000) | (i6<0?0:0xF000000) | (i7<0?0:0xF0000000);
    const int m2 = m1 & mz;

    // zeroing needed
    const bool dozero = ((i0|i1|i2|i3|i4|i5|i6|i7) & 0x80) != 0;

    // special case: all zero
    if (mz == 0) return  _mm512_setzero_epi32();

    // mask for elements not zeroed
    const __mmask16  z = __mmask16((i0>=0)<<0 | (i1>=0)<<1 | (i2>=0)<<2 | (i3>=0)<<3 | (i4>=0)<<4 | (i5>=0)<<5 | (i6>=0)<<6 | (i7>=0)<<7);
    // same with 2 bits for each element
    const __mmask16 zz = __mmask16((i0>=0?3:0) | (i1>=0?0xC:0) | (i2>=0?0x30:0) | (i3>=0?0xC0:0) | (i4>=0?0x300:0) | (i5>=0?0xC00:0) | (i6>=0?0x3000:0) | (i7>=0?0xC000:0));

    if (((m1 ^ 0x76543210) & mz) == 0) {
        // no shuffling
        if (dozero) {
            // zero some elements
            return _mm512_maskz_mov_epi64(z, a);
        }
        return a;                                 // do nothing
    }

    if (((m1 ^ 0x66442200) & 0x66666666 & mz) == 0) {
        // no exchange of data between the four 128-bit lanes
        const int pat = ((m2 | m2 >> 8 | m2 >> 16 | m2 >> 24) & 0x11) * 0x01010101;
        const int pmask = ((pat & 1) * 10 + 4) | ((((pat >> 4) & 1) * 10 + 4) << 4);
        if (((m1 ^ pat) & mz & 0x11111111) == 0) {
            // same permute pattern in all lanes
            if (dozero) {  // permute within lanes and zero
                return _mm512_maskz_shuffle_epi32(zz, a, (_MM_PERM_ENUM)pmask);
            }
            else {  // permute within lanes
                return _mm512_shuffle_epi32(a, (_MM_PERM_ENUM)pmask);
            }
        }
        // different permute patterns in each lane. It's faster to do a full permute than four masked permutes within lanes
    }
    if ((((m1 ^ 0x10101010) & 0x11111111 & mz) == 0) 
    &&  ((m1 ^ (m1 >> 4)) & 0x06060606 & mz & (mz >> 4)) == 0) {
        // permute lanes only. no permutation within each lane
        const int m3 = m2 | (m2 >> 4);
        const int s = ((m3 >> 1) & 3) | (((m3 >> 9) & 3) << 2) | (((m3 >> 17) & 3) << 4) | (((m3 >> 25) & 3) << 6);
        if (dozero) {
            // permute lanes and zero some 64-bit elements
            return  _mm512_maskz_shuffle_i64x2(z, a, a, (_MM_PERM_ENUM)s);
        }
        else {
            // permute lanes
            return _mm512_shuffle_i64x2(a, a, (_MM_PERM_ENUM)s);
        }
    }
    // full permute needed
    const __m512i pmask = constant16i<i0&7, 0, i1&7, 0, i2&7, 0, i3&7, 0, i4&7, 0, i5&7, 0, i6&7, 0, i7&7, 0>();
    if (dozero) {
        // full permute and zeroing
        return _mm512_maskz_permutexvar_epi64(z, pmask, a);
    }
    else {
        return _mm512_permutexvar_epi64(pmask, a);
    }
}

template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7>
static inline Vec8uq permute8(Vec8uq const & a) {
    return Vec8uq (permute8<i0,i1,i2,i3,i4,i5,i6,i7> (Vec8q(a)));
}


// Permute vector of 16 32-bit integers.
// Index -1 gives 0, index -256 means don't care.
template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7, int i8, int i9, int i10, int i11, int i12, int i13, int i14, int i15>
static inline Vec16i permute16(Vec16i const & a) {

    // Combine indexes into a single bitfield, with 4 bits for each
    const uint64_t m1 = (i0&15) | (i1&15)<<4 | (i2&15)<< 8 | (i3&15)<<12 | (i4&15)<<16 | (i5&15)<<20 | (i6&15)<<24 | (i7&15LL)<<28   // 15LL avoids sign extension of (int32_t | int64_t)
        | (i8&15LL)<<32 | (i9&15LL)<<36 | (i10&15LL)<<40 | (i11&15LL)<<44 | (i12&15LL)<<48 | (i13&15LL)<<52 | (i14&15LL)<<56 | (i15&15LL)<<60;

    // Mask to zero out negative indexes
    const uint64_t mz = (i0<0?0:0xF) | (i1<0?0:0xF0) | (i2<0?0:0xF00) | (i3<0?0:0xF000) | (i4<0?0:0xF0000) | (i5<0?0:0xF00000) | (i6<0?0:0xF000000) | (i7<0?0:0xF0000000ULL) | (i8<0?0:0xF00000000) 
        | (i9<0?0:0xF000000000) | (i10<0?0:0xF0000000000) | (i11<0?0:0xF00000000000) | (i12<0?0:0xF000000000000) | (i13<0?0:0xF0000000000000) | (i14<0?0:0xF00000000000000) | (i15<0?0:0xF000000000000000);

    const uint64_t m2 = m1 & mz;

    // zeroing needed
    const bool dozero = ((i0|i1|i2|i3|i4|i5|i6|i7|i8|i9|i10|i11|i12|i13|i14|i15) & 0x80) != 0;

    // special case: all zero
    if (mz == 0) return  _mm512_setzero_epi32();

    // mask for elements not zeroed
    const __mmask16 z = __mmask16((i0>=0)<<0 | (i1>=0)<<1 | (i2>=0)<<2 | (i3>=0)<<3 | (i4>=0)<<4 | (i5>=0)<<5 | (i6>=0)<<6 | (i7>=0)<<7
        | (i8>=0)<<8 | (i9>=0)<<9 | (i10>=0)<<10 | (i11>=0)<<11 | (i12>=0)<<12 | (i13>=0)<<13 | (i14>=0)<<14 | (i15>=0)<<15);

    if (((m1 ^ 0xFEDCBA9876543210) & mz) == 0) {
        // no shuffling
        if (dozero) {
            return _mm512_maskz_mov_epi32(z, a);
        }
        return a;                                 // do nothing
    }

    if (((m1 ^ 0xCCCC888844440000) & 0xCCCCCCCCCCCCCCCC & mz) == 0) {
        // no exchange of data between the four 128-bit lanes
        const uint64_t pat = ((m2 | (m2 >> 16) | (m2 >> 32) | (m2 >> 48)) & 0x3333) * 0x0001000100010001;
        const int pmask = (pat & 3) | (((pat >> 4) & 3) << 2) | (((pat >> 8) & 3) << 4) | (((pat >> 12) & 3) << 6);
        if (((m1 ^ pat) & 0x3333333333333333 & mz) == 0) {
            // same permute pattern in all lanes
            if (dozero) {  // permute within lanes and zero
                return _mm512_maskz_shuffle_epi32(z, a, (_MM_PERM_ENUM)pmask);
            }
            else {  // permute within lanes
                return _mm512_shuffle_epi32(a, (_MM_PERM_ENUM)pmask);
            }
        }
        // different permute patterns in each lane. It's faster to do a full permute than four masked permutes within lanes
    }
    const uint64_t lane = (m2 | m2 >> 4 | m2 >> 8 | m2 >> 12) & 0x000C000C000C000C;
    if ((((m1 ^ 0x3210321032103210) & 0x3333333333333333 & mz) == 0) 
    &&  ((m1 ^ (lane * 0x1111)) & 0xCCCCCCCCCCCCCCCC & mz) == 0) {
        // permute lanes only. no permutation within each lane
        const uint64_t s = ((lane >> 2) & 3) | (((lane >> 18) & 3) << 2) | (((lane >> 34) & 3) << 4) | (((lane >> 50) & 3) << 6);
        if (dozero) {
            // permute lanes and zero some 64-bit elements
            return  _mm512_maskz_shuffle_i32x4(z, a, a, (_MM_PERM_ENUM)s);
        }
        else {
            // permute lanes
            return _mm512_shuffle_i32x4(a, a, (_MM_PERM_ENUM)s);
        }
    }
    // full permute needed
    const __m512i pmask = constant16i<i0&15, i1&15, i2&15, i3&15, i4&15, i5&15, i6&15, i7&15, i8&15, i9&15, i10&15, i11&15, i12&15, i13&15, i14&15, i15&15>();
    if (dozero) {
        // full permute and zeroing
        return _mm512_maskz_permutexvar_epi32(z, pmask, a);
    }
    else {    
        return _mm512_permutexvar_epi32(pmask, a);
    }
}


template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7, int i8, int i9, int i10, int i11, int i12, int i13, int i14, int i15>
static inline Vec16ui permute4(Vec16ui const & a) {
    return Vec16ui (permute16<i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15> (Vec16i(a)));
}

// Support names prior to version 1.40 (now deprecated)
#define permute8q   permute8
#define permute8uq  permute8
#define permute16i  permute16
#define permute16ui permute16


/*****************************************************************************
*
*          Vector blend functions
*
******************************************************************************
*
* These blend functions can mix elements from two different vectors and
* optionally set some elements to zero. 
*
* The indexes are inserted as template parameters in <>. These indexes must be
* constants. Each template parameter is an index to the element you want to 
* select, where higher indexes indicate an element from the second source
* vector. For example, if each vector has 8 elements, then indexes 0 - 7
* will select an element from the first vector and indexes 8 - 15 will select 
* an element from the second vector. A negative index will generate zero.
*
* Example:
* Vec8q a(100,101,102,103,104,105,106,107); // a is (100, 101, 102, 103, 104, 105, 106, 107)
* Vec8q b(200,201,202,203,204,205,206,207); // b is (200, 201, 202, 203, 204, 205, 206, 207)
* Vec8q c;
* c = blend8<1,0,9,8,7,-1,15,15> (a,b);    // c is (101, 100, 201, 200, 107,   0, 207, 207)
*
* A lot of the code here is metaprogramming aiming to find the instructions
* that best fit the template parameters and instruction set. The metacode
* will be reduced out to leave only a few vector instructions in release
* mode with optimization on.
*****************************************************************************/


template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7> 
static inline Vec8q blend8(Vec8q const & a, Vec8q const & b) {  

    // Combine indexes into a single bitfield, with 4 bits for each
    const int m1 = (i0&0xF) | (i1&0xF)<<4 | (i2&0xF)<< 8 | (i3&0xF)<<12 | (i4&0xF)<<16 | (i5&0xF)<<20 | (i6&0xF)<<24 | (i7&0xF)<<28;

    // Mask to zero out negative indexes
    const int mz = (i0<0?0:0xF) | (i1<0?0:0xF0) | (i2<0?0:0xF00) | (i3<0?0:0xF000) | (i4<0?0:0xF0000) | (i5<0?0:0xF00000) | (i6<0?0:0xF000000) | (i7<0?0:0xF0000000);
    //const int m2 = m1 & mz;

    // zeroing needed
    const bool dozero = ((i0|i1|i2|i3|i4|i5|i6|i7) & 0x80) != 0;

    // mask for elements not zeroed
    const __mmask16 z = __mmask16((i0>=0)<<0 | (i1>=0)<<1 | (i2>=0)<<2 | (i3>=0)<<3 | (i4>=0)<<4 | (i5>=0)<<5 | (i6>=0)<<6 | (i7>=0)<<7);

    // special case: all zero
    if (mz == 0) return  _mm512_setzero_epi32();

    // special case: all from a
    if ((m1 & 0x88888888 & mz) == 0) {
        return permute8 <i0, i1, i2, i3, i4, i5, i6, i7> (a);
    }

    // special case: all from b
    if ((~m1 & 0x88888888 & mz) == 0) {
        return permute8 <i0^8, i1^8, i2^8, i3^8, i4^8, i5^8, i6^8, i7^8> (b);
    }

    // special case: blend without permute
    if (((m1 ^ 0x76543210) & 0x77777777 & mz) == 0) {
        __mmask16 blendmask = __mmask16((i0&8)>>3 | (i1&8)>>2 | (i2&8)>>1 | (i3&8)>>0 | (i4&8)<<1 | (i5&8)<<2 | (i6&8)<<3 | (i7&8)<<4 );
        __m512i t = _mm512_mask_blend_epi64(blendmask, a, b);
        if (dozero) {
            t = _mm512_maskz_mov_epi64(z, t);
        }
        return t;
    }
    // special case: all data stay within their lane
    if (((m1 ^ 0x66442200) & 0x66666666 & mz) == 0) {
        // mask for elements from a and b
        const uint32_t mb = ((i0&8)?0xF:0) | ((i1&8)?0xF0:0) | ((i2&8)?0xF00:0) | ((i3&8)?0xF000:0) | ((i4&8)?0xF0000:0) | ((i5&8)?0xF00000:0) | ((i6&8)?0xF000000:0) | ((i7&8)?0xF0000000:0);
        const uint32_t mbz = mb & mz;     // mask for nonzero elements from b
        const uint32_t maz = ~mb & mz;    // mask for nonzero elements from a
        const uint32_t m1a = m1 & maz;
        const uint32_t m1b = m1 & mbz;
        const uint32_t pata = ((m1a | m1a >> 8 | m1a >> 16 | m1a >> 24) & 0xFF) * 0x01010101;  // permute pattern for elements from a
        const uint32_t patb = ((m1b | m1b >> 8 | m1b >> 16 | m1b >> 24) & 0xFF) * 0x01010101;  // permute pattern for elements from b
        if (((m1 ^ pata) & 0x11111111 & maz) == 0 && ((m1 ^ patb) & 0x11111111 & mbz) == 0) {
            // Same permute pattern in all lanes:
            // This code generates two instructions instead of one, but we are avoiding the slow lane-crossing instruction,
            // and we are saving 64 bytes of data cache.
            // 1. Permute a, zero elements not from a (using _mm512_maskz_shuffle_epi32)
            __m512i ta = permute8< (maz&0xF)?i0&7:-1, (maz&0xF0)?i1&7:-1, (maz&0xF00)?i2&7:-1, (maz&0xF000)?i3&7:-1, 
                (maz&0xF0000)?i4&7:-1, (maz&0xF00000)?i5&7:-1, (maz&0xF000000)?i6&7:-1, (maz&0xF0000000)?i7&7:-1> (a);
            // write mask for elements from b
            const __mmask16 sb = ((mbz&0xF)?3:0) | ((mbz&0xF0)?0xC:0) | ((mbz&0xF00)?0x30:0) | ((mbz&0xF000)?0xC0:0) | ((mbz&0xF0000)?0x300:0) | ((mbz&0xF00000)?0xC00:0) | ((mbz&0xF000000)?0x3000:0) | ((mbz&0xF0000000)?0xC000:0);
            // permute index for elements from b
            const int pi = ((patb & 1) * 10 + 4) | ((((patb >> 4) & 1) * 10 + 4) << 4);
            // 2. Permute elements from b and combine with elements from a through write mask
            return _mm512_mask_shuffle_epi32(ta, sb, b, (_MM_PERM_ENUM)pi);
        }
        // not same permute pattern in all lanes. use full permute
    }
    // general case: full permute
    const __m512i pmask = constant16i<i0&0xF, 0, i1&0xF, 0, i2&0xF, 0, i3&0xF, 0, i4&0xF, 0, i5&0xF, 0, i6&0xF, 0, i7&0xF, 0>();
    if (dozero) {
        return _mm512_maskz_permutex2var_epi64(z, a, pmask, b);
    }
    else {
        return _mm512_permutex2var_epi64(a, pmask, b);
    }
}

template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7> 
static inline Vec8uq blend8(Vec8uq const & a, Vec8uq const & b) {
    return Vec8uq( blend8<i0,i1,i2,i3,i4,i5,i6,i7> (Vec8q(a),Vec8q(b)));
}


template <int i0,  int i1,  int i2,  int i3,  int i4,  int i5,  int i6,  int i7, 
          int i8,  int i9,  int i10, int i11, int i12, int i13, int i14, int i15 > 
static inline Vec16i blend16(Vec16i const & a, Vec16i const & b) {  

    // Combine indexes into a single bitfield, with 4 bits for each indicating shuffle, but not source
    const uint64_t m1 = (i0&0xF) | (i1&0xF)<<4 | (i2&0xF)<<8 | (i3&0xF)<<12 | (i4&0xF)<<16 | (i5&0xF)<<20 | (i6&0xF)<<24 | (i7&0xFLL)<<28
        | (i8&0xFLL)<<32 | (i9&0xFLL)<<36 | (i10&0xFLL)<<40 | (i11&0xFLL)<<44 | (i12&0xFLL)<<48 | (i13&0xFLL)<<52 | (i14&0xFLL)<<56 | (i15&0xFLL)<<60;

    // Mask to zero out negative indexes
    const uint64_t mz = (i0<0?0:0xF) | (i1<0?0:0xF0) | (i2<0?0:0xF00) | (i3<0?0:0xF000) | (i4<0?0:0xF0000) | (i5<0?0:0xF00000) | (i6<0?0:0xF000000) | (i7<0?0:0xF0000000ULL)
        | (i8<0?0:0xF00000000) | (i9<0?0:0xF000000000) | (i10<0?0:0xF0000000000) | (i11<0?0:0xF00000000000) | (i12<0?0:0xF000000000000) | (i13<0?0:0xF0000000000000) | (i14<0?0:0xF00000000000000) | (i15<0?0:0xF000000000000000);
    //const uint64_t m2 = m1 & mz;

    // collect bit 4 of each index = select source
    const uint64_t ms = ((i0&16)?0xF:0) | ((i1&16)?0xF0:0) | ((i2&16)?0xF00:0) | ((i3&16)?0xF000:0) | ((i4&16)?0xF0000:0) | ((i5&16)?0xF00000:0) | ((i6&16)?0xF000000:0) | ((i7&16)?0xF0000000ULL:0)
        | ((i8&16)?0xF00000000:0) | ((i9&16)?0xF000000000:0) | ((i10&16)?0xF0000000000:0) | ((i11&16)?0xF00000000000:0) | ((i12&16)?0xF000000000000:0) | ((i13&16)?0xF0000000000000:0) | ((i14&16)?0xF00000000000000:0) | ((i15&16)?0xF000000000000000:0);

    // zeroing needed
    const bool dozero = ((i0|i1|i2|i3|i4|i5|i6|i7|i8|i9|i10|i11|i12|i13|i14|i15) & 0x80) != 0;

    // mask for elements not zeroed
    const __mmask16 z = __mmask16((i0>=0)<<0 | (i1>=0)<<1 | (i2>=0)<<2 | (i3>=0)<<3 | (i4>=0)<<4 | (i5>=0)<<5 | (i6>=0)<<6 | (i7>=0)<<7 
        | (i8>=0)<<8 | (i9>=0)<<9 | (i10>=0)<<10 | (i11>=0)<<11 | (i12>=0)<<12 | (i13>=0)<<13 | (i14>=0)<<14 | (i15>=0)<<15);

    // special case: all zero
    if (mz == 0) return  _mm512_setzero_epi32();

    // special case: all from a
    if ((ms & mz) == 0) {
        return permute16<i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15> (a);
    }

    // special case: all from b
    if ((~ms & mz) == 0) {
        return permute16<i0^16,i1^16,i2^16,i3^16,i4^16,i5^16,i6^16,i7^16,i8^16,i9^16,i10^16,i11^16,i12^16,i13^16,i14^16,i15^16 > (b);
    }

    // special case: blend without permute
    if (((m1 ^ 0xFEDCBA9876543210) & mz) == 0) {
        __mmask16 blendmask = __mmask16((i0&16)>>4 | (i1&16)>>3 | (i2&16)>>2 | (i3&16)>>1 | (i4&16) | (i5&16)<<1 | (i6&16)<<2 | (i7&16)<<3
            | (i8&16)<<4 | (i9&16)<<5 | (i10&16)<<6 | (i11&16)<<7 | (i12&16)<<8 | (i13&16)<<9 | (i14&16)<<10 | (i15&16)<<11);
        __m512i t = _mm512_mask_blend_epi32(blendmask, a, b);
        if (dozero) {
            t = _mm512_maskz_mov_epi32(z, t);
        }
        return t;
    }

    // special case: all data stay within their lane
    if (((m1 ^ 0xCCCC888844440000) & 0xCCCCCCCCCCCCCCCC & mz) == 0) {

        // mask for elements from a and b
        const uint64_t mb  = ms;
        const uint64_t mbz = mb & mz;     // mask for nonzero elements from b
        const uint64_t maz = ~mb & mz;    // mask for nonzero elements from a
        const uint64_t m1a = m1 & maz;
        const uint64_t m1b = m1 & mbz;
        const uint64_t pata = ((m1a | m1a >> 16 | m1a >> 32 | m1a >> 48) & 0xFFFF) * 0x0001000100010001;  // permute pattern for elements from a
        const uint64_t patb = ((m1b | m1b >> 16 | m1b >> 32 | m1b >> 48) & 0xFFFF) * 0x0001000100010001;  // permute pattern for elements from b
        if (((m1 ^ pata) & 0x3333333333333333 & maz) == 0 && ((m1 ^ patb) & 0x3333333333333333 & mbz) == 0) {
            // Same permute pattern in all lanes:
            // This code generates two instructions instead of one, but we are avoiding the slow lane-crossing instruction,
            // and we are saving 64 bytes of data cache.
            // 1. Permute a, zero elements not from a (using _mm512_maskz_shuffle_epi32)
            __m512i ta = permute16< (maz&0xF)?i0&15:-1, (maz&0xF0)?i1&15:-1, (maz&0xF00)?i2&15:-1, (maz&0xF000)?i3&15:-1, 
                (maz&0xF0000)?i4&15:-1, (maz&0xF00000)?i5&15:-1, (maz&0xF000000)?i6&15:-1, (maz&0xF0000000)?i7&15:-1,
                (maz&0xF00000000)?i8&15:-1, (maz&0xF000000000)?i9&15:-1, (maz&0xF0000000000)?i10&15:-1, (maz&0xF00000000000)?i11&15:-1, 
                (maz&0xF000000000000)?i12&15:-1, (maz&0xF0000000000000)?i13&15:-1, (maz&0xF00000000000000)?i14&15:-1, (maz&0xF000000000000000)?i15&15:-1> (a);
            // write mask for elements from b
            const __mmask16 sb = ((mbz&0xF)?1:0) | ((mbz&0xF0)?0x2:0) | ((mbz&0xF00)?0x4:0) | ((mbz&0xF000)?0x8:0) | ((mbz&0xF0000)?0x10:0) | ((mbz&0xF00000)?0x20:0) | ((mbz&0xF000000)?0x40:0) | ((mbz&0xF0000000)?0x80:0) 
                | ((mbz&0xF00000000)?0x100:0) | ((mbz&0xF000000000)?0x200:0) | ((mbz&0xF0000000000)?0x400:0) | ((mbz&0xF00000000000)?0x800:0) | ((mbz&0xF000000000000)?0x1000:0) | ((mbz&0xF0000000000000)?0x2000:0) | ((mbz&0xF00000000000000)?0x4000:0) | ((mbz&0xF000000000000000)?0x8000:0);
            // permute index for elements from b
            const int pi = (patb & 3) | (((patb >> 4) & 3) << 2) | (((patb >> 8) & 3) << 4) | (((patb >> 12) & 3) << 6);
            // 2. Permute elements from b and combine with elements from a through write mask
            return _mm512_mask_shuffle_epi32(ta, sb, b, (_MM_PERM_ENUM)pi);
        }
        // not same permute pattern in all lanes. use full permute
    }

    // general case: full permute
    const __m512i pmask = constant16i<i0&0x1F, i1&0x1F, i2&0x1F, i3&0x1F, i4&0x1F, i5&0x1F, i6&0x1F, i7&0x1F, 
        i8&0x1F, i9&0x1F, i10&0x1F, i11&0x1F, i12&0x1F, i13&0x1F, i14&0x1F, i15&0x1F>();
    if (dozero) {
        return _mm512_maskz_permutex2var_epi32(z, a, pmask, b);        
    }
    else {
        return _mm512_permutex2var_epi32(a, pmask, b);
    }
}

template <int i0,  int i1,  int i2,  int i3,  int i4,  int i5,  int i6,  int i7, 
          int i8,  int i9,  int i10, int i11, int i12, int i13, int i14, int i15 > 
static inline Vec16ui blend16(Vec16ui const & a, Vec16ui const & b) {
    return Vec16ui( blend16<i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15> (Vec16i(a),Vec16i(b)));
}

// Support names prior to version 1.40 (now deprecated)
#define blend8q   blend8
#define blend8uq  blend8
#define blend16i  blend16
#define blend16ui blend16


/*****************************************************************************
*
*          Vector lookup functions
*
******************************************************************************
*
* These functions use vector elements as indexes into a table.
* The table is given as one or more vectors or as an array.
*
* This can be used for several purposes:
*  - table lookup
*  - permute or blend with variable indexes
*  - blend from more than two sources
*  - gather non-contiguous data
*
* An index out of range may produce any value - the actual value produced is
* implementation dependent and may be different for different instruction
* sets. An index out of range does not produce an error message or exception.
*
* Example:
* Vec8q a(2,0,0,6,4,3,5,0);                 // index a is (  2,   0,   0,   6,   4,   3,   5,   0)
* Vec8q b(100,101,102,103,104,105,106,107); // table b is (100, 101, 102, 103, 104, 105, 106, 107)
* Vec8q c = lookup8 (a,b);                  // c is       (102, 100, 100, 106, 104, 103, 105, 100)
*
*****************************************************************************/

static inline Vec16i lookup16(Vec16i const & index, Vec16i const & table) {
    return _mm512_permutexvar_epi32(index, table);
}

static inline Vec16i lookup32(Vec16i const & index, Vec16i const & table1, Vec16i const & table2) {
    return _mm512_permutex2var_epi32(table1, index, table2);
}

static inline Vec16i lookup64(Vec16i const & index, Vec16i const & table1, Vec16i const & table2, Vec16i const & table3, Vec16i const & table4) {
    Vec16i d12 = _mm512_permutex2var_epi32(table1, index, table2);
    Vec16i d34 = _mm512_permutex2var_epi32(table3, index, table4);
    return select((index >> 5) != 0, d34, d12);
}

template <int n>
static inline Vec16i lookup(Vec16i const & index, void const * table) {
    if (n <= 0) return 0;
    if (n <= 16) {
        Vec16i table1 = Vec16i().load(table);
        return lookup16(index, table1);
    }
    if (n <= 32) {
        Vec16i table1 = Vec16i().load(table);
        Vec16i table2 = Vec16i().load((int8_t*)table + 64);
        return _mm512_permutex2var_epi32(table1, index, table2);
    }
    // n > 32. Limit index
    Vec16ui index1;
    if ((n & (n-1)) == 0) {
        // n is a power of 2, make index modulo n
        index1 = Vec16ui(index) & (n-1);
    }
    else {
        // n is not a power of 2, limit to n-1
        index1 = min(Vec16ui(index), uint32_t(n-1));
    }
    return _mm512_i32gather_epi32(index1, (const int*)table, 4);
    // return  _mm512_i32gather_epi32(index1, table, _MM_UPCONV_EPI32_NONE, 4, 0);
}


static inline Vec8q lookup8(Vec8q const & index, Vec8q const & table) {
    return _mm512_permutexvar_epi64(index, table);
}

template <int n>
static inline Vec8q lookup(Vec8q const & index, void const * table) {
    if (n <= 0) return 0;
    if (n <= 8) {
        Vec8q table1 = Vec8q().load(table);
        return lookup8(index, table1);
    }
    if (n <= 16) {
        Vec8q table1 = Vec8q().load(table);
        Vec8q table2 = Vec8q().load((int8_t*)table + 64);
        return _mm512_permutex2var_epi64(table1, index, table2);
    }
    // n > 16. Limit index
    Vec8uq index1;
    if ((n & (n-1)) == 0) {
        // n is a power of 2, make index modulo n
        index1 = Vec8uq(index) & (n-1);
    }
    else {
        // n is not a power of 2, limit to n-1
        index1 = min(Vec8uq(index), uint32_t(n-1));
    }
    return _mm512_i64gather_epi64(index1, (const long long*)table, 8);
}


/*****************************************************************************
*
*          Gather functions with fixed indexes
*
*****************************************************************************/
// Load elements from array a with indices i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15
template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7, 
int i8, int i9, int i10, int i11, int i12, int i13, int i14, int i15>
static inline Vec16i gather16i(void const * a) {
    Static_error_check<(i0|i1|i2|i3|i4|i5|i6|i7|i8|i9|i10|i11|i12|i13|i14|i15)>=0> Negative_array_index;  // Error message if index is negative
    // find smallest and biggest index, using only compile-time constant expressions
    const int i01min   = i0  < i1  ? i0  : i1;
    const int i23min   = i2  < i3  ? i2  : i3;
    const int i45min   = i4  < i5  ? i4  : i5;
    const int i67min   = i6  < i7  ? i6  : i7;
    const int i89min   = i8  < i9  ? i8  : i9;
    const int i1011min = i10 < i11 ? i10 : i11;
    const int i1213min = i12 < i13 ? i12 : i13;
    const int i1415min = i14 < i15 ? i14 : i15;
    const int i0_3min   = i01min   < i23min    ? i01min   : i23min;
    const int i4_7min   = i45min   < i67min    ? i45min   : i67min;
    const int i8_11min  = i89min   < i1011min  ? i89min   : i1011min;
    const int i12_15min = i1213min < i1415min  ? i1213min : i1415min;
    const int i0_7min   = i0_3min  < i4_7min   ? i0_3min  : i4_7min;
    const int i8_15min  = i8_11min < i12_15min ? i8_11min : i12_15min;
    const int imin      = i0_7min  < i8_15min  ? i0_7min  : i8_15min;
    const int i01max   = i0  > i1  ? i0  : i1;
    const int i23max   = i2  > i3  ? i2  : i3;
    const int i45max   = i4  > i5  ? i4  : i5;
    const int i67max   = i6  > i7  ? i6  : i7;
    const int i89max   = i8  > i9  ? i8  : i9;
    const int i1011max = i10 > i11 ? i10 : i11;
    const int i1213max = i12 > i13 ? i12 : i13;
    const int i1415max = i14 > i15 ? i14 : i15;
    const int i0_3max   = i01max   > i23max    ? i01max   : i23max;
    const int i4_7max   = i45max   > i67max    ? i45max   : i67max;
    const int i8_11max  = i89max   > i1011max  ? i89max   : i1011max;
    const int i12_15max = i1213max > i1415max  ? i1213max : i1415max;
    const int i0_7max   = i0_3max  > i4_7max   ? i0_3max  : i4_7max;
    const int i8_15max  = i8_11max > i12_15max ? i8_11max : i12_15max;
    const int imax      = i0_7max  > i8_15max  ? i0_7max  : i8_15max;
    if (imax - imin <= 15) {
        // load one contiguous block and permute
        if (imax > 15) {
            // make sure we don't read past the end of the array
            Vec16i b = Vec16i().load((int32_t const *)a + imax-15);
            return permute16<i0-imax+15, i1-imax+15, i2-imax+15, i3-imax+15, i4-imax+15, i5-imax+15, i6-imax+15, i7-imax+15,
                i8-imax+15, i9-imax+15, i10-imax+15, i11-imax+15, i12-imax+15, i13-imax+15, i14-imax+15, i15-imax+15> (b);
        }
        else {
            Vec16i b = Vec16i().load((int32_t const *)a + imin);
            return permute16<i0-imin, i1-imin, i2-imin, i3-imin, i4-imin, i5-imin, i6-imin, i7-imin,
                i8-imin, i9-imin, i10-imin, i11-imin, i12-imin, i13-imin, i14-imin, i15-imin> (b);
        }
    }
    if ((i0<imin+16  || i0>imax-16)  && (i1<imin+16  || i1>imax-16)  && (i2<imin+16  || i2>imax-16)  && (i3<imin+16  || i3>imax-16)
    &&  (i4<imin+16  || i4>imax-16)  && (i5<imin+16  || i5>imax-16)  && (i6<imin+16  || i6>imax-16)  && (i7<imin+16  || i7>imax-16)    
    &&  (i8<imin+16  || i8>imax-16)  && (i9<imin+16  || i9>imax-16)  && (i10<imin+16 || i10>imax-16) && (i11<imin+16 || i11>imax-16)
    &&  (i12<imin+16 || i12>imax-16) && (i13<imin+16 || i13>imax-16) && (i14<imin+16 || i14>imax-16) && (i15<imin+16 || i15>imax-16) ) {
        // load two contiguous blocks and blend
        Vec16i b = Vec16i().load((int32_t const *)a + imin);
        Vec16i c = Vec16i().load((int32_t const *)a + imax-15);
        const int j0  = i0 <imin+16 ? i0 -imin : 31-imax+i0;
        const int j1  = i1 <imin+16 ? i1 -imin : 31-imax+i1;
        const int j2  = i2 <imin+16 ? i2 -imin : 31-imax+i2;
        const int j3  = i3 <imin+16 ? i3 -imin : 31-imax+i3;
        const int j4  = i4 <imin+16 ? i4 -imin : 31-imax+i4;
        const int j5  = i5 <imin+16 ? i5 -imin : 31-imax+i5;
        const int j6  = i6 <imin+16 ? i6 -imin : 31-imax+i6;
        const int j7  = i7 <imin+16 ? i7 -imin : 31-imax+i7;
        const int j8  = i8 <imin+16 ? i8 -imin : 31-imax+i8;
        const int j9  = i9 <imin+16 ? i9 -imin : 31-imax+i9;
        const int j10 = i10<imin+16 ? i10-imin : 31-imax+i10;
        const int j11 = i11<imin+16 ? i11-imin : 31-imax+i11;
        const int j12 = i12<imin+16 ? i12-imin : 31-imax+i12;
        const int j13 = i13<imin+16 ? i13-imin : 31-imax+i13;
        const int j14 = i14<imin+16 ? i14-imin : 31-imax+i14;
        const int j15 = i15<imin+16 ? i15-imin : 31-imax+i15;
        return blend16<j0,j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14,j15>(b, c);
    }
    // use gather instruction
    return _mm512_i32gather_epi32(Vec16i(i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15), (const int *)a, 4);
}


template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7>
static inline Vec8q gather8q(void const * a) {
    Static_error_check<(i0|i1|i2|i3|i4|i5|i6|i7)>=0> Negative_array_index;  // Error message if index is negative

    const int i01min = i0 < i1 ? i0 : i1;
    const int i23min = i2 < i3 ? i2 : i3;
    const int i45min = i4 < i5 ? i4 : i5;
    const int i67min = i6 < i7 ? i6 : i7;
    const int i0123min = i01min < i23min ? i01min : i23min;
    const int i4567min = i45min < i67min ? i45min : i67min;
    const int imin = i0123min < i4567min ? i0123min : i4567min;
    const int i01max = i0 > i1 ? i0 : i1;
    const int i23max = i2 > i3 ? i2 : i3;
    const int i45max = i4 > i5 ? i4 : i5;
    const int i67max = i6 > i7 ? i6 : i7;
    const int i0123max = i01max > i23max ? i01max : i23max;
    const int i4567max = i45max > i67max ? i45max : i67max;
    const int imax = i0123max > i4567max ? i0123max : i4567max;
    if (imax - imin <= 7) {
        // load one contiguous block and permute
        if (imax > 7) {
            // make sure we don't read past the end of the array
            Vec8q b = Vec8q().load((int64_t const *)a + imax-7);
            return permute8<i0-imax+7, i1-imax+7, i2-imax+7, i3-imax+7, i4-imax+7, i5-imax+7, i6-imax+7, i7-imax+7> (b);
        }
        else {
            Vec8q b = Vec8q().load((int64_t const *)a + imin);
            return permute8<i0-imin, i1-imin, i2-imin, i3-imin, i4-imin, i5-imin, i6-imin, i7-imin> (b);
        }
    }
    if ((i0<imin+8 || i0>imax-8) && (i1<imin+8 || i1>imax-8) && (i2<imin+8 || i2>imax-8) && (i3<imin+8 || i3>imax-8)
    &&  (i4<imin+8 || i4>imax-8) && (i5<imin+8 || i5>imax-8) && (i6<imin+8 || i6>imax-8) && (i7<imin+8 || i7>imax-8)) {
        // load two contiguous blocks and blend
        Vec8q b = Vec8q().load((int64_t const *)a + imin);
        Vec8q c = Vec8q().load((int64_t const *)a + imax-7);
        const int j0 = i0<imin+8 ? i0-imin : 15-imax+i0;
        const int j1 = i1<imin+8 ? i1-imin : 15-imax+i1;
        const int j2 = i2<imin+8 ? i2-imin : 15-imax+i2;
        const int j3 = i3<imin+8 ? i3-imin : 15-imax+i3;
        const int j4 = i4<imin+8 ? i4-imin : 15-imax+i4;
        const int j5 = i5<imin+8 ? i5-imin : 15-imax+i5;
        const int j6 = i6<imin+8 ? i6-imin : 15-imax+i6;
        const int j7 = i7<imin+8 ? i7-imin : 15-imax+i7;
        return blend8<j0, j1, j2, j3, j4, j5, j6, j7>(b, c);
    }
    // use gather instruction
    return _mm512_i64gather_epi64(Vec8q(i0,i1,i2,i3,i4,i5,i6,i7), (const long long *)a, 8);
}

/*****************************************************************************
*
*          Vector scatter functions
*
******************************************************************************
*
* These functions write the elements of a vector to arbitrary positions in an
* array in memory. Each vector element is written to an array position 
* determined by an index. An element is not written if the corresponding
* index is out of range.
* The indexes can be specified as constant template parameters or as an
* integer vector.
* 
* The scatter functions are useful if the data are distributed in a sparce
* manner into the array. If the array is dense then it is more efficient
* to permute the data into the right positions and then write the whole
* permuted vector into the array.
*
* Example:
* Vec8q a(10,11,12,13,14,15,16,17);
* int64_t b[16] = {0};
* scatter<0,2,14,10,1,-1,5,9>(a,b); 
* // Now, b = {10,14,11,0,0,16,0,0,0,17,13,0,0,0,12,0}
*
*****************************************************************************/

template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7,
int i8, int i9, int i10, int i11, int i12, int i13, int i14, int i15>
    static inline void scatter(Vec16i const & data, void * array) {
    __m512i indx = constant16i<i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15>();
    Vec16ib mask(i0>=0, i1>=0, i2>=0, i3>=0, i4>=0, i5>=0, i6>=0, i7>=0,
        i8>=0, i9>=0, i10>=0, i11>=0, i12>=0, i13>=0, i14>=0, i15>=0);
    _mm512_mask_i32scatter_epi32((int*)array, mask, indx, data, 4);
}

template <int i0, int i1, int i2, int i3, int i4, int i5, int i6, int i7>
static inline void scatter(Vec8q const & data, void * array) {
    __m256i indx = constant8i<i0,i1,i2,i3,i4,i5,i6,i7>();
    Vec8qb mask(i0>=0, i1>=0, i2>=0, i3>=0, i4>=0, i5>=0, i6>=0, i7>=0);
    _mm512_mask_i32scatter_epi64((long long *)array, mask, indx, data, 8);
}


/*****************************************************************************
*
*          Scatter functions with variable indexes
*
*****************************************************************************/

static inline void scatter(Vec16i const & index, uint32_t limit, Vec16i const & data, void * destination) {
    Vec16ib mask = Vec16ui(index) < limit;
    _mm512_mask_i32scatter_epi32((int*)destination, mask, index, data, 4);
}

static inline void scatter(Vec8q const & index, uint32_t limit, Vec8q const & data, void * destination) {
    Vec8qb mask = Vec8uq(index) < uint64_t(limit);
    _mm512_mask_i64scatter_epi64((long long *)destination, mask, index, data, 8);
}

static inline void scatter(Vec8i const & index, uint32_t limit, Vec8q const & data, void * destination) {
#if INSTRSET >= 10 //  __AVX512VL__
    __mmask16 mask = _mm256_cmplt_epu32_mask(index, Vec8ui(limit));
#else
    __mmask16 mask = _mm512_mask_cmplt_epu32_mask(0xFFu, _mm512_castsi256_si512(index), _mm512_castsi256_si512(Vec8ui(limit)));
#endif
    _mm512_mask_i32scatter_epi64((long long *)destination, (__mmask8)mask, index, data, 8);
}


/*****************************************************************************
*
*          Functions for conversion between integer sizes
*
*****************************************************************************/

// Extend 16-bit integers to 32-bit integers, signed and unsigned

// Function extend_to_int : extends Vec16s to Vec16i with sign extension
static inline Vec16i extend_to_int (Vec16s const & a) {
    return _mm512_cvtepi16_epi32(a);
}

// Function extend_to_int : extends Vec16us to Vec16ui with zero extension
static inline Vec16ui extend_to_int (Vec16us const & a) {
    return _mm512_cvtepu16_epi32(a);
}

// Function extend_to_int : extends Vec16c to Vec16i with sign extension
static inline Vec16i extend_to_int (Vec16c const & a) {
    return _mm512_cvtepi8_epi32(a);
}

// Function extend_to_int : extends Vec16uc to Vec16ui with zero extension
static inline Vec16ui extend_to_int (Vec16uc const & a) {
    return _mm512_cvtepu8_epi32(a);
}


// Extend 32-bit integers to 64-bit integers, signed and unsigned

// Function extend_low : extends the low 8 elements to 64 bits with sign extension
static inline Vec8q extend_low (Vec16i const & a) {
    return _mm512_cvtepi32_epi64(a.get_low());
}

// Function extend_high : extends the high 8 elements to 64 bits with sign extension
static inline Vec8q extend_high (Vec16i const & a) {
    return _mm512_cvtepi32_epi64(a.get_high());
}

// Function extend_low : extends the low 8 elements to 64 bits with zero extension
static inline Vec8uq extend_low (Vec16ui const & a) {
    return _mm512_cvtepu32_epi64(a.get_low());
}

// Function extend_high : extends the high 8 elements to 64 bits with zero extension
static inline Vec8uq extend_high (Vec16ui const & a) {
    return _mm512_cvtepu32_epi64(a.get_high());
}


// Compress 32-bit integers to 8-bit integers, signed and unsigned, with and without saturation

// Function compress : packs two vectors of 16-bit integers into one vector of 8-bit integers
// Overflow wraps around
static inline Vec16c compress_to_int8 (Vec16i const & a) {
    return _mm512_cvtepi32_epi8(a);
}

static inline Vec16s compress_to_int16 (Vec16i const & a) {
    return _mm512_cvtepi32_epi16(a);
}

// with signed saturation
static inline Vec16c compress_to_int8_saturated (Vec16i const & a) {
    return _mm512_cvtsepi32_epi8(a);
}

static inline Vec16s compress_to_int16_saturated (Vec16i const & a) {
    return _mm512_cvtsepi32_epi16(a);
}

// with unsigned saturation
static inline Vec16uc compress_to_int8_saturated (Vec16ui const & a) {
    return _mm512_cvtusepi32_epi8(a);
}

static inline Vec16us compress_to_int16_saturated (Vec16ui const & a) {
    return _mm512_cvtusepi32_epi16(a);
}

// Compress 64-bit integers to 32-bit integers, signed and unsigned, with and without saturation

// Function compress : packs two vectors of 64-bit integers into one vector of 32-bit integers
// Overflow wraps around
static inline Vec16i compress (Vec8q const & low, Vec8q const & high) {
    Vec8i low2   = _mm512_cvtepi64_epi32(low);
    Vec8i high2  = _mm512_cvtepi64_epi32(high);
    return Vec16i(low2, high2);
}
static inline Vec16ui compress (Vec8uq const & low, Vec8uq const & high) {
    return Vec16ui(compress(Vec8q(low), Vec8q(high)));
}

// Function compress_saturated : packs two vectors of 64-bit integers into one vector of 32-bit integers
// Signed, with saturation
static inline Vec16i compress_saturated (Vec8q const & low, Vec8q const & high) {
    Vec8i low2   = _mm512_cvtsepi64_epi32(low);
    Vec8i high2  = _mm512_cvtsepi64_epi32(high);
    return Vec16i(low2, high2);
}

// Function compress_saturated : packs two vectors of 64-bit integers into one vector of 32-bit integers
// Unsigned, with saturation
static inline Vec16ui compress_saturated (Vec8uq const & low, Vec8uq const & high) {
    Vec8ui low2   = _mm512_cvtusepi64_epi32(low);
    Vec8ui high2  = _mm512_cvtusepi64_epi32(high);
    return Vec16ui(low2, high2);
}


/*****************************************************************************
*
*          Integer division operators
*
*          Please see the file vectori128.h for explanation.
*
*****************************************************************************/

// vector operator / : divide each element by divisor

static inline __m512i broadcast128_512 (__m128i x) {
#if defined(_MSC_VER) && !defined(__INTEL_COMPILER)   // fix bug in VS2019, 32 bit mode, release
    Vec4i m1 = x;
    Vec8i m2 = Vec8i(m1,m1);
    Vec16i m = Vec16i(m2,m2);
    return m;
#else
    return _mm512_broadcast_i32x4(x);
#endif 
}

// vector of 16 32-bit signed integers
static inline Vec16i operator / (Vec16i const & a, Divisor_i const & d) {
//    __m512i m   = _mm512_broadcast_i32x4(d.getm());      // broadcast multiplier
//    __m512i sgn = _mm512_broadcast_i32x4(d.getsign());   // broadcast sign of d

    __m512i m   = broadcast128_512(d.getm());              // broadcast multiplier
    __m512i sgn = broadcast128_512(d.getsign());           // broadcast sign of d

    __m512i t1  = _mm512_mul_epi32(a,m);                   // 32x32->64 bit signed multiplication of even elements of a
    __m512i t3  = _mm512_srli_epi64(a,32);                 // get odd elements of a into position for multiplication
    __m512i t4  = _mm512_mul_epi32(t3,m);                  // 32x32->64 bit signed multiplication of odd elements
    __m512i t2  = _mm512_srli_epi64(t1,32);                // dword of even index results
    __m512i t7  = _mm512_mask_mov_epi32(t2, __mmask16(0xAAAA), t4);  // blend two results
    __m512i t8  = _mm512_add_epi32(t7,a);                  // add
    __m512i t9  = _mm512_sra_epi32(t8,d.gets1());          // shift right artihmetic
    __m512i t10 = _mm512_srai_epi32(a,31);                 // sign of a
    __m512i t11 = _mm512_sub_epi32(t10,sgn);               // sign of a - sign of d
    __m512i t12 = _mm512_sub_epi32(t9,t11);                // + 1 if a < 0, -1 if d < 0
    return        _mm512_xor_si512(t12,sgn);               // change sign if divisor negative
}

// vector of 16 32-bit unsigned integers
static inline Vec16ui operator / (Vec16ui const & a, Divisor_ui const & d) {
    __m512i m   = broadcast128_512(d.getm());              // broadcast multiplier
    __m512i t1  = _mm512_mul_epu32(a,m);                   // 32x32->64 bit unsigned multiplication of even elements of a
    __m512i t3  = _mm512_srli_epi64(a,32);                 // get odd elements of a into position for multiplication
    __m512i t4  = _mm512_mul_epu32(t3,m);                  // 32x32->64 bit unsigned multiplication of odd elements
    __m512i t2  = _mm512_srli_epi64(t1,32);                // high dword of even index results
    __m512i t7  = _mm512_mask_mov_epi32(t2, __mmask16(0xAAAA), t4);  // blend two results
    __m512i t8  = _mm512_sub_epi32(a,t7);                  // subtract
    __m512i t9  = _mm512_srl_epi32(t8,d.gets1());          // shift right logical
    __m512i t10 = _mm512_add_epi32(t7,t9);                 // add
    return        _mm512_srl_epi32(t10,d.gets2());         // shift right logical 
}

// vector operator /= : divide
static inline Vec16i & operator /= (Vec16i & a, Divisor_i const & d) {
    a = a / d;
    return a;
}

// vector operator /= : divide
static inline Vec16ui & operator /= (Vec16ui & a, Divisor_ui const & d) {
    a = a / d;
    return a;
}


/*****************************************************************************
*
*          Integer division 2: divisor is a compile-time constant
*
*****************************************************************************/

// Divide Vec16i by compile-time constant
template <int32_t d>
static inline Vec16i divide_by_i(Vec16i const & x) {
    Static_error_check<(d!=0)> Dividing_by_zero;                     // Error message if dividing by zero
    if (d ==  1) return  x;
    if (d == -1) return -x;
    if (uint32_t(d) == 0x80000000u) {
        return _mm512_maskz_set1_epi32(x == Vec16i(0x80000000), 1);  // avoid overflow of abs(d). return (x == 0x80000000) ? 1 : 0;
    }
    const uint32_t d1 = d > 0 ? uint32_t(d) : -uint32_t(d);          // compile-time abs(d). (force GCC compiler to treat d as 32 bits, not 64 bits)
    if ((d1 & (d1-1)) == 0) {
        // d1 is a power of 2. use shift
        const int k = bit_scan_reverse_const(d1);
        __m512i sign;
        if (k > 1) sign = _mm512_srai_epi32(x, k-1); else sign = x;  // k copies of sign bit
        __m512i bias    = _mm512_srli_epi32(sign, 32-k);             // bias = x >= 0 ? 0 : k-1
        __m512i xpbias  = _mm512_add_epi32 (x, bias);                // x + bias
        __m512i q       = _mm512_srai_epi32(xpbias, k);              // (x + bias) >> k
        if (d > 0)      return q;                                    // d > 0: return  q
        return _mm512_sub_epi32(_mm512_setzero_epi32(), q);          // d < 0: return -q

    }
    // general case
    const int32_t sh = bit_scan_reverse_const(uint32_t(d1)-1);       // ceil(log2(d1)) - 1. (d1 < 2 handled by power of 2 case)
    const int32_t mult = int(1 + (uint64_t(1) << (32+sh)) / uint32_t(d1) - (int64_t(1) << 32));   // multiplier
    const Divisor_i div(mult, sh, d < 0 ? -1 : 0);
    return x / div;
}

// define Vec8i a / const_int(d)
template <int32_t d>
static inline Vec16i operator / (Vec16i const & a, Const_int_t<d>) {
    return divide_by_i<d>(a);
}

// define Vec16i a / const_uint(d)
template <uint32_t d>
static inline Vec16i operator / (Vec16i const & a, Const_uint_t<d>) {
    Static_error_check< (d<0x80000000u) > Error_overflow_dividing_signed_by_unsigned; // Error: dividing signed by overflowing unsigned
    return divide_by_i<int32_t(d)>(a);                               // signed divide
}

// vector operator /= : divide
template <int32_t d>
static inline Vec16i & operator /= (Vec16i & a, Const_int_t<d> b) {
    a = a / b;
    return a;
}

// vector operator /= : divide
template <uint32_t d>
static inline Vec16i & operator /= (Vec16i & a, Const_uint_t<d> b) {
    a = a / b;
    return a;
}


// Divide Vec16ui by compile-time constant
template <uint32_t d>
static inline Vec16ui divide_by_ui(Vec16ui const & x) {
    Static_error_check<(d!=0)> Dividing_by_zero;                     // Error message if dividing by zero
    if (d == 1) return x;                                            // divide by 1
    const int b = bit_scan_reverse_const(d);                         // floor(log2(d))
    if ((uint32_t(d) & (uint32_t(d)-1)) == 0) {
        // d is a power of 2. use shift
        return  _mm512_srli_epi32(x, b);                             // x >> b
    }
    // general case (d > 2)
    uint32_t mult = uint32_t((uint64_t(1) << (b+32)) / d);           // multiplier = 2^(32+b) / d
    const uint64_t rem = (uint64_t(1) << (b+32)) - uint64_t(d)*mult; // remainder 2^(32+b) % d
    const bool round_down = (2*rem < d);                             // check if fraction is less than 0.5
    if (!round_down) {
        mult = mult + 1;                                             // round up mult
    }
    // do 32*32->64 bit unsigned multiplication and get high part of result
    const __m512i multv = _mm512_maskz_set1_epi32(0x5555, mult);     // zero-extend mult and broadcast
    __m512i t1 = _mm512_mul_epu32(x,multv);                          // 32x32->64 bit unsigned multiplication of even elements
    if (round_down) {
        t1      = _mm512_add_epi64(t1,multv);                        // compensate for rounding error. (x+1)*m replaced by x*m+m to avoid overflow
    }
    __m512i t2 = _mm512_srli_epi64(t1,32);                           // high dword of result 0 and 2
    __m512i t3 = _mm512_srli_epi64(x,32);                            // get odd elements into position for multiplication
    __m512i t4 = _mm512_mul_epu32(t3,multv);                         // 32x32->64 bit unsigned multiplication of x[1] and x[3]
    if (round_down) {
        t4      = _mm512_add_epi64(t4,multv);                        // compensate for rounding error. (x+1)*m replaced by x*m+m to avoid overflow
    }
    __m512i t7 = _mm512_mask_mov_epi32(t2, __mmask16(0xAAAA), t4);   // blend two results
    Vec16ui q  = _mm512_srli_epi32(t7, b);                           // shift right by b
    return q;                                                        // no overflow possible
}

// define Vec8ui a / const_uint(d)
template <uint32_t d>
static inline Vec16ui operator / (Vec16ui const & a, Const_uint_t<d>) {
    return divide_by_ui<d>(a);
}

// define Vec8ui a / const_int(d)
template <int32_t d>
static inline Vec16ui operator / (Vec16ui const & a, Const_int_t<d>) {
    Static_error_check< (d>=0) > Error_dividing_unsigned_by_negative;// Error: dividing unsigned by negative is ambiguous
    return divide_by_ui<d>(a);                                       // unsigned divide
}

// vector operator /= : divide
template <uint32_t d>
static inline Vec16ui & operator /= (Vec16ui & a, Const_uint_t<d> b) {
    a = a / b;
    return a;
}

// vector operator /= : divide
template <int32_t d>
static inline Vec16ui & operator /= (Vec16ui & a, Const_int_t<d> b) {
    a = a / b;
    return a;
}

/*****************************************************************************
*
*          Horizontal scan functions
*
*****************************************************************************/

// Get index to the first element that is true. Return -1 if all are false
static inline int horizontal_find_first(Vec16ib const & x) {
    uint32_t b = uint16_t(__mmask16(x));
    if (b) {
        return (int)bit_scan_forward(b);
    }
    else {
        return -1;
    }
}

static inline int horizontal_find_first(Vec8qb const & x) {
    uint32_t b = uint8_t(__mmask16(x));
    if (b) {
        return (int)bit_scan_forward(b);
    }
    else {
        return -1;
    }
}

// Count number of true elements
static inline uint32_t horizontal_count(Vec16ib const & x) {
    return vml_popcnt(uint32_t(uint16_t(__mmask16(x))));
}

static inline uint32_t horizontal_count(Vec8qb const & x) {
    return vml_popcnt(uint32_t(uint8_t(__mmask16(x))));
}


/*****************************************************************************
*
*          Boolean <-> bitfield conversion functions
*
*****************************************************************************/

// to_bits: convert boolean vector to integer bitfield
static inline uint8_t to_bits(Vec4ib x) {
    __m512i a = _mm512_castsi128_si512(x);
    __mmask16 b = _mm512_mask_testn_epi32_mask(0xF, a, a);
    return uint8_t(b) ^ 0xFu;
}

// to_Vec16c: convert integer bitfield to boolean vector
// DEPRECATED. REPLACED BY load_bits
static inline Vec4ib to_Vec4ib(uint8_t x) {
    return _mm512_castsi512_si128(_mm512_maskz_set1_epi32(__mmask16(x), -1));
}

// to_bits: convert boolean vector to integer bitfield
static inline uint8_t to_bits(Vec2qb x) {
    __m512i a = _mm512_castsi128_si512(x);
    __mmask16 b = _mm512_mask_testn_epi64_mask(0x3, a, a);
    return uint8_t(b) ^ 0x3u;
}

// to_Vec16c: convert integer bitfield to boolean vector
// DEPRECATED. REPLACED BY load_bits
static inline Vec2qb to_Vec2qb(uint8_t x) {
    return _mm512_castsi512_si128(_mm512_maskz_set1_epi64(__mmask16(x), -1LL));
}

// to_bits: convert boolean vector to integer bitfield
static inline uint8_t to_bits(Vec8ib x) {
    __m512i a = _mm512_castsi256_si512(x);
    __mmask16 b = _mm512_mask_testn_epi32_mask(0xFF, a, a);
    return uint8_t(~ uint8_t(b));
}

static inline Vec8ib to_Vec8ib(uint8_t x) {
    // DEPRECATED. REPLACED BY load_bits
    return _mm512_castsi512_si256(_mm512_maskz_set1_epi32(__mmask16(x), -1));
}

// to_bits: convert boolean vector to integer bitfield
static inline uint8_t to_bits(Vec4qb x) {
    __m512i a = _mm512_castsi256_si512(x);
    __mmask16 b = _mm512_mask_testn_epi64_mask(0xF, a, a);
    return uint8_t(b) ^ 0xFu;
}

// to_Vec16c: convert integer bitfield to boolean vector
// DEPRECATED. REPLACED BY load_bits
static inline Vec4qb to_Vec4qb(uint8_t x) {
    return _mm512_castsi512_si256(_mm512_maskz_set1_epi64(__mmask16(x), -1LL));
}


// to_bits: convert to integer bitfield
static inline uint16_t to_bits(Vec16b a) {
    return (uint16_t)(__mmask16)a;
}

// to_Vec16b: convert integer bitfield to boolean vector
// DEPRECATED. REPLACED BY load_bits
static inline Vec16b to_Vec16b(uint16_t x) {
    return (__mmask16)x;
}

// to_Vec16ib: convert integer bitfield to boolean vector
// DEPRECATED. REPLACED BY load_bits
static inline Vec16ib to_Vec16ib(uint16_t x) {
    return to_Vec16b(x);
}

// to_Vec8b: convert integer bitfield to boolean vector
// DEPRECATED. REPLACED BY load_bits
static inline Vec8qb to_Vec8qb(uint8_t x) {
    return (__mmask16)x;
}

#ifdef VCL_NAMESPACE
}
#endif

#endif // VECTORI512_H