Почему скаляр SSE sqrt(x) медленнее, чем rsqrt(x) * x?

Question

Я профилировал некоторые наши основные математические расчеты на Intel Core Duo, и, глядя на различные подходы к квадратному корню, я заметил кое-что странное: используя скалярные операции SSE, быстрее взять взаимный квадратный корень и умножить его получить sqrt, чем использовать собственный код операции sqrt!

Я проверяю это с помощью цикла что-то вроде:

inline float TestSqrtFunction( float in );

void TestFunc()
{
  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    {
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    }
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );
}

Я пробовал это с несколькими различными телами для TestSqrtFunction, и у меня есть некоторые моменты, которые действительно царапают мою голову. Хуже всего было использовать встроенную функцию sqrt() и позволить "умному" компилятору "оптимизировать". На 24ns/float, используя x87 FPU, это было патетически плохо:

inline float TestSqrtFunction( float in )
{  return sqrt(in); }

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

inline void SSESqrt( float * restrict pOut, float * restrict pIn )
{
   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss
}

Это было лучше, на 11,9 нс / плавать. Я также попробовал причудливую методику аппроксимации Ньютона-Рафсона Кармака, которая работала даже лучше, чем аппаратная, на скорости 4,3 нс / с плавающей точкой, хотя с ошибкой 1 в 2¹⁰ (что слишком много для моих целей).

Ошеломление было, когда я попробовал SSE op для обратного квадратного корня, а затем использовал умножение, чтобы получить квадратный корень ( x * 1/√x = √x). Несмотря на то, что для этого требуются две зависимые операции, на сегодняшний день это было самое быстрое решение с 1,24 нс / с плавающей запятой и точностью до 2^-14:

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )
{
   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss
}

Мой вопрос в основном, что дает? Почему встроенный в аппаратный код квадратного корня SSE медленнее, чем синтезировать его из двух других математических операций?

Я уверен, что это действительно стоимость самой операции, потому что я проверил:

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

(edit: stephentyrone правильно указывает, что операции с длинными строками чисел должны использовать векторизацию SIMD-упакованных операций, например, rsqrtps - но структура данных массива здесь только для целей тестирования: я действительно пытаюсь измерить скалярную производительность для использования в коде, который нельзя векторизовать.)

Answer 1

sqrtss дает правильно округленный результат. rsqrtss дает приближение к обратной, с точностью до 11 бит.

sqrtss дает гораздо более точный результат, когда требуется точность. rsqrtss существует для случаев, когда достаточно приближения, но требуется скорость. Если вы прочитаете документацию Intel, вы также найдете последовательность инструкций (обратное приближение квадратного корня, за которым следует один шаг Ньютона-Рафсона), которая дает почти полную точность (~23 бита, если я правильно помню), и все еще в некоторой степени быстрее, чем sqrtss,

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

Answer 2

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

• Инструкции rsqrt* вычисляют приближение к обратному квадратному корню, что составляет примерно 11–12 бит.
• Это реализовано с помощью таблицы поиска (то есть ПЗУ), индексированной мантиссой. (Фактически, это сжатая справочная таблица, похожая на старые математические таблицы, в которой используются корректировки младших битов для экономии на транзисторах.)
• Причина, по которой он доступен, заключается в том, что это начальная оценка, используемая FPU для "реального" алгоритма извлечения квадратного корня.
• Также есть примерная обратная инструкция rcp. Обе эти инструкции являются ключом к пониманию того, как FPU реализует извлечение квадратного корня и деление.

Вот что ошиблось в консенсусе:

• FPU эпохи SSE не используют метод Ньютона-Рафсона для вычисления квадратных корней. Это отличный программный метод, но было бы ошибкой реализовывать его таким образом на оборудовании.

Алгоритм NR для вычисления обратного квадратного корня имеет этот этап обновления, как отмечали другие:

x' = 0.5 * x * (3 - n*x*x);

Это много умножений, зависящих от данных, и одного вычитания.

Далее следует алгоритм, который фактически используют современные FPU.

Данный b[0] = n, предположим, что мы можем найти ряд чисел Y[i] такой, что b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2 подходы 1. Затем рассмотрите:

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

Ясно x[n] подходы sqrt(n) а также y[n] подходы 1/sqrt(n).

Мы можем использовать шаг обновления Ньютона-Рафсона для получения обратного квадратного корня, чтобы получить хорошее Y[i]:

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

Потом:

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

а также:

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

Следующее ключевое наблюдение: b[i] = x[i-1] * y[i-1]. Так:

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

Потом:

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

То есть, учитывая начальные x и y, мы можем использовать следующий шаг обновления:

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

Или, что еще интереснее, мы можем установить h = 0.5 * y. Это инициализация:

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

И это шаг обновления:

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

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

В 1999 году FPU уже нуждались в конвейерной схеме сложения / вычитания и конвейерной схеме умножения, иначе SSE не был бы очень "потоковым". В 1999 году потребовалась только одна из каждой схемы, чтобы реализовать этот внутренний цикл полностью конвейерным способом, не тратя много оборудования только на извлечение квадратного корня.

Сегодня, конечно, мы представили программисту слияние умножения-сложения. Опять же, внутренний цикл - это три конвейерных FMA, которые (опять же) обычно полезны, даже если вы не вычисляете квадратные корни.

Answer 3

Это также верно для разделения. MULSS(a,RCPSS(b)) намного быстрее, чем DIVSS(a,b). На самом деле это все еще быстрее, даже когда вы увеличиваете его точность с итерацией Ньютона-Рафсона.

Intel и AMD рекомендуют эту технику в своих руководствах по оптимизации. В приложениях, которые не требуют соответствия стандарту IEEE-754, единственной причиной для использования div/sqrt является удобочитаемость кода.

Answer 4

Вместо предоставления ответа, который на самом деле может быть неправильным (я также не собираюсь проверять или спорить о кеше и других вещах, скажем, они идентичны), я постараюсь указать вам источник, который может ответить на ваш вопрос.
Разница может заключаться в том, как вычисляются sqrt и rsqrt. Вы можете прочитать больше здесь http://www.intel.com/products/processor/manuals/. Я бы предложил начать с чтения о функциях процессора, которые вы используете, есть некоторая информация, особенно о rsqrt (cpu использует внутреннюю таблицу поиска с огромной аппроксимацией, что значительно упрощает получение результата). Может показаться, что rsqrt намного быстрее, чем sqrt, что одна дополнительная операция mul (не дорогостоящая) может не изменить ситуацию здесь.

Редактировать: Несколько фактов, о которых стоит упомянуть:
1. Однажды я выполнял некоторые микрооптимизации для моей графической библиотеки и использовал rsqrt для вычисления длины векторов. (вместо sqrt я умножил свою сумму в квадрате на rsqrt, что в точности соответствовало вашим тестам), и она работала лучше.
2. Вычисление rsqrt с использованием простой таблицы поиска может быть проще, как для rsqrt, когда x обращается в бесконечность, 1/sqrt(x) обращается в 0, поэтому для маленьких x значения функции не изменяются (много), тогда как для sqrt - уходит в бесконечность, так что это простой случай;).

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

Answer 5

Ньютон-Рафсон сходится к нулю f(x) используя приращения равны -f/f' где f' является производной.

За x=sqrt(y)можно попытаться решить f(x) = 0 за x с помощью f(x) = x^2 - y;

Тогда приращение: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / xкоторый имеет медленный разрыв в этом.

Вы можете попробовать другие функции (например, f(x) = 1/y - 1/x^2) но они будут одинаково сложными.

Давайте посмотрим на 1/sqrt(y) сейчас. Ты можешь попробовать f(x) = x^2 - 1/y, но это будет одинаково сложно: dx = 2xy / (y*x^2 - 1) например. Один неочевидный альтернативный выбор для f(x) является: f(x) = y - 1/x^2

Затем: dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

Ах! Это не тривиальное выражение, но у вас есть только умножения, без деления. => Быстрее!

И: полный шаг обновления new_x = x + dx затем читает:

x *= 3/2 - y/2 * x * x что тоже легко.

Answer 6

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