What Approximation Techniques Exist for Computing the Square Root?

Question

I have very limited resources as I'm working with a microcontroller. Is there a taylor-series expansion, common lookup table, or recursive approach?

I'd prefer to do something without using math.h's sqrt()

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

Answer 1

if you want a cheap and dirty optimized power-series expansion (the coefficients for Taylor series converge slowly) for sqrt() and a bunch of other trancendentals, i have some code from long ago. i used to sell this code, but no one has paid me for it for nearly a decade. so i think i'll release it for public consumption. this particular file was for an application where the processor had floating point (IEEE-754 single precision) and they had a C compiler and dev system, but they did not have (or they didn't want to link in) the stdlib that would have had the standard math functions. they did not need perfect precision, but they wanted things to be fast. you can pretty easily reverse engineer the code to see what the power series coefficients are and write your own code. this code assumes IEEE-754 and masked off the bits for mantissa and exponent.

file: __functions.h

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

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

The math that sorta explains this

Answer 2

If you haven't seen it, the "Quake square root" is simply mystifying. It uses some bit-level magic to give you a very good first approximation, and then uses a round or two of Newton's approximation to revise. It might help you if you're working with limited resources.

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

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

Answer 3

You could also approximate the square root function by using Newton's Method. Newton's Method is a way of approximating where the roots of a function are. It is also an iterative method where the result from the previous iteration is used in the next iteration until convergence. The equation for Newton's method to guess where the root is of a function $f(x)$ given an initial guess $x_0$ is defined as:

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

$x_1$ is the first guess of where the root is located. We keep recycling the equation and using results from previous iterations until the answer doesn't change. In general, to determine the guess of the root at the $(n+1)$ iteration, given the guess at the $n$ iteration is defined as:

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

To use Newton's method to approximate for the square root, suppose that we are given a number $a$. As such, to compute the square root, we need to compute $\sqrt{a}$ As such, we seek to find an answer such that $x = \sqrt{a}$. Square both sides, and moving $a$ to the other side of the equation yields $x^2 - a = 0$. As such, the answer to this equation is $\sqrt{a}$ and is thus the root of the function. As such, let $f(x) = x^2 - a$ be the equation we want to find the root of. By substituting this into Newton's method, $f'(x) = 2x$, and therefore:

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

Therefore, to compute the square root of $a$, we simply need to compute Newton's method until we converge. However, as noted by @robertbristow-johnson, division is a very expensive operation - especially for microcontrollers / DSPs with limited resources. In addition, there may be a case where a guess may be 0, which would result in a divide by 0 error due to the division operation. As such, what we can do is use Newton's method and solve for the reciprocal function instead, i.e. $\frac{1}{x} = \sqrt{a}$. This also avoids any division, as we will see later. Squaring both sides, and moving $\sqrt{a}$ over to the left hand side thus gives $\frac{1}{x^2} - a = 0$. Therefore, the solution to this would be $\frac{1}{\sqrt{a}}$. By multiplying by $a$, we would get our intended result. Again, using Newton's method, we thus have:

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

However, there is a warning that we should consider when looking at the above equation. For square roots, the solution should be positive and so in order for the iterations (and the result) to be positive, the following condition must be satisfied:

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

Therefore:

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

Therefore, given the number we wish to compute the square root of, the initial guess must satisfy the above condition. As this is going to ultimately be placed on a microcontroller, we could start with any value of $x_0$ (say 1), then keep looping and decreasing the value of $x_0$ until the above condition is satisfied. Note that I avoided doing division to directly compute what the value of $x_0$ should be as division is an expensive operation. Once we have our initial guess, iterate through Newton's method. Note that depending on the initial guess, it may take a shorter or longer amount of time to converge. It all depends on how close you are to the actual answer. You can either cap the number of iterations, or wait until the relative difference between the two roots is less than some threshold (like $10^{-6}$ or so).

As your tag is looking for an algorithm in C, let's write one very quickly:

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

This is a pretty basic implementation of Newton's method. Note that I keep decreasing the initial guess by half until the condition we talked about earlier is satisfied. I am also trying to find the square root of 5. We know that this is roughly equal to 2.236 or so. Using the above code gives the following output:

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!

Note that Newton's method is finding the solution of the reciprocal solution and we multiply by $a$ at the end to get our final answer. Also, take note that the initial guess was changed to ensure that the criteria I talked about above is satisfied. Just for fun, let's try finding the square root of 9876.

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!

As you can see, the only thing that is different is how many iterations are required to compute the square root. The higher the number of what you want to compute, the more iterations it will take.

I know that this method has already been suggested in an earlier post, but I figured I'd derive the method as well as provide some code!

Answer 4

because the code markup for SE seems to work like shit, i'll try to answer this more directly, specifically for the $\sqrt{x}$ function.

yes, a power series can quickly and efficiently approximate the square root function, and only over a limited domain. the wider the domain, the more terms you will need in your power series to keep the error sufficiently low.

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

where

$a_1$ = 0.49959804148061

$a_2$ = -0.12047308243453

$a_3$ = 0.04585425015501

$a_4$ = -0.01076564682800

these coefficients were determined using a modified Remez exchange algorithm such that equality is at $x=1$ and $x=2$ and the maximum relative error in between is minimized.

so, if your implementation is fixed point, then you need to shift your bits (counting the number of bit positions shifted) until your values is between 1 and 2 using the scaling of your fixed-point scheme. if you shift right[left] the argument by $2n$ bits to get the argument to be between 1 and 2, the result should be shift left[right] by $n$ bits. if the number of shift bits is odd, the extra bit shift is compensated with an additional multiplication of $\sqrt{2}$, which can be stored as a constant in your code.

if it's floating point, you need to separate the exponent and mantissa like my C code does in the other answer.

Answer 5

One square root expansion that puzzled me in the past is the one for complex magnitude (or diagonal in a rectangle); if $a>b$, then:

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

Within 4% precision, if I remember well. It was used by engineers, before logarithmic rulers and calculators. I learned it in Notes et formules de l'ingénieur, De Laharpe, 1923

Answer 6

Actually it is done by solving an quadratic equation using Newton Method - Wikipedia: Methods of Computing Square Roots.

For numbers greater than one you can use the following Taylor Expansion - Planet Math - Taylor Expansion for $\sqrt{x + 1}$.