Раздел 32 · Системное программирование: Zig, ассемблер, Verilog

Развёртывание циклов и параллелизм на уровне инструкций

lead~130 мин

открытый урокЭтот раздел читается без входа. Войди, чтобы отмечать прогресс, вести заметки и решать задачи в редакторе. войти

Развёртывание циклов и параллелизм на уровне инструкций

В прошлом уроке ты собрал модель процессора, который выбирает по несколько инструкций за такт, переименовывает регистры и исполняет всё, что готово, не дожидаясь очереди. Из таблицы задержек и графа потока данных вышли две границы CPE: граница задержки на цепочке зависимостей и граница пропускной способности на числе блоков. И там же обнаружилось неприятное: combine4 сидит ровно на границе задержки, а до границы пропускной способности ему ещё далеко: на дробном умножении больше чем в десять раз. Сегодня мы эту дистанцию пройдём. Сначала развернём цикл и увидим, что само по себе это почти ничего не даёт. Потом разорвём цепочку зависимостей на несколько независимых, и CPE упадёт вдвое, втрое, вчетверо. Разберём, почему компилятор не сделал этого сам, где предел, в который упрётся любая скалярная версия, и как тот же приём работает на скалярном произведении.

Цели урока

Понять, что развёртывание цикла k×1 снимает накладные расходы, но не трогает критический путь, и увидеть это в цифрах на combine5.
Написать combine6 с двумя аккумуляторами и обобщить его до k×k через inline for, зная, где живут аккумуляторы и почему массив из k штук раскладывается по регистрам.
Написать combine7 с переставленными скобками и объяснить по графу потока данных, почему одна перестановка даёт тот же выигрыш, что и второй аккумулятор.
Понимать, почему компилятор не переставит скобки у дробных чисел, и почему у целых он это делает по своему усмотрению и в любую сторону.
Считать три границы для любой развёртки k×a: задержка, делённая на число цепочек, пропускная способность блоков и загрузок, и вытеснение аккумуляторов в стек.
Читать ассемблер развёрнутого цикла и находить в нём цепочку, независимо от того, как компилятор перетасовал регистры.
Применить всё то же к скалярному произведению inner4 и объяснить, почему оно упирается в задержку сложения, а не в сумму задержек.

Идея: цикл ограничен не работой, а цепочкой

Начнём с того, где мы остановились. Все числа этого урока сняты на Apple M4 Max (aarch64, macOS), Zig 0.16.0, ReleaseFast, 2026-09-09, измерителем из урока про CPE. Рядом для сравнения таблица книги для Intel Haswell.

граница, тактов на операцию	i64 плюс	i64 умножить	f64 плюс	f64 умножить
задержка (M4 Max)	1.00	2.99	2.49	3.42
пропускная способность (M4)	0.16	0.33	0.25	0.25
задержка (Haswell, книга)	1.00	3.00	3.00	5.00
пропускная способность (Haswell)	0.50	1.00	1.00	0.50
`combine4` (M4 Max)	1.06	2.99	2.46	3.40

Смотри на последнюю строку и на первую. combine4 лежит на границе задержки с точностью до шума. Это значит, что процессор тратит на каждый элемент ровно столько тактов, сколько длится одна операция, и ни один из его блоков в это время не занят ничем другим: умножитель на M4 может начинать операцию каждые 0.25 такта, а начинает раз в 3.4 такта. Больше девяноста процентов его времени это ожидание.

Причина в графе потока данных из прошлого урока. Каждая итерация combine4 берёт аккумулятор, умножает его на элемент и кладёт обратно в аккумулятор. Следующая итерация не может начать своё умножение, пока не закончилось предыдущее, потому что ей нужен его результат. Получается одна

цепочка зависимостей

длиной в n умножений, и её длина в тактах равна n, умноженному на задержку.

Всё остальное в цикле, загрузки элементов, инкремент индекса, сравнение и переход, в цепочку не входит. Процессор исполняет это впереди, пока умножитель ждёт, и на времени это не сказывается. Поэтому есть ровно два способа сделать цикл быстрее.

Первый: сократить работу, которая всё-таки стоит тактов. У combine4 для целого сложения CPE равен 1.06 при задержке 1.00, и лишние сотые это накладные расходы цикла. Способ называется развёртывание, и сегодня мы увидим, что он даёт мало.

Второй: разорвать цепочку. Если вместо одной цепочки длиной n сделать две по n / 2, они пойдут параллельно на разных блоках, и цикл закончится вдвое раньше. Этот способ и есть

параллелизм на уровне инструкций

в чистом виде. Процессор умеет его использовать, но в combine4 использовать нечего: цепочка одна.

Файл combine.zig: что у нас уже есть

Мы продолжаем тот же файл, который тянется через все уроки о производительности. Чтобы урок читался без заглядывания назад, вот его шапка целиком: тип операции, нейтральный элемент, сама операция, вектор с интерфейсом и combine4 из урока про циклы и память. Версии с первой по третью я убрал, они нам сегодня не нужны.

//! Свёртка вектора одной операцией: файл, который мы тянем через все уроки о производительности.
//! Сегодня в нём появляются версии 5, 6 и 7 и их обобщения на любое k.

const std = @import("std");

pub const Op = enum { add, mul };

/// Нейтральный элемент операции: 0 для сложения, 1 для умножения.
pub fn ident(comptime T: type, comptime op: Op) T {
    return switch (op) {
        .add => 0,
        .mul => 1,
    };
}

/// Сама операция. Для целых с заворачиванием, для дробных обычная.
pub inline fn apply(comptime T: type, comptime op: Op, a: T, b: T) T {
    return switch (@typeInfo(T)) {
        .int => switch (op) {
            .add => a +% b,
            .mul => a *% b,
        },
        else => switch (op) {
            .add => a + b,
            .mul => a * b,
        },
    };
}

/// Вектор из книги: данные плюс интерфейс к ним, который компилятор
/// не имеет права встроить.
pub fn Vec(comptime T: type) type {
    return struct {
        const Self = @This();

        data: []T,

        pub fn init(data: []T) Self {
            return .{ .data = data };
        }

        pub noinline fn len(self: *const Self) usize {
            return self.data.len;
        }

        pub noinline fn get(self: *const Self, index: usize, dest: *T) bool {
            if (index >= self.data.len) return false;
            dest.* = self.data[index];
            return true;
        }
    };
}

/// Версия 4: аккумулятор в локальной переменной, то есть в регистре.
pub fn combine4(comptime T: type, comptime op: Op, v: *const Vec(T), dest: *T) void {
    const length = v.len();
    const data = v.data;
    var acc = ident(T, op);
    var i: usize = 0;
    while (i < length) : (i += 1) {
        acc = apply(T, op, acc, data[i]);
    }
    dest.* = acc;
}

Заворачивающие +% и *% для целых здесь не прихоть: книга считает в long, где переполнение молча отбрасывает старшие биты. Нам это ещё пригодится: заворачивание делает целое умножение ассоциативным при любых данных, а дробное умножение ассоциативным не бывает.

Развёртывание k×1: combine5

Развёртывание цикла

это первое, что приходит в голову, когда говорят «сократить накладные расходы». Одна итерация делает два элемента, проверок условия вдвое меньше, инкрементов вдвое меньше. Вот версия 5, развёртка 2×1: два элемента за итерацию, один аккумулятор.

/// Версия 5: развёртка 2×1. Два элемента за итерацию, один аккумулятор.
/// Меньше накладных расходов цикла, но цепочка зависимостей та же.
pub fn combine5(comptime T: type, comptime op: Op, v: *const Vec(T), dest: *T) void {
    const length = v.len();
    const data = v.data;
    var acc = ident(T, op);
    var i: usize = 0;
    while (i + 1 < length) : (i += 2) {
        acc = apply(T, op, apply(T, op, acc, data[i]), data[i + 1]);
    }
    while (i < length) : (i += 1) {
        acc = apply(T, op, acc, data[i]);
    }
    dest.* = acc;
}

Две детали, на которые стоит посмотреть до того, как мерить.

Условие основного цикла записано как i + 1 < length, а не как i < length - 1. В книге ровно на этом месте разобрана ошибка: length беззнаковое, и при пустом векторе length - 1 заворачивается в огромное число, цикл стартует и читает мимо массива. В Zig ситуация ещё жёстче: вычитание из usize нуля это переполнение, в Debug программа упадёт с паникой, а в ReleaseFast это неопределённое поведение. Форма i + 1 < length безопасна при любой длине. Общее правило для всех развёрток: сравнивай i + k <= length, никогда не вычитай из длины.

Второй цикл это хвост. Развёртка на два элемента обрабатывает пары, и если длина нечётная, последний элемент остаётся. Хвост досчитывает его по одному. При k равном двум это один элемент, при k равном десяти до девяти, и забыть хвост это классическая ошибка: тесты на длинах, кратных k, её не поймают. Поэтому тесты ниже гоняют длины 0, 1, k минус 1, k и 2k плюс 1.

Теперь граф потока данных одной итерации. Загрузки двух элементов независимы и идут впереди. А вот два умножения стоят одно за другим: второе берёт результат первого.

      acc ──▶ mul ──▶ mul ──▶ acc'
               ▲       ▲
           data[i]  data[i+1]

Цепочка из n умножений никуда не делась, только теперь на две операции цепочки приходится одна проверка условия. Критический путь на элемент прежний, задержка операции. Предсказание: CPE не изменится. Вот что намерил эталонный измеритель на M4 Max.

версия	i64 плюс	i64 умножить	f64 плюс	f64 умножить	приём
`combine4`	1.06	2.99	2.46	3.40	аккумулятор в регистре
`combine5`	1.12	3.00	2.15	3.31	развёртка 2×1

Целые остались на месте, в пределах шума. Дробное сложение чуть ниже границы задержки 2.49, дробное умножение на уровне шума. У книги на Haswell: 1.01, 3.01, 3.01, 5.01 против 1.27, 3.01, 3.01, 5.01 у combine4. Там развёртка помогла только целому сложению, и только потому, что у него задержка один такт, и накладные расходы цикла были заметны на его фоне. Там, где операция длиннее, накладные расходы и так прятались в тени цепочки.

Небольшое отступление про дробное сложение: 2.15 при границе 2.49 выглядит как нарушение нижней границы. Это не так. Граница 2.49 сама измерена как цепочка сложений, и на этом ядре часть сложений в цепочке идёт с меньшей задержкой, чем остальные: два значения из трёх в тестовых данных таковы, что сложение с ними обходится дешевле. Не увлекайся третьим знаком, шум серии на этой машине до четырнадцати процентов.

Задачка из книги: развёртка на пять элементов

Первое практическое упражнение главы просит переписать combine5 под k равное пяти. Попробуй сначала руками: пять вложенных apply, условие i + 4 < length, хвост по одному. Проверь себя на длинах 0, 4, 5 и 11. А потом посмотри, как это делается один раз для всех k.

/// Развёртка k×1 для любого k: обобщение combine5.
pub fn combineKx1(comptime T: type, comptime op: Op, comptime k: usize, v: *const Vec(T), dest: *T) void {
    const length = v.len();
    const data = v.data;
    var acc = ident(T, op);
    var i: usize = 0;
    while (i + k <= length) : (i += k) {
        inline for (0..k) |j| acc = apply(T, op, acc, data[i + j]);
    }
    while (i < length) : (i += 1) acc = apply(T, op, acc, data[i]);
    dest.* = acc;
}

k объявлен comptime, и inline for по диапазону от нуля до k разворачивается на этапе компиляции: в машинном коде будет ровно k обращений data[i + 0], data[i + 1] и так далее с константными смещениями, без всякого внутреннего цикла. Это то же самое, что ты написал бы руками для k равного пяти, только без риска ошибиться в смещении. Условие i + k <= length тот же безопасный трюк, что и раньше.

Что видно в ассемблере

Собери combine5 для f64-умножения как экспортируемую функцию и посмотри, что сделал компилятор. Ниже ассемблер для x86-64, снятый кросс-компиляцией с -target x86_64-linux -O ReleaseFast -femit-asm, основной цикл без директив.

// asm_demo.zig: тела трёх версий как отдельные символы, чтобы они не встроились.
const c = @import("combine.zig");

export fn mul5(data: [*]f64, len: usize) f64 {
    var v = c.Vec(f64).init(data[0..len]);
    var r: f64 = undefined;
    c.combine5(f64, .mul, &v, &r);
    return r;
}

; combine5, f64 умножение, основной цикл. x86-64, Zig 0.16.0, ReleaseFast.
.LBB2_13:
        mulsd   xmm0, qword ptr [rdi + 8*rax]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 8]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 16]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 24]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 32]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 40]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 48]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 56]
        add     rax, 8
        add     rdx, -4
        jne     .LBB2_13

Мы развернули цикл на два, а LLVM развернул его ещё на четыре: восемь элементов за итерацию, восемь mulsd подряд. Загрузка спрятана в операнд памяти каждой инструкции, поэтому отдельных movsd нет. Но смотри на регистры: все восемь умножений пишут в xmm0 и читают из xmm0. Это одна цепочка длиной в восемь, и на ней держится весь цикл. Развёртывание компилятор делает сам и лучше нас, накладные расходы у него и так были почти нулевые. Поэтому ручная развёртка k×1 на современном компиляторе почти ничего не даёт: она чинит то, что уже починено.

Несколько аккумуляторов: combine6 и k×k

Раз накладные расходы не проблема, остаётся цепочка. Единственный способ её укоротить это перестать складывать всё в одну переменную. Версия 6 держит два аккумулятора: чётные элементы копятся в одном, нечётные в другом, и только в самом конце два результата сводятся в один.

/// Версия 6: развёртка 2×2. Два аккумулятора, две независимые цепочки,
/// критический путь вдвое короче.
pub fn combine6(comptime T: type, comptime op: Op, v: *const Vec(T), dest: *T) void {
    const length = v.len();
    const data = v.data;
    var acc0 = ident(T, op);
    var acc1 = ident(T, op);
    var i: usize = 0;
    while (i + 1 < length) : (i += 2) {
        acc0 = apply(T, op, acc0, data[i]);
        acc1 = apply(T, op, acc1, data[i + 1]);
    }
    while (i < length) : (i += 1) {
        acc0 = apply(T, op, acc0, data[i]);
    }
    dest.* = apply(T, op, acc0, acc1);
}

Граф одной итерации теперь другой. Два умножения не связаны между собой: одно берёт acc0, другое acc1, и ни одно не ждёт другого.

      acc0 ──▶ mul ──▶ acc0'
                ▲
            data[i]

      acc1 ──▶ mul ──▶ acc1'
                ▲
           data[i+1]

Раскладываем n / 2 итераций подряд: получаются две цепочки по n / 2 умножений каждая. Внеочередное ядро видит, что они независимы, и запускает их вперемешку на свободных блоках. Критический путь на элемент теперь половина задержки. На M4 это предсказывает 0.50, 1.50, 1.25 и 1.71, и вот измерение.

версия	i64 плюс	i64 умножить	f64 плюс	f64 умножить	приём
`combine5`	1.12	3.00	2.15	3.31	развёртка 2×1
`combine6`	0.73	1.50	1.26	1.75	два аккумулятора

Три колонки из четырёх легли на предсказание с точностью до шума, целое сложение отстало: 0.73 вместо 0.50. Здесь одна цепочка стоит один такт на элемент, и её половина это полтакта, а на таком масштабе уже заметны сами загрузки и обвязка цикла, которые в остальных колонках прятались в тени. Книга на Haswell получает 0.81, 1.51, 1.51 и 2.51: те же половины задержек 3 и 5 у дробных, и то же отставание на целом сложении.

Как это выглядит в машинном коде

Тот же трюк с экспортом для combine6. Здесь компилятор сделал кое-что, чего ты, возможно, не ждёшь.

; combine6, f64 умножение, основной цикл. x86-64, Zig 0.16.0, ReleaseFast.
.LBB1_14:
        movupd  xmm1, xmmword ptr [rdi + 8*rax]
        mulpd   xmm1, xmm0
        movupd  xmm0, xmmword ptr [rdi + 8*rax + 16]
        mulpd   xmm0, xmm1
        movupd  xmm1, xmmword ptr [rdi + 8*rax + 32]
        mulpd   xmm1, xmm0
        movupd  xmm0, xmmword ptr [rdi + 8*rax + 48]
        mulpd   xmm0, xmm1
        add     rax, 8
        add     rdx, -4
        jne     .LBB1_14
...
        unpckhpd xmm0, xmm0
        mulsd   xmm0, xmm1

Вместо двух скалярных mulsd на два регистра здесь одна mulpd, упакованное умножение двух f64 в одном 128-битном регистре. Компилятор заметил, что acc0 и acc1 обновляются одинаково по соседним элементам, и сложил оба аккумулятора в две половины xmm0: нижняя это acc0, верхняя acc1. Загрузка movupd тоже берёт два элемента разом. После цикла unpckhpd достаёт верхнюю половину, и последний mulsd сводит аккумуляторы в один результат, это наш apply(acc0, acc1). Цепочка по-прежнему одна, xmm0 в xmm1 и обратно, но каждое её звено теперь обрабатывает два элемента, и критический путь на элемент ровно половина задержки. Две независимые цепочки и одна цепочка «в две полосы» по времени одно и то же, и процессор с компилятором выбрали вторую форму. Векторные регистры как самостоятельный приём мы разберём в следующем уроке.

Заметь, чего компилятор не сделал: combine4 он так не упаковал. Там был один аккумулятор и одна цепочка, и упаковать в две полосы её нельзя без перестановки скобок. Почему нельзя, чуть ниже.

Обобщение на k аккумуляторов

Два аккумулятора дают половину задержки, три дадут треть, k дадут k-ю часть, пока хватает блоков и регистров. Вот развёртка k×k: k элементов за итерацию, k аккумуляторов.

/// Развёртка k×k: k аккумуляторов, обобщение combine6. При k = 2
/// совпадает с ним.
pub fn combineKxK(comptime T: type, comptime op: Op, comptime k: usize, v: *const Vec(T), dest: *T) void {
    const length = v.len();
    const data = v.data;
    var acc: [k]T = @splat(ident(T, op));
    var i: usize = 0;
    while (i + k <= length) : (i += k) {
        inline for (0..k) |j| acc[j] = apply(T, op, acc[j], data[i + j]);
    }
    while (i < length) : (i += 1) acc[0] = apply(T, op, acc[0], data[i]);
    var result = acc[0];
    inline for (1..k) |j| result = apply(T, op, result, acc[j]);
    dest.* = result;
}

Аккумуляторы лежат в массиве [k]T, заполненном нейтральным элементом через @splat. Массив в локальной переменной с константными индексами компилятор разложит по регистрам так же, как разложил бы acc0, acc1, acc2, объявленные по отдельности: индексы acc[j] внутри inline for известны на этапе компиляции, обращения в память по ним не нужны. Хвост уходит в нулевой аккумулятор, а сведение это k минус 1 операций, цепочка, но короткая и вне цикла.

Таблица книги для k×k на Haswell при k от одного до десяти сходится к 0.55, 1.00, 1.01 и 0.52, и это уже не задержки, а граница пропускной способности. К ней мы вернёмся через раздел, а пока поиграй с моделью.

Виджет считает CPE как наибольшее из трёх чисел: задержка, делённая на число аккумуляторов, граница пропускной способности и штраф за аккумуляторы, не поместившиеся в регистры. Пресет M4 берёт границы, снятые нашим измерителем, и рядом показывает измеренный CPE для тех форм, которые есть в эталонном замере. Кривая по числу аккумуляторов падает как гипербола и упирается в горизонталь пропускной способности: после этого добавлять аккумуляторы бесполезно. А если уменьшить число свободных регистров, кривая ломается в другую сторону: лишние аккумуляторы уезжают в стек, и цепочка удлиняется на запись и чтение. Этот предел подробно разберёт следующий урок.

Почему компилятор не сделал этого сам

Вопрос, который здесь обязан возникнуть: если два аккумулятора это две строки кода, почему combine4 не превращается в combine6 автоматически? Компилятор развернул цикл на восемь без спроса, а тут остановился.

Потому что это разные преобразования. Развёртка меняет порядок проверок, но не порядок операций над данными: результат совпадает бит в бит. Два аккумулятора меняют порядок: combine4 считает ((((x0 · x1) · x2) · x3) · x4), а combine6 считает ((x0 · x2) · x4) · (x1 · x3). Это одно и то же число, только если операция

ассоциативна

и коммутативна. Для целых с заворачиванием это так при любых данных: арифметика по модулю 2 в шестьдесят четвёртой степени это кольцо, и там скобки можно ставить как угодно. Для чисел с плавающей точкой это неправда, и вот доказательство в виде теста.

const std = @import("std");

fn sumOneChain(data: []const f64) f64 {
    var acc: f64 = 0;
    for (data) |x| acc += x;
    return acc;
}

fn sumTwoChains(data: []const f64) f64 {
    var acc0: f64 = 0;
    var acc1: f64 = 0;
    var i: usize = 0;
    while (i + 1 < data.len) : (i += 2) {
        acc0 += data[i];
        acc1 += data[i + 1];
    }
    while (i < data.len) : (i += 1) acc0 += data[i];
    return acc0 + acc1;
}

test "две цепочки дают другой ответ, чем одна" {
    // Единица рядом с 1e16 теряется: в f64 после 1e16 следующее
    // представимое число это 1e16 + 2. В одной цепочке единицы пропадают
    // по очереди, доживает только последняя. В двух цепочках все четыре
    // единицы копятся отдельно и доживают до ответа.
    const data = [_]f64{ 1e16, 1, -1e16, 1, 1e16, 1, -1e16, 1 };
    try std.testing.expectEqual(@as(f64, 1), sumOneChain(&data));
    try std.testing.expectEqual(@as(f64, 4), sumTwoChains(&data));
}

1/1 test.две цепочки дают другой ответ, чем одна...OK
All 1 tests passed.

Одна цепочка отвечает 1, две цепочки отвечают 4. Это не ошибка ни в одной из них, это разные вычисления. Компилятор обязан сохранить результат программы бит в бит, поэтому для дробных он не имеет права ни на два аккумулятора, ни на перестановку скобок, и всё это остаётся тебе. Ты знаешь свои данные и знаешь, что для суммы зарплат разница в последнем бите неважна, а для итерационного метода, который сходится на границе точности, может быть важна. Компилятор этого не знает.

У этого правила есть выключатель: @setFloatMode(.optimized) разрешает компилятору считать дробную арифметику ассоциативной. Я проверил его на x86-64 в Zig 0.16.0: сумма f64 в цикле с этим режимом собралась в тот же самый код, что и без него, восемь addsd в одну цепочку через xmm0. Разрешение переставлять скобки не обязывает компилятор ими пользоваться, и в этой версии он не воспользовался. Домашнее задание просит повторить опыт на твоей платформе.

А теперь про целые, у которых разрешение есть всегда. Вернись к таблице: combine4 для целого сложения показал 1.06, то есть одну цепочку. Компилятор имел право разложить её на четыре и не стал. Зато для combineKx1 с k равным пяти на том же целом сложении мы ниже увидим 0.25: там он право использовал, потому что развёрнутое тело подсказало ему форму. Вывод неприятный, но честный: у целых компилятор переставляет скобки по своему усмотрению, иногда в твою пользу, иногда нет, и предсказать это по исходному коду нельзя, только по ассемблеру или по замеру.

Переупорядочение операций: combine7 и k×1a

Есть третий способ разорвать цепочку, и он самый неожиданный: не добавлять аккумуляторов, а переставить скобки в одном выражении. Версия 7 отличается от версии 5 расположением одной скобки.

/// Версия 7: развёртка 2×1a. Один аккумулятор, но скобки переставлены:
/// сначала перемножаются соседние элементы, и только произведение
/// попадает в цепочку аккумулятора.
pub fn combine7(comptime T: type, comptime op: Op, v: *const Vec(T), dest: *T) void {
    const length = v.len();
    const data = v.data;
    var acc = ident(T, op);
    var i: usize = 0;
    while (i + 1 < length) : (i += 2) {
        acc = apply(T, op, acc, apply(T, op, data[i], data[i + 1]));
    }
    while (i < length) : (i += 1) {
        acc = apply(T, op, acc, data[i]);
    }
    dest.* = acc;
}

В combine5 было (acc · x) · y, в combine7 стало acc · (x · y). Граф одной итерации:

   data[i] ──▶ mul ──▶ mul ──▶ acc'
   data[i+1] ─▶  ▲       ▲
                 │      acc

Произведение x · y не зависит от аккумулятора: оба его операнда только что загружены из памяти. Процессор может посчитать его заранее, ещё до того, как предыдущая итерация отдаст acc. В цепочку аккумулятора входит только второе умножение. Итого на два элемента одна операция в цепочке, критический путь на элемент половина задержки, как у двух аккумуляторов, а аккумулятор один.

версия	i64 плюс	i64 умножить	f64 плюс	f64 умножить	приём
`combine6`	0.73	1.50	1.26	1.75	два аккумулятора
`combine7`	1.12	3.00	1.28	1.76	переупорядочение 2×1a

Дробные колонки сошлись с combine6 до сотых. А целые не сдвинулись вовсе: 1.12 и 3.00, ровно как у combine5. Для целого сложения это ожидаемо, его задержка один такт, и половина такта на элемент упирается в загрузки и обвязку так же, как в combine6. Но целое умножение с задержкой три должно было дать полтора, а дало три. Книга на Haswell получает для 2×1a ровно 1.51 на целом умножении, у нас нет.

Ассемблер отвечает на вопрос без гаданий. Вот основной цикл combine7 для целого умножения на x86-64.

; combine7, i64 умножение, основной цикл. x86-64, Zig 0.16.0, ReleaseFast.
        imul    rax, qword ptr [rdi + 8*rcx]
        imul    rax, qword ptr [rdi + 8*rcx + 8]
        imul    rax, qword ptr [rdi + 8*rcx + 16]
        imul    rax, qword ptr [rdi + 8*rcx + 24]
        imul    rax, qword ptr [rdi + 8*rcx + 32]
        imul    rax, qword ptr [rdi + 8*rcx + 40]
        imul    rax, qword ptr [rdi + 8*rcx + 48]
        imul    rax, qword ptr [rdi + 8*rcx + 56]

Восемь imul в одну цепочку через rax. Компилятор взял наше acc · (x · y), вспомнил, что целое умножение ассоциативно, и переставил скобки обратно: (acc · x) · y, потому что так одна инструкция меньше и операнд берётся прямо из памяти. Он не знает про задержки, он оптимизирует число инструкций. У дробных он трогать скобки не имел права, и там наша перестановка дожила до машинного кода. Вот тот же цикл для f64:

; combine7, f64 умножение, основной цикл. x86-64, Zig 0.16.0, ReleaseFast.
.LBB0_13:
        movsd   xmm1, qword ptr [rdi + 8*rax]
        movsd   xmm2, qword ptr [rdi + 8*rax + 16]
        mulsd   xmm1, qword ptr [rdi + 8*rax + 8]
        mulsd   xmm2, qword ptr [rdi + 8*rax + 24]
        mulsd   xmm1, xmm0
        movsd   xmm3, qword ptr [rdi + 8*rax + 32]
        mulsd   xmm3, qword ptr [rdi + 8*rax + 40]
        mulsd   xmm2, xmm1
        movsd   xmm0, qword ptr [rdi + 8*rax + 48]
        mulsd   xmm0, qword ptr [rdi + 8*rax + 56]
        mulsd   xmm3, xmm2
        mulsd   xmm0, xmm3
        add     rax, 8
        add     rdx, -4
        jne     .LBB0_13

Восемь элементов, восемь умножений, но в цепочке только четыре: xmm0 в xmm1, xmm1 в xmm2, xmm2 в xmm3, xmm3 в xmm0. Остальные четыре перемножают пары из памяти и от аккумулятора не зависят. Ровно то, что нарисовано в графе, только компилятор ещё раз развернул на четыре.

Обобщение на любое k: k элементов сворачиваются между собой, и только результат присоединяется к аккумулятору.

/// Развёртка k×1a: обобщение combine7, k элементов свёрнуты между собой
/// и только потом присоединены к аккумулятору.
pub fn combineKx1a(comptime T: type, comptime op: Op, comptime k: usize, v: *const Vec(T), dest: *T) void {
    const length = v.len();
    const data = v.data;
    var acc = ident(T, op);
    var i: usize = 0;
    while (i + k <= length) : (i += k) {
        var group = data[i];
        inline for (1..k) |j| group = apply(T, op, group, data[i + j]);
        acc = apply(T, op, acc, group);
    }
    while (i < length) : (i += 1) acc = apply(T, op, acc, data[i]);
    dest.* = acc;
}

Внутри группы это цепочка длиной k минус 1, но она не переходит из итерации в итерацию, поэтому процессор считает группы следующих итераций впереди. В цепочку аккумулятора входит одна операция на k элементов.

Задачка из книги: пять способов расставить скобки

Второе практическое упражнение главы: развёртка 3×1 для умножения, три элемента x, y, z за итерацию, и пять способов расставить скобки в r · x · y · z. Для каждого нужно предсказать CPE на дробном умножении. Возьми задержку M4, 3.42 такта, и считай, сколько умножений на три элемента лежат в цепочке аккумулятора r.

//! Пять способов расставить скобки в развёртке 3×1 умножения.
//! Все считают одно и то же, но критический путь у них разной длины.

const std = @import("std");

fn Tail(comptime T: type) type {
    return struct {
        /// Хвост общий: досчитать элементы, не вошедшие в тройки.
        fn finish(acc: T, data: []const T, from: usize) T {
            var r = acc;
            var i = from;
            while (i < data.len) : (i += 1) r = r * data[i];
            return r;
        }
    };
}

/// A1: ((r * x) * y) * z. Три умножения подряд в цепочке аккумулятора.
pub fn mul3A1(data: []const f64) f64 {
    var r: f64 = 1;
    var i: usize = 0;
    while (i + 2 < data.len) : (i += 3) {
        r = ((r * data[i]) * data[i + 1]) * data[i + 2];
    }
    return Tail(f64).finish(r, data, i);
}

/// A2: (r * (x * y)) * z. Произведение x*y независимо, в цепочке два умножения.
pub fn mul3A2(data: []const f64) f64 {
    var r: f64 = 1;
    var i: usize = 0;
    while (i + 2 < data.len) : (i += 3) {
        r = (r * (data[i] * data[i + 1])) * data[i + 2];
    }
    return Tail(f64).finish(r, data, i);
}

/// A3: r * ((x * y) * z). Тройка перемножена отдельно, в цепочке одно умножение.
pub fn mul3A3(data: []const f64) f64 {
    var r: f64 = 1;
    var i: usize = 0;
    while (i + 2 < data.len) : (i += 3) {
        r = r * ((data[i] * data[i + 1]) * data[i + 2]);
    }
    return Tail(f64).finish(r, data, i);
}

/// A4: r * (x * (y * z)). То же, что A3, только тройка свёрнута с другого конца.
pub fn mul3A4(data: []const f64) f64 {
    var r: f64 = 1;
    var i: usize = 0;
    while (i + 2 < data.len) : (i += 3) {
        r = r * (data[i] * (data[i + 1] * data[i + 2]));
    }
    return Tail(f64).finish(r, data, i);
}

/// A5: (r * x) * (y * z). Два умножения в цепочке: r*x, потом умножить на y*z.
pub fn mul3A5(data: []const f64) f64 {
    var r: f64 = 1;
    var i: usize = 0;
    while (i + 2 < data.len) : (i += 3) {
        r = (r * data[i]) * (data[i + 1] * data[i + 2]);
    }
    return Tail(f64).finish(r, data, i);
}

test "все пять расстановок дают одно произведение" {
    const data = [_]f64{ 1.5, 2, 0.5, 4, 1.25, 2, 0.8, 3, 1, 2, 0.25 };
    const expected = mul3A1(&data);
    try std.testing.expectApproxEqRel(expected, mul3A2(&data), 1e-12);
    try std.testing.expectApproxEqRel(expected, mul3A3(&data), 1e-12);
    try std.testing.expectApproxEqRel(expected, mul3A4(&data), 1e-12);
    try std.testing.expectApproxEqRel(expected, mul3A5(&data), 1e-12);
    try std.testing.expectApproxEqRel(@as(f64, 18), expected, 1e-12);
}

Разбор. Считаем умножения, через которые проходит r от начала итерации до конца.

A1, ((r · x) · y) · z: три подряд, каждое ждёт предыдущего. Три умножения на три элемента, CPE равен задержке, 3.42.
A2, (r · (x · y)) · z: произведение x · y считается заранее, r проходит через два умножения. Две трети задержки, 2.28.
A3, r · ((x · y) · z): вся тройка перемножена отдельно, r проходит через одно умножение. Треть задержки, 1.14.
A4, r · (x · (y · z)): то же самое, одно умножение в цепочке, 1.14.
A5, (r · x) · (y · z): r · x первое звено, умножение на y · z второе. Две трети, 2.28.

На Haswell с задержкой пять это 5.00, 3.33, 1.67, 1.67 и 3.33. Проверим на M4 Max, замер ниже: A1 дал 3.19, A2 2.23, A3 1.20, A4 1.23, A5 2.25. Порядок и пропорции ровно те, что предсказаны, а сами числа чуть ниже задержки по той же причине, что и у combine5: тестовые данные дружелюбны к умножителю.

Ограничение пропускной способности

Аккумуляторов можно добавлять много, но кривая перестаёт падать. Причина в том, что кроме цепочки у цикла есть ещё и просто работа: на каждый элемент одна операция и одна загрузка, и у каждого блока есть предел, сколько операций он начинает за такт.

Граница пропускной способности

это наибольшее из двух чисел: операции и загрузки. У Haswell из книги дробное умножение умеют два блока, дробное сложение один, целое умножение один, целое сложение четыре, а портов загрузки два. Отсюда книжные 0.50, 1.00, 1.00 и 0.50: целое сложение упирается не в свои четыре сумматора, а в два порта загрузки, а дробное умножение в два своих блока. На M4 портов больше, и наш измеритель дал 0.16, 0.33, 0.25 и 0.25 для самих операций, при трёх портах загрузки на ядре производительности это 0.33 на элемент для любой скалярной свёртки.

Книга приводит таблицу k×k при k от одного до десяти для Haswell. Целое сложение доходит до 0.55 уже при k равном трём и дальше не падает; целое умножение и дробное сложение упираются в 1.00 при k около трёх; дробное умножение доходит до 0.52 при k равном десяти. Вот наш замер на M4 Max (aarch64, macOS, Zig 0.16.0, ReleaseFast, 2026-09-13), снятый маленьким измерителем, который приведён в конце раздела.

форма	i64 плюс	i64 умножить	f64 плюс	f64 умножить
1×1	0.99	2.97	2.45	3.29
2×2	0.62	1.49	1.15	1.68
3×3	0.47	0.99	0.82	1.13
4×4	0.49	0.76	0.66	0.84
6×6	0.35	0.60	0.46	0.61
8×8	0.25	0.40	0.34	0.46
10×10	0.20	0.33	0.29	0.37

Целое умножение на десяти аккумуляторах дошло до 0.33, это его граница по блокам. Остальные три колонки ещё падают: задержка, делённая на десять, у них 0.10, 0.25 и 0.34, а граница по загрузкам 0.33, и десяти цепочек мало, чтобы её увидеть. Дальше начинаются регистры: на aarch64 их тридцать один целочисленный и тридцать два векторных, поэтому M4 терпит больше аккумуляторов, чем Haswell с его шестнадцатью, но и у него предел есть. Что происходит, когда аккумуляторы перестают помещаться, покажет следующий урок.

Для сравнения k×1 и k×1a на тех же формах.

форма	i64 плюс	i64 умножить	f64 плюс	f64 умножить
2×1	0.99	2.97	2.13	3.10
3×1	0.66	2.96	2.20	3.16
5×1	0.25	2.97	2.15	3.12
2×1a	1.01	2.97	1.25	1.72
3×1a	0.69	2.97	0.86	1.20
4×1a	0.35	2.95	0.76	1.00
6×1a	0.22	2.97	0.45	0.60

Три вещи стоит увидеть. Дробные колонки k×1 стоят на задержке при любом k: развёртка без второй цепочки бесполезна. Дробные колонки k×1a падают как задержка, делённая на k: 3.42 на шесть это 0.57, замер 0.60. И целое умножение в k×1a не сдвигается ни на сотую при любом k: компилятор каждый раз переставляет скобки обратно, как мы видели в ассемблере. А вот целое сложение в k×1 при k равном пяти внезапно даёт 0.25, вчетверо лучше единицы: здесь компилятор, наоборот, воспользовался ассоциативностью и разложил развёрнутое тело на несколько цепочек. Тот же компилятор, та же операция, другой исходный код, другое решение.

Маленький измеритель

Числа выше сняты этой программой. Она короче измерителя из урока про CPE: частота калибруется той же цепочкой сложений, для каждой формы берётся минимум из семи прогонов на двух размерах, 2048 и 4096 элементов, а CPE это наклон между двумя точками. Оба размера сидят в L1, накладные расходы вызова одинаковы и из разности выпадают. Файл лежит рядом с combine.zig, inner.zig из следующего раздела и paren.zig с пятью расстановками скобок.

//! unroll_bench.zig. Собирать только так: zig build-exe unroll_bench.zig -O ReleaseFast

const std = @import("std");
const c = @import("combine.zig");
const inner = @import("inner.zig");
const paren = @import("paren.zig");

const Vec = c.Vec;
const Op = c.Op;

/// Пропускает значение через пустой asm, чтобы компилятор считал его неизвестным.
inline fn blackBox(value: anytype) @TypeOf(value) {
    var v = value;
    asm volatile (""
        : [v] "+r" (v),
    );
    return v;
}

/// Цепочка зависимых сложений: ровно по такту на сложение.
fn addChain(n: usize) u64 {
    var acc: u64 = 1;
    const step: u64 = blackBox(@as(u64, 7));
    var i: usize = 0;
    while (i < n) : (i += 1) acc = blackBox(acc +% step);
    return acc;
}

fn elapsedNs(io: std.Io, started: std.Io.Timestamp) f64 {
    return @floatFromInt(started.durationTo(std.Io.Timestamp.now(io, .awake)).nanoseconds);
}

/// Частота ядра в тактах на наносекунду: лучший из пяти замеров цепочки.
fn calibrate(io: std.Io) f64 {
    const adds: usize = 200_000_000;
    var best: f64 = std.math.inf(f64);
    var attempt: usize = 0;
    while (attempt < 5) : (attempt += 1) {
        const started = std.Io.Timestamp.now(io, .awake);
        std.mem.doNotOptimizeAway(addChain(adds));
        best = @min(best, elapsedNs(io, started));
    }
    return @as(f64, @floatFromInt(adds)) / best;
}

/// Минимальное время одного вызова над n элементами, наносекунды.
fn minNs(io: std.Io, ctx: anytype, n: usize) f64 {
    const reps = 1_000_000 / n;
    var warm: usize = 0;
    while (warm < 2) : (warm += 1) ctx.run(n, reps);
    var best: f64 = std.math.inf(f64);
    var run: usize = 0;
    while (run < 7) : (run += 1) {
        const started = std.Io.Timestamp.now(io, .awake);
        ctx.run(n, reps);
        best = @min(best, elapsedNs(io, started) / @as(f64, @floatFromInt(reps)));
    }
    return best;
}

/// CPE как наклон между n = 2048 и n = 4096: накладные расходы вызова
/// одинаковы в обеих точках и из разности выпадают.
fn cpe(io: std.Io, ghz: f64, ctx: anytype) f64 {
    const small = minNs(io, ctx, 2048);
    const large = minNs(io, ctx, 4096);
    return (large - small) * ghz / 2048.0;
}

/// Какую развёртку гонять: форма и k.
const Shape = enum { plain, kx1, kxk, kx1a };

fn Runner(comptime T: type, comptime op: Op, comptime shape: Shape, comptime k: usize) type {
    return struct {
        data: []T,

        pub fn run(self: @This(), n: usize, reps: usize) void {
            var rep: usize = 0;
            while (rep < reps) : (rep += 1) {
                var v = Vec(T).init(self.data[0..n]);
                var result: T = undefined;
                switch (shape) {
                    .plain => c.combine4(T, op, &v, &result),
                    .kx1 => c.combineKx1(T, op, k, &v, &result),
                    .kxk => c.combineKxK(T, op, k, &v, &result),
                    .kx1a => c.combineKx1a(T, op, k, &v, &result),
                }
                std.mem.doNotOptimizeAway(result);
            }
        }
    };
}

fn InnerRunner(comptime T: type, comptime shape: Shape, comptime k: usize) type {
    return struct {
        u: []T,
        v: []T,

        pub fn run(self: @This(), n: usize, reps: usize) void {
            var rep: usize = 0;
            while (rep < reps) : (rep += 1) {
                var result: T = undefined;
                switch (shape) {
                    .plain => inner.inner4(T, self.u[0..n], self.v[0..n], &result),
                    .kx1 => inner.innerKx1(T, k, self.u[0..n], self.v[0..n], &result),
                    .kxk => inner.innerKxK(T, k, self.u[0..n], self.v[0..n], &result),
                    .kx1a => inner.innerKx1a(T, k, self.u[0..n], self.v[0..n], &result),
                }
                std.mem.doNotOptimizeAway(result);
            }
        }
    };
}

fn ParenRunner(comptime f: fn ([]const f64) f64) type {
    return struct {
        data: []f64,

        pub fn run(self: @This(), n: usize, reps: usize) void {
            var rep: usize = 0;
            while (rep < reps) : (rep += 1) std.mem.doNotOptimizeAway(f(self.data[0..n]));
        }
    };
}

fn fill(comptime T: type, data: []T) void {
    for (data, 0..) |*cell, i| {
        cell.* = switch (@typeInfo(T)) {
            .float => 1.0 + @as(T, @floatFromInt(i % 3)) / 1024.0,
            else => @intCast(1 + i % 3),
        };
    }
}

const Column = struct { key: []const u8, T: type, op: Op };
const columns = [_]Column{
    .{ .key = "i64+", .T = i64, .op = .add },
    .{ .key = "i64*", .T = i64, .op = .mul },
    .{ .key = "f64+", .T = f64, .op = .add },
    .{ .key = "f64*", .T = f64, .op = .mul },
};

const Row = struct { name: []const u8, shape: Shape, k: usize };
const rows = [_]Row{
    .{ .name = "combine4  1x1", .shape = .plain, .k = 1 },
    .{ .name = "combine5  2x1", .shape = .kx1, .k = 2 },
    .{ .name = "combine6  2x2", .shape = .kxk, .k = 2 },
    .{ .name = "combine7  2x1a", .shape = .kx1a, .k = 2 },
    .{ .name = "k x 1     3x1", .shape = .kx1, .k = 3 },
    .{ .name = "k x 1     5x1", .shape = .kx1, .k = 5 },
    .{ .name = "k x k     3x3", .shape = .kxk, .k = 3 },
    .{ .name = "k x k     4x4", .shape = .kxk, .k = 4 },
    .{ .name = "k x k     6x6", .shape = .kxk, .k = 6 },
    .{ .name = "k x k     8x8", .shape = .kxk, .k = 8 },
    .{ .name = "k x k   10x10", .shape = .kxk, .k = 10 },
    .{ .name = "k x 1a    3x1a", .shape = .kx1a, .k = 3 },
    .{ .name = "k x 1a    4x1a", .shape = .kx1a, .k = 4 },
    .{ .name = "k x 1a    6x1a", .shape = .kx1a, .k = 6 },
};

pub fn main(init: std.process.Init) !void {
    const gpa = init.arena.allocator();
    var buf: [4096]u8 = undefined;
    var writer = std.Io.File.stdout().writer(init.io, &buf);
    const out = &writer.interface;

    const ghz = calibrate(init.io);
    try out.print("частота ядра по калибровке: {d:.2} ГГц\n\n", .{ghz});
    try out.writeAll("CPE свёртки, тактов на элемент\n");
    try out.writeAll("версия              i64+    i64*    f64+    f64*\n");
    inline for (rows) |row| {
        try out.print("{s:<16}", .{row.name});
        inline for (columns) |column| {
            const data = try gpa.alloc(column.T, 4096);
            fill(column.T, data);
            const runner = Runner(column.T, column.op, row.shape, row.k){ .data = data };
            try out.print(" {d:>7.2}", .{cpe(init.io, ghz, runner)});
        }
        try out.writeAll("\n");
    }

    try out.writeAll("\nпять расстановок скобок, f64 умножение, 3x1\n");
    const fns = .{ paren.mul3A1, paren.mul3A2, paren.mul3A3, paren.mul3A4, paren.mul3A5 };
    const names = [_][]const u8{ "A1 ((r*x)*y)*z", "A2 (r*(x*y))*z", "A3 r*((x*y)*z)", "A4 r*(x*(y*z))", "A5 (r*x)*(y*z)" };
    inline for (fns, names) |f, name| {
        const data = try gpa.alloc(f64, 4096);
        fill(f64, data);
        try out.print("{s:<16} {d:>7.2}\n", .{ name, cpe(init.io, ghz, ParenRunner(f){ .data = data }) });
    }

    try out.writeAll("\nскалярное произведение\n");
    try out.writeAll("версия              i64     f64\n");
    const inner_rows = [_]Row{
        .{ .name = "inner4    1x1", .shape = .plain, .k = 1 },
        .{ .name = "inner     6x1", .shape = .kx1, .k = 6 },
        .{ .name = "inner     6x6", .shape = .kxk, .k = 6 },
        .{ .name = "inner     6x1a", .shape = .kx1a, .k = 6 },
    };
    inline for (inner_rows) |row| {
        try out.print("{s:<16}", .{row.name});
        inline for (.{ i64, f64 }) |T| {
            const u = try gpa.alloc(T, 4096);
            const v = try gpa.alloc(T, 4096);
            fill(T, u);
            fill(T, v);
            const runner = InnerRunner(T, row.shape, row.k){ .u = u, .v = v };
            try out.print(" {d:>7.2}", .{cpe(init.io, ghz, runner)});
        }
        try out.writeAll("\n");
    }
    try out.flush();
}

Вывод целиком, тот самый прогон, из которого взяты таблицы выше.

частота ядра по калибровке: 4.46 ГГц

CPE свёртки, тактов на элемент
версия              i64+    i64*    f64+    f64*
combine4  1x1       0.99    2.97    2.45    3.29
combine5  2x1       0.99    2.97    2.13    3.10
combine6  2x2       0.62    1.49    1.15    1.68
combine7  2x1a      1.01    2.97    1.25    1.72
k x 1     3x1       0.66    2.96    2.20    3.16
k x 1     5x1       0.25    2.97    2.15    3.12
k x k     3x3       0.47    0.99    0.82    1.13
k x k     4x4       0.49    0.76    0.66    0.84
k x k     6x6       0.35    0.60    0.46    0.61
k x k     8x8       0.25    0.40    0.34    0.46
k x k   10x10       0.20    0.33    0.29    0.37
k x 1a    3x1a      0.69    2.97    0.86    1.20
k x 1a    4x1a      0.35    2.95    0.76    1.00
k x 1a    6x1a      0.22    2.97    0.45    0.60

пять расстановок скобок, f64 умножение, 3x1
A1 ((r*x)*y)*z      3.19
A2 (r*(x*y))*z      2.23
A3 r*((x*y)*z)      1.20
A4 r*(x*(y*z))      1.23
A5 (r*x)*(y*z)      2.25

скалярное произведение
версия              i64     f64
inner4    1x1       1.28    2.81
inner     6x1       1.07    2.33
inner     6x6       0.61    1.16
inner     6x1a      1.08    0.56

Одно предупреждение из практики. Первый раз я снял эти числа, пока на машине шла сборка в десять параллельных процессов, и калибровка выдала 0.93 ГГц вместо 4.46: часы считают реальное время, а не время нашего процесса, и всё, что процессор отдал соседям, попадает в замер. Все CPE вышли вдвое-втрое ниже правды и при этом не в одной пропорции. Правило прежнее, из урока про CPE: мерить на тихой машине, смотреть на калибровку, минимум из серии, и не верить числу, которое не воспроизвелось во втором прогоне.

Скалярное произведение: inner4 и его развёртки

Домашние задачи главы крутятся вокруг одного цикла: скалярное произведение двух векторов, sum += u[i] * v[i]. Он интереснее свёртки, потому что в итерации две операции, умножение и сложение, и две загрузки. Вот файл целиком, с развёртками и тестами.

//! inner.zig. Скалярное произведение: домашние задачи книги про границы CPE.
//!
//! `sum += u[i] * v[i]`. На каждом шаге умножение и сложение, но в цепочку
//! аккумулятора входит только сложение: произведения независимы и
//! считаются впереди. Поэтому нижняя граница по задержке это задержка
//! сложения, а не сумма задержек. Развёртка k×k снимает и её, и упирается
//! в пропускную способность загрузок: два чтения на элемент.

const std = @import("std");

inline fn mulAdd(comptime T: type, acc: T, x: T, y: T) T {
    return if (@typeInfo(T) == .int) acc +% (x *% y) else acc + x * y;
}

inline fn add(comptime T: type, a: T, b: T) T {
    return if (@typeInfo(T) == .int) a +% b else a + b;
}

/// Версия из книги: один аккумулятор, один элемент за итерацию.
pub fn inner4(comptime T: type, u: []const T, v: []const T, dest: *T) void {
    var sum: T = 0;
    var i: usize = 0;
    while (i < u.len) : (i += 1) {
        sum = mulAdd(T, sum, u[i], v[i]);
    }
    dest.* = sum;
}

/// Развёртка k×1: один аккумулятор, k элементов за итерацию.
pub fn innerKx1(comptime T: type, comptime k: usize, u: []const T, v: []const T, dest: *T) void {
    var sum: T = 0;
    var i: usize = 0;
    while (i + k <= u.len) : (i += k) {
        inline for (0..k) |j| sum = mulAdd(T, sum, u[i + j], v[i + j]);
    }
    while (i < u.len) : (i += 1) sum = mulAdd(T, sum, u[i], v[i]);
    dest.* = sum;
}

/// Развёртка k×k: k аккумуляторов.
pub fn innerKxK(comptime T: type, comptime k: usize, u: []const T, v: []const T, dest: *T) void {
    var sums: [k]T = @splat(0);
    var i: usize = 0;
    while (i + k <= u.len) : (i += k) {
        inline for (0..k) |j| sums[j] = mulAdd(T, sums[j], u[i + j], v[i + j]);
    }
    while (i < u.len) : (i += 1) sums[0] = mulAdd(T, sums[0], u[i], v[i]);
    var total = sums[0];
    inline for (1..k) |j| total = add(T, total, sums[j]);
    dest.* = total;
}

/// Развёртка k×1a: k произведений складываются между собой и только
/// сумма группы присоединяется к аккумулятору.
pub fn innerKx1a(comptime T: type, comptime k: usize, u: []const T, v: []const T, dest: *T) void {
    var sum: T = 0;
    var i: usize = 0;
    while (i + k <= u.len) : (i += k) {
        var group = mulAdd(T, 0, u[i], v[i]);
        inline for (1..k) |j| group = mulAdd(T, group, u[i + j], v[i + j]);
        sum = add(T, sum, group);
    }
    while (i < u.len) : (i += 1) sum = mulAdd(T, sum, u[i], v[i]);
    dest.* = sum;
}

/// Эталон для тестов: наивно, без всяких развёрток.
fn naive(comptime T: type, u: []const T, v: []const T) T {
    var sum: T = 0;
    for (u, v) |x, y| sum = mulAdd(T, sum, x, y);
    return sum;
}

test "развёртки дают то же скалярное произведение" {
    var u: [13]i64 = undefined;
    var v: [13]i64 = undefined;
    for (&u, &v, 0..) |*x, *y, i| {
        x.* = @intCast(i + 1);
        y.* = @as(i64, @intCast(i)) - 5;
    }
    const lengths = [_]usize{ 0, 1, 5, 6, 13 };
    for (lengths) |n| {
        const expected = naive(i64, u[0..n], v[0..n]);
        var got: i64 = undefined;
        inner4(i64, u[0..n], v[0..n], &got);
        try std.testing.expectEqual(expected, got);
        innerKx1(i64, 6, u[0..n], v[0..n], &got);
        try std.testing.expectEqual(expected, got);
        innerKxK(i64, 6, u[0..n], v[0..n], &got);
        try std.testing.expectEqual(expected, got);
        innerKx1a(i64, 6, u[0..n], v[0..n], &got);
        try std.testing.expectEqual(expected, got);
    }
}

test "дробные сходятся с допуском" {
    var u: [20]f64 = undefined;
    var v: [20]f64 = undefined;
    for (&u, &v, 0..) |*x, *y, i| {
        x.* = 0.5 + @as(f64, @floatFromInt(i));
        y.* = 1.0 / @as(f64, @floatFromInt(i + 1));
    }
    const expected = naive(f64, &u, &v);
    var got: f64 = undefined;
    innerKxK(f64, 6, &u, &v, &got);
    try std.testing.expectApproxEqRel(expected, got, 1e-12);
    innerKx1a(f64, 6, &u, &v, &got);
    try std.testing.expectApproxEqRel(expected, got, 1e-12);
}

1/2 inner.test.развёртки дают то же скалярное произведение...OK
2/2 inner.test.дробные сходятся с допуском...OK
All 2 tests passed.

Граф одной итерации inner4:

   u[i] ──▶ mul ──▶ add ──▶ sum'
   v[i] ──▶  ▲       ▲
             │      sum

Умножение берёт два только что загруженных значения и от аккумулятора не зависит. Процессор считает произведения впереди, на сколько хватит окна, и в цепочку sum входит только сложение. Поэтому нижняя граница inner4 по задержке это задержка сложения: 3.00 на Haswell для дробных и 1.00 для целых, а не восемь и не четыре. Книга приводит такой же замер, 3.00 для double, и спрашивает, почему не больше: ответ в том, что умножение лежит вне критического пути. На M4 граница 2.49 для дробных, замер 2.81, и 1.00 для целых, замер 1.28: чуть выше границ, потому что двух загрузок и умножения на элемент здесь уже достаточно, чтобы обвязка стала заметна.

Вторая граница, по пропускной способности, у скалярного произведения жёстче, чем у свёртки: две загрузки на элемент. На Haswell с двумя портами загрузки это ровно 1.00 такт на элемент для любой скалярной версии, целой или дробной. Отсюда книжные ответы на домашние задачи: 6×1 даёт 1.00 для целых (задержка сложения) и 3.00 для дробных (та же задержка, одна цепочка), 6×6 даёт 1.00 для обоих (загрузки), 6×1a тоже 1.00 для обоих (в цепочке одно сложение на шесть элементов, значит 0.5, но загрузки не пускают ниже 1.00). На M4 три порта загрузки дают границу 0.67, и наши замеры 6×6 и 6×1a на дробных, 1.16 и 0.56, показывают, что до неё ещё есть куда идти и что форма 6×1a здесь оказалась выгоднее шести аккумуляторов. Обе формы ждут тебя в упражнениях.

Мораль всего урока в одной таблице из книги: для double на Haswell combine4 даёт 5.01, 2×2 даёт 2.51, 10×10 даёт 0.52. Десятикратное ускорение без единого нового блока в процессоре, без векторных инструкций, без изменения алгоритма. Только форма цикла. Процессор всё это время был готов, ему не хватало независимых операций.

Практика

Задача повторяет combine6, но в другой сигнатуре: операция передаётся как comptime-функция fn (T, T) T вместе с нейтральным элементом, а число аккумуляторов k это comptime-параметр. В заготовке уже есть addOp(T), mulOp(T) и последовательная combine4, по образцу которой нужно написать combine6(T, op, identity, k, v).

Скрытые тесты проверяют три вещи. Первое: для i64 результат совпадает с combine4 при k из набора 1, 2, 3, 4, 8 и на длинах от нуля до ста, включая длины меньше k и не кратные ему, для сложения и умножения. Второе: для f64 на обычных данных результат совпадает с последовательной суммой с относительной точностью 1e-9, а на векторе 1e16, 1, -1e16, 1, 1e16, 1, -1e16, 1 два аккумулятора дают 4, а не 1, и тест на этом настаивает: если ответы совпали, цепочка осталась одна. Третье, самое строгое: тест передаёт в качестве op функцию со счётчиком и ждёт ровно длина + k - 1 вызовов. Одна операция на элемент и k - 1 на склейку, ни одной лишней. Склеивать аккумуляторы внутри цикла или пересчитывать хвост дважды не получится.

Упражнения

CHK-1★

Нарисуй граф потока данных одной итерации inner4 для f64 (форма домашней задачи 5.13): две загрузки, умножение, сложение, инкремент индекса, сравнение и переход. Обведи рёбра, которые уходят в следующую итерацию, и выпиши критический путь. Дай нижнюю границу CPE по задержке для дробных и для целых на Haswell из книги и на M4 из таблицы урока, и объясни одним предложением, почему задержка умножения в ответ не входит.

CHK-2★★

Напиши inner в форме 6×1 руками, без inline for (форма домашней задачи 5.14): шесть пар u[i + j] * v[i + j], один аккумулятор, хвост. Сверь с innerKx1(T, 6, ...) на длинах 0, 5, 6 и 13. Затем объясни: почему ни одна скалярная версия скалярного произведения на Haswell не может дать CPE ниже 1.00 независимо от типа данных, почему целая 6×1 достигает ровно 1.00, а дробная остаётся на 3.00, хотя работы в ней ровно столько же.

CHK-3★★

Форма домашней задачи 5.15: для innerKxK(T, 6, ...) предскажи CPE на Haswell для целых и дробных до того, как смотреть в книгу, через три границы из виджета: задержка сложения, делённая на шесть, пропускная способность сумматора и умножителя, и две загрузки на элемент при двух портах. Какая из границ побеждает? Потом сними на своей машине inner4, 6×1 и 6×6 маленьким измерителем из урока и объясни каждое число.

CHK-4★★

Форма домашней задачи 5.16: разбери innerKx1a(T, 6, ...). Нарисуй граф одной итерации и покажи, что в цепочку аккумулятора входит одно сложение на шесть элементов. Предскажи CPE на Haswell для целых и дробных и объясни, почему обе формы, шесть аккумуляторов и одна группа из шести произведений, дают на этой машине один и тот же ответ. Затем найди в замере урока строку 6×1a для дробных и объясни, почему на M4 она оказалась ниже строки 6×6.

CHK-5★★★

Собери combineKx1a для целого умножения при k равном четырём как экспортируемую функцию, кросс-компилируй под x86-64 с -femit-asm и найди основной цикл. Сосчитай imul, которые пишут в один и тот же регистр подряд. Потом сделай то же для f64. Объясни разницу через ассоциативность и через то, что оптимизирует компилятор. Наконец, придумай, как заставить компилятор оставить твои скобки у целых (подсказка: посмотри на blackBox в измерителе), и сними CPE до и после.

CHK-6★★★

Сними кривую k×k на своей машине для k от 1 до 16 по всем четырём колонкам и нарисуй её рядом с предсказанием виджета для пресета своей машины. Найди k, на котором каждая колонка легла на границу пропускной способности, и k, после которого CPE снова растёт. Для второго случая посмотри в ассемблер основного цикла и найди там обращения к стеку. Сколько аккумуляторов поместилось в регистры на твоей архитектуре, и как это число связано с тем, что в цикле кроме аккумуляторов живут ещё индекс, длина и указатель?

Итоги

combine4 сидит на границе задержки: одна цепочка зависимостей длиной в вектор, и каждый элемент стоит ровно одну задержку операции, пока блоки простаивают.
Развёртывание k×1 (combine5) сокращает проверки и инкременты, но не трогает цепочку. На всех колонках, кроме целого сложения с задержкой в один такт, CPE не меняется. Компилятор и так разворачивает цикл сам, у нас на восемь.
Условие развёрнутого цикла пишется как i + k <= length, никогда как i < length - k: беззнаковая длина при пустом векторе заворачивается, в Zig это паника в Debug и неопределённое поведение в ReleaseFast. Хвост досчитывается отдельным циклом, и тесты обязаны проверять длины 0, 1, k минус 1, k и 2k плюс 1.
Несколько аккумуляторов (combine6, k×k) делят цепочку на независимые части, и CPE падает до задержки, делённой на число аккумуляторов: 1.50, 1.26 и 1.75 против 2.99, 2.46 и 3.40 на M4. Компилятор упаковал два аккумулятора f64 в две половины одного xmm-регистра, это та же самая идея в другой форме.
Переупорядочение (combine7, k×1a) даёт тот же выигрыш с одним аккумулятором: в цепочку входит одна операция на группу, остальные считаются впереди. На дробных 1.28 и 1.76 против 2.15 и 3.31 у combine5.
Компилятор не переставляет скобки у дробных, потому что результат изменится: одна цепочка на векторе с 1e16 отвечает 1, две цепочки 4. У целых с заворачиванием он переставляет их когда хочет и в любую сторону: наше acc · (x · y) он вернул в (acc · x) · y, и целое умножение в combine7 осталось на 3.00.
Граница пропускной способности это наибольшее из двух: операция на своих блоках и одна загрузка на элемент через порты загрузки. К ней сходится кривая k×k: на Haswell 0.55, 1.00, 1.01, 0.52, на M4 целое умножение дошло до 0.33 при десяти аккумуляторах.
Скалярное произведение упирается в задержку сложения, а не в сумму задержек: умножение вне цепочки. Две загрузки на элемент делают его границу по пропускной способности вдвое жёстче, чем у свёртки: 1.00 на Haswell для любой скалярной формы.
Мерить нужно на тихой машине и с проверкой калибровки: под нагрузкой часы считают чужое время, и все CPE выходят ниже правды, не в одной пропорции.

Дальше

Сегодня мы ускорили один цикл в шесть раз, ничего не меняя в алгоритме, и упёрлись в границу пропускной способности. Но по дороге в трёх местах появились оговорки «если регистров хватает». В следующем уроке мы посмотрим, что происходит, когда их не хватает: аккумуляторы уезжают в стек, и цепочка удлиняется на запись и чтение. Потом два других предела, которых сегодня не было видно, потому что в свёртке нет условий и нет записей в память: предсказатель переходов, который делает непредсказуемое ветвление в горячем цикле дороже любого умножения, и зависимость записи от чтения, из-за которой копирование с перекрытием стоит в семь раз дороже копирования без него. А в конце пойдём дальше книги: раз компилятор сам сложил два аккумулятора в один xmm-регистр, значит регистр шире элемента, и @Vector в Zig позволяет сделать это явно и на четыре, и на восемь элементов сразу. Там же станет ясно, откуда в эталонной таблице строки simd1 и simd4 с CPE 0.19 и 0.25.

домашка

Домашка

Сними на своей машине кривую CPE для k×k при k от 1 до 12 по всем четырём колонкам и отметь, на каком k каждая колонка легла на границу пропускной способности. Сравни с предсказанием виджета для пресета своей машины.L
Собери combine4 для f64 в режиме @setFloatMode(.optimized) на своей целевой платформе и по ассемблеру ответь, разорвал ли компилятор цепочку. Если нет, объясни, почему разрешение переставлять скобки само по себе ничего не гарантирует.M
Для скалярного произведения на своей машине посчитай обе границы (задержка сложения и загрузки на элемент), сними inner4, 6×1, 6×6 и 6×1a и объясни каждое число одним предложением.L