/* The dictionary half of an inverted index: five layouts, one question.
 * Companion to retrieval.html.
 *
 *   cc -O2 -o search search.c && ./search
 *
 * Every variant answers exactly the same query -- "is this 32-bit key in the
 * set?" -- over the same n keys, and every variant except the scan is
 * O(log n) or O(1). What differs is where the bytes are and in what order they
 * are touched, which is the only thing the sweep is trying to isolate:
 *
 *   scan          n/2 sequential reads. The prefetcher's best case.
 *   binary        branchy binary search over the sorted array. One mispredict
 *                 per level until the range fits in cache.
 *   branchless    the same search with the compare folded into the index
 *                 arithmetic, so the pipeline never has to guess.
 *   eytzinger     the same keys in breadth-first order, so the hot top of the
 *                 tree is contiguous instead of spread over the whole array.
 *   btree         a 17-way static B-tree, one 64-byte node per cache line, so
 *                 a level costs one line rather than one line per comparison.
 *   hash          open addressing, load factor 0.5. One probe, but the table
 *                 is twice the size of the sorted array, so it leaves cache
 *                 an octave earlier.
 *
 * "Branchless" is a claim about the emitted code, not about the source. GCC 13
 * and Clang 18 do not agree on whether `base += cond ? half : 0` deserves a
 * cmov, so the sweep is worth running under both; the page that reports it
 * gives each compiler its own column rather than picking one.
 *
 * Two deliberate simplifications, both of which make the comparison *about
 * layout* rather than about anything else. Keys are 4i+1, evenly spaced: no
 * variant here is distribution-sensitive except the hash, and multiply-shift
 * scatters an arithmetic sequence perfectly well. And queries are generated by
 * an inline xorshift rather than read from an array, so the query stream costs
 * no cache of its own; the `control` row measures what that generator costs
 * and it is included, identically, in every other row.
 *
 * Every query hits. A membership test that misses is a different measurement
 * (it is a measurement of the *absence* structure -- the last level, or the
 * probe sequence) and mixing the two would report neither.
 */

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

#define B         16      /* keys per B-tree node: 16 x 4 bytes = one line */
#define QUERIES   (1 << 20)
#define TRIALS    7
#define INF       0xffffffffu
#define ROUND(b)  (((b) + 63) & ~(size_t)63)

static double now(void)
{
    struct timespec t;
    clock_gettime(CLOCK_MONOTONIC, &t);
    return t.tv_sec + t.tv_nsec * 1e-9;
}

static unsigned n, mask;
static unsigned *sorted;      /* sorted[i] = 4i+1                            */
static unsigned *eyt;         /* 1-indexed, breadth-first                    */
static unsigned *btree;       /* nblocks nodes of B keys, 17-ary heap order  */
static unsigned nblocks;
static unsigned *table;       /* open addressing, 0 = empty                  */
static unsigned tmask, tcap;

/* ------------------------------------------------------------------ scan */
static inline int f_scan(unsigned x)
{
    for (unsigned i = 0; i < n; i++)
        if (sorted[i] == x) return 1;
    return 0;
}

/* ---------------------------------------------------------------- binary */
static inline int f_binary(unsigned x)
{
    unsigned lo = 0, hi = n;
    while (lo < hi) {
        unsigned mid = lo + (hi - lo) / 2;
        if (sorted[mid] < x) lo = mid + 1;
        else hi = mid;
    }
    return lo < n && sorted[lo] == x;
}

/* The compare becomes an addend rather than a jump. The candidate range is
   halved unconditionally, so the loop trip count depends only on n and the
   branch predictor has nothing left to get wrong. */
static inline int f_branchless(unsigned x)
{
    const unsigned *base = sorted;
    unsigned len = n;
    while (len > 1) {
        unsigned half = len >> 1;
        base += (base[half - 1] < x) ? half : 0;
        len -= half;
    }
    return base[0] == x;
}

/* ------------------------------------------------------------- eytzinger */
static unsigned et;

static void build_eyt(unsigned k)
{
    if (k <= n) {
        build_eyt(2 * k);
        eyt[k] = sorted[et++];
        build_eyt(2 * k + 1);
    }
}

/* No equality test inside the loop: the descent is unconditional and the
   position of the answer is recovered afterwards from the trailing ones of
   the final index, which is where the search "turned left" last. */
static inline int f_eyt(unsigned x)
{
    unsigned k = 1;
    while (k <= n) k = 2 * k + (eyt[k] < x);
    k >>= __builtin_ffs((int)~k);
    return k != 0 && eyt[k] == x;
}

/* ---------------------------------------------------------------- btree */
static unsigned bt;

static unsigned go(unsigned k, unsigned i) { return k * (B + 1) + i + 1; }

static void build_btree(unsigned k)
{
    if (k < nblocks) {
        for (unsigned i = 0; i < B; i++) {
            build_btree(go(k, i));
            btree[(size_t)k * B + i] = bt < n ? sorted[bt++] : INF;
        }
        build_btree(go(k, B));
    }
}

/* Rank within one node, counted rather than searched: sixteen independent
   compares over one cache line, with no branch to mispredict. */
static inline unsigned rank(const unsigned *node, unsigned x)
{
    unsigned r = 0;
    for (unsigned i = 0; i < B; i++) r += node[i] < x;
    return r;
}

static inline int f_btree(unsigned x)
{
    unsigned k = 0, res = INF;
    while (k < nblocks) {
        unsigned i = rank(btree + (size_t)k * B, x);
        if (i < B) res = btree[(size_t)k * B + i];
        k = go(k, i);
    }
    return res == x;
}

/* ----------------------------------------------------------------- hash */
static inline int f_hash(unsigned x)
{
    unsigned h = (unsigned)(x * 2654435761u) & tmask;
    while (table[h]) {
        if (table[h] == x) return 1;
        h = (h + 1) & tmask;
    }
    return 0;
}

/* ---------------------------------------------------------------- timing
 *
 * One timed loop per variant, written out by a macro rather than reached
 * through a function pointer. An indirect call would be the same few
 * nanoseconds in every row, but it would also stop the compiler inlining the
 * probe, and at n = 16 the whole answer is a few nanoseconds: a measurement
 * whose overhead is the size of its subject is not a measurement.
 */
static volatile unsigned long long sink;

#define SEED 0x1234567u

/* Every name the macro introduces is suffixed, including the repeat count:
   the scan is called with a repeat count the caller computed into a variable,
   and an unsuffixed `reps` here would shadow it inside its own initialiser. */
#define TIME(DEST, REPS, EXPR)                                                \
    do {                                                                      \
        double best_ = 1e30;                                                  \
        long reps_ = (REPS);                                                  \
        for (int t_ = 0; t_ < TRIALS; t_++) {                                 \
            unsigned r_ = SEED;                                               \
            unsigned long long hits_ = 0;                                     \
            double a_ = now();                                                \
            for (long q_ = 0; q_ < reps_; q_++) {                             \
                r_ ^= r_ << 13; r_ ^= r_ >> 17; r_ ^= r_ << 5;                \
                unsigned x = 4 * (r_ & mask) + 1;                             \
                hits_ += (EXPR);                                              \
            }                                                                 \
            double b_ = now();                                                \
            sink += hits_;                                                    \
            double ns_ = (b_ - a_) / (double)reps_ * 1e9;                     \
            if (ns_ < best_) best_ = ns_;                                     \
        }                                                                     \
        DEST = best_;                                                         \
    } while (0)

int main(void)
{
    printf("%10s %10s %10s %10s %10s %10s %10s %10s\n",
           "n", "control", "scan", "binary", "branchless", "eytzinger",
           "btree", "hash");
    for (unsigned lg = 4; lg <= 24; lg++) {
        n = 1u << lg;
        mask = n - 1;
        nblocks = (n + B - 1) / B;
        tcap = 1u << (lg + 1);
        tmask = tcap - 1;

        /* aligned_alloc requires a size that is a multiple of the alignment,
           so every request is rounded up to a whole cache line. */
        sorted = aligned_alloc(64, ROUND((size_t)n * 4));
        eyt = aligned_alloc(64, ROUND(((size_t)n + 1) * 4));
        btree = aligned_alloc(64, ROUND((size_t)nblocks * B * 4));
        table = aligned_alloc(64, ROUND((size_t)tcap * 4));
        if (!sorted || !eyt || !btree || !table) { fprintf(stderr, "oom\n"); return 1; }

        for (unsigned i = 0; i < n; i++) sorted[i] = 4 * i + 1;
        et = 0; build_eyt(1);
        bt = 0; memset(btree, 0xff, (size_t)nblocks * B * 4); build_btree(0);
        memset(table, 0, (size_t)tcap * 4);
        for (unsigned i = 0; i < n; i++) {
            unsigned x = sorted[i], h = (unsigned)(x * 2654435761u) & tmask;
            while (table[h]) h = (h + 1) & tmask;
            table[h] = x;
        }

        /* Correctness before speed: a fast wrong answer is not a result. The
           scan is checked only at the sizes where it is also timed: it is
           O(n) per probe, so checking it at 16 Mi keys would cost more than
           the whole rest of the sweep. */
        int scan_here = n <= (1u << 16);
        for (unsigned i = 0; i < n; i += (n / 97) + 1) {
            unsigned x = sorted[i];
            if (!f_binary(x) || !f_branchless(x) || !f_eyt(x) || !f_btree(x) || !f_hash(x)
                || (scan_here && !f_scan(x))) {
                fprintf(stderr, "miss on present key %u at n=%u\n", x, n); return 1;
            }
            if (f_binary(x + 1) || f_branchless(x + 1) || f_eyt(x + 1)
                || f_btree(x + 1) || f_hash(x + 1) || (scan_here && f_scan(x + 1))) {
                fprintf(stderr, "hit on absent key %u at n=%u\n", x + 1, n); return 1;
            }
        }

        double c, b, l, e, t, h, s = -1;
        TIME(c, QUERIES, x);   /* the query generator, alone */
        /* The scan costs O(n) per query, so a fixed query count would make it
           the whole runtime of the sweep. It gets a budget proportional to
           1/n instead, floored so the timer still has something to measure;
           past 64Ki keys its shape is established and its value is an order
           of magnitude off the chart. */
        if (n <= (1u << 16)) {
            long reps = QUERIES / (long)n * 64;
            TIME(s, reps < 512 ? 512 : reps, f_scan(x));
        }
        TIME(b, QUERIES, f_binary(x));
        TIME(l, QUERIES, f_branchless(x));
        TIME(e, QUERIES, f_eyt(x));
        TIME(t, QUERIES, f_btree(x));
        TIME(h, QUERIES, f_hash(x));
        printf("%10u %10.2f", n, c);
        if (s < 0) printf(" %10s", "-"); else printf(" %10.1f", s);
        printf(" %10.2f %10.2f %10.2f %10.2f %10.2f\n", b, l, e, t, h);
        fflush(stdout);

        free(sorted); free(eyt); free(btree); free(table);
    }
    return 0;
}