/* The unit of transfer is a cache line, not a byte. Companion to retrieval.html.
 *
 *   cc -O2 -o stride stride.c && ./stride
 *
 * Three measurements over one 64 MiB array of 32-bit integers -- a working set
 * far larger than any cache here, so every part of it must come from DRAM at
 * least once:
 *
 *   1. A stride sweep. The loop touches every Nth element and nothing else, so
 *      the number of elements read falls by N while the number of *lines*
 *      touched stays constant until N reaches 16 (16 x 4 bytes = one 64-byte
 *      line). Time per element touched should therefore rise roughly linearly
 *      to N = 16 and flatten after it: past one element per line, a larger
 *      stride buys nothing because the line was already the smallest thing the
 *      machine would fetch.
 *
 *   2. Sequential versus random over the identical footprint. Same array, same
 *      element count, same arithmetic; only the order changes. The difference
 *      is entirely locality and prefetching.
 *
 *   3. The same random order, but with the loads made independent of each
 *      other (index read from a precomputed array rather than chained). This
 *      separates *latency* from *throughput*: the core can keep many
 *      independent misses in flight, so a random gather is far quicker than a
 *      random pointer chase over the same addresses.
 *
 * All three report nanoseconds per element touched, and (1) additionally
 * reports nanoseconds per 64-byte line, which is the figure that stays flat.
 */

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

#define MIB       64
#define N         ((size_t)MIB * 1024 * 1024 / sizeof(unsigned))
#define TRIALS    7
#define LINE      64

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

static unsigned long long rng_state = 0x243f6a8885a308d3ull;

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);
}

static volatile unsigned long long sink;

int main(void)
{
    unsigned *a = malloc(N * sizeof *a);
    unsigned *idx = malloc((N / 16) * sizeof *idx);
    if (!a || !idx) { fprintf(stderr, "out of memory\n"); return 1; }
    for (size_t i = 0; i < N; i++) a[i] = (unsigned)i;

    printf("# array %d MiB, %zu elements of 4 bytes\n", MIB, N);
    printf("%8s %12s %12s %12s\n", "stride", "elements", "ns/element", "ns/line");
    for (size_t s = 1; s <= 64; s *= 2) {
        size_t touched = N / s;
        double best = 1e30;
        for (int t = 0; t < TRIALS; t++) {
            unsigned long long acc = 0;
            double t0 = now();
            for (size_t i = 0; i < N; i += s) acc += a[i];
            double t1 = now();
            sink = acc;
            double ns = (t1 - t0) / (double)touched * 1e9;
            if (ns < best) best = ns;
        }
        /* Lines touched: one per stride step until the stride exceeds a line,
           after which it is one line per element. */
        double per_line = s < LINE / sizeof(unsigned)
                        ? best * (LINE / sizeof(unsigned)) / (double)s : best;
        printf("%8zu %12zu %12.3f %12.3f\n", s, touched, best, per_line);
    }

    /* Sequential, random-dependent and random-independent over one sixteenth
       of the array's elements -- one per cache line, so all three touch the
       same 64 MiB and the same one million lines, and only the order and the
       dependency structure differ. */
    size_t m = N / 16;
    for (size_t i = 0; i < m; i++) idx[i] = (unsigned)(i * 16);
    for (size_t i = m - 1; i > 0; i--) {
        size_t j = (size_t)(rnd() % (i + 1));
        unsigned t = idx[i]; idx[i] = idx[j]; idx[j] = t;
    }
    /* The dependent walk needs a cycle laid down in the array itself. It
       overwrites a[], which is why the strided sweep above runs first. */
    for (size_t i = 0; i < m; i++) a[idx[i]] = idx[(i + 1) % m];

    double seq = 1e30, dep = 1e30, ind = 1e30;
    for (int t = 0; t < TRIALS; t++) {
        unsigned long long acc = 0;
        double t0 = now();
        for (size_t i = 0; i < N; i += 16) acc += a[i];
        double t1 = now();
        unsigned p = idx[0];
        for (size_t i = 0; i < m; i++) p = a[p];
        double t2 = now();
        unsigned long long acc2 = 0;
        for (size_t i = 0; i < m; i++) acc2 += a[idx[i]];
        double t3 = now();
        sink = acc + p + acc2;
        double s1 = (t1 - t0) / (double)m * 1e9;
        double s2 = (t2 - t1) / (double)m * 1e9;
        double s3 = (t3 - t2) / (double)m * 1e9;
        if (s1 < seq) seq = s1;
        if (s2 < dep) dep = s2;
        if (s3 < ind) ind = s3;
    }
    printf("\n%28s %12s %10s\n", "order", "ns/element", "vs seq");
    printf("%28s %12.3f %9.1fx\n", "sequential", seq, seq / seq);
    printf("%28s %12.3f %9.1fx\n", "random, independent loads", ind, ind / seq);
    printf("%28s %12.3f %9.1fx\n", "random, dependent chase", dep, dep / seq);
    free(idx);
    free(a);
    return 0;
}