整数压缩：varint → PForDelta → SIMD-BP128

// delta decode = prefix sum
void delta_decode(uint32_t *data, size_t n) {
    for (size_t i = 1; i < n; i++) {
        data[i] += data[i - 1];
    }
}

// zigzag encode: 将有符号整数映射到无符号整数
static inline uint32_t zigzag_encode(int32_t n) {
    return (uint32_t)((n << 1) ^ (n >> 31));
}

// zigzag decode: 逆映射
static inline int32_t zigzag_decode(uint32_t n) {
    return (int32_t)((n >> 1) ^ -(n & 1));
}

#include <stdint.h>
#include <stddef.h>

// 编码一个 uint32，返回写入的字节数
size_t varint_encode_u32(uint8_t *buf, uint32_t value) {
    size_t i = 0;
    while (value >= 0x80) {
        buf[i++] = (uint8_t)(value | 0x80);
        value >>= 7;
    }
    buf[i++] = (uint8_t)value;
    return i;
}

// 解码一个 uint32，返回消耗的字节数
size_t varint_decode_u32(const uint8_t *buf, uint32_t *out) {
    uint32_t result = 0;
    size_t i = 0;
    uint32_t shift = 0;
    do {
        result |= (uint32_t)(buf[i] & 0x7F) << shift;
        shift += 7;
    } while (buf[i++] & 0x80);
    *out = result;
    return i;
}

#include <string.h>

// 预计算表：tag -> 每个整数的长度
static uint8_t gv_len_table[256][4];

void gv_init_tables(void) {
    for (int tag = 0; tag < 256; tag++) {
        gv_len_table[tag][0] = ((tag >> 0) & 3) + 1;
        gv_len_table[tag][1] = ((tag >> 2) & 3) + 1;
        gv_len_table[tag][2] = ((tag >> 4) & 3) + 1;
        gv_len_table[tag][3] = ((tag >> 6) & 3) + 1;
    }
}

// 解码 4 个整数，返回消耗的总字节数
size_t gv_decode4(const uint8_t *in, uint32_t out[4]) {
    uint8_t tag = in[0];
    const uint8_t *p = in + 1;
    for (int i = 0; i < 4; i++) {
        uint8_t len = gv_len_table[tag][i];
        uint32_t val = 0;
        memcpy(&val, p, len);  // 小端机器上直接 memcpy
        out[i] = val;
        p += len;
    }
    return (size_t)(p - in);
}

#include <stdint.h>

// Simple-8b selector 表
static const struct {
    int count;
    int bits;
} s8b_table[15] = {
    {240, 0}, {60, 1}, {30, 2}, {20, 3}, {15, 4},
    {12, 5},  {10, 6}, {8, 7},  {7, 8},  {6, 10},
    {5, 12},  {4, 15}, {3, 20}, {2, 30}, {1, 60}
};

// 尝试用 selector sel 编码 data 中的前若干个整数
// 成功返回编码的整数个数，失败返回 0
int s8b_try_pack(const uint32_t *data, int n, int sel, uint64_t *out) {
    int count = s8b_table[sel].count;
    int bits  = s8b_table[sel].bits;
    if (count > n) count = n;
    uint64_t mask = bits ? ((1ULL << bits) - 1) : 0;

    uint64_t word = (uint64_t)sel << 60;
    for (int i = 0; i < count; i++) {
        if (bits > 0 && data[i] > mask) return 0;
        word |= (uint64_t)data[i] << (i * bits);
    }
    *out = word;
    return count;
}

// 编码：选择能打包最多整数的 selector
int s8b_encode(const uint32_t *data, int n, uint64_t *out) {
    for (int sel = 0; sel < 15; sel++) {
        int packed = s8b_try_pack(data, n, sel, out);
        if (packed > 0) return packed;
    }
    return 0;  // should not happen for valid input
}

#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#define BLOCK_SIZE 128

typedef struct {
    uint8_t  bit_width;
    uint8_t  num_exceptions;
    uint32_t packed[BLOCK_SIZE];    // bit-packed base values
    uint8_t  exc_positions[BLOCK_SIZE];
    uint32_t exc_values[BLOCK_SIZE];
} pfd_block_t;

// 找到能覆盖 90% 值的最小位宽
static uint8_t find_bit_width(const uint32_t *deltas, int n) {
    // 计算每个位宽能覆盖的值的比例
    int counts[33] = {0};
    for (int i = 0; i < n; i++) {
        uint32_t v = deltas[i];
        int bits = 0;
        while (v > 0) { bits++; v >>= 1; }
        counts[bits]++;
    }

    int cumulative = 0;
    int threshold = (int)(n * 0.9);
    for (int b = 0; b <= 32; b++) {
        cumulative += counts[b];
        if (cumulative >= threshold) return (uint8_t)b;
    }
    return 32;
}

// PForDelta 编码一个块
void pfd_encode_block(const uint32_t *deltas, int n, pfd_block_t *blk) {
    blk->bit_width = find_bit_width(deltas, n);
    uint32_t mask = blk->bit_width < 32
                    ? (1U << blk->bit_width) - 1 : 0xFFFFFFFF;
    blk->num_exceptions = 0;

    for (int i = 0; i < n; i++) {
        blk->packed[i] = deltas[i] & mask;
        if (deltas[i] > mask) {
            blk->exc_positions[blk->num_exceptions] = (uint8_t)i;
            blk->exc_values[blk->num_exceptions]    = deltas[i];
            blk->num_exceptions++;
        }
    }
}

// PForDelta 解码一个块
void pfd_decode_block(const pfd_block_t *blk, uint32_t *out, int n) {
    // 第一步：复制 bit-packed 基础值
    memcpy(out, blk->packed, n * sizeof(uint32_t));

    // 第二步：用 exception 表修补
    for (int i = 0; i < blk->num_exceptions; i++) {
        out[blk->exc_positions[i]] = blk->exc_values[i];
    }
}

/*
 * simd_bitpack.c
 *
 * SIMD-friendly bit-packing / unpacking for 128-integer blocks.
 * Uses SSE2 intrinsics on x86; falls back to scalar on other platforms.
 *
 * Compile: gcc -O2 -msse2 -o simd_bitpack simd_bitpack.c
 * Run:     ./simd_bitpack
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <time.h>

#ifdef __SSE2__
#include <emmintrin.h>
#define USE_SIMD 1
#else
#define USE_SIMD 0
#endif

#define BLOCK_SIZE   128
#define MINI_BLOCK   32   /* SIMD processes 32 values at a time (8 x 128-bit) */
#define MAX_BITS     32

/* ------------------------------------------------------------------ */
/*  Scalar bit-pack: pack BLOCK_SIZE uint32s with uniform bit-width   */
/* ------------------------------------------------------------------ */
void scalar_pack(const uint32_t *in, uint8_t *out, int bit_width) {
    if (bit_width == 0) return;
    uint64_t buffer = 0;
    int bits_in_buf = 0;
    int out_pos = 0;
    uint32_t mask = (bit_width == 32) ? 0xFFFFFFFF : (1U << bit_width) - 1;

    for (int i = 0; i < BLOCK_SIZE; i++) {
        buffer |= ((uint64_t)(in[i] & mask)) << bits_in_buf;
        bits_in_buf += bit_width;
        while (bits_in_buf >= 8) {
            out[out_pos++] = (uint8_t)(buffer & 0xFF);
            buffer >>= 8;
            bits_in_buf -= 8;
        }
    }
    if (bits_in_buf > 0) {
        out[out_pos++] = (uint8_t)(buffer & 0xFF);
    }
}

/* ------------------------------------------------------------------ */
/*  Scalar bit-unpack                                                  */
/* ------------------------------------------------------------------ */
void scalar_unpack(const uint8_t *in, uint32_t *out, int bit_width) {
    if (bit_width == 0) {
        memset(out, 0, BLOCK_SIZE * sizeof(uint32_t));
        return;
    }
    uint64_t buffer = 0;
    int bits_in_buf = 0;
    int in_pos = 0;
    uint32_t mask = (bit_width == 32) ? 0xFFFFFFFF : (1U << bit_width) - 1;

    for (int i = 0; i < BLOCK_SIZE; i++) {
        while (bits_in_buf < bit_width) {
            buffer |= ((uint64_t)in[in_pos++]) << bits_in_buf;
            bits_in_buf += 8;
        }
        out[i] = (uint32_t)(buffer & mask);
        buffer >>= bit_width;
        bits_in_buf -= bit_width;
    }
}

/* ------------------------------------------------------------------ */
/*  SIMD-accelerated unpack (SSE2)                                     */
/*  Processes 4 integers at a time using 128-bit registers             */
/* ------------------------------------------------------------------ */
#if USE_SIMD
void simd_unpack(const uint8_t *in, uint32_t *out, int bit_width) {
    if (bit_width == 0) {
        __m128i zero = _mm_setzero_si128();
        for (int i = 0; i < BLOCK_SIZE; i += 4) {
            _mm_storeu_si128((__m128i *)(out + i), zero);
        }
        return;
    }
    if (bit_width == 32) {
        memcpy(out, in, BLOCK_SIZE * sizeof(uint32_t));
        return;
    }

    uint32_t mask = (1U << bit_width) - 1;
    __m128i vmask = _mm_set1_epi32((int)mask);
    uint64_t buffer = 0;
    int bits_in_buf = 0;
    int in_pos = 0;

    for (int i = 0; i < BLOCK_SIZE; i += 4) {
        uint32_t tmp[4];
        for (int j = 0; j < 4; j++) {
            while (bits_in_buf < bit_width) {
                buffer |= ((uint64_t)in[in_pos++]) << bits_in_buf;
                bits_in_buf += 8;
            }
            tmp[j] = (uint32_t)(buffer & mask);
            buffer >>= bit_width;
            bits_in_buf -= bit_width;
        }
        __m128i v = _mm_loadu_si128((const __m128i *)tmp);
        v = _mm_and_si128(v, vmask);
        _mm_storeu_si128((__m128i *)(out + i), v);
    }
}
#endif

/* ------------------------------------------------------------------ */
/*  Delta encode / decode                                              */
/* ------------------------------------------------------------------ */
void delta_encode(const uint32_t *in, uint32_t *out, int n) {
    out[0] = in[0];
    for (int i = 1; i < n; i++) {
        out[i] = in[i] - in[i - 1];
    }
}

void delta_decode(uint32_t *data, int n) {
    for (int i = 1; i < n; i++) {
        data[i] += data[i - 1];
    }
}

#if USE_SIMD
void delta_decode_simd(uint32_t *data, int n) {
    /* 利用 SIMD 做前缀和：分阶段累加 */
    /* 阶段1：4-wide 局部前缀和 */
    for (int i = 0; i < n; i += 4) {
        __m128i v = _mm_loadu_si128((const __m128i *)(data + i));
        /* shift-and-add pattern for prefix sum within 128-bit lane */
        v = _mm_add_epi32(v, _mm_slli_si128(v, 4));
        v = _mm_add_epi32(v, _mm_slli_si128(v, 8));
        _mm_storeu_si128((__m128i *)(data + i), v);
    }
    /* 阶段2：块间传播 */
    for (int i = 4; i < n; i += 4) {
        __m128i prev = _mm_set1_epi32((int)data[i - 1]);
        __m128i curr = _mm_loadu_si128((const __m128i *)(data + i));
        curr = _mm_add_epi32(curr, prev);
        _mm_storeu_si128((__m128i *)(data + i), curr);
    }
}
#endif

/* ------------------------------------------------------------------ */
/*  Compute required bit-width for a block                             */
/* ------------------------------------------------------------------ */
int compute_bit_width(const uint32_t *data, int n) {
    uint32_t max_val = 0;
    for (int i = 0; i < n; i++) {
        if (data[i] > max_val) max_val = data[i];
    }
    if (max_val == 0) return 0;
    int bits = 0;
    while (max_val > 0) {
        bits++;
        max_val >>= 1;
    }
    return bits;
}

/* ------------------------------------------------------------------ */
/*  Compressed block: header + packed data                             */
/* ------------------------------------------------------------------ */
typedef struct {
    uint8_t  bit_width;
    uint32_t first_value;
    uint8_t  packed_data[BLOCK_SIZE * 4]; /* max 128 * 32 / 8 = 512 */
    int      packed_bytes;
} compressed_block_t;

void compress_block(const uint32_t *in, compressed_block_t *blk) {
    uint32_t deltas[BLOCK_SIZE];
    delta_encode(in, deltas, BLOCK_SIZE);

    blk->first_value = deltas[0];
    blk->bit_width = (uint8_t)compute_bit_width(deltas + 1, BLOCK_SIZE - 1);
    blk->packed_bytes = (int)(((BLOCK_SIZE - 1) * blk->bit_width + 7) / 8);

    scalar_pack(deltas + 1, blk->packed_data, blk->bit_width);
}

void decompress_block(const compressed_block_t *blk, uint32_t *out) {
    out[0] = blk->first_value;
#if USE_SIMD
    simd_unpack(blk->packed_data, out + 1, blk->bit_width);
#else
    scalar_unpack(blk->packed_data, out + 1, blk->bit_width);
#endif
    delta_decode(out, BLOCK_SIZE);
}

/* ------------------------------------------------------------------ */
/*  Benchmark helper                                                   */
/* ------------------------------------------------------------------ */
static double time_diff_ms(struct timespec *start, struct timespec *end) {
    return (end->tv_sec - start->tv_sec) * 1000.0
         + (end->tv_nsec - start->tv_nsec) / 1e6;
}

/* ------------------------------------------------------------------ */
/*  Main: generate test data, compress, decompress, verify, benchmark  */
/* ------------------------------------------------------------------ */
int main(void) {
    printf("SIMD-friendly bit-packing demo (BLOCK_SIZE = %d)\n", BLOCK_SIZE);
    printf("SIMD support: %s\n\n", USE_SIMD ? "SSE2" : "scalar");

    /* 生成有序整数序列，模拟 posting list */
    uint32_t original[BLOCK_SIZE];
    uint32_t decoded[BLOCK_SIZE];
    srand(42);
    original[0] = (uint32_t)(rand() % 100);
    for (int i = 1; i < BLOCK_SIZE; i++) {
        original[i] = original[i - 1] + 1 + (uint32_t)(rand() % 20);
    }

    printf("Original (first 16): ");
    for (int i = 0; i < 16; i++) printf("%u ", original[i]);
    printf("...\n");

    /* 压缩 */
    compressed_block_t blk;
    compress_block(original, &blk);

    int uncompressed_bytes = BLOCK_SIZE * 4;
    int compressed_bytes = 1 + 4 + blk.packed_bytes;  /* header + data */
    printf("Bit-width: %d\n", blk.bit_width);
    printf("Compressed: %d bytes (from %d bytes, ratio: %.2fx)\n",
           compressed_bytes, uncompressed_bytes,
           (double)uncompressed_bytes / compressed_bytes);
    printf("Bits per integer: %.2f\n\n",
           (double)compressed_bytes * 8 / BLOCK_SIZE);

    /* 解压并验证 */
    decompress_block(&blk, decoded);
    int ok = 1;
    for (int i = 0; i < BLOCK_SIZE; i++) {
        if (original[i] != decoded[i]) { ok = 0; break; }
    }
    printf("Correctness: %s\n\n", ok ? "PASS" : "FAIL");

    /* 基准测试：解压速度 */
    const int ROUNDS = 1000000;
    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);
    for (int r = 0; r < ROUNDS; r++) {
        decompress_block(&blk, decoded);
    }
    clock_gettime(CLOCK_MONOTONIC, &t1);
    double ms = time_diff_ms(&t0, &t1);
    double total_ints = (double)ROUNDS * BLOCK_SIZE;
    double gints_per_sec = total_ints / (ms / 1000.0) / 1e9;
    printf("Benchmark: %d rounds, %.2f ms\n", ROUNDS, ms);
    printf("Throughput: %.2f Gints/s (%.2f billion integers per second)\n",
           gints_per_sec, gints_per_sec);

    return 0;
}

gcc -O2 -msse2 -o simd_bitpack simd_bitpack.c -lrt
./simd_bitpack

typedef struct {
    uint16_t high_bits;      // 高 16 位前缀
    uint8_t  container_type; // 0=array, 1=bitmap, 2=run
    union {
        struct {
            uint16_t *values;
            int       count;
        } array;
        struct {
            uint64_t bits[1024]; // 65536 / 64 = 1024
        } bitmap;
        struct {
            uint16_t *runs; // [start0, len0, start1, len1, ...]
            int       num_runs;
        } run;
    } data;
} roaring_container_t;

// 检查一个值是否存在
int roaring_contains(const roaring_container_t *c, uint16_t low) {
    switch (c->container_type) {
    case 0: { // array: 二分查找
        int lo = 0, hi = c->data.array.count - 1;
        while (lo <= hi) {
            int mid = (lo + hi) / 2;
            if (c->data.array.values[mid] == low) return 1;
            if (c->data.array.values[mid] < low) lo = mid + 1;
            else hi = mid - 1;
        }
        return 0;
    }
    case 1: // bitmap: 直接位查询
        return (c->data.bitmap.bits[low >> 6] >> (low & 63)) & 1;
    case 2: { // run: 二分查找区间
        for (int i = 0; i < c->data.run.num_runs; i++) {
            uint16_t start = c->data.run.runs[i * 2];
            uint16_t len   = c->data.run.runs[i * 2 + 1];
            if (low >= start && low <= start + len) return 1;
            if (low < start) break;
        }
        return 0;
    }
    }
    return 0;
}