static void encode_component (float x, uint8_t *code, int i) {

code[i] = (int)(255 * x);

}

static float decode_component (const uint8_t *code, int i) {

return (code[i] + 0.5f) / 255.0f;

}

#ifdef USE_AVX

static __m256 decode_8_components (const uint8_t *code, int i) {

uint64_t c8 = *(uint64_t*)(code + i);

__m128i c4lo = _mm_cvtepu8_epi32 (_mm_set1_epi32(c8));

__m128i c4hi = _mm_cvtepu8_epi32 (_mm_set1_epi32(c8 >> 32));

// __m256i i8 = _mm256_set_m128i(c4lo, c4hi);

__m256i i8 = _mm256_castsi128_si256 (c4lo);

i8 = _mm256_insertf128_si256 (i8, c4hi, 1);

__m256 f8 = _mm256_cvtepi32_ps (i8);

__m256 half = _mm256_set1_ps (0.5f);

f8 += half;

__m256 one_255 = _mm256_set1_ps (1.f / 255.f);

return f8 * one_255;

}

#endif

static void encode_component (float x, uint8_t *code, int i) {

code [i / 2] |= (int)(x * 15.0) << ((i & 1) << 2);

}

static float decode_component (const uint8_t *code, int i) {

return (((code[i / 2] >> ((i & 1) << 2)) & 0xf) + 0.5f) / 15.0f;

}

#ifdef USE_AVX

static __m256 decode_8_components (const uint8_t *code, int i) {

uint32_t c4 = *(uint32_t*)(code + (i >> 1));

uint32_t mask = 0x0f0f0f0f;

uint32_t c4ev = c4 & mask;

uint32_t c4od = (c4 >> 4) & mask;

// the 8 lower bytes of c8 contain the values

__m128i c8 = _mm_unpacklo_epi8 (_mm_set1_epi32(c4ev),

_mm_set1_epi32(c4od));

__m128i c4lo = _mm_cvtepu8_epi32 (c8);

__m128i c4hi = _mm_cvtepu8_epi32 (_mm_srli_si128(c8, 4));

__m256i i8 = _mm256_castsi128_si256 (c4lo);

i8 = _mm256_insertf128_si256 (i8, c4hi, 1);

__m256 f8 = _mm256_cvtepi32_ps (i8);

__m256 half = _mm256_set1_ps (0.5f);

f8 += half;

__m256 one_255 = _mm256_set1_ps (1.f / 15.f);

return f8 * one_255;

}

#endif

static void encode_component (float x, uint8_t *code, int i) {

int bits = (int)(x * 63.0);

code += (i >> 2) * 3;

switch(i & 3) {

case 0:

code[0] |= bits;

break;

case 1:

code[0] |= bits << 6;

code[1] |= bits >> 2;

break;

case 2:

code[1] |= bits << 4;

code[2] |= bits >> 4;

break;

case 3:

code[2] |= bits << 2;

break;

}

}

static float decode_component (const uint8_t *code, int i) {

uint8_t bits;

code += (i >> 2) * 3;

switch(i & 3) {

case 0:

bits = code[0] & 0x3f;

break;

case 1:

bits = code[0] >> 6;

bits |= (code[1] & 0xf) << 2;

break;

case 2:

bits = code[1] >> 4;

bits |= (code[2] & 3) << 4;

break;

case 3:

bits = code[2] >> 2;

break;

}

return (bits + 0.5f) / 63.0f;

}

#ifdef USE_AVX

static __m256 decode_8_components (const uint8_t *code, int i) {

return _mm256_set_ps

(decode_component(code, i + 7),

decode_component(code, i + 6),

decode_component(code, i + 5),

decode_component(code, i + 4),

decode_component(code, i + 3),

decode_component(code, i + 2),

decode_component(code, i + 1),

decode_component(code, i + 0));

}

#endif

__m128 xf = _mm_set1_ps (x);

__m128i xi = _mm_cvtps_ph (

xf, _MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC);

return _mm_cvtsi128_si32 (xi) & 0xffff;

__m128i xi = _mm_set1_epi16 (x);

__m128 xf = _mm_cvtph_ps (xi);

return _mm_cvtss_f32 (xf);

// via Fabian "ryg" Giesen.

// https://gist.github.com/2156668

uint32_t sign_mask = 0x80000000u;

int32_t o;

uint32_t fint = intbits(f);

uint32_t sign = fint & sign_mask;

fint ^= sign;

// NOTE all the integer compares in this function can be safely

// compiled into signed compares since all operands are below

// 0x80000000. Important if you want fast straight SSE2 code (since

// there's no unsigned PCMPGTD).

// Inf or NaN (all exponent bits set)

// NaN->qNaN and Inf->Inf

// unconditional assignment here, will override with right value for

// the regular case below.

uint32_t f32infty = 255u << 23;

o = (fint > f32infty) ? 0x7e00u : 0x7c00u;

// (De)normalized number or zero

// update fint unconditionally to save the blending; we don't need it

// anymore for the Inf/NaN case anyway.

const uint32_t round_mask = ~0xfffu;

const uint32_t magic = 15u << 23;

// Shift exponent down, denormalize if necessary.

// NOTE This represents half-float denormals using single

// precision denormals. The main reason to do this is that

// there's no shift with per-lane variable shifts in SSE*, which

// we'd otherwise need. It has some funky side effects though:

// - This conversion will actually respect the FTZ (Flush To Zero)

// flag in MXCSR - if it's set, no half-float denormals will be

// generated. I'm honestly not sure whether this is good or

// bad. It's definitely interesting.

// - If the underlying HW doesn't support denormals (not an issue

// with Intel CPUs, but might be a problem on GPUs or PS3 SPUs),

// you will always get flush-to-zero behavior. This is bad,

// unless you're on a CPU where you don't care.

// - Denormals tend to be slow. FP32 denormals are rare in

// practice outside of things like recursive filters in DSP -

// not a typical half-float application. Whether FP16 denormals

// are rare in practice, I don't know. Whatever slow path your

// HW may or may not have for denormals, this may well hit it.

float fscale = floatbits(fint & round_mask) * floatbits(magic);

fscale = std::min(fscale, floatbits((31u << 23) - 0x1000u));

int32_t fint2 = intbits(fscale) - round_mask;

if (fint < f32infty)

o = fint2 >> 13; // Take the bits!

return (o | (sign >> 16));

// https://gist.github.com/2144712

// Fabian "ryg" Giesen.

const uint32_t shifted_exp = 0x7c00u << 13; // exponent mask after shift

int32_t o = ((int32_t)(h & 0x7fffu)) << 13; // exponent/mantissa bits

int32_t exp = shifted_exp & o; // just the exponent

o += (int32_t)(127 - 15) << 23; // exponent adjust

int32_t infnan_val = o + ((int32_t)(128 - 16) << 23);

int32_t zerodenorm_val = intbits(

floatbits(o + (1u<<23)) - floatbits(113u << 23));

int32_t reg_val = (exp == 0) ? zerodenorm_val : o;

int32_t sign_bit = ((int32_t)(h & 0x8000u)) << 16;

return floatbits(((exp == shifted_exp) ? infnan_val : reg_val) | sign_bit);

for (size_t i = 0; i < d; i++) {

float xi = (x[i] - vmin) / vdiff;

if (xi < 0) {

xi = 0;

}

if (xi > 1.0) {

xi = 1.0;

}

Codec::encode_component(xi, code, i);

}

for (size_t i = 0; i < d; i++) {

float xi = Codec::decode_component(code, i);

x[i] = vmin + xi * vdiff;

}

{

float xi = Codec::decode_component (code, i);

return vmin + xi * vdiff;

{

__m256 xi = Codec::decode_8_components (code, i);

return _mm256_set1_ps(this->vmin) + xi * _mm256_set1_ps (this->vdiff);

for (size_t i = 0; i < d; i++) {

float xi = (x[i] - vmin[i]) / vdiff[i];

if (xi < 0)

xi = 0;

if (xi > 1.0)

xi = 1.0;

Codec::encode_component(xi, code, i);

}

for (size_t i = 0; i < d; i++) {

float xi = Codec::decode_component(code, i);

x[i] = vmin[i] + xi * vdiff[i];

}

{

float xi = Codec::decode_component (code, i);

return vmin[i] + xi * vdiff[i];

{

__m256 xi = Codec::decode_8_components (code, i);

return _mm256_loadu_ps (this->vmin + i) + xi * _mm256_loadu_ps (this->vdiff + i);

for (size_t i = 0; i < d; i++) {

((uint16_t*)code)[i] = encode_fp16(x[i]);

}

for (size_t i = 0; i < d; i++) {

x[i] = decode_fp16(((uint16_t*)code)[i]);

}

{

__m128i codei = _mm_loadu_si128 ((const __m128i*)(code + 2 * i));

return _mm256_cvtph_ps (codei);

for (size_t i = 0; i < d; i++) {

code[i] = (uint8_t)x[i];

}

for (size_t i = 0; i < d; i++) {

x[i] = code[i];

}

{

__m128i x8 = _mm_loadl_epi64((__m128i*)(code + i)); // 8 * int8

__m256i y8 = _mm256_cvtepu8_epi32 (x8); // 8 * int32

return _mm256_cvtepi32_ps (y8); // 8 * float32

idx_t n, int k, const float *x,

std::vector<float> & trained)

{

trained.resize (2);

float & vmin = trained[0];

float & vmax = trained[1];

if (rs == ScalarQuantizer::RS_minmax) {

vmin = HUGE_VAL; vmax = -HUGE_VAL;

for (size_t i = 0; i < n; i++) {

if (x[i] < vmin) vmin = x[i];

if (x[i] > vmax) vmax = x[i];

}

float vexp = (vmax - vmin) * rs_arg;

vmin -= vexp;

vmax += vexp;

} else if (rs == ScalarQuantizer::RS_meanstd) {

double sum = 0, sum2 = 0;

for (size_t i = 0; i < n; i++) {

sum += x[i];

sum2 += x[i] * x[i];

}

float mean = sum / n;

float var = sum2 / n - mean * mean;

float std = var <= 0 ? 1.0 : sqrt(var);

vmin = mean - std * rs_arg ;

vmax = mean + std * rs_arg ;

} else if (rs == ScalarQuantizer::RS_quantiles) {

std::vector<float> x_copy(n);

memcpy(x_copy.data(), x, n * sizeof(*x));

// TODO just do a qucikselect

std::sort(x_copy.begin(), x_copy.end());

int o = int(rs_arg * n);

if (o < 0) o = 0;

if (o > n - o) o = n / 2;

vmin = x_copy[o];

vmax = x_copy[n - 1 - o];

} else if (rs == ScalarQuantizer::RS_optim) {

float a, b;

float sx = 0;

{

vmin = HUGE_VAL, vmax = -HUGE_VAL;

for (size_t i = 0; i < n; i++) {

if (x[i] < vmin) vmin = x[i];

if (x[i] > vmax) vmax = x[i];

sx += x[i];

}

b = vmin;

a = (vmax - vmin) / (k - 1);

}

int verbose = false;

int niter = 2000;

float last_err = -1;

int iter_last_err = 0;

for (int it = 0; it < niter; it++) {

float sn = 0, sn2 = 0, sxn = 0, err1 = 0;

for (idx_t i = 0; i < n; i++) {

float xi = x[i];

float ni = floor ((xi - b) / a + 0.5);

if (ni < 0) ni = 0;

if (ni >= k) ni = k - 1;

err1 += sqr (xi - (ni * a + b));

sn += ni;

sn2 += ni * ni;

sxn += ni * xi;

}

if (err1 == last_err) {

iter_last_err ++;

if (iter_last_err == 16) break;

} else {

last_err = err1;

iter_last_err = 0;

}

float det = sqr (sn) - sn2 * n;

b = (sn * sxn - sn2 * sx) / det;

a = (sn * sx - n * sxn) / det;

if (verbose) {

printf ("it %d, err1=%g \r", it, err1);

fflush(stdout);

}

}

if (verbose) printf("\n");

vmin = b;

vmax = b + a * (k - 1);

} else {

FAISS_THROW_MSG ("Invalid qtype");

}

vmax -= vmin;

idx_t n, int d, int k, const float *x,

std::vector<float> & trained)

{

trained.resize (2 * d);

float * vmin = trained.data();

float * vmax = trained.data() + d;

if (rs == ScalarQuantizer::RS_minmax) {

memcpy (vmin, x, sizeof(*x) * d);

memcpy (vmax, x, sizeof(*x) * d);

for (size_t i = 1; i < n; i++) {

const float *xi = x + i * d;

for (size_t j = 0; j < d; j++) {

if (xi[j] < vmin[j]) vmin[j] = xi[j];

if (xi[j] > vmax[j]) vmax[j] = xi[j];

}

}

float *vdiff = vmax;

for (size_t j = 0; j < d; j++) {

float vexp = (vmax[j] - vmin[j]) * rs_arg;

vmin[j] -= vexp;

vmax[j] += vexp;

vdiff [j] = vmax[j] - vmin[j];

}

} else {

// transpose

std::vector<float> xt(n * d);

for (size_t i = 1; i < n; i++) {

const float *xi = x + i * d;

for (size_t j = 0; j < d; j++) {

xt[j * n + i] = xi[j];

}

}

std::vector<float> trained_d(2);

#pragma omp parallel for

for (size_t j = 0; j < d; j++) {

train_Uniform(rs, rs_arg,

n, k, xt.data() + j * n,

trained_d);

vmin[j] = trained_d[0];

vmax[j] = trained_d[1];

}

}

__m256 yiv = _mm256_loadu_ps (yi);

yi += 8;

__m256 tmp = yiv - x;

accu8 += tmp * tmp;

__m256 sum = _mm256_hadd_ps(accu8, accu8);

__m256 sum2 = _mm256_hadd_ps(sum, sum);

// now add the 0th and 4th component

return

_mm_cvtss_f32 (_mm256_castps256_ps128(sum2)) +

_mm_cvtss_f32 (_mm256_extractf128_ps(sum2, 1));

__m256 yiv = _mm256_loadu_ps (yi);

yi += 8;

accu8 += yiv * x;

__m256 sum = _mm256_hadd_ps(accu8, accu8);

__m256 sum2 = _mm256_hadd_ps(sum, sum);

// now add the 0th and 4th component

return

_mm_cvtss_f32 (_mm256_castps256_ps128(sum2)) +

_mm_cvtss_f32 (_mm256_extractf128_ps(sum2, 1));

Similarity sim(x);

sim.begin();

for (size_t i = 0; i < quant.d; i++) {

float xi = quant.reconstruct_component(code, i);

sim.add_component(xi);

}

return sim.result();

const {

Similarity sim(nullptr);

sim.begin();

for (size_t i = 0; i < quant.d; i++) {

float x1 = quant.reconstruct_component(code1, i);

float x2 = quant.reconstruct_component(code2, i);

sim.add_component_2(x1, x2);

}

return sim.result();

Similarity sim(x);

sim.begin_8();

for (size_t i = 0; i < quant.d; i += 8) {

__m256 xi = quant.reconstruct_8_components(code, i);

sim.add_8_components(xi);

}

return sim.result_8();

const {

Similarity sim(nullptr);

sim.begin_8();

for (size_t i = 0; i < quant.d; i += 8) {

__m256 x1 = quant.reconstruct_8_components(code1, i);

__m256 x2 = quant.reconstruct_8_components(code2, i);

sim.add_8_components_2(x1, x2);

}

return sim.result_8();