Que técnicas de aproximação existem para calcular a raiz quadrada?

Eu tenho recursos muito limitados, pois estou trabalhando com um microcontrolador. Existe uma expansão da série Taylor, uma tabela de pesquisa comum ou uma abordagem recursiva?

Eu preferiria fazer algo sem usar o sqrt () do math.h

http://www.cplusplus.com/reference/cmath/sqrt/

se você quer uma expansão de séries de potência otimizada barata e suja (os coeficientes da série de Taylor convergem lentamente) para sqrt()um monte de outros trancendentais, eu tenho algum código de muito tempo atrás. Eu costumava vender esse código, mas ninguém me pagou por quase uma década. então eu acho que vou lançá-lo para consumo público. esse arquivo específico era para um aplicativo em que o processador tinha ponto flutuante (precisão única IEEE-754) e eles tinham um compilador C e um sistema dev, mas não tinhampossui (ou não queria vincular) o stdlib que teria as funções matemáticas padrão. eles não precisavam de precisão perfeita, mas queriam que as coisas fossem rápidas. você pode facilmente fazer engenharia reversa do código para ver quais são os coeficientes da série de potência e escrever seu próprio código. esse código assume o IEEE-754 e mascara os bits para mantissa e expoente.

parece que a "marcação de código" que o SE possui não é amigável com caracteres angulares (você sabe ">" ou "<"); portanto, você provavelmente precisará clicar em "editar" para ver tudo.

//
//    FILE: __functions.h
//
//    fast and approximate transcendental functions
//
//    copyright (c) 2004  Robert Bristow-Johnson
//
//    [email protected]
//


#ifndef __FUNCTIONS_H
#define __FUNCTIONS_H

#define TINY 1.0e-8
#define HUGE 1.0e8

#define PI              (3.1415926535897932384626433832795028841972)        /* pi */
#define ONE_OVER_PI     (0.3183098861837906661338147750939)
#define TWOPI           (6.2831853071795864769252867665590057683943)        /* 2*pi */
#define ONE_OVER_TWOPI  (0.15915494309189535682609381638)
#define PI_2            (1.5707963267948966192313216916397514420986)        /* pi/2 */
#define TWO_OVER_PI     (0.636619772367581332267629550188)
#define LN2             (0.6931471805599453094172321214581765680755)        /* ln(2) */
#define ONE_OVER_LN2    (1.44269504088896333066907387547)
#define LN10            (2.3025850929940456840179914546843642076011)        /* ln(10) */
#define ONE_OVER_LN10   (0.43429448190325177635683940025)
#define ROOT2           (1.4142135623730950488016887242096980785697)        /* sqrt(2) */
#define ONE_OVER_ROOT2  (0.707106781186547438494264988549)

#define DB_LOG2_ENERGY          (3.01029995663981154631945610163)           /* dB = DB_LOG2_ENERGY*__log2(energy) */
#define DB_LOG2_AMPL            (6.02059991327962309263891220326)           /* dB = DB_LOG2_AMPL*__log2(amplitude) */
#define ONE_OVER_DB_LOG2_AMPL   (0.16609640474436811218256075335)           /* amplitude = __exp2(ONE_OVER_DB_LOG2_AMPL*dB) */

#define LONG_OFFSET     4096L
#define FLOAT_OFFSET    4096.0



float   __sqrt(float x);

float   __log2(float x);
float   __exp2(float x);

float   __log(float x);
float   __exp(float x);

float   __pow(float x, float y);

float   __sin_pi(float x);
float   __cos_pi(float x);

float   __sin(float x);
float   __cos(float x);
float   __tan(float x);

float   __atan(float x);
float   __asin(float x);
float   __acos(float x);

float   __arg(float Imag, float Real);

float   __poly(float *a, int order, float x);
float   __map(float *f, float scaler, float x);
float   __discreteMap(float *f, float scaler, float x);

unsigned long __random();

#endif




//
//    FILE: __functions.c
//
//    fast and approximate transcendental functions
//
//    copyright (c) 2004  Robert Bristow-Johnson
//
//    [email protected]
//

#define STD_MATH_LIB 0

#include "__functions.h"

#if STD_MATH_LIB
#include "math.h"   // angle brackets don't work with SE markup
#endif




float   __sqrt(register float x)
    {
#if STD_MATH_LIB
    return (float) sqrt((double)x);
#else
    if (x > 5.877471754e-39)
        {
        register float accumulator, xPower;
        register long intPart;
        register union {float f; long i;} xBits;

        xBits.f = x;

        intPart = ((xBits.i)>>23);                  /* get biased exponent */
        intPart -= 127;                             /* unbias it */

        x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
        x *= 1.192092895507812e-07;                 /* divide by 0x800000 */

        accumulator =  1.0 + 0.49959804148061*x;
        xPower = x*x;
        accumulator += -0.12047308243453*xPower;
        xPower *= x;
        accumulator += 0.04585425015501*xPower;
        xPower *= x;
        accumulator += -0.01076564682800*xPower;

        if (intPart & 0x00000001)
            {
            accumulator *= ROOT2;                   /* an odd input exponent means an extra sqrt(2) in the output */
            }

        xBits.i = intPart >> 1;                     /* divide exponent by 2, lose LSB */
        xBits.i += 127;                             /* rebias exponent */
        xBits.i <<= 23;                             /* move biased exponent into exponent bits */

        return accumulator * xBits.f;
        }
     else
        {
        return 0.0;
        }
#endif
    }




float   __log2(register float x)
    {
#if STD_MATH_LIB
    return (float) (ONE_OVER_LN2*log((double)x));
#else
    if (x > 5.877471754e-39)
        {
        register float accumulator, xPower;
        register long intPart;

        register union {float f; long i;} xBits;

        xBits.f = x;

        intPart = ((xBits.i)>>23);                  /* get biased exponent */
        intPart -= 127;                             /* unbias it */

        x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
        x *= 1.192092895507812e-07;                 /* divide by 0x800000 */

        accumulator = 1.44254494359510*x;
        xPower = x*x;
        accumulator += -0.71814525675041*xPower;
        xPower *= x;
        accumulator += 0.45754919692582*xPower;
        xPower *= x;
        accumulator += -0.27790534462866*xPower;
        xPower *= x;
        accumulator += 0.12179791068782*xPower;
        xPower *= x;
        accumulator += -0.02584144982967*xPower;

        return accumulator + (float)intPart;
        }
     else
        {
        return -HUGE;
        }
#endif
    }


float   __exp2(register float x)
    {
#if STD_MATH_LIB
    return (float) exp(LN2*(double)x);
#else
    if (x >= -127.0)
        {
        register float accumulator, xPower;
        register union {float f; long i;} xBits;

        xBits.i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;       /* integer part */
        x -= (float)(xBits.i);                                  /* fractional part */

        accumulator = 1.0 + 0.69303212081966*x;
        xPower = x*x;
        accumulator += 0.24137976293709*xPower;
        xPower *= x;
        accumulator += 0.05203236900844*xPower;
        xPower *= x;
        accumulator += 0.01355574723481*xPower;

        xBits.i += 127;                                         /* bias integer part */
        xBits.i <<= 23;                                         /* move biased int part into exponent bits */

        return accumulator * xBits.f;
        }
     else
        {
        return 0.0;
        }
#endif
    }


float   __log(register float x)
    {
#if STD_MATH_LIB
    return (float) log((double)x);
#else
    return LN2*__log2(x);
#endif
    }

float   __exp(register float x)
    {
#if STD_MATH_LIB
    return (float) exp((double)x);
#else
    return __exp2(ONE_OVER_LN2*x);
#endif
    }

float   __pow(float x, float y)
    {
#if STD_MATH_LIB
    return (float) pow((double)x, (double)y);
#else
    return __exp2(y*__log2(x));
#endif
    }




float   __sin_pi(register float x)
    {
#if STD_MATH_LIB
    return (float) sin(PI*(double)x);
#else
    register float accumulator, xPower, xSquared;

    register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024)<<1;
    x -= (float)evenIntPart;

    xSquared = x*x;
    accumulator = 3.14159265358979*x;
    xPower = xSquared*x;
    accumulator += -5.16731953364340*xPower;
    xPower *= xSquared;
    accumulator += 2.54620566822659*xPower;
    xPower *= xSquared;
    accumulator += -0.586027023087261*xPower;
    xPower *= xSquared;
    accumulator += 0.06554823491427*xPower;

    return accumulator;
#endif
    }


float   __cos_pi(register float x)
    {
#if STD_MATH_LIB
    return (float) cos(PI*(double)x);
#else
    register float accumulator, xPower, xSquared;

    register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024)<<1;
    x -= (float)evenIntPart;

    xSquared = x*x;
    accumulator = 1.57079632679490*x;                       /* series for sin(PI/2*x) */
    xPower = xSquared*x;
    accumulator += -0.64596406188166*xPower;
    xPower *= xSquared;
    accumulator += 0.07969158490912*xPower;
    xPower *= xSquared;
    accumulator += -0.00467687997706*xPower;
    xPower *= xSquared;
    accumulator += 0.00015303015470*xPower;

    return 1.0 - 2.0*accumulator*accumulator;               /* cos(w) = 1 - 2*(sin(w/2))^2 */
#endif
    }


float   __sin(register float x)
    {
#if STD_MATH_LIB
    return (float) sin((double)x);
#else
    x *= ONE_OVER_PI;
    return __sin_pi(x);
#endif
    }

float   __cos(register float x)
    {
#if STD_MATH_LIB
    return (float) cos((double)x);
#else
    x *= ONE_OVER_PI;
    return __cos_pi(x);
#endif
    }

float   __tan(register float x)
    {
#if STD_MATH_LIB
    return (float) tan((double)x);
#else
    x *= ONE_OVER_PI;
    return __sin_pi(x)/__cos_pi(x);
#endif
    }




float   __atan(register float x)
    {
#if STD_MATH_LIB
    return (float) atan((double)x);
#else
    register float accumulator, xPower, xSquared, offset;

    offset = 0.0;

    if (x < -1.0)
        {
        offset = -PI_2;
        x = -1.0/x;
        }
     else if (x > 1.0)
        {
        offset = PI_2;
        x = -1.0/x;
        }
    xSquared = x*x;
    accumulator = 1.0;
    xPower = xSquared;
    accumulator += 0.33288950512027*xPower;
    xPower *= xSquared;
    accumulator += -0.08467922817644*xPower;
    xPower *= xSquared;
    accumulator += 0.03252232640125*xPower;
    xPower *= xSquared;
    accumulator += -0.00749305860992*xPower;

    return offset + x/accumulator;
#endif
    }


float   __asin(register float x)
    {
#if STD_MATH_LIB
    return (float) asin((double)x);
#else
    return __atan(x/__sqrt(1.0 - x*x));
#endif
    }

float   __acos(register float x)
    {
#if STD_MATH_LIB
    return (float) acos((double)x);
#else
    return __atan(__sqrt(1.0 - x*x)/x);
#endif
    }


float   __arg(float Imag, float Real)
    {
#if STD_MATH_LIB
    return (float) atan2((double)Imag, (double)Real);
#else
    register float accumulator, xPower, xSquared, offset, x;

    if (Imag > 0.0)
        {
        if (Imag <= -Real)
            {
            offset = PI;
            x = Imag/Real;
            }
         else if (Imag > Real)
            {
            offset = PI_2;
            x = -Real/Imag;
            }
         else
            {
            offset = 0.0;
            x = Imag/Real;
            }
        }
     else
        {
        if (Imag >= Real)
            {
            offset = -PI;
            x = Imag/Real;
            }
         else if (Imag < -Real)
            {
            offset = -PI_2;
            x = -Real/Imag;
            }
         else
            {
            offset = 0.0;
            x = Imag/Real;
            }
        }

    xSquared = x*x;
    accumulator = 1.0;
    xPower = xSquared;
    accumulator += 0.33288950512027*xPower;
    xPower *= xSquared;
    accumulator += -0.08467922817644*xPower;
    xPower *= xSquared;
    accumulator += 0.03252232640125*xPower;
    xPower *= xSquared;
    accumulator += -0.00749305860992*xPower;

    return offset + x/accumulator;
#endif
    }




float   __poly(float *a, int order, float x)
    {
    register float accumulator = 0.0, xPower;
    register int n;

    accumulator = a[0];
    xPower = x;
    for (n=1; n<=order; n++)
        {
        accumulator += a[n]*xPower;
        xPower *= x;
        }

    return accumulator;
    }


float   __map(float *f, float scaler, float x)
    {
    register long i;

    x *= scaler;

    i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;         /* round down without floor() */

    return f[i] + (f[i+1] - f[i])*(x - (float)i);       /* linear interpolate between points */
    }


float   __discreteMap(float *f, float scaler, float x)
    {
    register long i;

    x *= scaler;

    i = (long)(x + (FLOAT_OFFSET+0.5)) - LONG_OFFSET;   /* round to nearest */

    return f[i];
    }


unsigned long __random()
    {
    static unsigned long seed0 = 0x5B7A2775, seed1 = 0x80C7169F;

    seed0 += seed1;
    seed1 += seed0;

    return seed1;
    }

Se você ainda não viu, a "raiz quadrada do Quake" é simplesmente confusa. Ele usa alguma mágica no nível de bit para fornecer uma boa primeira aproximação e, em seguida, usa uma ou duas rodadas da aproximação de Newton para revisar. Pode ser útil se você estiver trabalhando com recursos limitados.

https://en.wikipedia.org/wiki/Fast_inverse_square_root

http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/

Você também pode aproximar a função de raiz quadrada usando o Método de Newton . O método de Newton é uma maneira de aproximar onde estão as raízes de uma função. É também um método iterativo em que o resultado da iteração anterior é usado na próxima iteração até a convergência. A equação do método de Newton para adivinhar onde a raiz tem uma função dada uma suposição inicial é definida como:$f(x)$$x_0$

$$x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}$$

$x_1$ é o primeiro palpite de onde a raiz está localizada. Continuamos reciclando a equação e usando resultados de iterações anteriores até a resposta não mudar. Em geral, para determinar a suposição da raiz na iteração , dada a suposição na iteração, é definido como:$(n+1)$$n$

$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$

Para usar o método de Newton para aproximar a raiz quadrada, suponha que recebamos um número . Como tal, para calcular a raiz quadrada, precisamos calcular Dessa forma, procuramos encontrar uma resposta tal que . Esquadre os dois lados e movendo para o outro lado da equação produz . Como tal, a resposta para esta equação é e, portanto, é a raiz da função. Como tal, seja a equação da qual queremos encontrar a raiz. Substituindo isso no método de Newton, e, portanto:$a$$\sqrt{a}$$x = \sqrt{a}$$a$$x^2 - a = 0$$\sqrt{a}$$f(x) = x^2 - a$$f'(x) = 2x$

$$x_{n+1} = x_n - \frac{x_{n}^2 - a}{2x_{n}}$$ $$x_{n+1} = \frac{1}{2}\left(x_n + \frac{a}{x_n}\right)$$

Portanto, para calcular a raiz quadrada de , precisamos simplesmente calcular o método de Newton até convergir. No entanto, como observado por @ robertbristow-johnson, a divisão é uma operação muito cara - especialmente para microcontroladores / DSPs com recursos limitados. Além disso, pode haver um caso em que um palpite pode ser 0, o que resultaria em um erro de divisão por 0 devido à operação de divisão. Como tal, o que podemos fazer é usar o método de Newton e resolver a função recíproca , por exemplo, . Isso também evita qualquer divisão, como veremos mais adiante. Esquadrar os dois lados e mover para o lado esquerdo resulta em . Portanto, a solução para isso seria$a$$\frac{1}{x} = \sqrt{a}$$\sqrt{a}$$\frac{1}{x^2} - a = 0$$\frac{1}{\sqrt{a}}$ . Ao multiplicar por , obteríamos o resultado pretendido. Novamente, usando o método de Newton, temos assim:$a$

$$x_{n+1} = x_{n} - \frac{f(x_n)}{f'(x_n)}$$ $$x_{n+1} = x_{n} - \frac{\frac{1}{(x_n)^2} - a}{\frac{-2}{(x_n)^3}}$$ $$x_{n+1} = \frac{1}{2}\left(3x_n - (x_n)^3a\right)$$

No entanto, há um aviso que devemos considerar ao examinar a equação acima. Para raízes quadradas, a solução deve ser positiva e, para que as iterações (e o resultado) sejam positivas, a seguinte condição deve ser atendida:

3 x n > ( x n ) 3 a ( x n ) 2 a <

$$3x_n - (x_n)^3a > 0$$ $$3x_n > (x_n)^3a$$ $$(x_n)^2a < 3$$

Portanto:

$$(x_0)^2a < 3$$

Portanto, dado o número do qual desejamos calcular a raiz quadrada, a suposição inicial deve satisfazer a condição acima. Como isso finalmente será colocado em um microcontrolador, poderíamos começar com qualquer valor de (por exemplo, 1), e continuar repetindo e diminuindo o valor de até que a condição acima seja satisfeita. Observe que evitei fazer divisão para calcular diretamente qual o valor dex 0 x 0 10 - $x_0$$x_0$$x_0$deve ser como divisão é uma operação cara. Quando tivermos o nosso palpite inicial, repita o método de Newton. Observe que, dependendo da estimativa inicial, pode levar um tempo menor ou maior para convergir. Tudo depende de quão perto você está da resposta real. Você pode limitar o número de iterações ou aguardar até que a diferença relativa entre as duas raízes seja menor que algum limite (como ou mais).$10^{-6}$

Como sua tag está procurando um algoritmo C, vamos escrever um muito rapidamente:

#include <stdio.h> // For printf
#include <math.h> // For fabs
void main() 
{
   float a = 5.0; // Number we want to take the square root of
   float x = 1.0; // Initial guess
   float xprev; // Root for previous iteration
   int count; // Counter for iterations

   // Find a better initial guess
   // Half at each step until condition is satisfied
   while (x*x*a >= 3.0)
       x *= 0.5;

   printf("Initial guess: %f\n", x);

   count = 1; 
   do { 
       xprev = x; // Save for previous iteration
       printf("Iteration #%d: %f\n", count++, x);                   
       x = 0.5*(3*xprev - (xprev*xprev*xprev)*a); // Find square root of the reciprocal
   } while (fabs(x - xprev) > 1e-6); 

   x *= a; // Actual answer - Multiply by a
   printf("Square root is: %f\n", x);
   printf("Done!");
}

Esta é uma implementação bastante básica do método de Newton. Observe que eu continuo diminuindo a estimativa inicial pela metade até que a condição da qual falamos anteriormente seja satisfeita. Também estou tentando encontrar a raiz quadrada de 5. Sabemos que isso é aproximadamente igual a 2.236 ou mais. O uso do código acima fornece a seguinte saída:

Initial guess: 0.500000
Iteration #1: 0.500000
Iteration #2: 0.437500
Iteration #3: 0.446899
Iteration #4: 0.447213
Square root is: 2.236068
Done!

Observe que o método de Newton é encontrar a solução da solução recíproca e multiplicamos por no final para obter nossa resposta final. Além disso, observe que o palpite inicial foi alterado para garantir que os critérios mencionados acima sejam atendidos. Apenas por diversão, vamos tentar encontrar a raiz quadrada de 9876.$a$

Initial guess: 0.015625
Iteration #1: 0.015625
Iteration #2: 0.004601
Iteration #3: 0.006420
Iteration #4: 0.008323
Iteration #5: 0.009638
Iteration #6: 0.010036
Iteration #7: 0.010062
Square root is: 99.378067
Done!

Como você pode ver, a única coisa diferente é quantas iterações são necessárias para calcular a raiz quadrada. Quanto maior o número do que você deseja calcular, mais iterações serão necessárias.

Eu sei que esse método já foi sugerido em uma postagem anterior, mas achei que derivaria o método e forneceria algum código!

$\sqrt{x}$

Sim, uma série de potências pode aproximar rápida e eficientemente a função de raiz quadrada e apenas em um domínio limitado. quanto maior o domínio, mais termos serão necessários em sua série de potências para manter o erro suficientemente baixo.

$1 \le x \le 2$

$$\begin{align} \sqrt{x} \ &\approx \ 1 + a_1 (x-1) + a_2 (x-1)^2 + a_3 (x-1)^3 + a_4 (x-1)^4 \\ \\ & = 1 + (x-1)\bigg( a_1 + (x-1)\Big( a_2 + (x-1)\big( a_3 + (x-1)a_4 \big) \Big) \bigg) \\ \end{align}$$

Onde

$a_1$

$a_2$

$a_3$

$a_4$

$x=1$$x=2$

portanto, se sua implementação for um ponto fixo, você precisará alterar seus bits (contando o número de posições de bits alteradas) até que seus valores estejam entre 1 e 2 usando a escala do seu esquema de ponto fixo. se você deslocar para a direita [esquerda] o argumento em bits para obter o argumento entre 1 e 2, o resultado deverá ser deslocado para a esquerda [direita] em bits. se o número de bits de deslocamento for ímpar, o deslocamento extra de bits será compensado com uma multiplicação adicional de , que pode ser armazenada como uma constante no seu código.n $2n$$n$$\sqrt{2}$

se for ponto flutuante, você precisará separar o expoente e a mantissa como meu código C faz na outra resposta.

Na verdade, isso é feito resolvendo uma equação quadrática usando o Método Newton:

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots

Para números maiores que um, você pode usar a seguinte expansão de Taylor:

http://planetmath.org/taylorexpansionofsqrt1x

$a>b$

$$\sqrt{a^2+b^2}\approx 0.96a+0.4 b.$$

Dentro de 4% de precisão, se bem me lembro. Foi usado por engenheiros, antes de réguas e calculadoras logarítmicas. Eu aprendi isso em Notes et formules de l'ingénieur, De Laharpe , 1923