В этом материале речь пойдёт об оптимизациях, которые включает опция -ffast-math при компиляции кода, написанного на C или C++, с использованием GCC 11 для x86_64 Linux (при применении других языков, операционных систем, процессоров могут использоваться немного другие оптимизации).



Опция -ffast-math


Большинство оптимизаций, включаемых с помощью опции -ffast-math, можно включать и отключать в индивидуальном порядке. Эта опция включает их все. Вот их список:

  • -ffinite-math-only
  • -fno-signed-zeros
  • -fno-trapping-math
  • -fassociative-math
  • -fno-math-errno
  • -freciprocal-math
  • -funsafe-math-optimizations
  • -fcx-limited-range

При компиляции кода, написанного на стандартном C (то есть — при применении опции -std=c99 или чего-то в этом роде, что приводит к использованию стандартного C вместо диалекта «GNU C», используемого по умолчанию) опция -ffast-math включает ещё и оптимизацию -ffp-contract=fast, что позволяет компилятору комбинировать операции умножения и сложения с применением инструкции FMA (пример на godbolt.org). Эта оптимизация не влияет на C++ и GNU C, так как для них оптимизация -ffp-contract=fast применяется по умолчанию.

Опция -ffast-math, кроме того, включает оптимизации -fno-signaling-nans, -fno-rounding-math и -fexcess-precision=fast, которые включены по умолчанию при компиляции C- или C++-кода для x86_64 Linux, поэтому тут я об этих оптимизациях не рассказываю.

У линковки программ с применением gcc с опцией -ffast-math есть и ещё одно последствие: благодаря этому процессор воспринимает субнормальные числа как 0.0.

Восприятие субнормальных чисел как 0.0


В формате описания чисел с плавающей точкой имеется специальное представление для значений, близких к 0.0. Появление этих «субнормальных» чисел (их ещё называют «денормализованными») в некоторых случаях приводит к очень большим дополнительным затратам вычислительных ресурсов, так как процессор обрабатывает субнормальные результаты с использованием микрокода для обработки исключений.

Агнер Фог пишет, что, при использовании процессорной микроархитектуры Broadwell, если в результате выполнения операции над нормальными числами получается субнормальное число, на выполнение операции требуется примерно 124 лишних штрафных тактовых цикла.

У x86_64-процессоров есть возможность рассматривать субнормальные числа, как 0.0 и сбрасывать эти числа в значение 0.0, что позволяет устранять ситуации, ведущие к снижению производительности. Включить это можно так:

#define MXCSR_DAZ (1<<6)    /* Включить режим приведения денормализованных чисел к нулю */
#define MXCSR_FTZ (1<<15)   /* Включить режим сброса значений в ноль */

unsigned int mxcsr = __builtin_ia32_stmxcsr();
mxcsr |= MXCSR_DAZ | MXCSR_FTZ;
__builtin_ia32_ldmxcsr(mxcsr);

Линковка программы с применением gcc с опцией -ffast-math приводит к отключению использования субнормальных чисел в этой программе через добавление вышеприведённого кода в глобальный конструктор, который выполняется до функции main.

Опции -ffinite-math-only и -fno-signed-zeros


Проведению многих оптимизаций препятствуют особенности таких значений с плавающей точкой, как NaN, Inf и -0.0. Например:

  • x+0.0 нельзя свести к x, так как это даст неверный результат в том случае, когда x равно -0.0.
  • x-x нельзя свести к 0.0, так как это даст неверный результат в том случае, когда x равно NaN или Inf.
  • x*0.0 нельзя свести к 0.0, так как это даст неверный результат в том случае, когда x равно NaN или Inf.
  • Компилятор не может трансформировать такую конструкцию:

if (x > y) {
  do_something();
} else {
  do_something_else();
}

в такую:

if (x <= y) {
  do_something_else();
} else {
  do_something();
}

Дело в том, что такая трансформация некорректна в том случае, если x или y равняются NaN. (Это преобразование можно счесть удачной оптимизацией при упрощении управления потоком выполнения программы и для обеспечения того, чтобы при обработке общего случая не использовались ветвления).

Опции -ffinite-math-only и -fno-signed-zeros сообщают компилятору о том, что в программе нет вычислений, в результате выполнения которых получаются значения NaN, Inf или -0.0. Поэтому компилятор может применить к коду вышеописанные оптимизации.

Обратите внимание на то, что, если эти флаги установлены, программа может вести себя странно (например — не определять true — или false-части инструкции if) в тех случаях, когда вычисления приводят к появлению значений Inf, NaN или -0.0.

Опция -fno-trapping-math


Имеется возможность включить перехват исключений, имеющих отношение к вычислениям с плавающей точкой, используя функцию feenableexcept GNU libc для генерирования сигнала SIGFPE. Делается это в тех случаях, когда при выполнении некоей инструкции для работы с числами с плавающей точкой происходит переполнение, исчезновение значащих разрядов, появление NaN, и в прочих подобных ситуациях. Например, нижеприведённая функция compute выполняет вычисления, которые приводят к переполнению с появлением Inf, что, при включённом флаге FE_OVERFLOW, приводит к завершению программы с сигналом SIGFPE.

// Это компилируется так: "gcc example.c -D_GNU_SOURCE -O2 -lm"
#include <stdio.h>
#include <fenv.h>

void compute(void) {
  float f = 2.0;
  for (int i = 0; i < 7; ++i) {
    f = f * f;
    printf("%d: f = %f\n", i, f);
  }
}

int main(void) {
  compute();

  printf("\nWith overflow exceptions:\n");
  feenableexcept(FE_OVERFLOW);
  compute();

  return 0;
}

Это значит, что компилятор не может запланировать спекулятивное выполнение инструкций, выполняющих операции с числами с плавающей точкой, так как подобные инструкции могут после этого сгенерировать сигнал в случаях, в которых исходная программа этого не сделает. Например, при вычислении x/y в следующем цикле используются константы, но компилятор не может вынести эту операцию за пределы цикла, так как это может привести к тому, что функция выполнит операцию x/y для случаев, когда все элементы arr будут больше, чем 0.0, что, в результате, приведёт к сбою в случаях, когда в исходной программе сбоев не будет (пример).

double arr[1024];

void foo(int n, double x, double y) {
  for (int i = 0; i < n; ++i) {
    if (arr[i] > 0.0)
      arr[i] = x / y;
  }
}

Ещё один интересный особый случай связан с операциями, в которых используются атомарные числа с плавающей точкой. Стандарт C требует, чтобы исключения, связанные с такими операциями, отбрасывались бы при выполнении комбинированного присваивания (смотрите раздел «Compound assignment» в стандарте C). В результате, если не включена опция -fno-trapping-math, компилятору необходимо добавлять в готовую к выполнению программу дополнительный код (пример).

Передача компилятору опции -fno-trapping-math сообщает ему о том, что программа не допускает возникновения исключений при выполнении вычислений с плавающей точкой, что позволяет компилятору выполнить эти оптимизации.

Опция -fassociative-math


Опция -fassociative-math позволяет выполнять реассоциацию операндов в последовательностях операций с плавающей точкой (а также — выполнять ещё несколько более общих оптимизаций, связанных с переупорядочением выражений). Большинство таких оптимизаций нуждается также в опциях -fno-trapping-math, -ffinite-math-only и -fno-signed-zeros.

Вот некоторые примеры оптимизаций, возможных благодаря -fassociative-math:

Исходное выражение Оптимизированное выражение
(X + Y) - X
Y
(X * Z) + (Y * Z)
(X + Y) * Z
(X * C) + X
X * (C + 1.0) когда C — это константа
(C1 / X) * C2
(C1 * C2) / X когда C1 и C2 — это константы
(C1 - X) < C2
(C1 — C2) > X когда C1 и C2 — это константы

Реассоциация выражений особенно полезна для векторизации вычислений. Рассмотрим, скажем, следующий цикл (пример):

float a[1024];

float foo(void) {
  float sum = 0.0f;
  for (int i = 0; i < 1024; ++i) {
    sum += a[i];
  }
  return sum;
}

Все операции прибавления значений к sum выполняются последовательно, в результате в обычных условиях эти вычисления векторизовать нельзя. Но опция -fassociative-math позволяет компилятору изменить порядок операций:

sum = (a[0] + a[4] + ... + a[1020]) + (a[1] + a[5] + ... + a[1021]) + (a[2] + a[6] + ... + a[1022]) + (a[3] + a[7] + ... + a[1023]);

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

float a[1024];

float foo(void) {
  float sum0 = sum1 = sum2 = sum3 = 0.0f;
  for (int i = 0; i < 1024; i += 4) {
    sum0 += a[i    ];
    sum1 += a[i + 1];
    sum2 += a[i + 2];
    sum3 += a[i + 3];
  }
  return sum0 + sum1 + sum2 + sum3;
}

А такая конструкция легко поддаётся векторизации.

Опция -fno-math-errno


Математические функции C могут устанавливать errno в том случае, если вызываются с некорректными аргументами. Этот возможный побочный эффект вычислений означает, что компилятор должен вызывать функции libc вместо использования инструкций, которые могут напрямую вычислить результат.

Компилятор может, в некоторых случаях, смягчить эту проблему, не вызывая функцию тогда, когда ему известно, что операция будет выполнена успешно (пример):

double foo(double x) {
  return sqrt(x);
}

Эта конструкция компилируется в исполняемый код с использованием инструкции sqrtsd в том случае, если входное значение находится в допустимом диапазоне, а sqrt вызывается лишь для входных значений, вычисления над которыми приводят к появлению NaN:

foo:
        pxor    xmm1, xmm1
        ucomisd xmm1, xmm0
        ja      .L10
        sqrtsd  xmm0, xmm0
        ret
.L10:
        jmp     sqrt

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

Опция -fno-math-errno заставляет GCC оптимизировать все математические функции так, как будто они не устанавливают errno (то есть — компилятору не нужно вызывать функции libc в том случае, если процессорная архитектура, для которой компилируется код, обладает подходящими инструкциями).

▍Нематематические функции


Одним неожиданным эффектом применения опции -fno-math-errno является тот факт, что она заставляет GCC считать, что функции libc, ответственные за выделение памяти (вроде malloc и strdup) не устанавливают errno. Это можно видеть в следующей функции, где использование -fno-math-errno приводит к тому, что GCC убирает вызовы perror (пример):

void *foo(size_t size) {
  errno = 0;
  void *p = malloc(size);
  if (p == NULL) {
    if (errno)
      perror("error");
    exit(1);
  }
  return p;
}

По этому поводу было сделано сообщение об ошибке GCC 88576.

Опция -freciprocal-math


Опция -freciprocal-math позволяет компилятору вычислять выражения вида x/y как x*(1/y). Это полезно для тех случаев, когда приходится иметь дело с такими конструкциями (пример):

float length = sqrtf(x*x + y*y + z*z);
x = x / length;
y = y / length;
z = z / length;

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

float t = 1.0f / sqrtf(x*x + y*y + z*z);
x = x * t;
y = y * t;
z = z * t;

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

Опция -funsafe-math-optimizations


Опция -funsafe-math-optimizations включает различные «математически корректные» оптимизации, в результате применения которых могут изменяться результаты вычислений, что вызвано особенностями устройства чисел с плавающей точкой. Вот несколько примеров:

Исходное выражение Оптимизированное выражение
sqrt(x)*sqrt(x)
x
sqrt(x)*sqrt(y)
sqrt(x*y)
exp(x)*exp(y)
exp(x+y)
x/exp(y)
x*exp(-y)
x*pow(x,c)
pow(x,c+1)
pow(x,0.5)
sqrt(x)
(int)log(d)
ilog(d)
sin(x)/cos(x)
tan(x)

Обратите внимание на то, что для применения многих из этих оптимизаций нужны дополнительные флаги, например — такие, как -ffinite-math-only и -fno-math-errno.

Опция -funsafe-math-optimizations, кроме того, включает следующие оптимизации:

  • -fno-signed-zeros
  • -fno-trapping-math
  • -fassociative-math
  • -freciprocal-math

Включает она (для стандартного C) и -ffp-contract=fast, а так же — восприятие субнормальных чисел как 0.0 (так же, как был описано в разделе о -ffast-math).

Опция -fcx-limited-range


Математические формулы для умножения и деления комплексных чисел выглядят так:

a+ ibc+id=ac-bd+i(bc+ad)

a+ibc+id= ac+bdc2+d2+ibc-adc2+d2

Но это не очень хорошо подходит для значений с плавающей точкой.

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

double complex
mul(double a, double b, double c, double d) {
  double ac, bd, ad, bc, x, y;
  double complex res;

  ac = a * c;
  bd = b * d;
  ad = a * d;
  bc = b * c;

  x = ac - bd;
  y = ad + bc;

  if (isnan(x) && isnan(y)) {
    /* Восстановить бесконечности, вычисленные как NaN + iNaN.  */
    _Bool recalc = 0;
    if (isinf(a) || isinf(b)) {
      /* z - бесконечность.  "Упаковать" бесконечность и сделать значения NaN
       * в другом показателе равными 0.  */
      a = copysign(isinf(a) ? 1.0 : 0.0, a);
      b = copysign(isinf(b) ? 1.0 : 0.0, b);
      if (isnan(c)) c = copysign(0.0, c);
      if (isnan(d)) d = copysign(0.0, d);
      recalc = 1;
    }
    if (isinf(c) || isinf(d)) {
      /* w - бесконечность. "Упаковать" бесконечность и сделать значения NaN
       * в другом показателе равными 0.  */
      c = copysign(isinf(c) ? 1.0 : 0.0, c);
      d = copysign(isinf(d) ? 1.0 : 0.0, d);
      if (isnan(a)) a = copysign(0.0, a);
      if (isnan(b)) b = copysign(0.0, b);
      recalc = 1;
    }
    if (!recalc
        && (isinf(ac) || isinf(bd)
            || isinf(ad) || isinf(bc))) {
      /* Восстановить бесконечности из переполнения путём установки
       * значений NaN в 0.  */
      if (isnan(a)) a = copysign(0.0, a);
      if (isnan(b)) b = copysign(0.0, b);
      if (isnan(c)) c = copysign(0.0, c);
      if (isnan(d)) d = copysign(0.0, d);
      recalc = 1;
    }
    if (recalc) {
      x = INFINITY * (a * c - b * d);
      y = INFINITY * (a * d + b * c);
    }
  }

  __real__ res = x;
  __imag__ res = y;
  return res;
}

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

Применение опции -fcx-limited_range приводит к тому, что компилятор использует обычные математические формулы для умножения и деления комплексных чисел.

Пользуетесь ли вы оптимизациями математических вычислений в своих С/С++-проектах?

Комментарии (5)


  1. YuryB
    01.11.2021 22:14

     for (int i = 0; i < 1024; i += 4) {
        sum0 += a[i    ];
        sum1 += a[i + 1];
        sum2 += a[i + 2];
        sum3 += a[i + 3];
      }
      return sum0 + sum1 + sum2 + sum3;

    А такая конструкция легко поддаётся векторизации.

    такую конструкцию даже векторизировать не надо, она прекрасно работает в 4 раза быстрее чем обычное инкрементирование, т.к. давно в ядрах процессора блоки суммирования и другие продублированы (и их как раз штуки 3-4 минимум) и процессор без проблем выполнит 4 сложения за условный такт если нет зависимости по данным. чтобы извлечь заметный плюс от векторизации для простого суммирования нужно складывать не 4 а 8 интов, или байты, их в AVX много влезет, а обычных сумматоров всё же меньше будет.


    1. tzlom
      02.11.2021 01:31

      4 сложения подряд действительно быстрее чем 4*время на одно сложение, но все таки медленнее чем одна SSE инструкция, ьак что векторизовать стоит.


    1. dudinroman
      02.11.2021 12:03
      +3

      Ну вот зачем писать, если сами не разбираетесь? В современных x64 ядрах по два функциональных устройства для сложения-умножения чисел с плавающей точкой, а не "минимум 3-4". Вы можете в этом убедиться, посмотрев диаграму ядра Skylake (https://en.wikichip.org/w/images/7/7e/skylake_block_diagram.svg), и поискав там прямоугольнички с текстом "FP FMA". Либо можете зайти на uops.info (https://uops.info/table.html), и увидеть что на всех современных микроархитектурах у сложения пропускная способность 0.5 (то есть они делают два сложения за такт).

      Автор статьи верно утверждает, что такая конструкция поддается векторизации, потому что понимает, что использование четырех независимых аккумуляторов никак не связано с векторизацией. Действительно, те самые "FP FMA" - это векторные устройства, которые могут складывать по 128/256/512-битному вектору за раз. Вот так будет выглядеть этот цикл, если его векторизовать (для примера используем SSE):

      __m128 sum0 = _mm_setzero_ps();
      __m128 sum1 = _mm_setzero_ps();
      __m128 sum2 = _mm_setzero_ps();
      __m128 sum3 = _mm_setzero_ps();
      
      for (int i = 0; i < 1024; i += 16) {
          _mm_add_ps(sum0, _mm_load_ps(a + i + 0));
          _mm_add_ps(sum1, _mm_load_ps(a + i + 4));
          _mm_add_ps(sum2, _mm_load_ps(a + i + 8));
          _mm_add_ps(sum3, _mm_load_ps(a + i + 12));
      }
      
      _mm_add_ps(sum0, sum1);
      _mm_add_ps(sum2, sum3);
      _mm_add_ps(sum0, sum2);
      
      // TODO: складываем четыре линии sum0, получаем результат

      Наконец, почему именно четыре аккумулятора? Я вообще не понял что вы написали (да вы и сами не поняли, скорее всего). Настоящий ответ тут такой: в Skylake (опять же, для примера) сложение имеет глубину конвейера, равную четырем (см uops.info), то есть одновременно "in-flight" может находится четыре сложения, и в таком случае конвейер будет полностью загружен. Однако, так как функциональных устройства два, правильнее будет использовать восемь аккумуляторов. У какого-нибудь Zen2, например, латентность уже три, а не четыре, а значит аккумуляторов нужно делать шесть.


      1. YuryB
        02.11.2021 13:58

        а ну да, FP 2, а целочисленных 4, спутал здесь.

        спасибо за интересный комментарий!


        1. dudinroman
          02.11.2021 14:21

          Почти. Скалярных целочисленных ALU действительно четыре. А вот векторных три (см. опять же диаграму на wikichips "INT Vect ALU" либо инструкцию padd на uops.info)