Быстрый обратный квадрат двойного в C/С++

Недавно я профилировал программу, в которой точка доступа определенно это

****** d = somevalue();
****** d2=d*d;
****** c = 1.0/d2 // HOT SPOT

Значение d2 не используется после того, как мне нужно только значение c. Некоторое время назад я читал о методе Carmack быстрого обратного квадратного корня, это, очевидно, не так, но мне интересно, могут ли подобные алгоритмы помочь мне вычислить 1/x ^ 2.

Мне нужна довольно точная точность, я проверил, что моя программа не дает правильных результатов с опцией gcc -ffast-math. (g++ - 4.5)

3 ответа

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

  • Вы уверены, что вам нужна точность ******? Вы можете пожертвовать точностью достаточно легко:

    ****** d = somevalue();
    float c = 1.0f / ((float) d * (float) d);

    В этом случае 1.0f является абсолютно обязательным, если вы используете 1.0, вместо этого вы получите точность ******.

  • Вы пытались включить "неряшливую" математику в своем компиляторе? В GCC вы можете использовать -ffast-math, для других компиляторов есть аналогичные варианты. Небрежная математика может быть более чем достаточно для вашего приложения. ( Изменить: Я не видел никакой разницы в результирующей сборке.)

  • Если вы используете GCC, считаете ли вы, что используете -mrecip? Существует функция "взаимной оценки", которая имеет только около 12 бит точности, но она намного быстрее. Вы можете использовать метод Ньютона-Рафсона для повышения точности результата. Опция -mrecip заставит компилятор автоматически генерировать взаимные оценки и шаги Newton-Raphson для вас, хотя вы всегда можете написать сборку самостоятельно, если хотите точно настроить компромисс производительности. (Newton-Raphson сходится очень быстро.) ( Изменить: Мне не удалось заставить GCC генерировать RCPSS. См. Ниже.)

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

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

Трюки, которые не работают

  • Метод Carmack: он устарел для современных процессоров, которые имеют взаимные коды операций оценки. Для обратных лучей лучшая версия, которую я видел, дает только один бит точности - ничего по сравнению с 12 бит RCPSS. Я думаю, что совпадение заключается в том, что трюк работает так хорошо для взаимных квадратных корней; совпадение, которое вряд ли повторится.

  • Перемещающие переменные. Что касается компилятора, разница между 1.0/(x*x) и ****** x2 = x*x; 1.0/x2 очень мала. Я был бы удивлен, если бы вы нашли компилятор, который генерирует другой код для двух версий с оптимизацией, включенной даже на самый низкий уровень.

  • Использование pow. Библиотечная функция pow - это полный монстр. Когда GCC -ffast-math выключен, вызов библиотеки довольно дорог. При включении GCC -ffast-math вы получите тот же код сборки для pow(x, -2), что и для 1.0/(x*x), поэтому нет никакой пользы.

Update

Ниже приведен пример аппроксимации Ньютона-Рафсона для обратного квадрата значения с плавающей запятой с двойной точностью.

static ****** invsq(****** x)
{
 ****** y;
 int i;
 __asm__ (
 "cvtpd2ps %1, %0\n\t"
 "rcpss %0, %0\n\t"
 "cvtps2pd %0, %0"
 : "=x"(y)
 : "x"(x));
 for (i = 0; i < RECIP_ITER; ++i)
 y *= 2 - x * y;
 return y * y;
}

К сожалению, при использовании тестов RECIP_ITER=1 на моем компьютере это немного медленнее (~ 5%), чем простая версия 1.0/(x*x). Это быстрее (2x так же быстро) с нулевыми итерациями, но тогда вы получаете только 12 бит точности. Я не знаю, достаточно ли вам 12 бит.

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

Например, вы сказали, что -ffast-math вызвало нежелательную потерю точности; это может указывать на проблему численной стабильности в используемом алгоритме. При правильном выборе алгоритма многие проблемы могут быть решены с помощью float вместо ******. (Конечно, вам может понадобиться больше 24 бит. Я не знаю.)

Я подозреваю, что метод RCPSS светит, если вы хотите вычислить несколько из них параллельно.


Да, конечно, вы можете попробовать что-то сделать. Позвольте мне просто дать вам некоторые общие идеи, вы можете заполнить детали.

Во-первых, давайте посмотрим, почему работает корень Carmack:

Мы пишем x = M & times; 2 E обычным способом. Напомним, что поплавок IEEE сохраняет смещение экспоненты смещением: если e обозначено поле экспоненты, мы имеем e = Bias + E & ge; 0. Перестраивая, мы получаем E = e & minus; Bias.

Теперь для обратного квадратного корня: x & minus; 1/2 = M -1/2 & times; 2 & минус; Е /2. Новое поле экспоненты:

        е '  = Смещение & минус; Е /2 = 3/2 Смещение & минус; е /2

С битборгом мы можем получить значение e/2 из e, сдвинув его, и 3/2 Bias - это просто константа.

Кроме того, мантисса М хранится как 1.0 x с x < 1, и мы можем аппроксимировать M -1/2 как 1 + x/2. Опять же, тот факт, что только х хранится в двоичном виде, означает, что мы делим на два путем простого смещения битов.

Теперь мы рассмотрим x & minus; 2: это равно M & minus; 2 & times; 2 & minus; 2 E и мы ищем поле экспоненты:

        е '  = Смещение & минус; 2 Е = 3 Смещение & минус; 2 е

Опять же, 3 Bias - это просто константа, и вы можете получить 2 e от e с помощью битвыбора. Что касается мантиссы, вы можете аппроксимировать (1 + x) & minus; 2 на 1   2 x, и поэтому проблема сводится к получению 2 x x x от x.

Обратите внимание, что игра Carmack magic с плавающей запятой фактически не вычисляет результат в правильном направлении: скорее, он дает замечательно точную оценку, которая используется в качестве отправной точки для традиционного итеративного вычисления. Но поскольку оценка настолько хороша, вам нужно всего лишь несколько раундов последующей итерации, чтобы получить приемлемый результат.


Для вашей текущей программы вы определили точку доступа - хорошо. В качестве альтернативы ускорению 1/d ^ 2 у вас есть возможность изменить программу, чтобы она не так часто вычисляла 1/d ^ 2. Вы можете вытащить его из внутреннего цикла? За сколько разных значений d вы вычисляете 1/d ^ 2? Не могли бы вы предварительно вычислить все значения, которые вам нужны, а затем искать результаты? Это немного громоздко для 1/d ^ 2, но если 1/d ^ 2 является частью некоторого большего фрагмента кода, возможно, стоит применить этот трюк к этому. Вы говорите, что если вы понижаете точность, вы не получите достаточно ответов. Есть ли способ перефразировать код, который может обеспечить лучшее поведение? Численный анализ достаточно тонкий, что, возможно, стоит попробовать несколько вещей и посмотреть, что произошло.

В идеале, конечно, вы найдете некоторую оптимизированную рутину, которая опирается на годы исследований - есть ли что-нибудь в лапаке или linpack, на которые вы могли бы ссылаться?

licensed under cc by-sa 3.0 with attribution.