Допустим, у нас есть массив чисел с плавающей запятой, и мы хотим их суммировать. Можно наивно подумать, что их достаточно просто сложить, например, на Rust:

fn naive_sum(arr: &[f32]) -> f32 {
    let mut out = 0.0;
    for x in arr {
        out += *x;
    }
    out
}

Однако это запросто может привести к произвольно большой накопленной погрешности. Давайте проверим:

naive_sum(&vec![1.0;     1_000_000]) =  1000000.0
naive_sum(&vec![1.0;    10_000_000]) = 10000000.0
naive_sum(&vec![1.0;   100_000_000]) = 16777216.0
naive_sum(&vec![1.0; 1_000_000_000]) = 16777216.0

Ой-ёй… Что произошло? Проблема в том .что следующее 32-битное число с плавающей запятой после 16777216 — это 16777218. Так что при вычислении 16777216 + 1, значение округляется до ближайшего числа с плавающей запятой, имеющей чётную мантиссу, то есть снова до 16777216. Мы зашли в тупик.

К счастью, есть более совершенные способы суммирования массива.

Попарное суммирование

Чуть более умный способ — это попарное суммирование. Вместо абсолютно линейной суммы с единственным накопителем он рекурсивно суммирует массив, разбивая его пополам, суммируя половины, а затем складывая суммы.

fn pairwise_sum(arr: &[f32]) -> f32 {
    if arr.len() == 0 { return 0.0; }
    if arr.len() == 1 { return arr[0]; }
    let (first, second) = arr.split_at(arr.len() / 2);
    pairwise_sum(first) + pairwise_sum(second)
}

Этот способ более точен:

pairwise_sum(&vec![1.0;     1_000_000]) =    1000000.0
pairwise_sum(&vec![1.0;    10_000_000]) =   10000000.0
pairwise_sum(&vec![1.0;   100_000_000]) =  100000000.0
pairwise_sum(&vec![1.0; 1_000_000_000]) = 1000000000.0

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

Выбрав точку разделения как кратное 256, мы гарантируем, что базовый случай в рекурсии всегда будет равен точно 256 элементам, за исключением самого последнего блока. Это гарантирует использование самого оптимального редуцирования и всегда корректно прогнозирует условие цикла. Эта небольшая деталь повысила производительность в случае больших массивов на 40%!

fn block_pairwise_sum(arr: &[f32]) -> f32 {
    if arr.len() > 256 {
        let split = (arr.len() / 2).next_multiple_of(256);
        let (first, second) = arr.split_at(split);
        block_pairwise_sum(first) + block_pairwise_sum(second)
    } else {
        naive_sum(arr)
    }
}

Суммирование Кэхэна

Наихудший случай погрешности округления наивного сумирования масштабируется при суммировании ? элементов как ?(??), где ? — это машинный ноль вашего типа с плавающей запятой (здесь это 2⁻²⁴). Показатели попарного суммирования улучшают этот результат до ?((log⁡?)?+??²). Суммирование Кэхэна ещё больше улучшает показатель (до ?(??²)), полностью избавляя нас от члена ? и оставляя только ?² , который несущественен, если только вы не суммируете очень большое количество чисел.

Все эти границы масштабируются как ∑?∣??∣, то есть для суммирования Кэхена граница абсолютной погрешности в наихудшем случае всё равно остаётся квадратичной относительно ?.

На практике показатели всех алгоритмов суммирования существенно лучше, чем границы их наихудших случаев, так как в большинстве ситуаций погрешности округляются не только строго вверх или вниз, а в среднем обнуляют друг друга.

pub fn kahan_sum(arr: &[f32]) -> f32 {
    let mut sum = 0.0;
    let mut c = 0.0;
    for x in arr {
        let y = *x - c;
        let t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
    sum
}

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

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

pub fn block_kahan_sum(arr: &[f32]) -> f32 {
    let mut sum = 0.0;
    let mut c = 0.0;
    for chunk in arr.chunks(256) {
        let x = naive_sum(chunk);
        let y = x - c;
        let t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
    sum
}

Точное суммирование

Я знаю по крайней мере два обобщённых способа получения корректно округлённой суммы последовательности чисел с плавающей запятой. Они логически вычисляют сумму с бесконечной точностью, а в конце округляют её до значения с плавающей запятой.

Первый способ основан на примитиве 2Sum, то есть на избавленном от погрешностей преобразовании из двух чисел ?,?  в ?,? так, что ?+?=?+?, где ? — это малая погрешность. Многократно применяя это преобразование, пока погрешность не пропадёт, мы можем получить корректно округлённую сумму. Сложность заключается в отслеживании того, что и в каком порядке складывать, а для исчезновения всех членов в наихудшем случае требуется ?(?²) сложений. Именно это реализовано в math.fsum языка Python и в крейте Rust fsum, которые используют дополнительную память для хранения частичных сумм. Крейт accurate тоже реализует это при помощи преобразования на месте в i_fast_sum_in_place.

Другой способ заключается в хранении большого буфера целых чисел, по одному на экспоненту. В случае прибавления числа с плавающей запятой оно раскладывается на экспоненту и мантиссу, а мантисса прибавляется к соответствующему целому числу в буфере. Если происходит переполнение целого числа buf[i], то выполняется инкремент целого числа в buf[i + w], где w — ширина целого числа.

На самом деле, этот способ позволяет вычислять абсолютно точную сумму без какого-либо округления; по сути, это более свободное представление числа с фиксированной запятой, оптимизированное для накопления чисел с плавающей запятой. Этот способ требует времени ?(?), но использует большой, хотя и константный объём памяти (≈≈ 1 КБ в случае f32, ≈≈ 16 КБ в случае f64). Преимущество этого способа ещё и в том, что это онлайновый алгоритм — и прибавление числа к сумме, и получение текущего общего значения занимают амортизированное ?(1).

Вариант этого способа реализован в крейте accurate как крейт OnlineExactSum, в котором для буфера используются числа с плавающей запятой, а не целые.

Раскрываем мощь компилятора

Кроме точности, у naive_sum есть ещё одна проблема. Компилятор Rust не может изменять порядок сложения чисел с плавающей запятой, потому что такое сложение не ассоциативно. Поэтому он не может автоматически векторизировать naive_sum, чтобы та использовала для вычисления суммы команды SIMD, и не может применять параллелизм на уровне команд.

Для решения этой проблемы в компиляторе Rust есть встроенные средства (intrinsics), выполняющие суммирование чисел с плавающей запятой, в то же время обеспечивая ассоциативность, например, std::intrinsics::fadd_fast. Однако такие команды невероятно опасны, потому что они предполагают, что и входные, и выходные значения являются конечными числами (не бесконечностями и не NaN), в противном случае они оказываются неопределённым поведением. Функционально это делает их неприменимыми, поскольку только в самых ограниченных ситуациях при вычислениях суммы мы знаем, что все входные числа являются конечными и что их сумма не может привести к переполнению.

Недавно я поделился своим разочарованием этими операторами с Беном Кимоком, и вместе мы предложили (а он реализовал) новое множество операторов: std::intrinsics::fadd_algebraic и ему подобные.

Я предложил, чтобы мы назвали эти операторы алгебраическими (algebraic), потому что они позволяют (теоретически) выполнять любые преобразования, допустимые в вещественной алгебре. Например, подстановки ?−?→0, ??+??→?(?+?) или ?⁶→(?²)³. В общем случае эти операторы обрабатываются, как будто они выполняются с использованием вещественных чисел, и их можно отобразить на любое множество команд с плавающей запятой, которые будут эквивалентом исходного выражения при условии, что команды с плавающей запятой будут точными.

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

Кроме того, они позволяют генерировать команды fused multiply-add , так как в вещественной алгебре fma⁡(?,?,?)=??+?.

Благодаря этим новым командам можно очень легко сгенерировать автоматически векторизированную сумму:

#![allow(internal_features)]
#![feature(core_intrinsics)]
use std::intrinsics::fadd_algebraic;

fn naive_sum_autovec(arr: &[f32]) -> f32 {
    let mut out = 0.0;
    for x in arr {
        out = fadd_algebraic(out, *x);
    }
    out
}

Если выполнить компиляцию с -C target-cpu=broadwell , то мы увидим, что компилятор автоматически сгенерировал следующий короткий цикл с четырьмя накопителями и командами AVX2:

.LBB0_5:
    vaddps  ymm0, ymm0, ymmword ptr [rdi + 4*r8]
    vaddps  ymm1, ymm1, ymmword ptr [rdi + 4*r8 + 32]
    vaddps  ymm2, ymm2, ymmword ptr [rdi + 4*r8 + 64]
    vaddps  ymm3, ymm3, ymmword ptr [rdi + 4*r8 + 96]
    add     r8, 32
    cmp     rdx, r8
    jne     .LBB0_5

Он обрабатывает 128 байтов данных с плавающей запятой (то есть 32 элемента) за 7 команд. Кроме того, все команды vaddps не зависят друг от друга, потому что накопление выполняется в разные регистры. Если проанализировать код при помощи uiCA, то можно увидеть. что по его оценке выполнение этого цикла занимает 4 такта, то есть обрабатывается по 32 байта на цикл. При частоте 4 ГГц это до 128 ГБ/с! Стоит отметить, что это гораздо больше, чем пропускная способность ОЗУ моей машины, так что этой скорости можно достичь, только суммируя данные, уже находящиеся в кэше.

Помня об этом, мы также может легко определить block_pairwise_sum_autovec и block_kahan_sum_autovec, заменив их вызовы к naive_sum  на naive_sum_autovec.

Точность и скорость

Давайте попробуем сравнить разные способы суммирования. Попробуем в качестве бенчмарка просуммировать 100000 случайных float в интервале от -100000 до +100,000. Это 400 КБ данных, так что всё равно уместится в кэш моего AMD Threadripper 2950x.

Весь код выложен на Github. Компиляция выполнена с RUSTFLAGS=-C target-cpu=native и --release. Я получил следующие результаты:

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

Во-первых, хотелось бы отметить, что разница производительности между самым быстрым и самым медленным способом составляет более чем 100 раз. И это для суммирования массива! Это может быть не совсем честно, потому что самые медленные способы вычисляют нечто существенно более сложное, но всё равно существует разница в двадцать раз между достаточно разумной наивной реализацией и самым быстрым способом.

Мы обнаружили, что в общем случае способы _autovec, использующие fadd_algebraic, быстрее и точнее, чем те, которые используют обычное сложение с плавающей запятой. Они точнее по той же причине, по которой точнее попарная сумма: любое изменение порядка сложений лучше, поскольку используемая по умолчанию длинная цепочка сложений уже является наихудшим случаем точности точности в сумме.

Ограничившись Парето-оптимальными вариантами, мы получим следующие четыре реализации:

Отметим, что разницы в реализации могут влиять достаточно сильно; к тому же, вероятно, существуют ещё десятки способов компенсированного суммирования, которые я не сравнивал.

Думаю, в большинстве случаев здесь выигрывает block_kahan_autovec, он имеет хорошую точность (не происходит вырождение при больших входных значениях) и почти максимальную скорость. В большинстве сценариев использования дополнительная точность, получаемая благодаря корректному округлению сумм, не требуется, а она в 50-100 раз медленнее.

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

#![allow(internal_features)]
#![feature(core_intrinsics)]
use std::intrinsics::fadd_algebraic;

fn sum_block(arr: &[f32]) -> f32 {
    arr.iter().fold(0.0, |x, y| fadd_algebraic(x, *y))
}

pub fn sum_orlp(arr: &[f32]) -> f32 {
    let mut chunks = arr.chunks_exact(256);
    let mut sum = 0.0;
    let mut c = 0.0;
    for chunk in &mut chunks {
        let y = sum_block(chunk) - c;
        let t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
    sum + (sum_block(chunks.remainder()) - c)
}

Разумеется, можно попробовать менять 256 на другие значения; я выяснил, что при 128 производительность снижается примерно на 20%, а при 512 производительность особо не повышается, но страдает точность.

Заключение

Я считаю, что fadd_algebraic и схожие алгебраические intrinsics очень полезны в достижении высокоскоростных процедур работы с плавающей запятой, и что их нужно добавить и в другие языки. Глобальный -ffast-math недостаточно хорош, так как мы уже видели, что наилучшая реализация оказалась гибридом между автоматически оптимизированными для скорости вычислениями и вручную реализованными неассоциативными компенсированными операциями.

Кроме того, при использовании LLVM следует опасаться -ffast-math. Когда этот флаг установлен в LLVM, получение NaN или бесконечности становится неопределённым поведением. Понятия не имею, почему разработчики выбрали такую жёсткую позицию, делающую ненадёжной практически все использующие его программы. Если вы используете в качестве целевой платформы для своего языка LLVM, то избегайте флагов fast-math nnan и ninf.