Single Instruction,
Multiple Data

/* ── PSHUFD: rearrange 32-bit integers (in-lane, immediate) ── */
// Immediate format: two 2-bit selectors per output lane (bits 0..7 → lanes 3..0)
// _MM_SHUFFLE(z, y, x, w) macro: builds the immediate
__m128i a = _mm_set_epi32(3, 2, 1, 0);
__m128i r = _mm_shuffle_epi32(a, _MM_SHUFFLE(0,1,2,3));
// r = [3,2,1,0] → [0,1,2,3] — reversed

/* ── SHUFPS: blend elements from two 128-bit float vectors ── */
__m128 fa = _mm_set_ps(3.f,2.f,1.f,0.f);
__m128 fb = _mm_set_ps(7.f,6.f,5.f,4.f);
__m128 fr = _mm_shuffle_ps(fa, fb, _MM_SHUFFLE(3,2,1,0));
// lower 2 lanes from fa, upper 2 from fb

/* ── PSHUFB (SSSE3): byte-granularity shuffle, in-lane, variable ── */
// Index byte 0..15: which source byte goes to each destination byte
// If index bit 7 = 1: output lane is zeroed
__m128i data  = _mm_set_epi8(15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0);
__m128i ctrl  = _mm_set_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
__m128i shuf  = _mm_shuffle_epi8(data, ctrl); // byte-reverse

/* ── AVX2 VPERMD: fully variable 32-bit element permute across 256b ── */
__m256i src = _mm256_set_epi32(7,6,5,4,3,2,1,0);
__m256i idx = _mm256_set_epi32(0,2,4,6,1,3,5,7); // pick any src lane
__m256i out = _mm256_permutevar8x32_epi32(src, idx);
// Crosses the 128-bit lane boundary! Indices 0..7 select from any lane in src.
// Very useful — this is the "full permute" instruction AVX2 was missing until Haswell.

/* ── Unpack / interleave ── */
__m128i lo = _mm_unpacklo_epi32(a, b); // [a0,b0,a1,b1]
__m128i hi = _mm_unpackhi_epi32(a, b); // [a2,b2,a3,b3]
// Unpack is the fundamental building block for transposing data (AoS ↔ SoA)

/* ── Blend (immediate) ── */
// BLENDPS: pick each float32 from a or b based on immediate bit
__m128 blended = _mm_blend_ps(fa, fb, 0b1010);
// bit 0=0: take fa[0], bit 1=1: take fb[1], etc.
// [fa[0], fb[1], fa[2], fb[3]]

// Input: 4 rows of 4 floats (row-major)
// Output: transposed in-place using xmm registers
static inline void transpose4x4(__m128 *r0, __m128 *r1, __m128 *r2, __m128 *r3) {
    // Step 1: interleave pairs of rows
    __m128 t0 = _mm_unpacklo_ps(*r0, *r1); // [r0[0],r1[0],r0[1],r1[1]]
    __m128 t1 = _mm_unpackhi_ps(*r0, *r1); // [r0[2],r1[2],r0[3],r1[3]]
    __m128 t2 = _mm_unpacklo_ps(*r2, *r3); // [r2[0],r3[0],r2[1],r3[1]]
    __m128 t3 = _mm_unpackhi_ps(*r2, *r3); // [r2[2],r3[2],r2[3],r3[3]]

    // Step 2: interleave pairs of pairs
    *r0 = _mm_movelh_ps(t0, t2);  // [r0[0],r1[0],r2[0],r3[0]] — column 0
    *r1 = _mm_movehl_ps(t2, t0);  // [r0[1],r1[1],r2[1],r3[1]] — column 1
    *r2 = _mm_movelh_ps(t1, t3);  // [r0[2],r1[2],r2[2],r3[2]] — column 2
    *r3 = _mm_movehl_ps(t3, t1);  // [r0[3],r1[3],r2[3],r3[3]] — column 3
}
// 8 instructions total. The macro _MM_TRANSPOSE4_PS does exactly this.

float dot_product_avx(const float* a, const float* b, int n) {
    __m256 acc0 = _mm256_setzero_ps();   // 4 accumulators to break FMA chain
    __m256 acc1 = _mm256_setzero_ps();   // dependency: each FMA waits ~4 cycles
    __m256 acc2 = _mm256_setzero_ps();   // using 4 hides that latency
    __m256 acc3 = _mm256_setzero_ps();

    int i = 0;
    for (; i <= n - 32; i += 32) {   // process 32 floats per iteration
        acc0 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i),
                                _mm256_loadu_ps(b+i), acc0);
        acc1 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+ 8),
                                _mm256_loadu_ps(b+i+ 8), acc1);
        acc2 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+16),
                                _mm256_loadu_ps(b+i+16), acc2);
        acc3 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+24),
                                _mm256_loadu_ps(b+i+24), acc3);
    }

    /* Combine accumulators */
    acc0 = _mm256_add_ps(acc0, acc1);
    acc2 = _mm256_add_ps(acc2, acc3);
    acc0 = _mm256_add_ps(acc0, acc2);

    /* ── Horizontal reduction: 8 floats → 1 float ── */
    // Step 1: add upper 128b half to lower half
    __m128 hi  = _mm256_extractf128_ps(acc0, 1); // lanes 4..7
    __m128 lo  = _mm256_castps256_ps128(acc0);     // lanes 0..3 (zero-cost cast)
    __m128 sum = _mm_add_ps(lo, hi);    // [0+4, 1+5, 2+6, 3+7]

    // Step 2: hadd pairs
    sum = _mm_hadd_ps(sum, sum);         // [(0+4)+(1+5), (2+6)+(3+7), ...]
    sum = _mm_hadd_ps(sum, sum);         // [total, ...]

    float result = _mm_cvtss_f32(sum);   // extract lowest lane as scalar

    /* Scalar tail for remaining elements */
    for (; i < n; i++) result += a[i] * b[i];

    return result;
}

// Convert ASCII uppercase to lowercase for 16 bytes simultaneously
void to_lower_sse(char* s, int n) {
    __m128i upper_A  = _mm_set1_epi8('A' - 1);   // 0x40
    __m128i upper_Z  = _mm_set1_epi8('Z');       // 0x5A
    __m128i lo_delta = _mm_set1_epi8(32);        // 'a' - 'A' = 32

    int i = 0;
    for (; i <= n - 16; i += 16) {
        __m128i chunk = _mm_loadu_si128((__m128i_u*)(s + i));

        // Compute mask: byte is uppercase if 'A'-1 < byte <= 'Z'
        __m128i above_A = _mm_cmpgt_epi8(chunk, upper_A);
        __m128i below_Z = _mm_cmpgt_epi8(upper_Z, chunk); // Z > chunk means chunk < Z

        // Note: _mm_cmple_epi8 doesn't exist, use Z >= chunk via cmpgt(Z+1, chunk)
        __m128i is_upper = _mm_and_si128(above_A, below_Z); // &: both conditions

        // Only add 32 to uppercase bytes (AND mask with delta)
        __m128i delta = _mm_and_si128(is_upper, lo_delta);
        chunk = _mm_add_epi8(chunk, delta);

        _mm_storeu_si128((__m128i_u*)(s + i), chunk);
    }
    for (; i < n; i++) s[i] |= (tolower((unsigned char)s[i]));
}

# GCC / Clang — enable auto-vectorization and target ISA
gcc -O2 -march=native -fvectorize -ftree-vectorize myfile.c
# -O2 enables basic vectorization; -O3 enables more aggressive
# -march=native: use all features of the current CPU
# -march=x86-64-v3: target AVX2+FMA baseline (common server minimum)
# -march=x86-64-v4: target AVX-512

# Get vectorization report (what was/wasn't vectorized and why)
gcc -O2 -march=native -fopt-info-vec-optimized -fopt-info-vec-missed myfile.c
clang -O2 -march=native -Rpass=loop-vectorize -Rpass-missed=loop-vectorize myfile.c

# Enable FMA and AVX2 explicitly
gcc -O2 -mavx2 -mfma myfile.c

# MSVC
cl /O2 /arch:AVX2 /Qvec-report:2 myfile.c

/* ── Use restrict to eliminate aliasing concerns ── */
void add_arrays(float *restrict c,
                const float *restrict a,
                const float *restrict b,
                int n) {
    for (int i = 0; i < n; i++) c[i] = a[i] + b[i];  // vectorizes cleanly
}

/* ── Alignment hints ── */
__attribute__((aligned(32))) float buf[1024];  // static aligned array
// or: void f(float *__restrict__ __attribute__((aligned(32))) p);

/* ── Loop that auto-vectorizes to AVX2 (with -mavx2 -mfma) ── */
void saxpy(float *restrict y, const float *restrict x, float a, int n) {
    for (int i = 0; i < n; i++) y[i] = a * x[i] + y[i];  // → VFMADD231PS
}

/* ── Loop that does NOT auto-vectorize (loop-carried dependency) ── */
void prefix_sum(float *a, int n) {
    for (int i = 1; i < n; i++) a[i] += a[i-1]; // NOT vectorizable as-is
    // Vectorizing prefix sums requires a different algorithm (parallel scan)
}

/* ── GCC pragma to force/prevent vectorization ── */
#pragma GCC ivdep              // "ignore vector dependencies" — trust me, no deps
for (int i = 0; i < n; i++) a[i] = b[i] + c[i];

#pragma GCC unroll 4           // manually unroll loop (helps expose ILP to vectorizer)
for (int i = 0; i < n; i++) a[i] *= 2.f;

/* Clang OpenMP SIMD (portable) */
#pragma omp simd
for (int i = 0; i < n; i++) c[i] = a[i] * b[i];

#include <immintrin.h>

/* ── Flush denormals to zero — critical for audio/ML performance ── */
static inline void disable_denormals(void) {
    // MXCSR bits: bit 15 = FZ (Flush-to-Zero), bit 6 = DAZ (Denormals-Are-Zero)
    _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);       // outputs: flush to zero
    _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); // inputs: treat as zero
    // Equivalent: _mm_setcsr(_mm_getcsr() | 0x8040);
}

/* Check if a float value is denormal */
static inline bool is_denormal(float x) {
    uint32_t bits;
    memcpy(&bits, &x, 4);
    return ((bits & 0x7F800000) == 0) && ((bits & 0x007FFFFF) != 0);
}

// WRONG: undefined behavior (strict aliasing violation)
float f = 1.0f;
uint32_t i = *(uint32_t*)&f;  // ← DO NOT DO THIS

// CORRECT: use memcpy
float f2 = 1.0f;
uint32_t i2;
memcpy(&i2, &f2, 4);  // compiles to a single register move

// ALSO CORRECT: C++ std::bit_cast (C++20)
uint32_t i3 = std::bit_cast<uint32_t>(f2);

// ALSO CORRECT: union (C only, not C++)
union { float f; uint32_t i; } pun;
pun.f = 1.0f;
uint32_t i4 = pun.i;  // valid in C99, undefined in C++ (but works in practice)

#include <arm_neon.h>

/* ── NEON types (analogous to Intel __m128) ── */
float32x4_t a, b, c;   // 4×float32 (128-bit)
int16x8_t   ia, ib;   // 8×int16 (128-bit)
uint8x16_t  ua;        // 16×uint8 (128-bit)
float32x2_t da;        // 2×float32 (64-bit, "doubleword")

/* ── Load / Store ── */
a = vld1q_f32(ptr);              // unaligned 128-bit load (aligned auto-selected)
vst1q_f32(ptr, a);               // store
a = vdupq_n_f32(1.0f);           // broadcast (dup = duplicate)
a = vld1q_dup_f32(ptr);          // load and broadcast

/* Multi-register load (de-interleave) — unique NEON feature */
float32x4x3_t rgb = vld3q_f32(ptr);  // loads 12 floats, de-interleaves RGB
// rgb.val[0] = R channel, rgb.val[1] = G, rgb.val[2] = B
// This is the NEON killer feature for AoS→SoA — no manual transpose needed!

/* ── Arithmetic ── */
c = vaddq_f32(a, b);             // c[i] = a[i] + b[i]
c = vmulq_f32(a, b);             // c[i] = a[i] * b[i]
c = vmlaq_f32(c, a, b);          // c[i] = c[i] + a[i]*b[i]  (mul-accumulate)
c = vmaxq_f32(a, b);             // c[i] = max(a[i], b[i])
c = vrecpeq_f32(a);             // 1/a[i] approximation

/* FMA (AArch64 only, not ARMv7 NEON) */
c = vfmaq_f32(c, a, b);          // c[i] = c[i] + a[i]*b[i] (true FMA)
c = vfmaq_lane_f32(c, a, da, 0); // multiply each a[i] by da[0], accumulate

/* ── Comparison ── */
uint32x4_t mask = vcgtq_f32(a, b);  // a > b, returns 0xFFFFFFFF or 0
c = vbslq_f32(mask, a, b);          // blend: c[i] = mask[i] ? a[i] : b[i]

/* ── Shuffle / permute ── */
uint8x16_t tbl = vtbl1q_u8(src, idx);  // table lookup (like x86 PSHUFB)
c = vextq_f32(a, b, 2);               // [a2,a3,b0,b1] — extract and shift
c = vrev64q_f32(a);                   // reverse pairs: [a1,a0,a3,a2]
c = vtrn1q_f32(a, b);                 // [a0,b0,a2,b2] — transpose even

/* Horizontal reduce */
float total = vaddvq_f32(a);  // sum all 4 lanes → scalar (AArch64)

#include <arm_sve.h>

void saxpy_sve(float* y, const float* x, float a, uint64_t n) {
    // svcntw() = number of float32 lanes — determined at RUNTIME by the hardware
    // On 128-bit SVE: svcntw() = 4
    // On 512-bit SVE: svcntw() = 16
    // Your code doesn't change.
    svfloat32_t va = svdup_n_f32(a);       // broadcast a to all SVE lanes

    uint64_t i = 0;
    for (; i < n; i += svcntw()) {
        // svwhilelt_b32: create predicate mask for i .. i+vl-1, limited to n
        svbool_t pg = svwhilelt_b32(i, n);   // handles loop tail automatically!

        svfloat32_t vx = svld1_f32(pg, x + i);  // masked load
        svfloat32_t vy = svld1_f32(pg, y + i);  // masked load

        vy = svmla_f32_z(pg, vy, va, vx);      // y = a*x + y (FMA), zeroed mask

        svst1_f32(pg, y + i, vy);               // masked store (safe at tail)
    }
    // No scalar tail needed! svwhilelt automatically limits the last iteration.
}
// This single function runs at optimal width on any SVE implementation:
// Cortex-A55 (128b), Fujitsu A64FX (512b), AWS Graviton, Apple M-series, etc.