Я профилировал некоторые из наших основных математических вычислений на 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 с использованием FPU x87 это было ужасно плохо:
inline float TestSqrtFunction( float in )
{ return sqrt(in); }
Следующее, что я пробовал, - это использовать встроенную функцию, чтобы заставить компилятор использовать скалярный код операции sqrt 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 10 (что слишком много для моих целей).
Неприятность была, когда я попробовал операцию SSE для получения обратного квадратного корня, а затем использовал умножение, чтобы получить квадратный корень (x * 1 / √x = √x). Даже если это занимает две зависимые операции, он был самым быстрым решением на сегодняшний день, в 1.24ns / поплавком и с точностью до 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
- но структура данных массива здесь предназначена только для целей тестирования: то, что я действительно пытаюсь измерить, - это скалярная производительность для использования в коде которые нельзя векторизовать.)
источник
inline float SSESqrt( float restrict fIn ) { float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; }
. Но это плохая идея, потому что это может легко вызвать остановку загрузки-попадания-сохранения, если ЦП записывает числа с плавающей запятой в стек, а затем сразу же считывает их обратно - в частности, перестановка из векторного регистра в регистр с плавающей точкой для возвращаемого значения плохие новости. Кроме того, коды операций базовой машины, которые встроенные функции SSE представляют, в любом случае принимают адресные операнды.eax
) очень плохое, в то время как круговой обход между xmm0 и стеком и обратно нет из-за переадресации магазина Intel. Вы можете рассчитать время сами, чтобы убедиться в этом. Как правило, самый простой способ увидеть потенциальную LHS - это посмотреть на выпущенную сборку и увидеть, где данные перебираются между наборами регистров; ваш компилятор может сделать умную вещь, а может и нет. Что касается нормализации векторов, я записал свои результаты здесь: bit.ly/9W5zoUОтветы:
sqrtss
дает правильно округленный результат.rsqrtss
дает приближение к обратной величине с точностью примерно до 11 бит.sqrtss
генерирует гораздо более точный результат, когда требуется точность.rsqrtss
существует для случаев, когда достаточно приближения, но требуется скорость. Если вы прочитаете документацию Intel, вы также найдете последовательность инструкций (аппроксимация обратного квадратного корня, за которой следует один шаг Ньютона-Рафсона), которая дает почти полную точность (~ 23 бита точности, если я правильно помню), и все же несколько быстрее чемsqrtss
.edit: Если скорость критична, и вы действительно вызываете это в цикле для многих значений, вы должны использовать векторизованные версии этих инструкций,
rsqrtps
илиsqrtps
оба из которых обрабатывают четыре числа с плавающей запятой на инструкцию.источник
sqrtss
это правильно закругленные , который требует ~ 50 бит до округления, и не может быть достигнуто с помощью простой N / R итерации в одинарной точности.Это верно и для деления. MULSS (a, RCPSS (b)) намного быстрее, чем DIVSS (a, b). Фактически, он все еще быстрее, даже если вы увеличите его точность с помощью итерации Ньютона-Рафсона.
Intel и AMD рекомендуют этот метод в своих руководствах по оптимизации. В приложениях, которые не требуют соответствия IEEE-754, единственной причиной использования div / sqrt является удобочитаемость кода.
источник
div
это не единственная операция, поэтому общая пропускная способность uop часто оказывается узким местом, даже если естьdivps
илиdivss
. См. Раздел Деление с плавающей запятой против умножения с плавающей запятой , где в моем ответе есть раздел о том, почемуrcpps
больше не выигрывает пропускная способность. (Или выигрыш в задержке) и цифры по разделению пропускной способности / задержки.a * rcpss(b)
может быть быстрее, но это все равно больше, чемa/b
!Вместо того, чтобы давать ответ, который на самом деле может быть неверным (я также не собираюсь проверять или спорить о кеше и других вещах, скажем, они идентичны), я попытаюсь указать вам источник, который может ответить на ваш вопрос.
Разница может заключаться в том, как вычисляются sqrt и rsqrt. Вы можете прочитать больше здесь http://www.intel.com/products/processor/manuals/ . Я предлагаю начать с чтения о функциях процессора, которые вы используете, есть некоторая информация, особенно о rsqrt (процессор использует внутреннюю таблицу поиска с огромным приближением, что значительно упрощает получение результата). Может показаться, что rsqrt настолько быстрее, чем sqrt, что 1 дополнительная операция mul (что не слишком дорого) может не изменить ситуацию.
Изменить: несколько фактов, которые стоит упомянуть:
1. Однажды я делал некоторые микрооптимизации для своей графической библиотеки и использовал rsqrt для вычисления длины векторов. (вместо sqrt я умножил свою сумму квадратов на rsqrt, что и было сделано в ваших тестах), и он работал лучше.
2. Вычисление rsqrt с использованием простой таблицы поиска может быть проще, как и для rsqrt, когда x переходит в бесконечность, 1 / sqrt (x) переходит в 0, поэтому для малых x значения функции не меняются (сильно), тогда как для sqrt - он уходит в бесконечность, так что это тот простой случай;).
Кроме того, пояснение: я не уверен, где я нашел его в книгах, которые я связал, но я почти уверен, что читал, что rsqrt использует некоторую таблицу поиска, и ее следует использовать только тогда, когда результат не обязательно быть точным, хотя - я тоже могу ошибаться, как это было некоторое время назад :).
источник
Ньютон-Рафсон сходится к нулю,
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
что тоже легко.источник
Уже несколько лет назад на этот вопрос есть ряд других ответов. Вот что было правильным при консенсусе:
Вот что ошиблось в консенсусе:
Алгоритм NR для вычисления обратного квадратного корня имеет этот шаг обновления, как отмечали другие:
Это много умножений, зависящих от данных, и одного вычитания.
Далее следует алгоритм, который фактически используют современные FPU.
Дано
b[0] = n
, предположим, что мы можем найти ряд чиселY[i]
,b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2
приближающихся к 1. Затем рассмотрим:Ясно
x[n]
подходитsqrt(n)
иy[n]
подходит1/sqrt(n)
.Мы можем использовать шаг обновления Ньютона-Рафсона для получения обратного квадратного корня, чтобы получить хорошее
Y[i]
:Затем:
и:
Следующее ключевое наблюдение таково
b[i] = x[i-1] * y[i-1]
. Так:Затем:
То есть, учитывая начальные x и y, мы можем использовать следующий шаг обновления:
Или, что еще интереснее, мы можем установить
h = 0.5 * y
. Это инициализация:И это шаг обновления:
Это алгоритм Гольдшмидта, и он имеет огромное преимущество, если вы реализуете его на аппаратном уровне: «внутренний цикл» состоит из трех операций умножения и сложения и ничего больше, а два из них независимы и могут быть конвейеризованы.
В 1999 году FPU уже нуждались в конвейерной схеме сложения / вычитания и конвейерной схеме умножения, иначе SSE не был бы очень «потоковым». В 1999 году потребовалась только одна из каждой схемы, чтобы реализовать этот внутренний цикл полностью конвейерным способом, не тратя много оборудования только на извлечение квадратного корня.
Сегодня, конечно, мы представили программисту слияние умножения-сложения. Опять же, внутренний цикл - это три конвейерных FMA, которые (опять же) обычно полезны, даже если вы не вычисляете квадратные корни.
источник
_mm256_rsqrt_ps
анализом производительности Haswell. Обычно это хорошая идея, если у вас нет другой работы в цикле, и это сильно затруднит пропускную способность делителя. HW sqrt - это один uop, поэтому его можно смешивать с другой работой.Это быстрее, потому что эти инструкции игнорируют режимы округления и не обрабатывают исключения с плавающей запятой или ненормализованные числа. По этим причинам намного проще конвейеризовать, спекулировать и выполнять другие инструкции fp вне очереди.
источник
rsqrt
«s намного ниже точность, что означает гораздо меньше работы (или вообще?) После табличной получить начальное предположение.