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

Ответ 1

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

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

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

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

Ответ 2

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

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

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

Ответ 3

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

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

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

Ответ 4

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

Ответ 5

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

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

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

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

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

Ответ 6

Вот быстрый способ с поисковой таблицей всего 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;
}