É possível rolar uma forma consideravelmente mais rápida versão do sqrt

Question

Em um aplicativo estou perfil, achei que em algumas situações esta função é capaz de demorar mais de 10% do tempo total de execução.

Eu vi o debate sobre os anos da mais rápido sqrt implementações usando sneaky de ponto flutuante de artifícios, mas não sei se essas coisas estão desatualizados em CPUs modernas.

MSVC++ 2008 compilador está sendo usado, como referência...apesar de eu assumir sqrt não vai adicionar muita sobrecarga embora.

Ver também aqui para a discussão similar sobre modf função.

EDITAR:para referência, este é amplamente utilizado o método, mas é, na verdade, muito mais rápido?Quantos ciclos é SQRT de qualquer maneira esses dias?

Answer 1

Sim, é possível, mesmo sem artifícios:

1) o sacrifício de precisão para a velocidade:o sqrt algoritmo é iterativo, re-implementar com menos iterações.

2) tabelas de pesquisa:apenas para o ponto de início da iteração, ou combinado com interpolação para você todo o caminho.

3) cache:você está sempre sqrting o mesmo conjunto limitado de valores?se assim for, o cache pode funcionar bem.Eu achei que este útil em aplicações gráficas, onde a mesma coisa está sendo calculado para muitas formas do mesmo tamanho, por isso os resultados podem ser utilmente em cache.

Answer 2

Há uma ótima tabela de comparação aqui:http://assemblyrequired.crashworks.org/timing-square-root/

Para encurtar a história, o SSQRTS do SSE2 é cerca de 2x mais rápido que o FPU FSQRT, e uma iteração de aproximação + é cerca de 4x mais rápida que isso (8x em geral).

Além disso, se você está tentando tomar um SQRT de precisão única, verifique se é realmente isso que você está recebendo. Ouvi falar de pelo menos um compilador que converteria o argumento do float em um SQRT duplo, chamado de precisão dupla e depois converterá de volta para flutuar.

Answer 3

É muito provável que você obtenha mais melhorias de velocidade alterando seu Algoritmos do que mudar seu implementações: Tente ligar sqrt() menos em vez de fazer chamadas mais rapidamente. (E se você acha que isso não é possível - as melhorias para sqrt() Você menciona é exatamente isso: melhorias do algoritmo usado para calcular uma raiz quadrada.)

Como é usado com muita frequência, é provável que a implementação de sua biblioteca padrão de sqrt() é quase ideal para o caso geral. A menos que você tenha um domínio restrito (por exemplo, se precisar de menos precisão), onde o algoritmo pode levar alguns atalhos, é muito improvável que alguém apresente uma implementação mais rápida.

Observe que, uma vez que essa função usa 10% do seu tempo de execução, mesmo que você consiga encontrar uma implementação que leva apenas 75% do tempo de std::sqrt(), isso ainda só trará seu tempo de execução para baixo 2,5%. Para a maioria dos aplicativos, os usuários nem perceberiam isso, exceto se usarem um relógio para medir.

Answer 4

Quão preciso você precisa de sua sqrt para ser?Você pode obter aproximações razoáveis muito rapidamente:ver Quake3 excelente inverso da raiz quadrada função para a inspiração (note que o código está sob a GPL, de modo que você pode não querer integrá-lo diretamente).

Answer 5

Não sei se você consertou isso, mas já li sobre isso antes, e parece que a coisa mais rápida a fazer é substituir o sqrt função com uma versão de montagem embutida;

Você pode ver uma descrição de uma carga de alternativas aqui.

O melhor é este trecho de magia:

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

É cerca de 4,7x mais rápido que o padrão sqrt Ligue com a mesma precisão.

Answer 6

Aqui está uma maneira rápida, com uma mesa de pesquisa de apenas 8kb. O erro é ~ 0,5% do resultado. Você pode aumentar facilmente a tabela, reduzindo assim o erro. Funciona cerca de 5 vezes mais rápido que o sqrt () regular

// 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;
}