Welche Approximationstechniken gibt es für die Berechnung der Quadratwurzel?

Ich habe sehr begrenzte Ressourcen, da ich mit einem Mikrocontroller arbeite. Gibt es eine Erweiterung der Taylor-Serie, eine gemeinsame Nachschlagetabelle oder einen rekursiven Ansatz?

Ich würde es vorziehen, etwas zu tun, ohne math.hs sqrt () zu verwenden.

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

Wenn Sie eine billige und schmutzig optimierte Power-Series-Erweiterung (die Koeffizienten für Taylor-Serien konvergieren langsam) für sqrt()und eine Reihe anderer Trancendentale wollen, habe ich einen Code von vor langer Zeit. Früher habe ich diesen Code verkauft, aber seit fast einem Jahrzehnt hat mich niemand mehr dafür bezahlt. Also denke ich, ich werde es für den öffentlichen Konsum freigeben. Diese spezielle Datei war für eine Anwendung gedacht, bei der der Prozessor Gleitkomma (IEEE-754 mit einfacher Genauigkeit) und einen C-Compiler und ein Dev-System hatte, dies jedoch nichthaben (oder wollten nicht verlinken) die stdlib, die die Standard-Mathematikfunktionen gehabt hätte. Sie brauchten keine perfekte Präzision, aber sie wollten, dass die Dinge schnell waren. Sie können den Code ziemlich einfach rückentwickeln, um die Leistungsreihenkoeffizienten zu ermitteln und Ihren eigenen Code zu schreiben. Dieser Code setzt IEEE-754 voraus und maskiert die Bits für Mantisse und Exponent.

Es scheint, dass das "Code-Markup", das SE hat, mit Winkelzeichen unfreundlich ist (Sie kennen ">" oder "<"), daher müssen Sie wahrscheinlich auf "Bearbeiten" klicken, um alles zu sehen.

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

Wenn Sie es nicht gesehen haben, ist die "Quake Quadratwurzel" einfach mystifizierend. Es verwendet etwas Magie auf Bitebene, um Ihnen eine sehr gute erste Annäherung zu geben, und verwendet dann ein oder zwei Runden von Newtons Näherung, um sie zu überarbeiten. Es kann Ihnen helfen, wenn Sie mit begrenzten Ressourcen arbeiten.

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

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

Sie können die Quadratwurzelfunktion auch mithilfe der Newtonschen Methode approximieren . Die Newtonsche Methode ist eine Methode zur Annäherung an die Wurzeln einer Funktion. Es ist auch eine iterative Methode, bei der das Ergebnis der vorherigen Iteration in der nächsten Iteration bis zur Konvergenz verwendet wird. Die Gleichung für Newtons Methode zum Erraten, wo die Wurzel einer Funktion wenn eine anfängliche Schätzung x 0 gegeben ist, ist definiert als:$f(x)$$x_0$

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

ist die erste Vermutung, wo sich die Wurzel befindet. Wir recyceln die Gleichung weiter und verwenden die Ergebnisse früherer Iterationen, bis sich die Antwort nicht mehr ändert. Um die Schätzung der Wurzel bei der ( n + 1 ) -Iterationzu bestimmen, wird im Allgemeinendie Schätzung bei der n- Iteration wie folgt definiert:$x_1$$(n+1)$$n$

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

Nehmen wir an, wir erhalten eine Zahl um die Newtonsche Methode zur Approximation der Quadratwurzel zu verwenden . Um die Quadratwurzel zu berechnen, müssen wir √ berechnen$a$ Als solches suchen wir nach einer Antwort, bei derx=$\sqrt{a}$ . Quadrieren Sie beide Seiten, und wenn Sieaauf die andere Seite der Gleichung verschieben, erhalten Siex2-a=0. Daher lautet die Antwort auf diese Gleichung$x = \sqrt{a}$$a$$x^2 - a = 0$ und ist somit dieWurzelder Funktion. Als solches seif(x)=x2-adie Gleichung, deren Wurzel wir finden wollen. Durch Einsetzen in Newtons Methode istf'(x)=2xund daher:$\sqrt{a}$$f(x) = x^2 - a$$f'(x) = 2x$

xn+1=

$$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)$$

Um die Quadratwurzel von zu berechnen, müssen wir einfach die Newtonsche Methode berechnen, bis wir konvergieren. Wie von @ robertbristow-johnson festgestellt, ist die Aufteilung jedoch eine sehr teure Operation - insbesondere für Mikrocontroller / DSPs mit begrenzten Ressourcen. Zusätzlich kann es einen Fall geben, in dem eine Schätzung 0 sein kann, was aufgrund der Teilungsoperation zu einem Fehler beim Teilen durch 0 führen würde. Als solches können wir die Newtonsche Methode verwenden und stattdessen nach der reziproken Funktion lösen , dh $a$ . Dies vermeidet auch jegliche Teilung, wie wir später sehen werden. Quadrieren Sie beide Seiten und bewegen Sie$\frac{1}{x} = \sqrt{a}$ Over auf der linken Seite ergibt also$\sqrt{a}$. Daher wäre die Lösung hierfür$\frac{1}{x^2} - a = 0$ . Durch Multiplikation mitawürden wir unser beabsichtigtes Ergebnis erhalten. Mit der Newtonschen Methode haben wir also wieder:$\frac{1}{\sqrt{a}}$$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)$$

Es gibt jedoch eine Warnung, die wir bei der Betrachtung der obigen Gleichung berücksichtigen sollten. Für Quadratwurzeln sollte die Lösung positiv sein. Damit die Iterationen (und das Ergebnis) positiv sind, muss die folgende Bedingung erfüllt sein:

3 x n > ( x n ) 3

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

Deshalb:

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

$x_0$$x_0$$x_0$$10^{-6}$

Da Ihr Tag nach einem Algorithmus sucht C, schreiben wir sehr schnell einen:

#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!");
}

Dies ist eine ziemlich grundlegende Implementierung der Newtonschen Methode. Beachten Sie, dass ich die anfängliche Schätzung immer wieder um die Hälfte reduziere, bis die Bedingung, über die wir zuvor gesprochen haben, erfüllt ist. Ich versuche auch, die Quadratwurzel von 5 zu finden. Wir wissen, dass dies ungefähr 2,236 oder so entspricht. Die Verwendung des obigen Codes ergibt die folgende Ausgabe:

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!

$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!

Wie Sie sehen, unterscheidet sich nur, wie viele Iterationen erforderlich sind, um die Quadratwurzel zu berechnen. Je höher die Anzahl der zu berechnenden Daten ist, desto mehr Iterationen sind erforderlich.

Ich weiß, dass diese Methode bereits in einem früheren Beitrag vorgeschlagen wurde, aber ich dachte, ich würde die Methode ableiten und Code bereitstellen!

Da das Code-Markup für SE wie Scheiße zu funktionieren scheint, werde ich versuchen, dies direkter zu beantworten, speziell für das $\sqrt{x}$ Funktion.

Ja, eine Potenzreihe kann schnell und effizient die Quadratwurzelfunktion und nur über einen begrenzten Bereich approximieren . Je breiter die Domäne, desto mehr Begriffe benötigen Sie in Ihrer Potenzreihe, um den Fehler ausreichend gering zu halten.

$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}$$

wo

$a_1$

$a_2$

$a_3$

$a_4$

$x=1$$x=2$

$2n$$n$$\sqrt{2}$

Wenn es sich um Gleitkomma handelt, müssen Sie den Exponenten und die Mantisse trennen, wie es mein C-Code in der anderen Antwort tut.

Tatsächlich erfolgt dies durch Lösen einer quadratischen Gleichung mit der Newton-Methode:

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

Für Zahlen größer als eins können Sie die folgende Taylor-Erweiterung verwenden:

http://planetmath.org/taylorexpansionofsqrt1x

$a>b$

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

Innerhalb von 4% Genauigkeit, wenn ich mich gut erinnere. Es wurde von Ingenieuren vor logarithmischen Linealen und Taschenrechnern verwendet. Ich habe es in Notes et formules de l'ingénieur, De Laharpe , 1923 gelernt