Оптимизация математических вычислений и опция -ffast-math в GCC 11

Опция -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 

(C1 / X) * C2 

(C1 - X) < 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 приводит к тому, что компилятор использует обычные математические формулы для умножения и деления комплексных чисел.

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