Почему GCC или Clang не оптимизируют обратную 1 инструкцию при использовании fast-math

Question

Кто-нибудь знает, почему GCC/Clang не будет оптимизировать функцию test1 в приведенном ниже примере кода, чтобы просто использовать только инструкцию RCPPS при использовании опции fast-math? Есть ли другой флаг компилятора, который будет генерировать этот код?

typedef float float4 __attribute__((vector_size(16)));

float4 test1(float4 v)
{
    return 1.0f / v;
}

Вы можете увидеть скомпилированный вывод здесь: https://goo.gl/jXsqat

Answer 1

Потому что точность RCPPS намного ниже, чем float разделение.

Опция, позволяющая включить эту оптимизацию, не подходит для -ffast-math,

Параметры цели x86 руководства gcc говорят, что на самом деле есть опция, которая (с -ffast-math) действительно заставляет gcc использовать их (с итерацией Ньютона-Рафсона):

• -mrecip Эта опция позволяет использовать команды RCPSS и RSQRTSS (и их векторизованные варианты RCPPS и RSQRTPS) с дополнительным шагом Ньютона-Рафсона для повышения точности вместо DIVSS и SQRTSS (и их векторизованных вариантов) для аргументов с плавающей точкой одинарной точности. Эти инструкции генерируются только тогда, когда -funsafe-math-optimizations включен вместе с -finite-math-only и -fno-trapping-math. Обратите внимание, что в то время как пропускная способность последовательности выше, чем пропускная способность невзаимной команды, точность последовательности может быть уменьшена до 2 ulp (то есть обратное значение 1,0 равно 0,999999994).

  Обратите внимание, что GCC реализует 1.0f/sqrtf(x) в терминах RSQRTSS (или RSQRTPS) уже с -ffast-math (или с вышеуказанной комбинацией параметров) и не нуждается в -mrecip.

  Также обратите внимание, что GCC испускает указанную выше последовательность с дополнительным шагом Ньютона-Рафсона для векторизованного деления с плавающей запятой и векторизованного sqrtf (x) уже с -ffast-math (или вышеуказанной комбинацией опций), и не нуждается в -mrecip.

• -mrecip=opt

Эта опция контролирует, какие инструкции обратной оценки могут быть использованы. opt - это список параметров, разделенных запятыми, которым может предшествовать '!' инвертировать опцию:

’all’
      Enable all estimate instructions.
‘default’
    Enable the default instructions, equivalent to -mrecip.
‘none’
    Disable all estimate instructions, equivalent to -mno-recip.
‘div’
    Enable the approximation for scalar division.
‘vec-div’
    Enable the approximation for vectorized division.
‘sqrt’
    Enable the approximation for scalar square root.
‘vec-sqrt’
    Enable the approximation for vectorized square root. 

Так, например, -mrecip=all,! Sqrt включает все обратные приближения, кроме квадратного корня.

Обратите внимание, что новый дизайн Intel Skylake дополнительно улучшает производительность деления FP, задержку до 8-11c, пропускную способность 1/3c. (Или один на 5c пропускной способности для 256b векторов, но такая же задержка для vdivps). Они расширили разделители, поэтому AVX vdivps ymm теперь та же задержка, что и для 128b векторов.

(SnB to Haswell сделал 256b div и sqrt с примерно вдвое большей задержкой / пропускной способностью, поэтому у них явно было только делители шириной 128b.) Skylake также конвейеризует обе операции больше, так что в полете может выполняться около 4 операций div. sqrt тоже быстрее.

Так что через несколько лет, когда Skylake получит широкое распространение, это будет стоить rcpps если вам нужно разделить на одну и ту же вещь несколько раз. rcpps и пара fma может иметь немного более высокую пропускную способность, но худшую задержку. Также, vdivps только один моп; поэтому больше ресурсов для выполнения будет доступно одновременно с делением.

Еще неизвестно, на что будет похожа первоначальная реализация AVX512. предположительно rcpps и пара FMA для итераций Ньютона-Рафсона будет победой, если производительность деления FP является узким местом. Если пропускная способность Uop является узким местом, и есть много другой работы, пока подразделения находятся в полете, vdivps zmm вероятно, все еще хорош (если, конечно, один и тот же делитель не используется повторно).

Answer 2

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

Вот набор инструкций, сгенерированных gcc, аннотированных мной выполненными вычислениями:

test1(float __vector): ; xmm0               = a
    rcpps   xmm1, xmm0 ; xmm1 = 1 / xmm0    = 1/a
    mulps   xmm0, xmm1 ; xmm0 = xmm0 * xmm1 = a * 1/a
    mulps   xmm0, xmm1 ; xmm0 = xmm0 * xmm1 = a * (1/a)^2
    addps   xmm1, xmm1 ; xmm1 = xmm1 + xmm1 = 2 * (1/a)
    subps   xmm1, xmm0 ; xmm1 = xmm1 - xmm0 = 2 * (1/a) - a * (1/a)^2
    movaps  xmm0, xmm1 ; xmm0 = xmm1        = 2 * (1/a) - a * (1/a)^2
    ret

Так что здесь происходит? Зачем тратить дополнительные 4 инструкции на вычисление выражения, которое математически эквивалентно просто обратной величине?

Ну а rcppsИнструкции являются лишь приблизительными. Остальные арифметические инструкции (mulps, addps, subps) точны с точностью до одинарной точности. Давайте напишем r(x)для приближенной обратной функции. Тогда финал становится y = 2*r(a) - a*r(a)^2. Если мы заменим r(x) = (1 + eps) * (1/x), с участием eps будучи относительной ошибкой, получаем:

y = 2 * (1 + eps) * (1/a) - a * (1 + eps)^2 * (1/a)^2
  = (2 + 2*eps - (1 + eps)^2) * (1/a)
  = (2 + 2*eps - (1 + 2*eps + eps^2)) * (1/a)
  = (1 - eps^2) * (1/a)

Относительная ошибка rcpps меньше чем 1.5 * 2^-12, так eps <= 1.5 * 2^-12, так:

eps^2 <= 2.25 * 2^-24
      <  1.5  * 2^-23

Таким образом, выполнив эти дополнительные инструкции, мы перешли с 12-битной точности на 23-битную. Обратите внимание, что число с плавающей запятой одинарной точности имеет 24 бита точности, поэтому здесь мы получаем почти полную точность.

Так это просто волшебная последовательность инструкций, которая помогает нам добиться дополнительной точности? Не совсем. Он основан на методе Ньютона (который, как я понимаю, люди, которые много работают со сборкой, называют его Ньютоном-Рафсоном).

Метод Ньютона - это метод поиска корней. Учитывая некоторую функцию f(x) находит приближенные решения f(x) = 0, начиная с приближенного решения x_0и многократно улучшая его. Итерация Ньютона определяется следующим образом:

x_n+1 = x_n - f(x_n) / f'(x_n)

В нашем случае мы можем переформулировать нахождение обратного 1/a из a как найти корень функции f(x) = a*x - 1, с производной f'(x) = a. Подставляя это в уравнение для итерации Ньютона, мы получаем:

x_n+1 = x_n - (a*x_n - 1) / a

Два наблюдения:

1. В этом случае итерация Ньютона фактически дает нам точный результат, а не просто лучшее приближение. Это имеет смысл, потому что метод Ньютона работает путем линейного приближения f вокруг x_n. В этом случае fполностью линейна, поэтому приближение идеальное. Тем не мение...

2. Для вычисления итерации Ньютона необходимо разделить на a, которое является точным вычислением, которое мы пытаемся приблизить. Это создает круговую проблему. Мы прерываем цикл, изменяя итерацию Ньютона, чтобы использовать нашу приближенную обратную x_n для деления на a:

   x_n+1 =  x_n - (a*x_n - 1) * x_n
         ~= x_n - (a*x_n - 1) / a
   

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

Чтобы этого избежать, можно переписать итерацию:

x_n+1 = x_n - (a*x_n - 1) * x_n
      = x_n - (a*x_n^2 - x_n)
      = 2*x_n - a*x_n^2

Теперь замените x_0 = r(a) и мы восстанавливаем наше выражение сверху:

y = x_1 = 2*r(a) - a*r(a)^2