/* The memory hierarchy, measured from C. Companion to retrieval.html.
 *
 *   cc -O2 -o latency latency.c && ./latency
 *
 * One dependent load at a time. The array is filled with a single random
 * permutation cycle, so `p = next[p]` cannot issue until the previous load has
 * returned: the loop measures load-to-use latency and nothing else. That is
 * the whole method, and every part of it is defensive.
 *
 *   - A *cycle*, not a random index per step, because an independent index
 *     stream lets the core keep a dozen loads in flight and reports memory
 *     level parallelism instead of latency.
 *   - A permutation, so every element is visited exactly once per pass and the
 *     working set is exactly the array. A random walk revisits and the
 *     effective footprint shrinks.
 *   - Stride at least one cache line between successive elements is NOT
 *     enforced; the permutation is over line-sized slots, so each hop lands on
 *     a distinct line and the line count equals the slot count.
 *   - The hardware prefetcher cannot follow a pointer chase, which is exactly
 *     why this shape is the standard one: it defeats the prefetcher rather
 *     than measuring it.
 *
 * The reported figure is the minimum over TRIALS passes. Noise on a shared
 * machine only ever adds time, so the minimum is the least contaminated
 * estimate rather than a cherry-picked best case.
 *
 * Run it twice:
 *
 *   ./latency          4 KiB pages
 *   ./latency huge     the same buffer under MADV_HUGEPAGE
 *
 * The pair is the point. A single curve cannot tell a cache capacity from a
 * TLB reach, because on 4 KiB pages the second-level TLB runs out first and
 * the resulting cliff looks exactly like a cache one. Re-running with 2 MiB
 * pages moves that cliff and leaves the real cache edges where they were.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/mman.h>
#include <time.h>

#define LINE     64          /* coherency_line_size on every target here */
#define TRIALS    5
#define MIN_STEPS 2000000L   /* per pass, independent of the working set */

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

/* splitmix64: a fixed generator so the permutation, and therefore the whole
   measurement, is reproducible from the source alone. */
static unsigned long long rng_state = 0x9e3779b97f4a7c15ull;

static unsigned long long rnd(void)
{
    unsigned long long z = (rng_state += 0x9e3779b97f4a7c15ull);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ull;
    z = (z ^ (z >> 27)) * 0x94d049bb133111ebull;
    return z ^ (z >> 31);
}

/* A single cycle through all n slots, built by shuffling 0..n-1 and then
   linking each element to its successor in the shuffled order. Sattolo's
   algorithm would give a cycle directly; this is the same thing said in two
   steps, and it keeps the "one cycle, every slot once" property explicit. */
static void build_cycle(size_t *next, size_t n, size_t *order)
{
    for (size_t i = 0; i < n; i++) order[i] = i;
    for (size_t i = n - 1; i > 0; i--) {
        size_t j = (size_t)(rnd() % (i + 1));
        size_t t = order[i]; order[i] = order[j]; order[j] = t;
    }
    for (size_t i = 0; i < n; i++)
        next[order[i] * (LINE / sizeof(size_t))] =
            order[(i + 1) % n] * (LINE / sizeof(size_t));
}

static volatile size_t sink;

/* What the kernel actually did with the advice. MADV_HUGEPAGE is a request,
   and a run that silently got 4 KiB pages anyway would look like a result. */
static void report_pages(void)
{
    FILE *f = fopen("/proc/self/smaps_rollup", "r");
    char line[256];
    if (!f) return;
    while (fgets(line, sizeof line, f))
        if (strncmp(line, "AnonHugePages:", 14) == 0) printf("# %s", line);
    fclose(f);
}

int main(int argc, char **argv)
{
    int huge = argc > 1 && strcmp(argv[1], "huge") == 0;
    printf("# pages: %s\n", huge ? "2 MiB (MADV_HUGEPAGE)" : "4 KiB (default)");
    printf("%10s %12s %10s %8s\n", "KiB", "lines", "ns/load", "trials");
    /* 4 KiB to 512 MiB. The interesting sizes are the ones straddling a cache
       capacity, so the sweep is by powers of two with a midpoint between each
       pair: a plateau edge that falls between 2^k and 2^(k+1) is otherwise
       reported as a single step of unknown shape. */
    static const size_t kib[] = {
        4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768,
        1024, 1536, 2048, 3072, 4096, 6144, 8192, 12288, 16384, 24576, 32768,
        49152, 65536, 98304, 131072, 262144, 524288,
    };

    size_t max_lines = kib[sizeof kib / sizeof *kib - 1] * 1024 / LINE;
    size_t bytes = max_lines * LINE;

    /* mmap rather than malloc, so the page size is this program's decision and
       not the allocator's. MADV_HUGEPAGE is advice: the kernel is free to
       refuse it, and the run prints AnonHugePages afterwards so a refusal is
       visible rather than assumed. */
    size_t *buf = mmap(NULL, bytes, PROT_READ | PROT_WRITE,
                       MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
    if (buf == MAP_FAILED) { fprintf(stderr, "mmap failed\n"); return 1; }
    if (huge && madvise(buf, bytes, MADV_HUGEPAGE) != 0)
        fprintf(stderr, "# MADV_HUGEPAGE refused\n");
    size_t *order = malloc(max_lines * sizeof *order);
    if (!order) { fprintf(stderr, "out of memory\n"); return 1; }

    for (size_t k = 0; k < sizeof kib / sizeof *kib; k++) {
        size_t lines = kib[k] * 1024 / LINE;
        rng_state = 0x9e3779b97f4a7c15ull;      /* same permutation shape per size */
        memset(buf, 0, lines * LINE);
        build_cycle(buf, lines, order);

        /* Enough steps that the timer resolution is irrelevant, and at least
           a few full laps of the cycle so a partial lap cannot bias a size. */
        long steps = MIN_STEPS;
        if ((size_t)steps < lines * 4) steps = (long)lines * 4;

        double best = 1e30;
        for (int t = 0; t < TRIALS; t++) {
            size_t p = 0;
            /* One untimed lap: the measurement is of a warm TLB and a settled
               page table, not of first-touch page faults. */
            for (size_t i = 0; i < lines; i++) p = buf[p];
            double t0 = now();
            for (long i = 0; i < steps; i++) p = buf[p];
            double t1 = now();
            sink = p;
            double ns = (t1 - t0) / (double)steps * 1e9;
            if (ns < best) best = ns;
        }
        printf("%10zu %12zu %10.2f %8d\n", kib[k], lines, best, TRIALS);
        fflush(stdout);
    }
    report_pages();
    free(order);
    munmap(buf, bytes);
    return 0;
}