/* cdc_sim.c -- content-defined chunking: size distribution, boundary shift,
 * degenerate input and relative throughput.
 *
 * gcc -O2 -Wall -Wextra -o cdc_sim cdc_sim.c
 * ./cdc_sim [hist.tsv]
 *
 * All chunkers target an 8 KiB average (13 mask bits); the bounded ones use the
 * LBFS limits min 2 KiB / max 64 KiB. Input is 64 MiB from splitmix64(seed 1).
 */
#define _POSIX_C_SOURCE 199309L
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#define KiB 1024u
#define MiB (1024u * 1024u)
#define WIN 48

static uint64_t sm_state;
static uint64_t splitmix64(void)
{
    uint64_t z = (sm_state += 0x9e3779b97f4a7c15ULL);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
    z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
    return z ^ (z >> 31);
}

/* ---------- Rabin fingerprint over GF(2), random irreducible P of degree 53 ---------- */
#define DEG 53
static uint64_t rb_pol, rb_out[256], rb_mod[256];

static uint64_t gf2_mulmod(uint64_t a, uint64_t b, uint64_t p)
{
    unsigned __int128 r = 0;
    for (int i = 0; i < DEG; i++)
        if (b >> i & 1) r ^= (unsigned __int128)a << i;
    for (int i = 2 * DEG - 2; i >= DEG; i--)
        if ((uint64_t)(r >> i) & 1) r ^= (unsigned __int128)p << (i - DEG);
    return (uint64_t)r;
}

/* Rabin's test for prime degree k: x^(2^k) = x mod P and gcd(x^2 - x, P) = 1,
 * i.e. P(0) = P(1) = 1. */
static int irreducible53(uint64_t p)
{
    if (!(p & 1) || !(__builtin_popcountll(p) & 1)) return 0;
    uint64_t x = 2;
    for (int i = 0; i < DEG; i++) x = gf2_mulmod(x, x, p);
    return x == 2;
}

static uint64_t pol_mod(uint64_t v, uint64_t p)   /* v may have degree up to 63 */
{
    for (int i = 63; i >= DEG; i--)
        if (v >> i & 1) v ^= p << (i - DEG);
    return v;
}

static void rabin_init(uint64_t seed)
{
    sm_state = seed;
    do {
        rb_pol = (splitmix64() & ((1ULL << DEG) - 1)) | (1ULL << DEG) | 1;
    } while (!irreducible53(rb_pol));
    for (int b = 0; b < 256; b++) {
        uint64_t h = (uint64_t)b;                     /* b * x^(8*(WIN-1)) mod P */
        for (int i = 0; i < WIN - 1; i++) h = pol_mod(h << 8, rb_pol);
        rb_out[b] = h;
        rb_mod[b] = pol_mod((uint64_t)b << DEG, rb_pol) | ((uint64_t)b << DEG);
    }
}

static inline uint64_t rabin_append(uint64_t d, uint8_t b)
{
    uint64_t idx = d >> (DEG - 8);
    return ((d << 8) | b) ^ rb_mod[idx];
}

/* ---------- 64-bit Buzhash, Gear ---------- */
static uint64_t bz_tab[256], gear_tab[256];
static inline uint64_t rotl(uint64_t x, unsigned r) { r &= 63; return r ? x << r | x >> (64 - r) : x; }

/* ---------- chunkers: return the length of the chunk starting at p[0] ---------- */
typedef size_t (*chunker_fn)(const uint8_t *p, size_t n, size_t mn, size_t mx);

static size_t cut_rabin(const uint8_t *p, size_t n, size_t mn, size_t mx)
{
    const uint64_t mask = (1u << 13) - 1;
    if (n <= mn) return n;
    if (n > mx) n = mx;
    size_t i = mn > WIN ? mn - WIN : 0;
    uint64_t d = 0;
    for (size_t j = i; j < i + WIN && j < n; j++) {   /* fill window, no removal */
        d = rabin_append(d, p[j]);
        if (j + 1 >= mn && j + 1 >= WIN && (d & mask) == 0) return j + 1;
    }
    for (i += WIN; i < n; i++) {
        d ^= rb_out[p[i - WIN]];
        d = rabin_append(d, p[i]);
        if (i + 1 >= mn && (d & mask) == 0) return i + 1;
    }
    return n;
}

static size_t cut_buz(const uint8_t *p, size_t n, size_t mn, size_t mx)
{
    const uint64_t mask = (1u << 13) - 1;
    if (n <= mn) return n;
    if (n > mx) n = mx;
    size_t i = mn > WIN ? mn - WIN : 0;
    uint64_t h = 0;
    for (size_t j = i; j < i + WIN && j < n; j++) {
        h = rotl(h, 1) ^ bz_tab[p[j]];
        if (j + 1 >= mn && j + 1 >= WIN && (h & mask) == 0) return j + 1;
    }
    for (i += WIN; i < n; i++) {
        h = rotl(h, 1) ^ rotl(bz_tab[p[i - WIN]], WIN) ^ bz_tab[p[i]];
        if (i + 1 >= mn && (h & mask) == 0) return i + 1;
    }
    return n;
}

static size_t gear_mask(const uint8_t *p, size_t n, size_t mn, size_t mx, uint64_t mask)
{
    if (n <= mn) return n;
    if (n > mx) n = mx;
    uint64_t fp = 0;
    for (size_t i = mn; i < n; i++) {
        fp = (fp << 1) + gear_tab[p[i]];
        if (!(fp & mask)) return i + 1;
    }
    return n;
}
static size_t cut_gear(const uint8_t *p, size_t n, size_t mn, size_t mx)
{ return gear_mask(p, n, mn, mx, (1u << 13) - 1); }
static size_t cut_fastcdc_nc0(const uint8_t *p, size_t n, size_t mn, size_t mx)
{ return gear_mask(p, n, mn, mx, 0x0000d90303530000ULL); }   /* MaskA, 13 bits */

/* FastCDC Algorithm 1 (ATC'16), normalization level 2: MaskS 15 bits, MaskL 11 bits. */
static size_t cut_fastcdc_nc2(const uint8_t *p, size_t n, size_t mn, size_t mx)
{
    const uint64_t maskS = 0x0003590703530000ULL, maskL = 0x0000d90003530000ULL;
    size_t normal = 8 * KiB;
    if (n <= mn) return n;
    if (n >= mx) n = mx;
    else if (n <= normal) normal = n;
    uint64_t fp = 0;
    size_t i = mn;
    for (; i < normal; i++) {
        fp = (fp << 1) + gear_tab[p[i]];
        if (!(fp & maskS)) return i + 1;
    }
    for (; i < n; i++) {
        fp = (fp << 1) + gear_tab[p[i]];
        if (!(fp & maskL)) return i + 1;
    }
    return n;
}

static size_t cut_fixed(const uint8_t *p, size_t n, size_t mn, size_t mx)
{ (void)p; (void)mn; (void)mx; return n < 8 * KiB ? n : 8 * KiB; }

typedef struct { const char *name; chunker_fn fn; size_t mn, mx; } Chunker;
static const Chunker CH[] = {
    {"rabin-unbounded", cut_rabin, 0, (size_t)-1},
    {"rabin",           cut_rabin, 2 * KiB, 64 * KiB},
    {"buzhash",         cut_buz,   2 * KiB, 64 * KiB},
    {"gear",            cut_gear,  2 * KiB, 64 * KiB},
    {"fastcdc-nc0",     cut_fastcdc_nc0, 2 * KiB, 64 * KiB},
    {"fastcdc-nc2",     cut_fastcdc_nc2, 2 * KiB, 64 * KiB},
    {"fixed-8k",        cut_fixed, 0, 0},
};
#define NCH (sizeof CH / sizeof CH[0])

/* Split buf into chunks; store lengths in *out (caller frees). */
static size_t split(const Chunker *c, const uint8_t *buf, size_t n, uint32_t **out)
{
    size_t cap = n / 1024 + 16, k = 0;
    uint32_t *len = malloc(cap * sizeof *len);
    if (!len) exit(1);
    for (size_t pos = 0; pos < n; ) {
        size_t l = c->fn(buf + pos, n - pos, c->mn, c->mx);
        if (k == cap) { cap *= 2; len = realloc(len, cap * sizeof *len); if (!len) exit(1); }
        len[k++] = (uint32_t)l;
        pos += l;
    }
    *out = len;
    return k;
}

/* ---------- chunk identity: 64-bit hash + memcmp confirmation ---------- */
typedef struct { uint64_t h; const uint8_t *p; uint32_t len; } Ent;

static uint64_t fnv1a64(const uint8_t *p, size_t n)
{
    uint64_t h = 0xcbf29ce484222325ULL;
    for (size_t i = 0; i < n; i++) { h ^= p[i]; h *= 0x100000001b3ULL; }
    return h;
}
static int cmp_ent(const void *a, const void *b)
{
    const Ent *x = a, *y = b;
    return x->h < y->h ? -1 : x->h > y->h;
}
static int present(const Ent *set, size_t k, const uint8_t *p, uint32_t len)
{
    uint64_t h = fnv1a64(p, len);
    size_t lo = 0, hi = k;
    while (lo < hi) { size_t m = (lo + hi) / 2; if (set[m].h < h) lo = m + 1; else hi = m; }
    for (; lo < k && set[lo].h == h; lo++)
        if (set[lo].len == len && memcmp(set[lo].p, p, len) == 0) return 1;
    return 0;
}

/* Bytes and chunks of b that are absent from a's chunk set. */
static void shift_cost(const Chunker *c, const uint8_t *a, size_t na,
                       const uint8_t *b, size_t nb, size_t *new_bytes, size_t *new_chunks,
                       size_t *total_chunks)
{
    uint32_t *la, *lb;
    size_t ka = split(c, a, na, &la), kb = split(c, b, nb, &lb);
    Ent *set = malloc(ka * sizeof *set);
    if (!set) exit(1);
    size_t pos = 0;
    for (size_t i = 0; i < ka; i++) {
        set[i] = (Ent){fnv1a64(a + pos, la[i]), a + pos, la[i]};
        pos += la[i];
    }
    qsort(set, ka, sizeof *set, cmp_ent);
    *new_bytes = *new_chunks = 0;
    pos = 0;
    for (size_t i = 0; i < kb; i++) {
        if (!present(set, ka, b + pos, lb[i])) { *new_bytes += lb[i]; (*new_chunks)++; }
        pos += lb[i];
    }
    *total_chunks = kb;
    free(set); free(la); free(lb);
}

static int cmp_u32(const void *a, const void *b)
{
    uint32_t x = *(const uint32_t *)a, y = *(const uint32_t *)b;
    return x < y ? -1 : x > y;
}

static double now_s(void)
{
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return ts.tv_sec + ts.tv_nsec * 1e-9;
}
static int cmp_d(const void *a, const void *b)
{
    double x = *(const double *)a, y = *(const double *)b;
    return x < y ? -1 : x > y;
}

int main(int argc, char **argv)
{
    const size_t N = 64 * MiB;
    uint8_t *a = malloc(N), *b = malloc(N + 256);
    if (!a || !b) return 1;
    rabin_init(0x5eed);
    sm_state = 0xb022;
    for (int i = 0; i < 256; i++) bz_tab[i] = splitmix64();
    sm_state = 0x6ea5;
    for (int i = 0; i < 256; i++) gear_tab[i] = splitmix64();
    sm_state = 1;
    for (size_t i = 0; i < N; i += 8) { uint64_t r = splitmix64(); memcpy(a + i, &r, 8); }
    printf("Rabin polynomial (degree 53): 0x%014llx\n\n", (unsigned long long)rb_pol);

    /* D1: size distribution on random data */
    FILE *hf = argc > 1 ? fopen(argv[1], "w") : NULL;
    printf("D1  chunk sizes on %zu MiB random data (bytes)\n", N / MiB);
    printf("    %-16s %7s %8s %8s %7s %7s %7s %7s  %6s %6s\n", "chunker", "chunks",
           "mean", "stdev", "p5", "p50", "p95", "max", "<4KiB", "=max");
    for (size_t c = 0; c < NCH - 1; c++) {
        uint32_t *len;
        size_t k = split(&CH[c], a, N, &len);
        double sum = 0, sq = 0;
        size_t small = 0, atmax = 0, bins[81] = {0};    /* 1 KiB bins, last is >= 80 KiB */
        for (size_t i = 0; i + 1 < k; i++) {            /* ignore the tail chunk */
            sum += len[i]; sq += (double)len[i] * len[i];
            if (len[i] < 4 * KiB) small++;
            if (len[i] == CH[c].mx) atmax++;
            bins[len[i] / KiB < 80 ? len[i] / KiB : 80]++;
        }
        for (int j = 0; hf && j <= 80; j++)
            fprintf(hf, "%s\t%d\t%zu\n", CH[c].name, j, bins[j]);
        size_t m = k - 1;
        double mean = sum / m, sd = sqrt(sq / m - mean * mean);
        qsort(len, m, sizeof *len, cmp_u32);
        printf("    %-16s %7zu %8.0f %8.0f %7u %7u %7u %7u  %5.1f%% %5.2f%%\n", CH[c].name, m,
               mean, sd, len[m * 5 / 100], len[m / 2], len[m * 95 / 100], len[m - 1],
               100.0 * small / m, 100.0 * atmax / m);
        free(len);
    }
    if (hf) fclose(hf);

    /* D2: boundary shift */
    printf("\nD2  chunks of B missing from A (A = 64 MiB random)\n");
    printf("    %-16s %24s %24s\n", "chunker", "B = A + 1 byte @1000000", "B = A + 100 inserts");
    size_t ins[100];
    sm_state = 99;
    for (int i = 0; i < 100; i++) ins[i] = splitmix64() % N;
    for (int i = 1; i < 100; i++)                      /* insertion sort, tiny */
        for (int j = i; j > 0 && ins[j - 1] > ins[j]; j--) { size_t t = ins[j]; ins[j] = ins[j - 1]; ins[j - 1] = t; }
    uint8_t *b100 = malloc(N + 100);
    if (!b100) return 1;
    {
        size_t src = 0, dst = 0;
        for (int i = 0; i < 100; i++) {
            memcpy(b100 + dst, a + src, ins[i] - src); dst += ins[i] - src; src = ins[i];
            b100[dst++] = (uint8_t)(0xA5 ^ i);
        }
        memcpy(b100 + dst, a + src, N - src);
    }
    memcpy(b, a, 1000000); b[1000000] = 0x5a; memcpy(b + 1000001, a + 1000000, N - 1000000);
    for (size_t c = 1; c < NCH; c++) {
        size_t nb1, nc1, t1, nb2, nc2, t2;
        shift_cost(&CH[c], a, N, b, N + 1, &nb1, &nc1, &t1);
        shift_cost(&CH[c], a, N, b100, N + 100, &nb2, &nc2, &t2);
        printf("    %-16s %9zu B in %2zu/%-5zu    %9zu B in %3zu/%-5zu\n", CH[c].name,
               nb1, nc1, t1, nb2, nc2, t2);
    }
    free(b100);

    /* D3: 8 MiB of zero bytes */
    printf("\nD3  8 MiB of zero bytes\n");
    uint8_t *z = calloc(8 * MiB, 1);
    if (!z) return 1;
    for (size_t c = 1; c < NCH - 1; c++) {
        uint32_t *len;
        size_t k = split(&CH[c], z, 8 * MiB, &len);
        printf("    %-16s %6zu chunks, first chunk %u B\n", CH[c].name, k, len[0]);
        free(len);
    }
    free(z);

    /* D4: relative throughput, median of 5 runs */
    printf("\nD4  throughput on 64 MiB random data, median of 5 runs (MiB/s)\n");
    for (size_t c = 1; c < NCH - 1; c++) {
        double t[5];
        size_t chunks = 0;
        for (int r = 0; r < 5; r++) {
            double t0 = now_s();
            chunks = 0;
            for (size_t pos = 0; pos < N; chunks++)
                pos += CH[c].fn(a + pos, N - pos, CH[c].mn, CH[c].mx);
            t[r] = now_s() - t0;
        }
        qsort(t, 5, sizeof t[0], cmp_d);
        printf("    %-16s %8.0f   (min %.0f, max %.0f, %zu chunks)\n", CH[c].name,
               N / MiB / t[2], N / MiB / t[4], N / MiB / t[0], chunks);
    }
    free(a); free(b);
    return 0;
}