Можно ли накатить значительно более быструю версию sqrt

Question

В приложении, которое я профилирую, я обнаружил, что в некоторых сценариях эта функция может занимать более 10% общего времени выполнения.

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

Компилятор MSVC ++ 2008 используется, для справки ... хотя я бы предположил, что sqrt не собирается добавлять много накладных расходов.

Смотрите также здесь для аналогичного обсуждения функции modf .

РЕДАКТИРОВАТЬ: для справки, этот является одним из широко используемых методов, но на самом ли деле он намного быстрее? Сколько циклов SQRT в любом случае в эти дни?

Answer 1

Да, это возможно даже без обмана:

1) жертвовать точностью ради скорости: алгоритм sqrt является итеративным, перестраивать с меньшим количеством итераций.

2) справочные таблицы: либо просто для начальной точки итерации, либо в сочетании с интерполяцией, чтобы вы могли пройти весь путь туда.

3) кеширование: всегда ли вы используете один и тот же ограниченный набор значений? если это так, кэширование может работать хорошо. Я нашел это полезным в графических приложениях, где для множества фигур одинакового размера вычисляется одна и та же вещь, поэтому результаты могут быть эффективно кэшированы.

Answer 2

Здесь есть отличная таблица сравнения: http://assemblyrequired.crashworks.org/timing-square-root/

Короче говоря, ssqrts SSE2 примерно в 2 раза быстрее, чем fsqrt FPU, а итерация приближения + примерно в 4 раза быстрее (в 8 раз больше).

Кроме того, если вы пытаетесь получить sqrt одинарной точности, убедитесь, что это именно то, что вы получаете. Я слышал, по крайней мере, об одном компиляторе, который преобразует аргумент float в double, вызывает sqrt двойной точности, а затем конвертирует обратно в float.

Answer 3

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

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

Обратите внимание, что, поскольку эта функция использует 10% вашего времени выполнения, даже если вам удастся придумать реализацию, которая занимает всего 75% времени std::sqrt(), это все равно приведет только к сокращению времени выполнения на 2,5% * * 1016. Для большинства приложений пользователи даже не заметят этого, если не использовать часы для измерения.

Answer 4

Насколько точным вам должен быть ваш sqrt?Вы можете получить разумные приближения очень быстро: см. Превосходную функцию обратного квадратного корня в Quake3 для вдохновения (обратите внимание, что код написан под GPL, поэтому вы можете не захотеть интегрировать его напрямую).

Answer 5

Вот быстрый способ с таблицей поиска только 8 КБ. Ошибка ~ 0,5% от результата. Вы можете легко увеличить стол, тем самым уменьшив ошибку. Работает примерно в 5 раз быстрее, чем обычный sqrt ()

// LUT for fast sqrt of floats. Table will be consist of 2 parts, half for sqrt(X) and half for sqrt(2X).
const int nBitsForSQRTprecision = 11;                       // Use only 11 most sagnificant bits from the 23 of float. We can use 15 bits instead. It will produce less error but take more place in a memory. 
const int nUnusedBits   = 23 - nBitsForSQRTprecision;       // Amount of bits we will disregard
const int tableSize     = (1 << (nBitsForSQRTprecision+1)); // 2^nBits*2 because we have 2 halves of the table.
static short sqrtTab[tableSize]; 
static unsigned char is_sqrttab_initialized = FALSE;        // Once initialized will be true

// Table of precalculated sqrt() for future fast calculation. Approximates the exact with an error of about 0.5%
// Note: To access the bits of a float in C quickly we must misuse pointers.
// More info in: http://en.wikipedia.org/wiki/Single_precision
void build_fsqrt_table(void){
    unsigned short i;
    float f;
    UINT32 *fi = (UINT32*)&f;

    if (is_sqrttab_initialized)
        return;

    const int halfTableSize = (tableSize>>1);
    for (i=0; i < halfTableSize; i++){
         *fi = 0;
         *fi = (i << nUnusedBits) | (127 << 23); // Build a float with the bit pattern i as mantissa, and an exponent of 0, stored as 127

         // Take the square root then strip the first 'nBitsForSQRTprecision' bits of the mantissa into the table
         f = sqrtf(f);
         sqrtTab[i] = (short)((*fi & 0x7fffff) >> nUnusedBits);

         // Repeat the process, this time with an exponent of 1, stored as 128
         *fi = 0;
         *fi = (i << nUnusedBits) | (128 << 23);
         f = sqrtf(f);
         sqrtTab[i+halfTableSize] = (short)((*fi & 0x7fffff) >> nUnusedBits);
    }
    is_sqrttab_initialized = TRUE;
}

// Calculation of a square root. Divide the exponent of float by 2 and sqrt() its mantissa using the precalculated table.
float fast_float_sqrt(float n){
    if (n <= 0.f) 
        return 0.f;                           // On 0 or negative return 0.
    UINT32 *num = (UINT32*)&n;
    short e;                                  // Exponent
    e = (*num >> 23) - 127;                   // In 'float' the exponent is stored with 127 added.
    *num &= 0x7fffff;                         // leave only the mantissa 

    // If the exponent is odd so we have to look it up in the second half of the lookup table, so we set the high bit.
    const int halfTableSize = (tableSize>>1);
    const int secondHalphTableIdBit = halfTableSize << nUnusedBits;
    if (e & 0x01) 
        *num |= secondHalphTableIdBit;  
    e >>= 1;                                  // Divide the exponent by two (note that in C the shift operators are sign preserving for signed operands

    // Do the table lookup, based on the quaternary mantissa, then reconstruct the result back into a float
    *num = ((sqrtTab[*num >> nUnusedBits]) << nUnusedBits) | ((e + 127) << 23);
    return n;
}

Answer 6

Не знаю, исправили ли вы это, но я читал об этом раньше, и кажется, что самое быстрое, что нужно сделать, это заменить функцию sqrt версией встроенной сборки;

Вы можете увидеть описание множества альтернатив здесь .

Лучше всего этот фрагмент магии:

double inline __declspec (naked) __fastcall sqrt(double n)
{
    _asm fld qword ptr [esp+4]
    _asm fsqrt
    _asm ret 8
} 

Это примерно в 4,7 раза быстрее, чем стандартный вызов sqrt с той же точностью.