Cache-SIMD-DSA

Problem:

For every index, return the product of all other elements without division.

Diagram:

naive:
  pass 1 prefix: [x|cold...][x|cold...][x|cold...] -> work[]
  pass 2 suffix: [x|cold...][x|cold...][x|cold...] -> work[]
  issue: two passes repeatedly drag cold record fields through cache

cache-aware:
  pass 1 prefix: x[] -> work[]
  pass 2 suffix: x[] -> work[]
  change: same math, dense input stream

SIMD:
  prefix/suffix are vector lanes, but each x comes from gather(AoS)

cache-aware + SIMD:
  prefix/suffix are vector lanes and each x comes from packed load(SoA)

Code:

// Naive: prefix/suffix over wide records. Each pass pulls cold fields into cache.
// Major difference: the algorithm scans the input twice, so AoS waste is paid twice.
NOINLINE NOVECTOR std::uint64_t product_naive_aos(const WideRecord<std::uint32_t>* aos,
                                                  std::size_t length) {
    std::vector<std::uint32_t> work(length);
    std::uint64_t total = 0;
    for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
        const WideRecord<std::uint32_t>* values = aos + lane * length;
        std::uint32_t prefix = 1U;
        for (std::size_t index = 0; index < length; ++index) {
            work[index] = prefix;
            // Cost center: values[index].value is useful, but the wide record
            // brings cold fields into the cache on the prefix pass.
            prefix *= values[index].value;
        }
        std::uint32_t suffix = 1U;
        for (std::size_t index = length; index-- > 0;) {
            work[index] *= suffix;
            // Same wide-record penalty again on the reverse suffix pass.
            suffix *= values[index].value;
        }
        total += sample_checksum(work.data(), length);
    }
    return total;
}

// Cache-aware: same prefix/suffix math over dense values.
// Major difference: only the representation changes; prefix/suffix semantics match naive.
NOINLINE NOVECTOR std::uint64_t product_cache_soa(const std::uint32_t* hot,
                                                  std::size_t length) {
    std::vector<std::uint32_t> work(length);
    std::uint64_t total = 0;
    for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
        const std::uint32_t* values = hot + lane * length;
        std::uint32_t prefix = 1U;
        for (std::size_t index = 0; index < length; ++index) {
            work[index] = prefix;
            // Dense linear load: each cache line contains many useful values.
            prefix *= values[index];
        }
        std::uint32_t suffix = 1U;
        for (std::size_t index = length; index-- > 0;) {
            work[index] *= suffix;
            // Reverse scan is still linear over compact data.
            suffix *= values[index];
        }
        total += sample_checksum(work.data(), length);
    }
    return total;
}

// SIMD: vectorized prefix/suffix, but values are gathered from wide records.
// Major difference from cache-aware + SIMD: gather keeps the AoS memory penalty.
NOINLINE std::uint64_t product_simd_aos(const WideRecord<std::uint32_t>* aos,
                                        std::size_t length,
                                        std::uint32_t* temp) {
    __m256i prefix = _mm256_set1_epi32(1);
    for (std::size_t index = 0; index < length; ++index) {
        _mm256_store_si256(reinterpret_cast<__m256i*>(temp + index * kBatchLanes), prefix);
        // Eight lanes advance together, but every lane fetches from a different wide record.
        prefix = _mm256_mullo_epi32(prefix, gather_hot_u32_values(aos, length, index));
    }

    __m256i suffix = _mm256_set1_epi32(1);
    for (std::size_t index = length; index-- > 0;) {
        const __m256i prefix_values =
            _mm256_load_si256(reinterpret_cast<const __m256i*>(temp + index * kBatchLanes));
        _mm256_store_si256(reinterpret_cast<__m256i*>(temp + index * kBatchLanes),
                           _mm256_mullo_epi32(prefix_values, suffix));
        // Reverse pass repeats the gather-heavy access pattern.
        suffix = _mm256_mullo_epi32(suffix, gather_hot_u32_values(aos, length, index));
    }
    return sample_checksum(temp, length * kBatchLanes);
}

// Cache-aware + SIMD: vectorized prefix/suffix over packed lanes.
// Major difference from plain SIMD: hot values are already interleaved for AVX2.
NOINLINE std::uint64_t product_cache_simd_soa(const std::uint32_t* hot,
                                              std::size_t length,
                                              std::uint32_t* temp) {
    __m256i prefix = _mm256_set1_epi32(1);
    for (std::size_t index = 0; index < length; ++index) {
        _mm256_store_si256(reinterpret_cast<__m256i*>(temp + index * kBatchLanes), prefix);
        // Packed load: one 32-byte read supplies the next value for all eight lanes.
        prefix = _mm256_mullo_epi32(
            prefix,
            _mm256_load_si256(reinterpret_cast<const __m256i*>(hot + index * kBatchLanes)));
    }

    __m256i suffix = _mm256_set1_epi32(1);
    for (std::size_t index = length; index-- > 0;) {
        const __m256i prefix_values =
            _mm256_load_si256(reinterpret_cast<const __m256i*>(temp + index * kBatchLanes));
        _mm256_store_si256(reinterpret_cast<__m256i*>(temp + index * kBatchLanes),
                           _mm256_mullo_epi32(prefix_values, suffix));
        // The suffix pass stays packed too, so SIMD and cache layout reinforce each other.
        suffix = _mm256_mullo_epi32(
            suffix,
            _mm256_load_si256(reinterpret_cast<const __m256i*>(hot + index * kBatchLanes)));
    }
    return sample_checksum(temp, length * kBatchLanes);
}

Result:

Evidence says cache-aware + SIMD wins. I expected the prefix/suffix passes to punish the wide AoS record because every pass rereads cold payload. Timing supports that. Perf supports the layout part through lower generic cache-miss pressure for cache-aware rows; the SIMD win is mostly throughput from doing more arithmetic per pass.

Problem:

Given bar heights, compute how much water is trapped.

Diagram:

naive:
  left pass:  [h|cold...][h|cold...][h|cold...] -> left_max[]
  right pass: [h|cold...][h|cold...][h|cold...] -> water sum
  issue: both passes reread wide records

cache-aware:
  left pass:  h[] -> left_max[]
  right pass: h[] -> water sum
  change: same two-pass method over dense heights

SIMD:
  vector left/right maxima, but h values come from gathers

cache-aware + SIMD:
  packed h vectors -> packed left_max[] -> packed water sum

Code:

// Naive: build left maxima and scan right maxima over wide records.
// Major difference: the algorithm touches heights twice, so AoS waste is repeated.
NOINLINE NOVECTOR float rain_naive_aos(const WideRecord<float>* aos, std::size_t length) {
    std::vector<float> left(length);
    float total = 0.0f;
    for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
        const WideRecord<float>* values = aos + lane * length;
        float left_max = values[0].value;
        for (std::size_t index = 0; index < length; ++index) {
            // Forward pass: each height load also pulls cold fields into cache.
            left_max = std::max(left_max, values[index].value);
            left[index] = left_max;
        }
        float right_max = values[length - 1].value;
        for (std::size_t index = length; index-- > 0;) {
            // Reverse pass: the same wide-record pattern is paid again.
            right_max = std::max(right_max, values[index].value);
            total += std::max(0.0f, std::min(left[index], right_max) - values[index].value);
        }
    }
    return total;
}

// Cache-aware: same two-pass idea over dense heights.
// Major difference from naive: left/right scans are linear reads of hot heights.
NOINLINE NOVECTOR float rain_cache_soa(const float* hot, std::size_t length) {
    std::vector<float> left(length);
    float total = 0.0f;
    for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
        const float* values = hot + lane * length;
        float left_max = values[0];
        for (std::size_t index = 0; index < length; ++index) {
            // Dense forward pass: left[] still costs memory, but input no longer wastes cache.
            left_max = std::max(left_max, values[index]);
            left[index] = left_max;
        }
        float right_max = values[length - 1];
        for (std::size_t index = length; index-- > 0;) {
            // Dense reverse pass: right_max and water use contiguous heights.
            right_max = std::max(right_max, values[index]);
            total += std::max(0.0f, std::min(left[index], right_max) - values[index]);
        }
    }
    return total;
}

// SIMD: vectorized left/right maxima, but heights are gathered from AoS.
// Major difference: vector math is wide, but both passes still gather from wide records.
NOINLINE float rain_simd_aos(const WideRecord<float>* aos, std::size_t length, float* left_temp) {
    __m256 left_max = _mm256_set1_ps(std::numeric_limits<float>::lowest());
    for (std::size_t index = 0; index < length; ++index) {
        // Gathered heights feed eight left-max streams.
        left_max = _mm256_max_ps(left_max, gather_hot_values(aos, length, index));
        _mm256_store_ps(left_temp + index * kBatchLanes, left_max);
    }

    __m256 right_max = _mm256_set1_ps(std::numeric_limits<float>::lowest());
    __m256 total = _mm256_setzero_ps();
    for (std::size_t index = length; index-- > 0;) {
        // Reverse pass gathers again, so SIMD cannot hide the poor layout.
        const __m256 values = gather_hot_values(aos, length, index);
        right_max = _mm256_max_ps(right_max, values);
        const __m256 left_values = _mm256_load_ps(left_temp + index * kBatchLanes);
        const __m256 water = _mm256_sub_ps(_mm256_min_ps(left_values, right_max), values);
        total = _mm256_add_ps(total, _mm256_max_ps(water, _mm256_setzero_ps()));
    }
    return horizontal_sum_ps(total);
}

// Cache-aware + SIMD: vectorized two-pass scan over packed heights.
// Major difference from SIMD AoS: every height vector is a packed load.
NOINLINE float rain_cache_simd_soa(const float* hot, std::size_t length, float* left_temp) {
    __m256 left_max = _mm256_set1_ps(std::numeric_limits<float>::lowest());
    for (std::size_t index = 0; index < length; ++index) {
        // Packed forward pass: one vector load updates eight left maxima.
        const __m256 values = _mm256_load_ps(hot + index * kBatchLanes);
        left_max = _mm256_max_ps(left_max, values);
        _mm256_store_ps(left_temp + index * kBatchLanes, left_max);
    }

    __m256 right_max = _mm256_set1_ps(std::numeric_limits<float>::lowest());
    __m256 total = _mm256_setzero_ps();
    for (std::size_t index = length; index-- > 0;) {
        // Packed reverse pass: this is where cache-aware layout and SIMD combine.
        const __m256 values = _mm256_load_ps(hot + index * kBatchLanes);
        right_max = _mm256_max_ps(right_max, values);
        const __m256 left_values = _mm256_load_ps(left_temp + index * kBatchLanes);
        const __m256 water = _mm256_sub_ps(_mm256_min_ps(left_values, right_max), values);
        total = _mm256_add_ps(total, _mm256_max_ps(water, _mm256_setzero_ps()));
    }
    return horizontal_sum_ps(total);
}

Result:

Evidence says cache-aware + SIMD wins. I expected multiple memory passes to reward contiguous hot arrays. Timing supports that. Perf supports the layout expectation for the scalar rows; the SIMD AoS timing shows that vectorizing gathers over wide records is not enough.

Problem:

Search target values in a sorted array that has been rotated.

Diagram:

naive:
  query q:
    left ---- mid ---- right
        branch decides which half survives
  issue: each probe reads values from wide records

cache-aware:
  same binary-search decisions, but probes read dense values[]

SIMD:
  lane0: left -> right half
  lane1: right -> left half
  lane2: found
  lane3: still active
  issue: lanes diverge, masks update per lane, gathers still happen

cache-aware + SIMD:
  dense/interleaved values reduce memory waste, but branchy lane divergence remains

Code:

// Scalar helper used by both scalar layouts.
// Major behavior: every iteration depends on branch decisions from previous probes.
template <typename LoadAt>
NOINLINE NOVECTOR int rotated_search_scalar_value(std::size_t length,
                                                  std::uint32_t target,
                                                  LoadAt&& load_at) {
    int left = 0;
    int right = static_cast<int>(length) - 1;
    while (left <= right) {
        // The branch state is sequential: mid, left, and right change after each decision.
        const int mid = left + (right - left) / 2;
        const std::uint32_t left_value = load_at(left);
        const std::uint32_t mid_value = load_at(mid);
        const std::uint32_t right_value = load_at(right);
        if (mid_value == target) {
            return mid;
        }
        if (left_value <= mid_value) {
            if (left_value <= target && target < mid_value) {
                right = mid - 1;
            } else {
                left = mid + 1;
            }
        } else if (mid_value < target && target <= right_value) {
            left = mid + 1;
        } else {
            right = mid - 1;
        }
    }
    return -1;
}

// Naive: binary search over wide records for each lane/query.
// Major difference: scalar control flow is already hard to predict, and AoS adds
// wasted memory traffic to every left/mid/right probe.
NOINLINE NOVECTOR std::uint64_t rotated_search_naive_aos(const WideRecord<std::uint32_t>* aos,
                                                         const std::uint32_t* queries,
                                                         std::size_t length,
                                                         std::size_t query_count) {
    std::uint64_t sum = 0;
    for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
        const WideRecord<std::uint32_t>* values = aos + lane * length;
        const std::uint32_t* lane_queries = queries + lane * query_count;
        for (std::size_t query = 0; query < query_count; ++query) {
            // The load lambda exposes the AoS penalty at every binary-search probe.
            sum += static_cast<std::uint32_t>(
                rotated_search_scalar_value(length, lane_queries[query],
                                            [values](int i) { return values[i].value; }) + 1);
        }
    }
    return sum;
}

// Cache-aware: same control flow, but array and queries are dense.
// Major difference from naive: same branches, lower memory waste.
NOINLINE NOVECTOR std::uint64_t rotated_search_cache_soa(const std::uint32_t* hot,
                                                         const std::uint32_t* queries,
                                                         std::size_t length,
                                                         std::size_t query_count) {
    std::uint64_t sum = 0;
    for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
        const std::uint32_t* values = hot + lane * length;
        const std::uint32_t* lane_queries = queries + lane * query_count;
        for (std::size_t query = 0; query < query_count; ++query) {
            // The measured winner: scalar branch behavior remains, but probes are compact.
            sum += static_cast<std::uint32_t>(
                rotated_search_scalar_value(length, lane_queries[query],
                                            [values](int i) { return values[i]; }) + 1);
        }
    }
    return sum;
}

// SIMD helper: eight searches advance together, but lanes diverge quickly.
// Major difference: SIMD is fighting control flow. Each lane may choose a different half.
template <typename GatherFn>
std::uint64_t rotated_search_simd_common(const std::uint32_t* interleaved_queries,
                                         std::size_t length,
                                         std::size_t query_count,
                                         GatherFn&& gather_values) {
    const __m256i zero = _mm256_setzero_si256();
    const __m256i one = _mm256_set1_epi32(1);
    const __m256i all_bits = _mm256_set1_epi32(-1);
    const __m256i max_index = _mm256_set1_epi32(static_cast<int>(length - 1));
    const std::size_t max_iterations = std::bit_width(length) + 2;
    std::uint64_t sum = 0;
    alignas(32) std::array<std::int32_t, kBatchLanes> result_lanes{};

    for (std::size_t query = 0; query < query_count; ++query) {
        const __m256i target =
            _mm256_load_si256(reinterpret_cast<const __m256i*>(interleaved_queries + query * kBatchLanes));
        __m256i left = zero;
        __m256i right = max_index;
        __m256i result = _mm256_set1_epi32(-1);
        __m256i active = all_bits;

        for (std::size_t iteration = 0; iteration < max_iterations && mask_bits(active) != 0; ++iteration) {
            // Inactive lanes are masked out, but the vector loop still carries their state.
            const __m256i safe_left = _mm256_blendv_epi8(zero, left, active);
            const __m256i safe_right = _mm256_blendv_epi8(zero, right, active);
            const __m256i mid = _mm256_srli_epi32(_mm256_add_epi32(safe_left, safe_right), 1);
            // This probe may be a gather in both SIMD variants because every lane has
            // a different mid after a few iterations.
            const __m256i mid_values = gather_values(mid);
            const __m256i equal_mask = _mm256_and_si256(active, _mm256_cmpeq_epi32(mid_values, target));
            result = _mm256_blendv_epi8(result, mid, equal_mask);
            active = _mm256_andnot_si256(equal_mask, active);
            if (mask_bits(active) == 0) {
                break;
            }

            const __m256i left_values = gather_values(safe_left);
            const __m256i right_values = gather_values(safe_right);
            const __m256i left_sorted = cmp_le_epi32(left_values, mid_values);
            // These masks implement the scalar if/else tree. This is the major
            // reason SIMD loses here: decisions are per lane, not one shared branch.
            const __m256i target_in_left =
                _mm256_and_si256(cmp_le_epi32(left_values, target), cmp_lt_epi32(target, mid_values));
            const __m256i target_in_right =
                _mm256_and_si256(cmp_lt_epi32(mid_values, target), cmp_le_epi32(target, right_values));
            const __m256i go_right = _mm256_and_si256(
                active,
                _mm256_or_si256(
                    _mm256_and_si256(left_sorted, _mm256_andnot_si256(target_in_left, all_bits)),
                    _mm256_andnot_si256(left_sorted, target_in_right)));
            const __m256i go_left = _mm256_and_si256(active, _mm256_andnot_si256(go_right, all_bits));

            left = _mm256_blendv_epi8(left, _mm256_add_epi32(mid, one), go_right);
            right = _mm256_blendv_epi8(right, _mm256_sub_epi32(mid, one), go_left);
            active = _mm256_and_si256(active, cmp_le_epi32(left, right));
        }

        _mm256_store_si256(reinterpret_cast<__m256i*>(result_lanes.data()), result);
        for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
            sum += static_cast<std::uint32_t>(result_lanes[lane] + 1);
        }
    }
    return sum;
}

// SIMD: vectorized search with AoS gathers.
// Major difference from scalar AoS: eight queries run together, but all probes are gathered.
NOINLINE std::uint64_t rotated_search_simd_aos(const WideRecord<std::uint32_t>* aos,
                                               const std::uint32_t* interleaved_queries,
                                               std::size_t length,
                                               std::size_t query_count) {
    return rotated_search_simd_common(
        interleaved_queries,
        length,
        query_count,
        [aos, length](__m256i logical_indices) {
            alignas(32) std::array<std::int32_t, kBatchLanes> indices{};
            _mm256_store_si256(reinterpret_cast<__m256i*>(indices.data()), logical_indices);
            for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
                // Convert logical lane indexes into AoS byte/word positions.
                indices[lane] = indices[lane] * 16 + static_cast<std::int32_t>(lane * length * 16);
            }
            // Gather from wide records. This combines branch divergence with scattered memory.
            return _mm256_i32gather_epi32(reinterpret_cast<const int*>(aos),
                                          _mm256_load_si256(reinterpret_cast<const __m256i*>(indices.data())),
                                          4);
        });
}

// Cache-aware + SIMD: vectorized search over interleaved dense lanes.
// Major difference from SIMD AoS: compact data, but still one scattered probe per lane.
NOINLINE std::uint64_t rotated_search_cache_simd_soa(const std::uint32_t* interleaved_hot,
                                                     const std::uint32_t* interleaved_queries,
                                                     std::size_t length,
                                                     std::size_t query_count) {
    return rotated_search_simd_common(
        interleaved_queries,
        length,
        query_count,
        [interleaved_hot](__m256i logical_indices) {
            alignas(32) std::array<std::int32_t, kBatchLanes> indices{};
            _mm256_store_si256(reinterpret_cast<__m256i*>(indices.data()), logical_indices);
            for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
                // Interleaving reduces stride versus AoS, but per-lane mids still differ.
                indices[lane] = indices[lane] * static_cast<std::int32_t>(kBatchLanes) + static_cast<std::int32_t>(lane);
            }
            // This is cache-friendlier than AoS gather, but not as good as scalar SoA here.
            return _mm256_i32gather_epi32(reinterpret_cast<const int*>(interleaved_hot),
                                          _mm256_load_si256(reinterpret_cast<const __m256i*>(indices.data())),
                                          4);
        });
}

Result:

Evidence says cache-aware scalar wins. I expected branchy per-query control flow to limit SIMD usefulness, and the timing verifies that. Perf supports the scalar SoA verdict: it has better IPC than naive and much lower generic cache-miss and dTLB-miss rates than naive.

Problem:

Traverse graph vertices in breadth-first order.

Diagram:

naive adjacency list:
  vertex u -> [(dst|cold cold cold), (dst|cold cold cold), ...]
  issue: each neighbor visit carries unused edge fields

cache-aware CSR:
  offsets: [0, 8, 16, 24, ...]
  edges:   [dst dst dst dst dst dst ...]
  change: neighbors are compact and contiguous

SIMD:
  gather dst values + gather visited[dst]
  issue: enqueue/visited updates still happen lane by lane

cache-aware + SIMD:
  load contiguous CSR dst values, but visited checks still serialize

Code:

// Naive: adjacency-list records carry cold fields with each destination.
// Major difference: each edge record has more data than BFS needs.
NOINLINE NOVECTOR std::uint64_t bfs_naive_adjlist(const GraphAdjList& graph, std::size_t passes) {
    const std::size_t vertices = graph.edges.size();
    std::vector<std::uint32_t> queue(vertices);
    std::vector<int> visited(vertices, 0);
    std::uint64_t total = 0;
    for (std::size_t pass = 0; pass < passes; ++pass) {
        std::fill(visited.begin(), visited.end(), 0);
        const std::uint32_t source = static_cast<std::uint32_t>((pass * 17) % vertices);
        std::size_t head = 0;
        std::size_t tail = 0;
        visited[source] = 1;
        queue[tail++] = source;
        while (head < tail) {
            const std::uint32_t node = queue[head++];
            total += node;
            for (const GraphEdgeRecord& edge : graph.edges[node]) {
                // BFS only needs edge.dst; cold0/cold1/cold2 still traveled with it.
                if (visited[edge.dst] == 0) {
                    // This visited check is data dependent and controls queue growth.
                    visited[edge.dst] = 1;
                    queue[tail++] = edge.dst;
                }
            }
        }
        total += tail;
    }
    return total;
}

// Cache-aware: CSR stores only destination IDs in one dense edge array.
// Major difference from naive: adjacency data is offsets[] + compact edges[].
NOINLINE NOVECTOR std::uint64_t bfs_cache_csr(const GraphCsr& graph, std::size_t passes) {
    const std::size_t vertices = graph.offsets.size() - 1;
    std::vector<std::uint32_t> queue(vertices);
    std::vector<int> visited(vertices, 0);
    std::uint64_t total = 0;
    for (std::size_t pass = 0; pass < passes; ++pass) {
        std::fill(visited.begin(), visited.end(), 0);
        const std::uint32_t source = static_cast<std::uint32_t>((pass * 17) % vertices);
        std::size_t head = 0;
        std::size_t tail = 0;
        visited[source] = 1;
        queue[tail++] = source;
        while (head < tail) {
            const std::uint32_t node = queue[head++];
            total += node;
            // The neighbor range is contiguous in graph.edges.
            for (std::uint32_t edge = graph.offsets[node]; edge < graph.offsets[node + 1]; ++edge) {
                const std::uint32_t next = graph.edges[edge];
                if (visited[next] == 0) {
                    // Still branchy, but now the edge stream is cache-friendly.
                    visited[next] = 1;
                    queue[tail++] = next;
                }
            }
        }
        total += tail;
    }
    return total;
}

// SIMD: gather destination and visited state for adjacency-list records.
// Major difference: tries to batch neighbor loads, but cannot batch queue mutation cleanly.
NOINLINE std::uint64_t bfs_simd_adjlist(const GraphAdjList& graph, std::size_t passes) {
    const std::size_t vertices = graph.edges.size();
    std::vector<std::uint32_t> queue(vertices);
    std::vector<int> visited(vertices, 0);
    std::uint64_t total = 0;
    alignas(32) std::array<std::uint32_t, kBatchLanes> dst_lanes{};
    constexpr std::int32_t kStrideWords = static_cast<std::int32_t>(sizeof(GraphEdgeRecord) / sizeof(std::uint32_t));
    const __m256i zero = _mm256_setzero_si256();

    for (std::size_t pass = 0; pass < passes; ++pass) {
        std::fill(visited.begin(), visited.end(), 0);
        const std::uint32_t source = static_cast<std::uint32_t>((pass * 17) % vertices);
        std::size_t head = 0;
        std::size_t tail = 0;
        visited[source] = 1;
        queue[tail++] = source;
        while (head < tail) {
            const std::uint32_t node = queue[head++];
            total += node;
            const auto& edges = graph.edges[node];
            const int* base = reinterpret_cast<const int*>(edges.data());
            std::size_t edge = 0;
            for (; edge + kBatchLanes <= edges.size(); edge += kBatchLanes) {
                // Gather destination fields from records that also contain cold fields.
                const __m256i edge_indices = _mm256_setr_epi32(static_cast<int>((edge + 0) * kStrideWords),
                                                               static_cast<int>((edge + 1) * kStrideWords),
                                                               static_cast<int>((edge + 2) * kStrideWords),
                                                               static_cast<int>((edge + 3) * kStrideWords),
                                                               static_cast<int>((edge + 4) * kStrideWords),
                                                               static_cast<int>((edge + 5) * kStrideWords),
                                                               static_cast<int>((edge + 6) * kStrideWords),
                                                               static_cast<int>((edge + 7) * kStrideWords));
                const __m256i next_nodes =
                    _mm256_i32gather_epi32(base, edge_indices, sizeof(std::uint32_t));
                // Second gather: visited[] is indexed by arbitrary destination IDs.
                const __m256i seen =
                    _mm256_i32gather_epi32(visited.data(), next_nodes, static_cast<int>(sizeof(int)));
                _mm256_store_si256(reinterpret_cast<__m256i*>(dst_lanes.data()), next_nodes);
                const int unseen_mask = mask_bits(_mm256_cmpeq_epi32(seen, zero));
                // Major limiter: queue insertion and duplicate checks remain lane-by-lane.
                for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
                    if (((unseen_mask >> lane) & 1) != 0 && visited[dst_lanes[lane]] == 0) {
                        visited[dst_lanes[lane]] = 1;
                        queue[tail++] = dst_lanes[lane];
                    }
                }
            }
            for (; edge < edges.size(); ++edge) {
                const std::uint32_t next = edges[edge].dst;
                if (visited[next] == 0) {
                    visited[next] = 1;
                    queue[tail++] = next;
                }
            }
        }
        total += tail;
    }
    return total;
}

// Cache-aware + SIMD: contiguous CSR edge loads, but visited checks still serialize.
// Major difference from SIMD adjacency list: destination load is dense, visited gather remains.
NOINLINE std::uint64_t bfs_cache_simd_csr(const GraphCsr& graph, std::size_t passes) {
    const std::size_t vertices = graph.offsets.size() - 1;
    std::vector<std::uint32_t> queue(vertices);
    std::vector<int> visited(vertices, 0);
    std::uint64_t total = 0;
    alignas(32) std::array<std::uint32_t, kBatchLanes> dst_lanes{};
    const __m256i zero = _mm256_setzero_si256();

    for (std::size_t pass = 0; pass < passes; ++pass) {
        std::fill(visited.begin(), visited.end(), 0);
        const std::uint32_t source = static_cast<std::uint32_t>((pass * 17) % vertices);
        std::size_t head = 0;
        std::size_t tail = 0;
        visited[source] = 1;
        queue[tail++] = source;
        while (head < tail) {
            const std::uint32_t node = queue[head++];
            total += node;
            std::uint32_t edge = graph.offsets[node];
            const std::uint32_t edge_limit = graph.offsets[node + 1];
            for (; edge + kBatchLanes <= edge_limit; edge += static_cast<std::uint32_t>(kBatchLanes)) {
                // CSR gives a real packed load for destinations.
                const __m256i next_nodes =
                    _mm256_loadu_si256(reinterpret_cast<const __m256i*>(graph.edges.data() + edge));
                // But visited state is still scattered by node id.
                const __m256i seen =
                    _mm256_i32gather_epi32(visited.data(), next_nodes, static_cast<int>(sizeof(int)));
                _mm256_store_si256(reinterpret_cast<__m256i*>(dst_lanes.data()), next_nodes);
                const int unseen_mask = mask_bits(_mm256_cmpeq_epi32(seen, zero));
                // This scalar enqueue loop is why cache-aware scalar can still beat SIMD.
                for (std::size_t lane = 0; lane < kBatchLanes; ++lane) {
                    if (((unseen_mask >> lane) & 1) != 0 && visited[dst_lanes[lane]] == 0) {
                        visited[dst_lanes[lane]] = 1;
                        queue[tail++] = dst_lanes[lane];
                    }
                }
            }
            for (; edge < edge_limit; ++edge) {
                const std::uint32_t next = graph.edges[edge];
                if (visited[next] == 0) {
                    visited[next] = 1;
                    queue[tail++] = next;
                }
            }
        }
        total += tail;
    }
    return total;
}