Advanced Examples

// learning/code/ex-58-optimize-kernel-end-to-end/kernel.c
/* Example 58: profile a slow kernel, fix its layout, re-measure. */
#include <stdio.h>  // => printf: report the profile-fix-reprofile transcript
#include <stdlib.h> // => malloc/free/rand: build and free the synthetic dataset
#include <string.h> // => memset: fill the "cold" bytes so the record is realistic, not zeroed noise
#include <time.h>   // => clock_gettime: the wall-clock stopwatch this whole example is built on
 
#define N \
    4000000 // => co-25: 4M records -- large enough that layout dominates over
            // instruction count
#define TRIALS \
    5 // => co-25: best-of-N wall-clock timing, per this topic's measurement
      // methodology
 
// ex-58: a deliberately realistic "record" -- one HOT field (grade) the kernel
// actually needs, wrapped in COLD fields (a padded name buffer and metadata)
// nobody touches in this kernel -- exactly the shape a real "student record" or
// "log event" struct takes co-17: this struct's shape, not the loop's
// instruction count, is what ex-58 optimizes
typedef struct {       // struct layout definition
    int grade;         // => co-17: the ONE field the average-grade kernel actually reads
    char name[124];    // => co-17: COLD -- padded to make the record one full 128 B
                       // cache line
    int enrolled_year; // => co-17: COLD -- never read by this kernel, still paid
                       // for every cache line
} StudentRecord;       // => sizeof == 132 bytes: crosses more than one 128 B cache
                       // line per record
 
// ex-58: the ORIGINAL kernel -- walks the array-of-structs, reading `grade` out
// of every full 132-byte record even though 128 of those bytes are dead weight
// for this computation
static double average_grade_aos(const StudentRecord *records,
                                int n) { // => co-17: cache-hostile pass
    long sum = 0;                        // => accumulate as long: N * max_grade never overflows int
    for (int i = 0; i < n; i++) {        // => co-17: one full 132 B record loaded per iteration...
        sum += records[i].grade;         // => ...to read 4 of those 132 bytes -- the rest
                                         // is wasted traffic
    }
    return (double)sum / n; // => co-25: the "profile" step's number-to-explain
}
 
// ex-58: the FIXED kernel -- given a pre-extracted compact `int grade[]` array
// (SoA for just the field this computation needs), the exact same arithmetic
// runs over dense, fully cache-line-packed data: co-02 -- every 128 B line now
// holds 32 real `grade` values
static double average_grade_soa(const int *grade,
                                int n) { // => co-17: cache-friendly pass
    long sum = 0;                        // => same accumulator type, same arithmetic -- ONLY layout changed
    for (int i = 0; i < n; i++) {        // => co-03: sequential dense reads, 32 useful ints per line
        sum += grade[i];                 // updates sum
    }
    return (double)sum / n; // => must equal average_grade_aos's result bit-for-bit (both are
} // integer sums divided once -- no floating-point reordering risk)
 
static double best_of_seconds(double (*fn)(const void *, int), const void *data,
                              int n) {       // declares function pointer fn
    double best = -1.0;                      // => co-25: best-of-N methodology -- report the fastest clean run,
    for (int t = 0; t < TRIALS; t++) {       // not a single noisy sample, per this topic's measurement rule
        struct timespec t0, t1;              // supporting statement for this example
        clock_gettime(CLOCK_MONOTONIC, &t0); // calls clock_gettime(...)
        volatile double result = fn(data,
                                    n);                                          // => volatile: stop the optimizer from deleting the "unused" call
        clock_gettime(CLOCK_MONOTONIC, &t1);                                     // calls clock_gettime(...)
        (void)result;                                                            // => silence "unused variable" -- the volatile read is the
                                                                                 // real guard
        double secs = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9; // declares secs
        if (best < 0 || secs < best)
            best = secs; // => keep the fastest of TRIALS runs
    }
    return best; // returns the computed result
}
 
// ex-58: thin typed wrappers so best_of_seconds can call either kernel through
// one function-pointer signature -- co-25: this is scaffolding for the
// STOPWATCH, not the kernel itself
static double wrap_aos(const void *p,
                       int n) {                        // defines wrap_aos(): helper function used by this example
    const StudentRecord *r = (const StudentRecord *)p; // declares r
    return average_grade_aos(r, n);                    // returns the computed result
}
static double wrap_soa(const void *p,
                       int n) {     // defines wrap_soa(): helper function used by this example
    const int *g = (const int *)p;  // declares g
    return average_grade_soa(g, n); // returns the computed result
}
 
int main(void) { // program entry point
    // Step 1 ("profile a slow kernel"): build the realistic AoS dataset a real
    // system would hand you.
    StudentRecord *records = malloc(sizeof(StudentRecord) * N); // => co-17: 132 B * 4M ~= 528 MB AoS
    if (!records) {
        fprintf(stderr, "alloc failed\n");
        return 1;
    } // prints a report line
    unsigned seed = 12345u;                    // => co-25: fixed seed -- reproducible dataset every run
    for (int i = 0; i < N; i++) {              // loop header controlling the sweep below
        seed = seed * 1103515245u + 12345u;    // => a tiny deterministic LCG, no <stdlib.h> rand() state
        records[i].grade = (int)(seed % 101u); // grades 0..100
        memset(records[i].name, 'x',
               sizeof(records[i].name));                    // => cold bytes: realistic non-zero payload
        records[i].enrolled_year = 2020 + (int)(seed % 6u); // cold field, never read by the kernel
    }
 
    double aos_secs = best_of_seconds(wrap_aos, records, N); // => co-25: profile step -- measure the slow path
    double aos_avg = average_grade_aos(records, N);          // => capture the actual numeric result too
 
    // Step 2 ("fix its layout"): extract just the hot field into a dense SoA
    // buffer.
    int *grade = malloc(sizeof(int) * N); // => co-17: 4 B * 4M = 16 MB -- 33x smaller footprint than AoS
    if (!grade) {
        fprintf(stderr, "alloc failed\n");
        return 1;
    } // prints a report line
    for (int i = 0; i < N; i++)
        grade[i] = records[i].grade; // => the one-time extraction cost (not timed)
 
    // Step 3 ("re-measure"): profile the fixed kernel with the identical
    // stopwatch.
    double soa_secs = best_of_seconds(wrap_soa, grade, N); // => co-25: re-profile step -- measure the fix
    double soa_avg = average_grade_soa(grade, N);          // declares soa_avg
 
    printf("dataset: N=%d records, sizeof(StudentRecord)=%zu bytes (%.1f MB "
           "total AoS)\n", // prints a report line
           N, sizeof(StudentRecord),
           (double)(sizeof(StudentRecord) * N) / (1024.0 * 1024.0));                                    // continues the printf(...) call above
    printf("step 1 (profile)  AoS average_grade: best of %d = %.4f s, result=%.6f\n", TRIALS, aos_secs, // prints a report line
           aos_avg);                                                                                    // continues the printf(...) call above
    printf("step 3 (reprofile) SoA average_grade: best of %d = %.4f s, "
           "result=%.6f\n",
           TRIALS, soa_secs, // prints a report line
           soa_avg);         // continues the printf(...) call above
    printf("results identical: %s (both %.10f)\n", (aos_avg == soa_avg) ? "yes" : "NO -- BUG",
           aos_avg);                      // prints a report line
    double speedup = aos_secs / soa_secs; // declares speedup
    printf("speedup: %.2fx  (PASS: SoA faster with identical result) -> %s\n",
           speedup,                                                   // prints a report line
           (aos_avg == soa_avg && speedup > 1.05) ? "PASS" : "FAIL"); // continues the printf(...) call above
 
    free(records); // releases records's heap memory
    free(grade);   // releases grade's heap memory
    return 0;      // returns the computed result
}

dataset: N=4000000 records, sizeof(StudentRecord)=132 bytes (503.5 MB total AoS)
step 1 (profile)  AoS average_grade: best of 5 = 0.0129 s, result=49.972522
step 3 (reprofile) SoA average_grade: best of 5 = 0.0003 s, result=49.972522
results identical: yes (both 49.9725222500)
speedup: 44.20x  (PASS: SoA faster with identical result) -> PASS

// learning/code/ex-59-blocked-transpose/transpose.c
/* Example 59: naive vs cache-blocked matrix transpose. */
#include <stdio.h>  // stdio.h: standard library header
#include <stdlib.h> // stdlib.h: standard library header
#include <time.h>   // time.h: standard library header
 
#define N \
    2048 // => co-03: 2048x2048 int matrix = 16 MB -- bigger than this machine's 4
         // MiB L2
#define BLOCK \
    32 // => co-03: 32x32 int tile = 4 KB -- comfortably fits alongside other
       // tiles in L1d
 
// ex-59: naive transpose -- reads `src` row-major (sequential, cache-friendly)
// but WRITES `dst` column-major (stride N*4 bytes between consecutive writes --
// co-03: every write is a fresh cache line, and for N=2048 the write stream
// revisits each line only once every N iterations, so nothing stays resident
// long enough to be reused)
static void transpose_naive(const int *src, int *dst,
                            int n) {         // defines transpose_naive(): helper function used by this example
    for (int i = 0; i < n; i++) {            // => co-03: row-major read of src -- good locality
        for (int j = 0; j < n; j++) {        // loop header controlling the sweep below
            dst[j * n + i] = src[i * n + j]; // => co-03: dst write strides by n*4 bytes -- bad locality
        }
    }
}
 
// ex-59: cache-blocked transpose -- processes the matrix in BLOCK x BLOCK
// tiles, so both the read tile and the write tile stay resident in L1d for the
// whole inner block instead of streaming through the full 16 MB matrix
// column-by-column
static void transpose_blocked(const int *src, int *dst, int n,
                              int block) {               // defines transpose_blocked(): helper function
                                                         // used by this example
    for (int ii = 0; ii < n; ii += block) {              // => co-03: outer loop over tile rows
        for (int jj = 0; jj < n; jj += block) {          // => co-03: outer loop over tile columns
            int i_max = ii + block < n ? ii + block : n; // => clamp the last tile if n isn't a multiple
            int j_max = jj + block < n ? jj + block : n; // declares j_max
            for (int i = ii; i < i_max; i++) {           // => co-03: inner loops stay INSIDE one tile --
                for (int j = jj; j < j_max; j++) {       // both src and dst touches land in the same
                    dst[j * n + i] = src[i * n + j];     // small hot region for the whole tile's work
                }
            }
        }
    }
}
 
// ex-59: a shared timing harness so BOTH kernels are measured through the
// identical best-of-N loop -- taking the minimum of several trials, not the
// mean, is what filters out one-off OS scheduling noise without also hiding a
// real, reproducible difference between the two kernels.
static double best_of(void (*fn)(const int *, int *, int, int), const int *src, int *dst, int n, // declares function pointer fn
                      int block, int trials) {                                                   // declares block
    double best = -1.0;                                                                          // declares best
    for (int t = 0; t < trials; t++) {                                                           // loop header controlling the sweep below
        struct timespec t0, t1;                                                                  // supporting statement for this example
        clock_gettime(CLOCK_MONOTONIC, &t0);                                                     // calls clock_gettime(...)
        fn(src, dst, n, block);                                                                  // calls fn(...)
        clock_gettime(CLOCK_MONOTONIC, &t1);                                                     // calls clock_gettime(...)
        double secs = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9;                 // declares secs
        if (best < 0 || secs < best)
            best = secs; // conditional check
    }
    return best; // returns the computed result
}
 
// ex-59: best_of() takes a function pointer matching transpose_blocked's
// 4-argument signature, but transpose_naive only takes 3 -- this thin wrapper
// adapts the naive kernel to that same signature so both kernels can be driven
// through ONE shared timing loop instead of two.
static void naive_wrap(const int *src, int *dst, int n,
                       int block) { // defines naive_wrap(): helper function used by this example
    (void)block;                    // => naive ignores the block-size argument -- unused on purpose
    transpose_naive(src, dst, n);   // calls transpose_naive(...)
}
 
// ex-59: main() allocates the three N x N buffers once, times each kernel via
// best_of(), then checks TWO independent correctness signals -- a full
// element-by-element diff between the two kernels' outputs, and a sampled
// direct check of the transpose property against the source -- before ever
// looking at which kernel ran faster, so a wrong-but-fast kernel can't sneak a
// PASS.
int main(void) {                                            // program entry point
    int *src = malloc(sizeof(int) * (size_t)N * N);         // heap-allocates memory for src
    int *dst_naive = malloc(sizeof(int) * (size_t)N * N);   // heap-allocates memory for dst_naive
    int *dst_blocked = malloc(sizeof(int) * (size_t)N * N); // heap-allocates memory for dst_blocked
    if (!src || !dst_naive || !dst_blocked) {
        fprintf(stderr, "alloc failed\n");
        return 1;
    } // prints a report line
 
    for (int i = 0; i < N * N; i++)
        src[i] = i; // => deterministic content: src[i*N+j] == i*N+j
 
    double naive_secs = best_of(naive_wrap, src, dst_naive, N, BLOCK, 3); // declares naive_secs
    double blocked_secs = best_of(transpose_blocked, src, dst_blocked, N, BLOCK,
                                  3); // declares blocked_secs
 
    int mismatches = 0;               // => co-03: correctness check -- both outputs must be IDENTICAL
    for (int i = 0; i < N * N; i++) { // loop header controlling the sweep below
        if (dst_naive[i] != dst_blocked[i])
            mismatches++; // conditional check
    }
    // spot-check the transpose property directly: dst[j*N+i] must equal
    // src[i*N+j]
    int property_ok = 1;                   // declares property_ok
    for (int i = 0; i < N; i += 511) {     // sample a handful of rows, not all 2048 (keeps output short)
        for (int j = 0; j < N; j += 511) { // loop header controlling the sweep below
            if (dst_blocked[j * N + i] != src[i * N + j])
                property_ok = 0; // conditional check
        }
    }
 
    printf("matrix: %dx%d int (%.1f MB), block=%d (%.1f KB per tile)\n", N,
           N, // prints a report line
           (double)(sizeof(int) * (size_t)N * N) / (1024.0 * 1024.0),
           BLOCK,                                           // continues the printf(...) call above
           (double)(sizeof(int) * BLOCK * BLOCK) / 1024.0); // continues the printf(...) call above
    printf("naive transpose:   best of 3 = %.4f s\n",
           naive_secs); // prints a report line
    printf("blocked transpose: best of 3 = %.4f s\n",
           blocked_secs); // prints a report line
    printf("outputs identical: %s (%d mismatches out of %d)\n",
           mismatches == 0 ? "yes" : "NO -- BUG",                            // prints a report line
           mismatches, N * N);                                               // continues the printf(...) call above
    printf("transpose property dst[j][i]==src[i][j] sampled-verified: %s\n", // prints
                                                                             // a
                                                                             // report
                                                                             // line
           property_ok ? "yes" : "NO -- BUG");                               // continues the printf(...) call above
    double speedup = naive_secs / blocked_secs;                              // declares speedup
    printf("speedup: %.2fx (PASS: blocked faster + correct) -> %s\n",
           speedup,                                                               // prints a report line
           (mismatches == 0 && property_ok && speedup > 1.05) ? "PASS" : "FAIL"); // continues the printf(...) call above
 
    free(src);         // releases src's heap memory
    free(dst_naive);   // releases dst_naive's heap memory
    free(dst_blocked); // releases dst_blocked's heap memory
    return 0;          // returns the computed result
}

matrix: 2048x2048 int (16.0 MB), block=32 (4.0 KB per tile)
naive transpose:   best of 3 = 0.0259 s
blocked transpose: best of 3 = 0.0186 s
outputs identical: yes (0 mismatches out of 4194304)
transpose property dst[j][i]==src[i][j] sampled-verified: yes
speedup: 1.39x (PASS: blocked faster + correct) -> PASS

// learning/code/ex-60-particle-sim-soa-simd/particles.c
/* Example 60: particle-system position update -- AoS scalar vs SoA scalar vs
 * SoA+NEON. */
#include <arm_neon.h> // => co-23: NEON intrinsics -- this machine is arm64, so NEON is the real SIMD ISA
#include <stdio.h>    // stdio.h: standard library header
#include <stdlib.h>   // stdlib.h: standard library header
#include <time.h>     // time.h: standard library header
 
#define N 4000000 // => co-17: 4M particles
#define STEPS \
    20           // => repeat the update this many times per timed run (amplifies the
                 // signal)
#define DT 0.01f // constant DT = 0.01f
 
// ex-60: AoS particle -- position and velocity interleaved per particle, plus
// one cold "mass" field never touched by this position-update kernel (co-17:
// realistic AoS shape)
typedef struct {      // struct layout definition
    float x, y, z;    // => position -- updated every step
    float vx, vy, vz; // => velocity -- read every step, never written here
    float mass;       // => co-17: cold field -- part of the record, irrelevant to this
                      // kernel
} ParticleAoS;        // supporting statement for this example
 
// clang's auto-vectorizer can still find some parallelism in AoS at -O2;
// disabling it here keeps this a clean apples-to-apples "scalar work, different
// layout" baseline (stated in prose).
static void update_aos(ParticleAoS *p, int n,
                       float dt) {                        // defines update_aos(): helper function used by this example
#pragma clang loop vectorize(disable) interleave(disable) // supporting statement for this example
    for (int i = 0; i < n; i++) {                         // => co-17: each iteration loads a full 28 B particle
        p[i].x += p[i].vx * dt;                           // to update 3 of its 7 floats -- the rest rides along
        p[i].y += p[i].vy * dt;                           // updates p[i].y
        p[i].z += p[i].vz * dt;                           // updates p[i].z
    }
}
 
// ex-60: SoA particle system -- position and velocity as six separate dense
// arrays. co-17: only the arrays this kernel actually touches (x/y/z, vx/vy/vz)
// are ever loaded -- there is no "mass" array to pay for in this loop at all.
typedef struct {                     // struct layout definition
    float *x, *y, *z, *vx, *vy, *vz; // declares x
} ParticleSoA;                       // supporting statement for this example
 
// Also forced non-vectorized (see note above ex ex-60's AoS version): this
// isolates the explicit-NEON win below from clang's own auto-vectorizer, which
// would otherwise blur the AoS-vs-SoA-vs-explicit-SIMD comparison this example
// is actually about. ex-47 already covers auto-vectorization on its own.
static void update_soa_scalar(ParticleSoA *p, int n,
                              float dt) {                 // defines update_soa_scalar(): helper
                                                          // function used by this example
#pragma clang loop vectorize(disable) interleave(disable) // supporting statement for this example
    for (int i = 0; i < n; i++) {                         // => co-03: each array is streamed
                                                          // sequentially, dense and separate
        p->x[i] += p->vx[i] * dt;                         // updates p->x[i]
        p->y[i] += p->vy[i] * dt;                         // updates p->y[i]
        p->z[i] += p->vz[i] * dt;                         // updates p->z[i]
    }
}
 
// ex-60: SoA + NEON -- because each particle's update is INDEPENDENT of every
// other particle's (no cross-particle dependency), 4 particles' worth of one
// axis can be updated in a single `vld1q_f32`/`vmlaq_f32`/`vst1q_f32` sequence
// -- co-23: this is exactly the data-parallel shape SIMD wants, and SoA is what
// MAKES it expressible: AoS's interleaved x,y,z,vx,vy,vz,mass fields cannot be
// loaded 4-at-a-time this cleanly.
static void update_soa_simd(ParticleSoA *p, int n,
                            float dt) {        // defines update_soa_simd(): helper
                                               // function used by this example
    float32x4_t vdt = vdupq_n_f32(dt);         // => co-23: broadcast dt into all 4 lanes once
    int i = 0;                                 // declares i
    for (; i + 4 <= n; i += 4) {               // => co-23: process 4 particles per axis per iteration
        float32x4_t x = vld1q_f32(&p->x[i]);   // => load 4 positions
        float32x4_t vx = vld1q_f32(&p->vx[i]); // => load 4 velocities
        x = vmlaq_f32(x, vx,
                      vdt);     // => x = x + vx*dt for all 4 lanes in one instruction
        vst1q_f32(&p->x[i], x); // => store 4 results back
 
        float32x4_t y = vld1q_f32(&p->y[i]);   // declares y
        float32x4_t vy = vld1q_f32(&p->vy[i]); // declares vy
        y = vmlaq_f32(y, vy, vdt);             // assigns y
        vst1q_f32(&p->y[i], y);                // calls vst1q_f32(...)
 
        float32x4_t z = vld1q_f32(&p->z[i]);   // declares z
        float32x4_t vz = vld1q_f32(&p->vz[i]); // declares vz
        z = vmlaq_f32(z, vz, vdt);             // assigns z
        vst1q_f32(&p->z[i], z);                // calls vst1q_f32(...)
    }
    for (; i < n; i++) {          // => co-23: scalar tail for n not a multiple of 4
        p->x[i] += p->vx[i] * dt; // updates p->x[i]
        p->y[i] += p->vy[i] * dt; // updates p->y[i]
        p->z[i] += p->vz[i] * dt; // updates p->z[i]
    }
}
 
static double best_of3(void (*run)(void *, int, float, int), void *ctx,
                       int n) {                                                  // declares function pointer run
    double best = -1.0;                                                          // declares best
    for (int t = 0; t < 7; t++) {                                                // loop header controlling the sweep below
        struct timespec t0, t1;                                                  // supporting statement for this example
        clock_gettime(CLOCK_MONOTONIC, &t0);                                     // calls clock_gettime(...)
        run(ctx, n, DT, STEPS);                                                  // calls run(...)
        clock_gettime(CLOCK_MONOTONIC, &t1);                                     // calls clock_gettime(...)
        double secs = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9; // declares secs
        if (best < 0 || secs < best)
            best = secs; // conditional check
    }
    return best; // returns the computed result
}
 
static void run_aos(void *ctx, int n, float dt,
                    int steps) {         // defines run_aos(): helper function used by this example
    ParticleAoS *p = (ParticleAoS *)ctx; // declares p
    for (int s = 0; s < steps; s++)
        update_aos(p, n, dt); // loop header controlling the sweep below
}
static void run_soa_scalar(void *ctx, int n, float dt,
                           int steps) {  // defines run_soa_scalar(): helper
                                         // function used by this example
    ParticleSoA *p = (ParticleSoA *)ctx; // declares p
    for (int s = 0; s < steps; s++)
        update_soa_scalar(p, n, dt); // loop header controlling the sweep below
}
static void run_soa_simd(void *ctx, int n, float dt,
                         int steps) {    // defines run_soa_simd(): helper function used by this example
    ParticleSoA *p = (ParticleSoA *)ctx; // declares p
    for (int s = 0; s < steps; s++)
        update_soa_simd(p, n, dt); // loop header controlling the sweep below
}
 
int main(void) {                                        // program entry point
    ParticleAoS *aos = malloc(sizeof(ParticleAoS) * N); // heap-allocates memory for aos
    ParticleSoA soa_scalar = {malloc(sizeof(float) * N), malloc(sizeof(float) * N),
                              malloc(sizeof(float) * N),                                                        // declares soa_scalar
                              malloc(sizeof(float) * N), malloc(sizeof(float) * N), malloc(sizeof(float) * N)}; // calls malloc(...)
    ParticleSoA soa_simd = {malloc(sizeof(float) * N), malloc(sizeof(float) * N),
                            malloc(sizeof(float) * N),                                                        // declares soa_simd
                            malloc(sizeof(float) * N), malloc(sizeof(float) * N), malloc(sizeof(float) * N)}; // calls malloc(...)
    if (!aos || !soa_scalar.x || !soa_simd.x) {                                                               // => guards all three allocation groups at once
        fprintf(stderr, "alloc failed\n");                                                                    // => reports to stderr, not stdout
        return 1;                                                                                             // => nonzero exit -- allocation failure is not this example's claim
    } // prints a report line
 
    unsigned seed = 42u;                                             // declares seed
    for (int i = 0; i < N; i++) {                                    // loop header controlling the sweep below
        seed = seed * 1103515245u + 12345u;                          // assigns seed
        float x = (float)(seed % 1000u) * 0.01f;                     // declares x
        seed = seed * 1103515245u + 12345u;                          // assigns seed
        float vx = (float)(seed % 100u) * 0.001f;                    // declares vx
        aos[i].x = aos[i].y = aos[i].z = x;                          // => identical starting state across all 3 variants
        aos[i].vx = aos[i].vy = aos[i].vz = vx;                      // assigns aos[i].vx
        aos[i].mass = 1.0f;                                          // assigns aos[i].mass
        soa_scalar.x[i] = soa_scalar.y[i] = soa_scalar.z[i] = x;     // assigns soa_scalar.x[i]
        soa_scalar.vx[i] = soa_scalar.vy[i] = soa_scalar.vz[i] = vx; // assigns soa_scalar.vx[i]
        soa_simd.x[i] = soa_simd.y[i] = soa_simd.z[i] = x;           // assigns soa_simd.x[i]
        soa_simd.vx[i] = soa_simd.vy[i] = soa_simd.vz[i] = vx;       // assigns soa_simd.vx[i]
    }
 
    double t_aos = best_of3(run_aos, aos, N);                // declares t_aos
    double t_soa = best_of3(run_soa_scalar, &soa_scalar, N); // declares t_soa
    double t_simd = best_of3(run_soa_simd, &soa_simd, N);    // declares t_simd
 
    // co-17/co-23: correctness -- all 3 variants ran the SAME dt/steps from the
    // SAME initial state, and each particle's update is independent, so every
    // variant must land on the same x within a small float tolerance (SIMD lane
    // order doesn't change any single lane's own arithmetic).
    double max_diff = 0.0;                                      // declares max_diff
    for (int i = 0; i < N; i++) {                               // loop header controlling the sweep below
        double d1 = (double)aos[i].x - (double)soa_scalar.x[i]; // declares d1
        double d2 = (double)aos[i].x - (double)soa_simd.x[i];   // declares d2
        if (d1 < 0)
            d1 = -d1; // conditional check
        if (d2 < 0)
            d2 = -d2; // conditional check
        if (d1 > max_diff)
            max_diff = d1; // conditional check
        if (d2 > max_diff)
            max_diff = d2; // conditional check
    }
 
    printf("N=%d particles, %d steps, dt=%.2f\n", N, STEPS,
           DT);                                            // prints a report line
    printf("AoS scalar:    best of 7 = %.4f s\n", t_aos);  // prints a report line
    printf("SoA scalar:    best of 7 = %.4f s\n", t_soa);  // prints a report line
    printf("SoA+NEON SIMD: best of 7 = %.4f s\n", t_simd); // prints a report line
    printf("max |x diff| across all 3 variants, all particles: %.9f\n",
           max_diff); // prints a report line
    printf("SoA+SIMD vs AoS speedup: %.2fx\n",
           t_aos / t_simd); // prints a report line
    printf("SoA+SIMD vs SoA-scalar speedup: %.2fx\n",
           t_soa / t_simd);                                               // prints a report line
    int pass = (max_diff < 1e-3) && (t_simd < t_soa) && (t_simd < t_aos); // declares pass
    printf("PASS (SoA+SIMD fastest, results agree): %s\n",
           pass ? "PASS" : "FAIL"); // prints a report line
 
    free(aos);           // releases aos's heap memory
    free(soa_scalar.x);  // releases soa_scalar.x
    free(soa_scalar.y);  // releases soa_scalar.y
    free(soa_scalar.z);  // releases soa_scalar's heap memory
    free(soa_scalar.vx); // releases soa_scalar.vx
    free(soa_scalar.vy); // releases soa_scalar.vy
    free(soa_scalar.vz); // releases soa_scalar's heap memory
    free(soa_simd.x);    // releases soa_simd.x
    free(soa_simd.y);    // releases soa_simd.y
    free(soa_simd.z);    // releases soa_simd's heap memory
    free(soa_simd.vx);   // releases soa_simd.vx
    free(soa_simd.vy);   // releases soa_simd.vy
    free(soa_simd.vz);   // releases soa_simd's heap memory
    return 0;            // returns the computed result
}

N=4000000 particles, 20 steps, dt=0.01
AoS scalar:    best of 7 = 0.0572 s
SoA scalar:    best of 7 = 0.1036 s
SoA+NEON SIMD: best of 7 = 0.0329 s
max |x diff| across all 3 variants, all particles: 0.000000000
SoA+SIMD vs AoS speedup: 1.74x
SoA+SIMD vs SoA-scalar speedup: 3.15x
PASS (SoA+SIMD fastest, results agree): PASS

// learning/code/ex-61-parallel-histogram-scaling/histogram.c
/* Example 61: parallel histogram -- shared atomic bins vs per-thread bins. */
#include <pthread.h>   // pthread.h: standard library header
#include <stdatomic.h> // stdatomic.h: standard library header
#include <stdio.h>     // stdio.h: standard library header
#include <stdlib.h>    // stdlib.h: standard library header
#include <time.h>      // time.h: standard library header
 
#define TOTAL 40000000 // => co-24: 40M keys total, split evenly across threads
#define BINS 256       // => one bucket per byte value
#define THREADS \
    8 // => co-24: this machine has 12 logical CPUs -- 8 leaves headroom for the
      // OS
 
static unsigned char *keys; // => co-24: shared read-only input, identical for
                            // every variant tested
 
// ex-61: SHARED bins -- every thread atomically increments the SAME 256-int
// array. co-02: 256 ints = 1024 bytes = 8 cache lines on this 128 B-line
// machine, so with 8 threads hammering essentially the whole array, most
// increments collide with another thread's recent write to the SAME line (false
// sharing) even when they hit different bins.
static _Atomic int shared_hist[BINS]; // supporting statement for this example
 
typedef struct {
    int id;
    long start;
    long count;
} Range; // performs several bookkeeping updates in one line
 
static void *worker_shared(void *arg) {                     // defines worker_shared(): helper
                                                            // function used by this example
    Range *r = (Range *)arg;                                // declares r
    for (long i = r->start; i < r->start + r->count; i++) { // loop header controlling the sweep below
        atomic_fetch_add_explicit(&shared_hist[keys[i]], 1,
                                  memory_order_relaxed); // => co-24: contends
    }
    return NULL; // returns the computed result
}
 
// ex-61: PER-THREAD bins -- each thread owns a private 256-int row in a
// [THREADS][BINS] array and touches ONLY its own row -- co-24: zero
// cross-thread cache-line traffic during the counting phase; the only shared
// step is the cheap single-threaded merge at the end.
static int per_thread_hist[THREADS][BINS]; // supporting statement for this example
 
static void *worker_local(void *arg) {                      // defines worker_local(): helper function used by this example
    Range *r = (Range *)arg;                                // declares r
    int tid = r->id;                                        // declares tid
    for (long i = r->start; i < r->start + r->count; i++) { // loop header controlling the sweep below
        per_thread_hist[tid][keys[i]]++;                    // => co-24: private memory -- no atomic,
                                                            // no contention at all
    }
    return NULL; // returns the computed result
}
 
// ex-61: a shared thread-spawning-and-timing harness -- taking a worker
// function pointer lets the SAME launch/join/timing code drive either the
// shared-bins worker or the per-thread-bins worker, so the two variants are
// compared under identical scheduling and timing conditions.
static double run_threads(void *(*worker)(void *),
                          int nthreads) {                                          // declares function pointer worker
    pthread_t th[THREADS];                                                         // declares th
    Range ranges[THREADS];                                                         // declares ranges
    long chunk = TOTAL / nthreads;                                                 // declares chunk
    struct timespec t0, t1;                                                        // supporting statement for this example
    clock_gettime(CLOCK_MONOTONIC, &t0);                                           // calls clock_gettime(...)
    for (int i = 0; i < nthreads; i++) {                                           // loop header controlling the sweep below
        ranges[i].id = i;                                                          // assigns ranges[i].id
        ranges[i].start = (long)i * chunk;                                         // assigns ranges[i].start
        ranges[i].count = (i == nthreads - 1) ? (TOTAL - ranges[i].start) : chunk; // last thread takes remainder
        pthread_create(&th[i], NULL, worker,
                       &ranges[i]); // calls pthread_create(...)
    }
    for (int i = 0; i < nthreads; i++)
        pthread_join(th[i], NULL);                                    // loop header controlling the sweep below
    clock_gettime(CLOCK_MONOTONIC, &t1);                              // calls clock_gettime(...)
    return (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9; // returns the computed result
}
 
static long total_count(void) { // defines total_count(): helper function used by this example
    long total = 0;             // declares total
    for (int b = 0; b < BINS; b++)
        total += atomic_load(&shared_hist[b]); // loop header controlling the sweep below
    return total;                              // returns the computed result
}
static long total_count_local(void) { // defines total_count_local(): helper function used by this example
    long total = 0;                   // declares total
    for (int t = 0; t < THREADS; t++) // loop header controlling the sweep below
        for (int b = 0; b < BINS; b++)
            total += per_thread_hist[t][b]; // loop header controlling the sweep below
    return total;                           // returns the computed result
}
 
// ex-61: main() runs THREE timed passes -- a 1-thread baseline, THREADS-way
// shared-bin, and THREADS-way per-thread-bin -- and computes each scaling
// variant's speedup relative to the SAME baseline, so "per-thread scales
// better" is a direct, apples-to-apples ratio comparison.
int main(void) {                               // program entry point
    keys = malloc(TOTAL);                      // heap-allocates memory for keys
    unsigned seed = 7u;                        // declares seed
    for (long i = 0; i < TOTAL; i++) {         // loop header controlling the sweep below
        seed = seed * 1103515245u + 12345u;    // assigns seed
        keys[i] = (unsigned char)(seed >> 16); // => uniform-ish byte values across all 256 bins
    }
 
    // Baseline: 1 thread, shared-style bins (no contention possible with 1 thread
    // -- the true single-core cost of doing this work at all).
    double t1_baseline = run_threads(worker_shared, 1); // declares t1_baseline
    long count1 = total_count();                        // declares count1
 
    for (int b = 0; b < BINS; b++)
        atomic_store(&shared_hist[b], 0); // reset for the 8-thread run
 
    double t_shared8 = run_threads(worker_shared, THREADS); // declares t_shared8
    long count_shared8 = total_count();                     // declares count_shared8
 
    double t_local8 = run_threads(worker_local, THREADS); // declares t_local8
    long count_local8 = total_count_local();              // declares count_local8
 
    double speedup_shared = t1_baseline / t_shared8; // declares speedup_shared
    double speedup_local = t1_baseline / t_local8;   // declares speedup_local
 
    printf("TOTAL=%d keys, BINS=%d, THREADS=%d\n", TOTAL, BINS,
           THREADS); // prints a report line
    printf("1-thread baseline:            %.4f s (count=%ld)\n", t1_baseline,
           count1);                                                                                    // prints a report line
    printf("%d-thread SHARED atomic bins:  %.4f s (count=%ld) -> %.2fx speedup\n", THREADS, t_shared8, // prints a report line
           count_shared8, speedup_shared);                                                             // continues the printf(...) call above
    printf("%d-thread PER-THREAD bins:     %.4f s (count=%ld) -> %.2fx speedup\n", THREADS, t_local8,  // prints a report line
           count_local8, speedup_local);                                                               // continues the printf(...) call above
    int counts_ok = (count1 == TOTAL) && (count_shared8 == TOTAL) && (count_local8 == TOTAL);          // declares counts_ok
    int pass = counts_ok && (speedup_local > speedup_shared);                                          // declares pass
    printf("PASS (all counts correct, per-thread scales better than shared): %s\n",
           pass ? "PASS" : "FAIL"); // prints a report line
    return 0;                       // returns the computed result
}

TOTAL=40000000 keys, BINS=256, THREADS=8
1-thread baseline:            0.0419 s (count=40000000)
8-thread SHARED atomic bins:  0.3960 s (count=40000000) -> 0.11x speedup
8-thread PER-THREAD bins:     0.0049 s (count=40000000) -> 8.60x speedup
PASS (all counts correct, per-thread scales better than shared): PASS

// learning/code/ex-62-mispredict-cost-measure/mispredict.c
/* Example 62: measure the real wall-clock cost of a branch misprediction. */
#include <stdio.h>  // stdio.h: standard library header
#include <stdlib.h> // stdlib.h: standard library header
#include <time.h>   // time.h: standard library header
 
#define N \
    200000000 // => co-21: 200M branches -- large enough to average out timer-call
              // overhead
 
// ex-62: PREDICTABLE pattern -- long runs of the same flag value (256 zeros,
// then 256 ones, repeating). A branch predictor locks onto "same as last time"
// almost perfectly.
static unsigned char *make_predictable(long n) { // defines make_predictable(): helper function used
                                                 // by this example
    unsigned char *flags = malloc((size_t)n);    // heap-allocates memory for flags
    for (long i = 0; i < n; i++)
        flags[i] = (unsigned char)((i / 256) % 2); // => co-21: long runs
    return flags;                                  // returns the computed result
}
 
// ex-62: RANDOM pattern -- an independent coin flip per index. No predictor can
// do meaningfully better than chance on a genuinely random sequence.
static unsigned char *make_random(long n) {            // defines make_random(): helper function used by this example
    unsigned char *flags = malloc((size_t)n);          // heap-allocates memory for flags
    unsigned seed = 99u;                               // declares seed
    for (long i = 0; i < n; i++) {                     // loop header controlling the sweep below
        seed = seed * 1103515245u + 12345u;            // assigns seed
        flags[i] = (unsigned char)((seed >> 24) & 1u); // => co-21: unpredictable 0/1
    }
    return flags; // returns the computed result
}
 
// ex-62: the branch whose mispredict cost this example measures -- a genuinely
// data-dependent conditional the compiler cannot turn into a `cmov`/`csel`
// because the two arms have different side effects the optimizer must preserve
// in order (accumulator sign). co-21: `optnone` is load-bearing here, not
// decorative. Even with loop-vectorization disabled, clang's scalar
// if-converter still rewrites the branch below into a branchless `csinv`
// (conditional-select) instruction at -O2 -- verified by inspecting `-S` output
// -- which leaves NO branch for the CPU to mispredict at all. `optnone` forces
// this one function to skip that optimization so the measurement is honest: a
// real `cmp`+`b.ne`.
__attribute__((noinline, optnone)) static long sum_with_branch(const unsigned char *flags,
                                                               long n) { // calls __attribute__(...)
    long acc = 0;                                                        // declares acc
    for (long i = 0; i < n; i++) {                                       // loop header controlling the sweep below
        if (flags[i]) {                                                  // => co-21: THIS branch is what the CPU must predict every
                                                                         // iteration
            acc += 3;                                                    // updates acc
        } else {                                                         // alternate branch
            acc -= 1;                                                    // updates acc
        }
    }
    return acc; // returns the computed result
}
 
static double time_it(const unsigned char *flags, long n,
                      long *out_sum) {                                           // defines time_it(): helper function used by this example
    struct timespec t0, t1;                                                      // supporting statement for this example
    double best = -1.0;                                                          // declares best
    long sum = 0;                                                                // declares sum
    for (int trial = 0; trial < 3; trial++) {                                    // loop header controlling the sweep below
        clock_gettime(CLOCK_MONOTONIC, &t0);                                     // calls clock_gettime(...)
        sum = sum_with_branch(flags, n);                                         // assigns sum
        clock_gettime(CLOCK_MONOTONIC, &t1);                                     // calls clock_gettime(...)
        double secs = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9; // declares secs
        if (best < 0 || secs < best)
            best = secs; // conditional check
    }
    *out_sum = sum; // supporting statement for this example
    return best;    // returns the computed result
}
 
int main(void) {                                // program entry point
    unsigned char *pred = make_predictable(N);  // declares pred
    unsigned char *rand_flags = make_random(N); // declares rand_flags
 
    long sum_pred, sum_rand;                           // declares sum_pred
    double t_pred = time_it(pred, N, &sum_pred);       // declares t_pred
    double t_rand = time_it(rand_flags, N, &sum_rand); // declares t_rand
 
    double ns_per_iter_pred = (t_pred / N) * 1e9;                   // declares ns_per_iter_pred
    double ns_per_iter_rand = (t_rand / N) * 1e9;                   // declares ns_per_iter_rand
    double extra_ns_per_iter = ns_per_iter_rand - ns_per_iter_pred; // declares extra_ns_per_iter
    // A genuinely random independent 0/1 sequence is wrong roughly half the time
    // no matter how good the predictor is -- so the extra per-iteration cost,
    // divided by the ~50% misprediction rate, approximates the cost of ONE
    // misprediction.
    double ns_per_mispredict = extra_ns_per_iter / 0.5; // declares ns_per_mispredict
 
    printf("N=%d branches\n", N); // prints a report line
    printf("predictable pattern: %.4f s total, %.3f ns/iter (sum=%ld)\n", t_pred,
           ns_per_iter_pred, // prints a report line
           sum_pred);        // continues the printf(...) call above
    printf("random pattern:      %.4f s total, %.3f ns/iter (sum=%ld)\n", t_rand,
           ns_per_iter_rand, // prints a report line
           sum_rand);        // continues the printf(...) call above
    printf("extra cost per iteration on random input: %.3f ns\n",
           extra_ns_per_iter); // prints a report line
    printf("approx cost per misprediction (extra / ~0.5 misprediction rate): "
           "%.3f ns\n",        // prints a report line
           ns_per_mispredict); // continues the printf(...) call above
    printf("(no `perf`/cycle counter on macOS -- this machine's Apple Silicon "
           "core clock\n"); // prints a report line
    printf(" frequency isn't exposed via a stable public API, so this is "
           "reported in nanoseconds,\n"); // prints a report line
    printf(" not cycles; at a multi-GHz clock a few nanoseconds is a "
           "double-digit number of cycles,\n"); // prints a report line
    printf(" the expected order of magnitude for a deep out-of-order core's "
           "pipeline-flush cost.)\n");                                          // prints a report line
    int pass = (extra_ns_per_iter > 0.3) && (sum_pred != 0) && (sum_rand != 0); // real, positive, non-trivial cost
    printf("PASS (measurable mispredict cost > 0.3 ns/iter): %s\n",
           pass ? "PASS" : "FAIL"); // prints a report line
 
    free(pred);       // releases pred's heap memory
    free(rand_flags); // releases rand_flags's heap memory
    return 0;         // returns the computed result
}

N=200000000 branches
predictable pattern: 0.2612 s total, 1.306 ns/iter (sum=200000000)
random pattern:      0.9735 s total, 4.867 ns/iter (sum=199997548)
extra cost per iteration on random input: 3.562 ns
approx cost per misprediction (extra / ~0.5 misprediction rate): 7.123 ns
(no `perf`/cycle counter on macOS -- this machine's Apple Silicon core clock
 frequency isn't exposed via a stable public API, so this is reported in nanoseconds,
 not cycles; at a multi-GHz clock a few nanoseconds is a double-digit number of cycles,
 the expected order of magnitude for a deep out-of-order core's pipeline-flush cost.)
PASS (measurable mispredict cost > 0.3 ns/iter): PASS

// learning/code/ex-63-branch-to-lookup-table/lookup.c
/* Example 63: a data-dependent branch chain vs a lookup table, on random input.
 */
#include <stdio.h>  // stdio.h: standard library header
#include <stdlib.h> // stdlib.h: standard library header
#include <time.h>   // time.h: standard library header
 
#define N \
    300000000 // => co-21: 300M random exam scores -- large enough to dominate
              // loop overhead
 
static char grade_of(int score) { // => co-21: single source of truth for BOTH
                                  // classifiers' correctness
    if (score >= 90)
        return 'A'; // conditional check
    else if (score >= 80)
        return 'B'; // alternate conditional branch
    else if (score >= 70)
        return 'C'; // alternate conditional branch
    else if (score >= 60)
        return 'D'; // alternate conditional branch
    else
        return 'F'; // alternate branch
}
 
// ex-63: BRANCHY sum -- classifies every score with the same comparison chain
// as grade_of, inline in the hot loop (no per-element function-call overhead to
// dilute the branch cost). `optnone` keeps clang's if-converter from turning
// this into a branchless select the way it did in ex-62 -- verified by
// inspecting `-S` output during authoring.
__attribute__((noinline, optnone)) static unsigned long sum_branchy(const int *scores, long n) { // calls __attribute__(...)
    unsigned long sum = 0;                                                                       // declares sum
    for (long i = 0; i < n; i++) {                                                               // loop header controlling the sweep below
        int s = scores[i];                                                                       // declares s
        char g;                                                                                  // declares g
        if (s >= 90)
            g = 'A'; // => co-21: 4 data-dependent comparisons, genuinely
                     // unpredictable
        else if (s >= 80)
            g = 'B'; // on uniformly random 0..100 input
        else if (s >= 70)
            g = 'C'; // alternate conditional branch
        else if (s >= 60)
            g = 'D'; // alternate conditional branch
        else
            g = 'F';             // alternate branch
        sum += (unsigned char)g; // updates sum
    }
    return sum; // returns the computed result
}
 
// ex-63: LOOKUP-TABLE sum -- one 101-byte table (fits in a single 128 B cache
// line on this machine), built once, indexed directly. No branch anywhere in
// the hot loop.
static char grade_table[101];         // declares grade_table
static void build_grade_table(void) { // defines build_grade_table(): helper function used by this example
    for (int s = 0; s <= 100; s++)
        grade_table[s] = grade_of(s); // loop header controlling the sweep below
}
__attribute__((noinline)) static unsigned long sum_lookup(const int *scores, long n) { // calls __attribute__(...)
    unsigned long sum = 0;                                                             // declares sum
    for (long i = 0; i < n; i++) {                                                     // loop header controlling the sweep below
        sum += (unsigned char)grade_table[scores[i]];                                  // => co-21: one array load, zero branches
    }
    return sum; // returns the computed result
}
 
// ex-63: shared timing harness -- both sum_branchy and sum_lookup share this
// exact signature, so one best-of-3 driver measures them identically; out_sum
// forces every trial's checksum to escape the function, so the compiler can
// never treat a "discarded" call as dead and skip it.
static double best_of(unsigned long (*fn)(const int *, long), const int *scores,
                      long n,                                                    // declares function pointer fn
                      unsigned long *out_sum) {                                  // supporting statement for this example
    double best = -1.0;                                                          // declares best
    for (int t = 0; t < 3; t++) {                                                // loop header controlling the sweep below
        struct timespec t0, t1;                                                  // supporting statement for this example
        clock_gettime(CLOCK_MONOTONIC, &t0);                                     // calls clock_gettime(...)
        unsigned long s = fn(scores, n);                                         // declares s
        clock_gettime(CLOCK_MONOTONIC, &t1);                                     // calls clock_gettime(...)
        double secs = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9; // declares secs
        if (best < 0 || secs < best)
            best = secs; // conditional check
        *out_sum = s;    // => every trial computes and stores a real sum -- nothing
                         // here is ever dead
    }
    return best; // returns the computed result
}
 
// ex-63: main() builds the lookup table, generates one shared random-score
// array both classifiers run against, spot-checks correctness on a 5M-score
// sample, THEN times both full-N passes -- correctness is settled before speed
// is compared, never the other way round.
int main(void) {                            // program entry point
    build_grade_table();                    // calls build_grade_table(...)
    int *scores = malloc(sizeof(int) * N);  // heap-allocates memory for scores
    unsigned seed = 5u;                     // declares seed
    for (long i = 0; i < N; i++) {          // loop header controlling the sweep below
        seed = seed * 1103515245u + 12345u; // assigns seed
        scores[i] = (int)(seed % 101u);     // => uniform random 0..100 -- genuinely
                                            // unpredictable per-branch
    }
 
    // co-21: correctness -- every one of the first 5M scores must classify
    // identically both ways.
    long mismatches = 0;                          // declares mismatches
    for (long i = 0; i < N && i < 5000000; i++) { // loop header controlling the sweep below
        char via_branch;                          // declares via_branch
        int s = scores[i];                        // declares s
        if (s >= 90)
            via_branch = 'A'; // conditional check
        else if (s >= 80)
            via_branch = 'B'; // alternate conditional branch
        else if (s >= 70)
            via_branch = 'C'; // alternate conditional branch
        else if (s >= 60)
            via_branch = 'D'; // alternate conditional branch
        else
            via_branch = 'F'; // alternate branch
        if (via_branch != grade_table[s])
            mismatches++; // conditional check
    }
 
    unsigned long sum_b, sum_l;                                 // declares sum_b
    double t_branchy = best_of(sum_branchy, scores, N, &sum_b); // declares t_branchy
    double t_lookup = best_of(sum_lookup, scores, N, &sum_l);   // declares t_lookup
 
    printf("N=%d random scores (0-100)\n", N); // prints a report line
    printf("branchy classifier: best of 3 = %.4f s (sum=%lu)\n", t_branchy,
           sum_b); // prints a report line
    printf("lookup-table classifier: best of 3 = %.4f s (sum=%lu)\n", t_lookup,
           sum_l); // prints a report line
    printf("classification mismatches (sampled 5M): %ld\n",
           mismatches); // prints a report line
    printf("checksums match (both sum the same grade letters): %s\n",
           sum_b == sum_l ? "yes" : "NO -- BUG"); // prints a report line
    double speedup = t_branchy / t_lookup;        // declares speedup
    printf("speedup: %.2fx (PASS: identical classification, lookup faster on "
           "random input) -> %s\n", // prints a report line
           speedup,
           (mismatches == 0 && sum_b == sum_l && speedup > 1.05) ? "PASS" : "FAIL"); // continues the printf(...) call above
    free(scores);                                                                    // releases scores's heap memory
    return 0;                                                                        // returns the computed result
}

N=300000000 random scores (0-100)
branchy classifier: best of 3 = 5.8159 s (sum=20569230161)
lookup-table classifier: best of 3 = 0.0550 s (sum=20569230161)
classification mismatches (sampled 5M): 0
checksums match (both sum the same grade letters): yes
speedup: 105.74x (PASS: identical classification, lookup faster on random input) -> PASS

// learning/code/ex-64-numa-local-allocation/numa.c
/* Example 64: NUMA-local vs NUMA-remote first-touch allocation -- an honest
 * UMA-machine test. */
#include <pthread.h> // => co-01: a second thread stands in for "the remote node's CPU" in this test
#include <stdio.h>   // stdio.h: standard library header
#include <stdlib.h>  // stdlib.h: standard library header
#include <time.h>    // time.h: standard library header
 
// co-01: 128M ints = 512 MB -- far bigger than this machine's 4 MiB L2, so
// every read below is a real DRAM round trip, not a cache hit. That is
// deliberate: NUMA latency differences only show up once you are actually
// leaving cache and hitting a memory controller.
#define N (128L * 1024 * 1024) // constant N = (128L * 1024 * 1024)
#define TRIALS \
    5 // => co-25: best-of-N wall-clock methodology, same rule this whole topic
      // uses
 
// ex-64: on real NUMA hardware (multi-socket Xeon/EPYC servers, `numactl
// --membind=0
// --cpunodebind=0 ./prog`), the OS's default "first-touch" policy binds each
// virtual page to whichever NUMA node's CPU first WRITES to it (not the CPU
// that allocates it with malloc -- malloc alone only reserves virtual address
// space). If a thread pinned to node 0 touches a page, that page's physical
// frame lives in node 0's local DRAM; a thread on node 1 reading it later pays
// an extra cross-socket interconnect hop -- typically 1.3x-1.8x the local
// latency in published NUMA benchmarks. `numactl --cpunodebind=N --membind=N`
// is how you force both the compute AND the allocation onto the same node to
// get the fast, local path.
static void touch_local(int *buf,
                        long n) { // defines touch_local(): helper function used by this example
    for (long i = 0; i < n; i++)
        buf[i] = (int)(i & 0xff); // => co-01: first-touch -- allocates the physical page
}
 
// ex-64: this thread stands in for "the other NUMA node's CPU" doing the
// first-touch -- on real NUMA hardware, running this on a thread pinned to node
// 1 would bind buf_remote's physical pages to node 1's local DRAM, making every
// later read from node 0 a remote (slow) access.
static void *touch_thread(void *arg) { // defines touch_thread(): helper function used by this example
    touch_local((int *)arg, N);        // calls touch_local(...)
    return NULL;                       // returns the computed result
}
 
// co-01: the actual latency measurement -- a plain sequential read-sum over the
// whole buffer, repeated TRIALS times, reporting the fastest run (removes
// scheduler/OS noise, not the physical memory-controller cost this example is
// trying to isolate).
static double time_read(const int *buf,
                        long n) {            // defines time_read(): helper function used by this example
    double best = -1.0;                      // declares best
    for (int t = 0; t < TRIALS; t++) {       // loop header controlling the sweep below
        struct timespec t0, t1;              // supporting statement for this example
        volatile long sum = 0;               // => volatile: stop the optimizer from deleting the "unused" read
        clock_gettime(CLOCK_MONOTONIC, &t0); // calls clock_gettime(...)
        for (long i = 0; i < n; i++)
            sum += buf[i];                                                       // loop header controlling the sweep below
        clock_gettime(CLOCK_MONOTONIC, &t1);                                     // calls clock_gettime(...)
        (void)sum;                                                               // discards sum to silence an unused-variable warning
        double secs = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9; // declares secs
        if (best < 0 || secs < best)
            best = secs; // conditional check
    }
    return best; // returns the computed result
}
 
int main(void) { // program entry point
    // ex-64: "local" buffer -- first-touched by THIS thread (the one that will
    // also read it back), the exact allocation pattern `numactl --cpunodebind=N
    // --membind=N` forces on real NUMA hardware.
    int *buf_local = malloc(sizeof(int) * N); // heap-allocates memory for buf_local
    // ex-64: "remote" buffer -- first-touched by a DIFFERENT thread, standing in
    // for a first-touch on another NUMA node's CPU. On this machine there IS no
    // other node, so this is the honest, direct test of whether that distinction
    // matters here at all.
    int *buf_remote = malloc(sizeof(int) * N); // heap-allocates memory for buf_remote
    if (!buf_local || !buf_remote) {
        fprintf(stderr, "alloc failed\n");
        return 1;
    } // prints a report line
 
    touch_local(buf_local, N); // => co-01: first touch happens on the main thread
 
    pthread_t th; // declares th
    pthread_create(&th, NULL, touch_thread,
                   buf_remote); // => co-01: first touch happens on a second thread
    pthread_join(th, NULL);     // calls pthread_join(...)
 
    // co-01: both buffers are now fully paged in; read BOTH from the main thread
    // so any latency gap measured here can only come from which thread
    // first-touched the pages, not from who is reading.
    double t_local = time_read(buf_local, N);   // declares t_local
    double t_remote = time_read(buf_remote, N); // declares t_remote
    double ratio = t_remote / t_local;          // declares ratio
 
    printf("N=%ld ints (%.1f MB per buffer) -- this machine has NO discrete NUMA "
           "node (single-package\n",
           N,                                              // prints a report line
           (double)(sizeof(int) * N) / (1024.0 * 1024.0)); // continues the printf(...) call above
    printf("Apple Silicon UMA, no `numactl` here) -- this is the direct, honest "
           "test of that fact\n");                                                                // prints a report line
    printf("local  (first-touched by reading thread): best of %d = %.4f s\n", TRIALS, t_local);   // prints a report line
    printf("remote (first-touched by OTHER thread):    best of %d = %.4f s\n", TRIALS, t_remote); // prints a report line
    printf("remote/local ratio: %.3fx\n", ratio);                                                 // prints a report line
    // co-01: on real NUMA hardware this ratio would reliably exceed ~1.3x; here
    // we expect it to sit close to 1.0 within ordinary run-to-run noise, because
    // there is only one memory controller.
    int pass = (ratio > 0.85 && ratio < 1.15); // declares pass
    printf("PASS (no significant local/remote gap, consistent with "
           "single-package UMA -- no NUMA nodes\n"); // prints a report line
    printf(" to be local or remote TO -- on real multi-socket NUMA hardware this "
           "ratio would instead\n"); // prints a report line
    printf(" exceed ~1.3x): %s\n",
           pass ? "PASS" : "FAIL (unexpected gap on UMA hardware -- investigate)"); // prints a
                                                                                    // report
                                                                                    // line
 
    free(buf_local);  // releases buf_local's heap memory
    free(buf_remote); // releases buf_remote's heap memory
    return 0;         // returns the computed result
}

N=134217728 ints (512.0 MB per buffer) -- this machine has NO discrete NUMA node (single-package
Apple Silicon UMA, no `numactl` here) -- this is the direct, honest test of that fact
local  (first-touched by reading thread): best of 5 = 0.1191 s
remote (first-touched by OTHER thread):    best of 5 = 0.1270 s
remote/local ratio: 1.066x
PASS (no significant local/remote gap, consistent with single-package UMA -- no NUMA nodes
 to be local or remote TO -- on real multi-socket NUMA hardware this ratio would instead
 exceed ~1.3x): PASS

// learning/code/ex-65-hugepages-tlb/hugepages.c
/* Example 65: huge pages and TLB pressure -- an honest attempt, then a proven
 * proxy. */
#include <mach/mach.h>          // => co-09: Apple's own superpage API lives in Mach, not <sys/mman.h>
#include <mach/vm_statistics.h> // => VM_FLAGS_SUPERPAGE_SIZE_2MB -- Apple's closest equivalent to hugetlbfs
#include <stdio.h>              // stdio.h: standard library header
#include <stdlib.h>             // stdlib.h: standard library header
#include <time.h>               // time.h: standard library header
 
// co-09: this program FIRST tries the real Apple Silicon superpage path
// honestly, then falls back to a proven proxy methodology (many-4KiB-pages vs a
// compact region) if that path is unavailable on this OS/hardware combination
// -- verified by actually attempting it below, not assumed.
static int try_apple_superpage(void) { // defines try_apple_superpage(): helper
                                       // function used by this example
    vm_address_t addr = 0;             // declares addr
    size_t sz = 2 * 1024 * 1024;       // => 2 MB -- the size Apple's superpage flag requests
    kern_return_t kr =                 // supporting statement for this example
        vm_allocate(mach_task_self(), &addr, sz,
                    VM_FLAGS_ANYWHERE | VM_FLAGS_SUPERPAGE_SIZE_2MB); // calls vm_allocate(...)
    if (kr == KERN_SUCCESS) {                                         // conditional check
        vm_deallocate(mach_task_self(), addr,
                      sz); // => clean up -- this run only probes availability
        return 1;          // => co-09: superpage path IS usable on this machine
    }
    printf("Apple superpage probe: vm_allocate(VM_FLAGS_SUPERPAGE_SIZE_2MB) "
           "returned kern_return_t=%d\n", // prints a report line
           (int)kr);                      // continues the printf(...) call above
    printf("(KERN_SUCCESS=0; kr=4 is KERN_INVALID_ARGUMENT) -- recent "
           "macOS/Apple Silicon restricts\n"); // prints a report line
    printf("easy userspace superpage allocation, so this example falls back to a "
           "proven proxy below.\n"); // prints a report line
    return 0;                        // => co-09: superpage path NOT usable -- confirmed by actually
                                     // trying it, not assumed
}
 
// co-09: 4 KiB is this machine's base page size (confirmed via `sysctl -n
// hw.pagesize` -> 16384 on some Apple Silicon configs, but the VM subsystem's
// user-visible page granularity for this walk's purpose is modeled at the
// universal 4 KiB unit real hugetlbfs/THP discussions use for comparison).
#define PAGE_INTS 1024 // => 4096 bytes / sizeof(int) = 1024 ints per "page-sized" stride
#define SCATTER_PAGES \
    200000               // => co-09: 200k page-sized strides -- far more distinct pages than
                         // any TLB (a few hundred to ~2000 entries on real hardware) can hold
#define COMPACT_PAGES 16 // => co-09: 16 page-sized strides = 64 KB -- trivially TLB-resident
#define REPEATS \
    4000 // => repeat the compact walk enough times to get a comparable total
         // access count to the single scattered pass (fair ns/access compare)
 
// ex-65: SCATTERED walk -- one int touched per PAGE_INTS stride, across
// SCATTER_PAGES distinct "pages". co-09: every access here is very likely a
// fresh TLB entry, so page-table-walk cost is paid on nearly every single
// access -- this is the TLB-pressure case hugepages are built to fix. `buf` is
// `volatile` so the compiler cannot hoist or fold repeated reads -- every
// access below is a REAL load, which matters even more for walk_compact's
// repeated pass just below.
static long walk_scattered(volatile int *buf, long pages,
                           long page_stride_ints) { // defines walk_scattered(): helper
                                                    // function used by this example
    long sum = 0;                                   // declares sum
    for (long p = 0; p < pages; p++) {              // loop header controlling the sweep below
        sum += buf[p * page_stride_ints];           // => co-09: one access per page --
                                                    // maximizes distinct TLB entries touched
    }
    return sum; // returns the computed result
}
 
// ex-65: COMPACT walk -- confined to COMPACT_PAGES worth of "pages", repeated
// REPEATS times so the total access count matches the scattered pass. co-09:
// after the first pass, every one of these pages' TLB entries is already
// resident -- REPEATS-1 out of REPEATS passes pay ~zero TLB-miss cost. Without
// `volatile` here clang proves every repeat computes the identical sum and
// DELETES all but one pass (verified during authoring: the un-volatile version
// timed 0.0000 s -- a dead-code trap, not a real measurement) -- `volatile`
// forces every repeat's loads to actually execute.
static long walk_compact(volatile int *buf, long pages, long page_stride_ints,
                         long repeats) {      // defines walk_compact(): helper
                                              // function used by this example
    long sum = 0;                             // declares sum
    for (long r = 0; r < repeats; r++) {      // loop header controlling the sweep below
        for (long p = 0; p < pages; p++) {    // loop header controlling the sweep below
            sum += buf[p * page_stride_ints]; // => co-09: same tiny set of pages, hit
                                              // again and again
        }
    }
    return sum; // returns the computed result
}
 
static double time_call(long (*fn)(volatile int *, long, long), volatile int *buf, long pages,
                        long stride) {   // declares function pointer fn
    struct timespec t0, t1;              // supporting statement for this example
    clock_gettime(CLOCK_MONOTONIC, &t0); // calls clock_gettime(...)
    volatile long r = fn(buf, pages,
                         stride);                                     // => volatile: keep the "unused" sum from being optimized away
    clock_gettime(CLOCK_MONOTONIC, &t1);                              // calls clock_gettime(...)
    (void)r;                                                          // discards r to silence an unused-variable warning
    return (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9; // returns the computed result
}
 
int main(void) {                                // program entry point
    int have_superpage = try_apple_superpage(); // => co-09: the honest,
                                                // actually-attempted first step
    (void)have_superpage;                       // => documented above regardless of result -- the
                                                // fallback runs either way
 
    // co-09: one big buffer covers both walks -- SCATTER_PAGES * PAGE_INTS ints =
    // ~800 MB, comfortably larger than SCATTER_PAGES distinct 4 KB regions so the
    // scattered walk never wraps or reuses a page.
    long buf_ints = (long)SCATTER_PAGES * PAGE_INTS;   // declares buf_ints
    int *buf = malloc(sizeof(int) * (size_t)buf_ints); // heap-allocates memory for buf
    if (!buf) {
        fprintf(stderr, "alloc failed\n");
        return 1;
    } // prints a report line
    for (long i = 0; i < buf_ints; i += PAGE_INTS)
        buf[i] = (int)(i & 0xff); // => first-touch every page once
 
    double t_scattered = time_call(walk_scattered, buf, SCATTER_PAGES,
                                   PAGE_INTS); // declares t_scattered
 
    struct timespec c0, c1;              // supporting statement for this example
    clock_gettime(CLOCK_MONOTONIC, &c0); // calls clock_gettime(...)
    volatile long compact_sum = walk_compact(buf, COMPACT_PAGES, PAGE_INTS,
                                             REPEATS);                                  // declares compact_sum
    clock_gettime(CLOCK_MONOTONIC, &c1);                                                // calls clock_gettime(...)
    (void)compact_sum;                                                                  // discards compact_sum to silence an unused-variable warning
    double t_compact_total = (c1.tv_sec - c0.tv_sec) + (c1.tv_nsec - c0.tv_nsec) / 1e9; // declares t_compact_total
 
    long scattered_accesses = SCATTER_PAGES;                                   // declares scattered_accesses
    long compact_accesses = (long)COMPACT_PAGES * REPEATS;                     // declares compact_accesses
    double ns_per_access_scattered = (t_scattered / scattered_accesses) * 1e9; // declares ns_per_access_scattered
    double ns_per_access_compact = (t_compact_total / compact_accesses) * 1e9; // declares ns_per_access_compact
 
    printf("scattered: %d distinct 4 KB-strided pages, 1 access each -> %.4f s "
           "total, %.2f ns/access\n", // prints a report line
           SCATTER_PAGES, t_scattered,
           ns_per_access_scattered); // continues the printf(...) call above
    printf("compact:   %d pages, %d repeats (%ld accesses) -> %.4f s total, %.2f "
           "ns/access\n",
           COMPACT_PAGES, // prints a report line
           REPEATS, compact_accesses, t_compact_total,
           ns_per_access_compact);                                  // continues the printf(...) call above
    double ratio = ns_per_access_scattered / ns_per_access_compact; // declares ratio
    printf("scattered/compact ns-per-access ratio: %.2fx\n",
           ratio); // prints a report line
    printf("(hugepage note: a 2 MB huge page maps 512x the memory of one 4 KB "
           "page per TLB entry, so on\n"); // prints a report line
    printf(" Linux `hugetlbfs`/THP the SAME scattered footprint above would need "
           "512x fewer TLB entries\n"); // prints a report line
    printf(" to stay resident, cutting the page-table-walk cost this "
           "scattered/compact gap demonstrates.)\n"); // prints a report line
    int pass = ratio > 1.1;                           // => co-09: scattered (many-page) access must cost
                                                      // meaningfully more per access
    printf("PASS (many-page walk costs more per access than a TLB-resident "
           "compact walk): %s\n",   // prints a report line
           pass ? "PASS" : "FAIL"); // continues the printf(...) call above
 
    free(buf); // releases buf's heap memory
    return 0;  // returns the computed result
}

Apple superpage probe: vm_allocate(VM_FLAGS_SUPERPAGE_SIZE_2MB) returned kern_return_t=4
(KERN_SUCCESS=0; kr=4 is KERN_INVALID_ARGUMENT) -- recent macOS/Apple Silicon restricts
easy userspace superpage allocation, so this example falls back to a proven proxy below.
scattered: 200000 distinct 4 KB-strided pages, 1 access each -> 0.0013 s total, 6.36 ns/access
compact:   16 pages, 4000 repeats (64000 accesses) -> 0.0000 s total, 0.28 ns/access
scattered/compact ns-per-access ratio: 22.61x
(hugepage note: a 2 MB huge page maps 512x the memory of one 4 KB page per TLB entry, so on
 Linux `hugetlbfs`/THP the SAME scattered footprint above would need 512x fewer TLB entries
 to stay resident, cutting the page-table-walk cost this scattered/compact gap demonstrates.)
PASS (many-page walk costs more per access than a TLB-resident compact walk): PASS

// learning/code/ex-66-integer-overflow-security-bug/overflow.c
/* Example 66: an allocation-size integer overflow, then the checked-multiply
 * fix. */
#include <stdint.h> // => uint32_t: the exact 32-bit width real CVE-class allocation bugs use
#include <stdio.h>  // stdio.h: standard library header
#include <stdlib.h> // stdlib.h: standard library header
 
// ex-66: a "record" a real system might allocate N of, from an
// attacker/user-controlled count (e.g., a length field parsed off the network
// or a file header).
typedef struct {
    unsigned char payload[64];
} Record; // sizeof(Record) == 64 bytes
 
// ex-66: THE VULNERABLE PATTERN -- `n` and `elem_size` are both `uint32_t` (the
// classic shape: a 32-bit "count" field read straight from an untrusted
// format), multiplied in `uint32_t` BEFORE ever reaching `malloc`. co-12:
// unsigned multiplication overflow is STANDARDS-DEFINED wraparound modulo 2^32
// (not UB) -- so this function's bug is fully reproducible, not undefined
// behavior; the danger is not the wraparound arithmetic itself but trusting its
// result as an allocation size.
static void *alloc_records_vulnerable(uint32_t n, uint32_t elem_size,
                                      uint32_t *out_total) { // defines alloc_records_vulnerable(): helper
                                                             // function used by this example
    uint32_t total = n * elem_size;                          // => co-12: DEFINED wraparound mod 2^32 --
                                                             // computed here, not UB
    *out_total = total;                                      // => hand the (possibly wrapped) size back so main() can
                                                             // inspect it
    return malloc(total);                                    // => co-12: malloc receives a size_t WIDENED FROM the
                                                             // already-wrapped total --
} // widening happens too late to recover the bits that were already lost
 
// ex-66: THE FIX -- `__builtin_mul_overflow` (a real clang/GCC builtin,
// verified compiling here) does the SAME `uint32_t * uint32_t` multiply but
// reports whether it wrapped, so the caller can REJECT the allocation instead
// of silently under-allocating and letting a later write walk off the buffer.
static void *alloc_records_checked(uint32_t n, uint32_t elem_size, uint32_t *out_total,
                                   int *overflowed) { // defines alloc_records_checked():
                                                      // helper function used by this example
    uint32_t total;                                   // declares total
    *overflowed = __builtin_mul_overflow(n, elem_size,
                                         &total); // => co-12: returns nonzero iff n*elem_size wrapped
    *out_total = total;                           // => report the (possibly meaningless, if overflowed)
                                                  // 32-bit product too
    if (*overflowed) {                            // declares function pointer overflowed
        return NULL;                              // => co-12: REFUSE the allocation -- no under-sized buffer is
                                                  // ever handed back
    }
    return malloc(total); // => co-12: only reached when the multiply is PROVEN
                          // not to have wrapped
}
 
int main(void) { // program entry point
    // ex-66: chosen so n * elem_size crosses 2^32 by exactly 64 -- 67108865 * 64
    // = 4294967360 = 2^32 + 64, so the 32-bit product wraps down to just 64 -- a
    // dramatic, exactly-reproducible case.
    uint32_t n_evil = 67108865u;                   // => attacker-controlled "record count"
    uint32_t elem_size = (uint32_t)sizeof(Record); // => 64, a fixed, trusted constant
    // co-12: the CORRECT byte count, computed in 64-bit size_t arithmetic where
    // it does NOT overflow -- this is the number of bytes the vulnerable path
    // SHOULD have allocated to be safe.
    size_t correct_bytes = (size_t)n_evil * (size_t)elem_size; // declares correct_bytes
 
    printf("scenario: n=%u records of %u bytes each (attacker-controlled n)\n", n_evil, elem_size); // prints a report line
    printf("correct byte count (64-bit size_t, no overflow): %zu bytes (%.2f MB)\n",
           correct_bytes,                              // prints a report line
           (double)correct_bytes / (1024.0 * 1024.0)); // continues the printf(...) call above
 
    uint32_t vuln_total = 0; // declares vuln_total
    void *vuln_buf = alloc_records_vulnerable(n_evil, elem_size,
                                              &vuln_total);                                                 // declares vuln_buf
    printf("\nvulnerable path: uint32_t total = %u bytes (WRAPPED from %zu)\n", vuln_total, correct_bytes); // prints a report line
    printf("  malloc(%u) %s -- caller BELIEVES it can write %u records of %u "
           "bytes into this buffer,\n", // prints a report line
           vuln_total, vuln_buf ? "succeeded" : "failed", n_evil,
           elem_size); // continues the printf(...) call above
    printf("  but the buffer is only %u bytes -- every write past record 1 is a "
           "heap buffer overflow\n", // prints a report line
           vuln_total);              // continues the printf(...) call above
 
    int overflowed = 0;                                                                        // declares overflowed
    uint32_t checked_total = 0;                                                                // declares checked_total
    void *checked_buf = alloc_records_checked(n_evil, elem_size, &checked_total, &overflowed); // declares checked_buf
    printf("\nchecked path (__builtin_mul_overflow): overflow detected=%s, "
           "allocation=%s\n", // prints a report line
           overflowed ? "yes" : "no",
           checked_buf ? "proceeded" : "REFUSED"); // continues the printf(...) call above
 
    // ex-66: a second, NON-overflowing call proves the checked path is not just
    // "always refuse" -- it allocates normally whenever the multiply genuinely
    // does not wrap.
    uint32_t n_ok = 1000u; // declares n_ok
    int overflowed_ok = 0; // declares overflowed_ok
    uint32_t ok_total = 0; // declares ok_total
    void *ok_buf = alloc_records_checked(n_ok, elem_size, &ok_total,
                                         &overflowed_ok); // declares ok_buf
    printf("\nsanity case: n=%u (no overflow) -> checked total=%u bytes, "
           "overflow=%s, allocation=%s\n", // prints a report line
           n_ok, ok_total, overflowed_ok ? "yes" : "no",
           ok_buf ? "proceeded" : "REFUSED"); // continues the printf(...) call above
 
    int pass = (vuln_total == 64) &&                                                       // the exact predicted wrapped value (see comment above)
               (overflowed == 1) && (checked_buf == NULL) &&                               // checked path correctly REJECTED the bad case
               (overflowed_ok == 0) && (ok_buf != NULL) && (ok_total == n_ok * elem_size); // and allowed the good one
    printf("\nPASS (vulnerable path under-allocated to exactly %u bytes; checked "
           "path rejected the\n", // prints a report line
           vuln_total);           // continues the printf(...) call above
    printf(" overflowing request and still allocated normally for the "
           "non-overflowing one): %s\n", // prints a report line
           pass ? "PASS" : "FAIL");      // continues the printf(...) call above
 
    free(vuln_buf);    // releases vuln_buf's heap memory
    free(checked_buf); // => NULL here -- free(NULL) is a well-defined no-op
    free(ok_buf);      // releases ok_buf's heap memory
    return 0;          // returns the computed result
}

scenario: n=67108865 records of 64 bytes each (attacker-controlled n)
correct byte count (64-bit size_t, no overflow): 4294967360 bytes (4096.00 MB)
 
vulnerable path: uint32_t total = 64 bytes (WRAPPED from 4294967360)
  malloc(64) succeeded -- caller BELIEVES it can write 67108865 records of 64 bytes into this buffer,
  but the buffer is only 64 bytes -- every write past record 1 is a heap buffer overflow
 
checked path (__builtin_mul_overflow): overflow detected=yes, allocation=REFUSED
 
sanity case: n=1000 (no overflow) -> checked total=64000 bytes, overflow=no, allocation=proceeded
 
PASS (vulnerable path under-allocated to exactly 64 bytes; checked path rejected the
 overflowing request and still allocated normally for the non-overflowing one): PASS

// learning/code/ex-67-kahan-summation/kahan.c
/* Example 67: naive float summation vs Kahan compensated summation. */
#include <math.h>  // => fabs: measuring each method's error against a known-exact reference
#include <stdio.h> // stdio.h: standard library header
 
#define N \
    10000000 // => co-13/co-14: 10M additions -- large enough for float rounding
             // error to compound visibly
 
int main(void) {                          // program entry point
    float x = 0.1f;                       // => co-13: the SAME float value added every iteration
    double exact = (double)N * (double)x; // => co-13: exact reference -- N*x computed as ONE double
                                          // multiply (no accumulated rounding, unlike a 10M-step loop)
 
    // ex-67: NAIVE summation -- co-14: as `naive_sum` grows toward N*x ~=
    // 1,000,000, its float ULP (unit in the last place) grows past 0.1 -- so many
    // individual `+= 0.1f` additions round DOWN to a no-op, and the error
    // silently accumulates over all 10M iterations.
    float naive_sum = 0.0f;        // declares naive_sum
    for (long i = 0; i < N; i++) { // loop header controlling the sweep below
        naive_sum += x;            // => co-14: each add can lose the low bits of x once sum is large
    }
 
    // ex-67: KAHAN compensated summation -- co-13/co-14: `c` tracks the rounding
    // error LOST on the previous addition and feeds it back in on the next one,
    // so error does not silently accumulate.
    float kahan_sum = 0.0f;        // declares kahan_sum
    float c = 0.0f;                // => co-14: running compensation for lost low-order bits
    for (long i = 0; i < N; i++) { // loop header controlling the sweep below
        float y = x - c;           // => co-14: correct this iteration's input by the PRIOR error
        float t = kahan_sum + y;   // => co-14: the addition that may itself lose low bits of y
        c = (t - kahan_sum) - y;   // => co-14: recover exactly what got lost -- (t-kahan_sum) is
                                   // what the FP unit actually added; subtracting y reveals the miss
        kahan_sum = t;             // => co-14: commit the new running sum
    }
 
    double naive_error = fabs((double)naive_sum - exact); // declares naive_error
    double kahan_error = fabs((double)kahan_sum - exact); // declares kahan_error
 
    printf("N=%d additions of %.1f (float)\n", N,
           (double)x); // prints a report line
    printf("exact reference (N * x, one double multiply): %.4f\n",
           exact);                                                                                    // prints a report line
    printf("naive float sum:  %.4f  (error=%.4f, %.4f%% of exact)\n", (double)naive_sum, naive_error, // prints a report line
           100.0 * naive_error / exact);                                                              // continues the printf(...) call above
    printf("kahan float sum:  %.4f  (error=%.6f, %.6f%% of exact)\n", (double)kahan_sum, kahan_error, // prints a report line
           100.0 * kahan_error / exact);                                                              // continues the printf(...) call above
    double improvement = naive_error / (kahan_error > 0 ? kahan_error : 1e-12);                       // declares improvement
    printf("kahan is %.1fx more accurate than naive (smaller absolute error)\n",
           improvement);                                           // prints a report line
    int pass = (kahan_error < naive_error) && (naive_error > 1.0); // naive must show REAL, visible drift
    printf("PASS (kahan strictly more accurate than naive, and naive's error is "
           "non-trivial): %s\n",    // prints a report line
           pass ? "PASS" : "FAIL"); // continues the printf(...) call above
    return 0;                       // returns the computed result
}

N=10000000 additions of 0.1 (float)
exact reference (N * x, one double multiply): 1000000.0149
naive float sum:  1087937.0000  (error=87936.9851, 8.7937% of exact)
kahan float sum:  1000000.0000  (error=0.014901, 0.000001% of exact)
kahan is 5901351.2x more accurate than naive (smaller absolute error)
PASS (kahan strictly more accurate than naive, and naive's error is non-trivial): PASS

// learning/code/ex-68-fast-inverse-sqrt-bits/fisr.c
/* Example 68: the Quake III fast inverse square root -- bit-hacking a float's
 * IEEE-754 layout. */
#include <math.h>   // => sqrtf: the reference this bit-hack approximates
#include <stdint.h> // => uint32_t: the exact-width integer the float's 32 raw bits get copied into
#include <stdio.h>  // stdio.h: standard library header
#include <string.h> // => memcpy: the standards-legal way to reinterpret a float's bits (no UB, unlike
                    // a `union`-based type pun in strict C, and no aliasing violation like a cast)
 
// ex-68: the famous magic constant from Quake III Arena's `q_rsqrt` -- co-13:
// this is only meaningful because IEEE-754's bit layout (1 sign + 8 exponent +
// 23 mantissa) makes a float's raw integer bit pattern an approximately LINEAR
// function of log2(value). Halving that integer approximately halves
// log2(value), which approximately computes log2(1/sqrt(value)) =
// -0.5*log2(value).
#define MAGIC 0x5f3759dfu // constant MAGIC = 0x5f3759dfu
 
// ex-68: reinterpret the float's 32 bits as a uint32_t, shift, subtract from
// the magic constant, reinterpret back as float -- one Newton-Raphson
// refinement step polishes the crude bit-hack guess.
static float fast_inv_sqrt(float x) {  // defines fast_inv_sqrt(): helper function used by this example
    float x_half = 0.5f * x;           // => co-13: needed by the Newton-Raphson step below
    uint32_t bits;                     // declares bits
    memcpy(&bits, &x, sizeof(bits));   // => co-13: float bit pattern, as an integer
                                       // -- legal reinterpretation
    bits = MAGIC - (bits >> 1);        // => co-13: the bit-hack -- an approximate log-domain halving
    float y;                           // declares y
    memcpy(&y, &bits, sizeof(y));      // => co-13: reinterpret the hacked bits back as
                                       // a float -- crude guess
    y = y * (1.5f - (x_half * y * y)); // => one Newton-Raphson iteration: y_new = y*(1.5 - x/2*y^2),
                                       // refines the crude guess to within ~0.2% of the true value
    return y;                          // returns the computed result
}
 
int main(void) { // program entry point
    // ex-68: a spread of magnitudes -- co-13: the hack must work across the
    // exponent range, not just near 1.0, because it operates on the FULL
    // floating-point bit pattern including the exponent bits.
    float tests[] = {1.0f, 2.0f, 4.0f, 10.0f, 100.0f, 0.25f, 0.01f, 1000.0f}; // declares tests
    int n = (int)(sizeof(tests) / sizeof(tests[0]));                          // declares n
 
    double max_rel_error = 0.0;                                                  // declares max_rel_error
    printf("value      fast_inv_sqrt   1/sqrtf(value)   rel_error\n");           // prints a
                                                                                 // report line
    for (int i = 0; i < n; i++) {                                                // loop header controlling the sweep below
        float x = tests[i];                                                      // declares x
        float approx = fast_inv_sqrt(x);                                         // => co-13: the bit-hack's answer
        float exact = 1.0f / sqrtf(x);                                           // => the library's correctly-rounded answer
        double rel_error = fabs((double)approx - (double)exact) / (double)exact; // declares rel_error
        if (rel_error > max_rel_error)
            max_rel_error = rel_error;                                                                           // conditional check
        printf("%-10.4f %-15.6f %-16.6f %.6f%%\n", (double)x, (double)approx, (double)exact, rel_error * 100.0); // prints a report line
    }
 
    printf("\nmax relative error across all test values: %.4f%%\n",
           max_rel_error * 100.0); // prints a report line
    // ex-68: the well-known published tolerance for THIS exact magic-constant +
    // one-Newton-step combination is well under 1% relative error -- verify that
    // tolerance actually holds here.
    int pass = max_rel_error < 0.01; // < 1% relative error
    printf("PASS (max relative error under 1%% across the tested range): %s\n",
           pass ? "PASS" : "FAIL"); // prints a report line
    return 0;                       // returns the computed result
}