Apa Teknik Perkiraan yang Ada untuk Komputasi Akar Kuadrat?

Saya memiliki sumber daya yang sangat terbatas karena saya bekerja dengan mikrokontroler. Apakah ada ekspansi seri taylor, tabel pencarian umum, atau pendekatan rekursif?

Saya lebih suka melakukan sesuatu tanpa menggunakan math.h's sqrt ()

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

jika Anda ingin ekspansi seri daya murah dan kotor dioptimalkan (koefisien untuk seri Taylor bertemu perlahan) untuk sqrt()dan sekelompok trancendentals lainnya, saya memiliki beberapa kode dari dulu. saya dulu menjual kode ini, tetapi tidak ada yang membayar saya untuk itu selama hampir satu dekade. jadi saya pikir saya akan merilisnya untuk konsumsi publik. file khusus ini adalah untuk aplikasi di mana prosesor memiliki titik apung (presisi tunggal IEEE-754) dan mereka memiliki kompiler C dan sistem dev, tetapi mereka tidakmemiliki (atau mereka tidak ingin menautkan) stdlib yang akan memiliki fungsi matematika standar. mereka tidak membutuhkan ketelitian yang sempurna, tetapi mereka ingin segala sesuatunya menjadi cepat. Anda dapat dengan mudah merekayasa balik kode untuk melihat apa koefisien seri daya dan menulis kode Anda sendiri. kode ini mengasumsikan IEEE-754 dan menutup bit untuk mantissa dan eksponen.

tampaknya "markup kode" yang dimiliki SE tidak ramah dengan karakter sudut (Anda tahu ">" atau "<"), jadi Anda mungkin harus menekan "edit" untuk melihat semuanya.

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

Jika Anda belum melihatnya, "Quake square root" hanyalah membingungkan. Ia menggunakan beberapa bit-level magic untuk memberi Anda perkiraan pertama yang sangat baik, dan kemudian menggunakan satu atau dua putaran pendekatan Newton untuk merevisinya. Mungkin membantu Anda jika Anda bekerja dengan sumber daya terbatas.

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

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

Anda juga bisa memperkirakan fungsi akar kuadrat dengan menggunakan Metode Newton . Metode Newton adalah cara perkiraan di mana akar suatu fungsi berada. Ini juga merupakan metode iteratif di mana hasil dari iterasi sebelumnya digunakan dalam iterasi berikutnya hingga konvergensi. Persamaan untuk metode Newton untuk menebak di mana root adalah fungsi diberi tebakan awal x 0 didefinisikan sebagai:$f(x)$$x_0$

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

adalah tebakan pertama di mana root berada. Kami terus mendaur ulang persamaan dan menggunakan hasil dari iterasi sebelumnya hingga jawabannya tidak berubah. Secara umum, untuk menentukan menebak akar di ( n + 1 ) iterasi, mengingat menebak di n iterasi didefinisikan sebagai:$x_1$$(n+1)$$n$

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

Untuk menggunakan metode Newton untuk memperkirakan akar kuadrat, anggaplah kita diberi nomor . Karena itu, untuk menghitung akar kuadrat, kita perlu menghitung $a$ Karena itu, kami berusaha menemukan jawaban sedemikian rupa sehinggax=$\sqrt{a}$ . Kuadratkan kedua sisi, dan memindahkanake sisi lain dari persamaan menghasilkanx2-a=0. Dengan demikian, jawaban untuk persamaan ini adalah$x = \sqrt{a}$$a$$x^2 - a = 0$ dan dengan demikian merupakanakardari fungsi tersebut. Karena itu, misalkanf(x)=x2-aadalah persamaan yang ingin kita temukan akar. Dengan mengganti ini ke dalam metode Newton,f′(x)=2x, dan oleh karena itu:$\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)$$

$a$$\frac{1}{x} = \sqrt{a}$$\sqrt{a}$$\frac{1}{x^2} - a = 0$$\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)$$

Namun, ada peringatan yang harus kita pertimbangkan ketika melihat persamaan di atas. Untuk akar kuadrat, solusinya harus positif dan agar iterasi (dan hasilnya) menjadi positif, kondisi berikut harus dipenuhi:

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

Karena itu:

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

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

Saat tag Anda mencari algoritme C, mari kita tuliskan dengan sangat cepat:

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

Ini adalah implementasi yang cukup mendasar dari metode Newton. Perhatikan bahwa saya terus mengurangi setengah tebakan awal sampai kondisi yang kita bicarakan sebelumnya puas. Saya juga mencoba mencari akar kuadrat dari 5. Kita tahu bahwa ini kira-kira sama dengan 2,236 atau lebih. Menggunakan kode di atas memberikan hasil sebagai berikut:

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!

Seperti yang Anda lihat, satu-satunya hal yang berbeda adalah berapa banyak iterasi yang diperlukan untuk menghitung akar kuadrat. Semakin tinggi jumlah yang ingin Anda hitung, semakin banyak iterasi yang akan diambil.

Saya tahu bahwa metode ini telah disarankan di posting sebelumnya, tetapi saya pikir saya akan mendapatkan metode ini dan juga menyediakan beberapa kode!

$\sqrt{x}$

ya, rangkaian daya dapat dengan cepat dan efisien memperkirakan fungsi akar kuadrat, dan hanya melalui domain terbatas. semakin luas domain, semakin banyak istilah yang Anda butuhkan dalam seri daya Anda untuk menjaga kesalahan cukup rendah.

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

dimana

$a_1$

$a_2$

$a_3$

$a_4$

$x=1$$x=2$

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

jika itu floating point, Anda perlu memisahkan eksponen dan mantissa seperti kode C saya lakukan di jawaban yang lain.

Sebenarnya ini dilakukan dengan menyelesaikan persamaan kuadrat menggunakan Metode Newton:

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

Untuk angka yang lebih besar dari satu, Anda dapat menggunakan ekspansi Taylor berikut:

http://planetmath.org/taylorexpansionofsqrt1x

$a>b$

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

Dalam presisi 4%, jika saya ingat dengan baik. Itu digunakan oleh insinyur, sebelum penguasa dan kalkulator logaritmik. Saya mempelajarinya dalam Notes et formules de l'ingénieur, De Laharpe , 1923