Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

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/

share|improve this question
4  
Check out this link: codeproject.com/Articles/69941/… –  Matt L. Jul 9 at 17:29
1  
Except the fact that it is more programming question, why not to make it an answer Matt? –  jojek Jul 9 at 17:54
    
Floating-point or fixed-point input? For fixed-point an iterative method can be preferable, but I won't bother explaining it unless you really want it. –  Oscar Jul 11 at 13:44
    
@Oscar, I'd love to learn the fixed-point method as I try to not require the use of floats in my firmware :). –  tarabyte Jul 11 at 16:54

5 Answers 5

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.

it appears that the "code markup" that SE has is unfriendly with angle characters (you know ">" or "<"), so you'll probably have to hit "edit" to see all of it.

//
//    FILE: __functions.h
//
//    fast and approximate transcendental functions
//
//    copyright (c) 2004  Robert Bristow-Johnson
//
//    rbj@audioimagination.com
//


#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
//
//    rbj@audioimagination.com
//

#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;
    }
share|improve this answer
    
does anybody know how this code markup works with SE? if you hit "edit" you can see the code that i intended, but what we view here has many lines of code omitted, and not only on the end of the file. i am using the markup reference that we're directed to by the SE markup help. if anyone can figure it out, please edit the answer and tell us what you did. –  robert bristow-johnson Jul 9 at 20:35

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/

share|improve this answer

Actually it is done by solving an quadratic equation using Newton Method:

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

For numbers greater than one you can use the following Taylor expansion:

http://planetmath.org/taylorexpansionofsqrt1x

share|improve this answer

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 $

$$\sqrt{x} \ \approx \ 1 + a_1 (x-1) + a_2 (x-1)^2 + a_3 (x-1)^3 + a_4 (x-1)^4$$

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 it's floating point, you need to separate the exponent and mantissa like my C code does in the other answer.

share|improve this answer

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!

share|improve this answer
2  
ray, might i suggest that the function you aim for is $f(x)=\frac{1}{\sqrt{x}}$ instead. no division will be needed in the iteration and all you need to do is multiply the result with $x$ to get $\sqrt{x}$. that's why you see all this stuff about the reciprocal square root in the lit and in the real-world implementations. –  robert bristow-johnson Jul 12 at 6:15
    
@robertbristow-johnson - Good point! Every time I looked at the equation for Newton's method, I cringed when I saw the division operator. Not only would we have a divide by zero error, but multiplying by the reciprocal of a number is faster than dividing. I'll work on an edit that does this. Thanks for the suggestion! –  rayryeng Jul 12 at 6:20
2  
it's just that, for people coding DSPs and some other chips, that division is particularly expensive, whereas these chips can multiply numbers as fast as they can move the numbers. –  robert bristow-johnson Jul 12 at 6:23
    
@robertbristow-johnson - and another excellent point. I remember back when I worked with the Motorola 6811 that multiplication took a few cycles while division took several hundred. Wasn't pretty. –  rayryeng Jul 12 at 6:25
2  
ahh, the good ol' 68HC11. had some stuff from the 6809 (like a quick multiply) but more of a microcontroller. –  robert bristow-johnson Jul 12 at 6:29

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.