tracking square root of moving value

D

4

7

I have a control loop running at high frequency and need to compute a square root each cycle. Typical square root functions work fine but take excessive time. Since the value I'm taking the square root of doesn't change too much on each cycle, I would like to find an iterative square root that will converge and then track the correct result. This way I could do a single iteration at each time step, rather than many.

The problem is that all of the iterative square root methods I've seen will probably fail when the input is changing. In particular it looks like there will be problems when the input goes to zero and then increases again - the methods don't like to start with a guess of zero.

My input range is 0-4.5 and I need a precision of around 0.01 so using an increment/decrement of 0.01 could take far too long - I want it to mostly converge in 10 cycles or less.

FYI I'm using 16/32bit fixed point the input is 16bit q12. It's on a micro-controller so I'm not interested in using 1K for a lookup table. The code is also generated from a simulink model and their table lookup functions are rather full of overhead.

Is there a nice solution to this?

Desexualize answered 12/4, 2012 at 14:10 Comment(5)

One shot of Halley's method (mathpath.org/Algor/squareroot/algor.square.root.halley.htm) should do fine. If you want to avoid division, update 1 / sqrt(x) instead and use Newton or Halley. – Rudie 12/4, 2012 at 15:58

What do you mean, the value changes? Are you saying you want to find sqrt(x + epsilon) knowing x and sqrt(x) without having to calculate it directly? Or are you saying the register that contains x is volatile, and can change in the middle of the calculation(!?!)? – Austerity 12/4, 2012 at 20:50

Look at this FastSqrt function used in gaming gamedev.net/topic/278840-fast-sqrt – Fairway 13/4, 2012 at 19:31

@BlueRaja-DannyPflughoeft: It's a control system running at 10KHz. Each 100us I need to take the square root of a signal. That signal doesn't change much from one pass to the next. I can also tolerate a cycle or two of delay, so doing 1 iteration each pass would be great if it could be convergent. Most algorithms seem to require starting from scratch and then doing a lot of iterations just to get a single square root. So yes, the input may change from one iteration to the next, but not in the middle of a calculation. – Desexualize 16/4, 2012 at 13:52

@Alexandre C. :Please post Halley's method as an answer instead of a comment so I can accept it - it works well. – Desexualize 26/4, 2012 at 19:19

R

1

You can use one shot of Halley method. It has cubic convergence and therefore should be quite precise if the value moves slightly:

x_{n+1} = x_n * (x_n^2 + 3Q) / (3 x_n^2 + Q)

This converges cubcially to sqrt(Q).

Reference: http://www.mathpath.org/Algor/squareroot/algor.square.root.halley.htm

Rudie answered 26/4, 2012 at 20:15 Comment(4)

My simulations show this to work the best. It's also the only thing I've tried that works reasonably near zero. x_n needs to be clipped to a small value at least in the feedback to prevent division by zero. – Desexualize 27/4, 2012 at 14:30

If this is floating point, you can just halve the exponent of the input and use the resulting value as a first estimate to throw into this or any other iterative square root algorithm. – Testimony 21/5, 2012 at 19:50

@R.. Initial guess is the last updated value (see question). But yes, in a more generic setting this would work very well. – Rudie 21/5, 2012 at 20:11

I noticed OP was having trouble with zero not being a valid initial-guess for some algorithms, so I figured this approach might open up more algorithms. – Testimony 21/5, 2012 at 20:56

W

5

the range 0-4.5 is fairly small. With precision of 0.01, that's only 450 possible calculations. You could compute them all at compile time as constants and just do a look up during runtime.

Willywillynilly answered 12/4, 2012 at 14:15 Comment(1)

Have you tried posting this question or something similiar to the math stack exchange site? – Willywillynilly 12/4, 2012 at 15:22

I

2

I would suggest you use a lookup table, if you know in advanced the ranges you are dealing with. Generate an array or hash table (depending on the language you're working in) to the level of precision you need and refer to this when you need your roots.

Immensurable answered 12/4, 2012 at 14:18 Comment(0)

F

2

I tried a second order Taylor expansion on sqrt(x) and go the following result

if y=sqrt(x) and you know y_c = sqrt(x_c) already then:

t = x-3*x_c;
y = (12*x_c*x_c-t*t)/(8*y_c*y_c*y_c);

The larger x is the better the approximation. For the worst case with x_c=0.01 and x=0.02 the result comes out 0.1375 vs. the real result of sqrt(0.02)=0.1414 or a difference of 0.0039 which is under 0.01.

I tested the code with C# and saw a steady 33% speedup vs Math.Sqrt().

Fairway answered 12/4, 2012 at 14:37 Comment(6)

I'll look into this. The division is undesirable but may be OK. What does it do when the input is zero for a while? – Desexualize 12/4, 2012 at 15:11

@phkahler: Shortcut the zero to return 0 duh! if y_c is zero you have to evaluate the sqrt(). – Fairway 12/4, 2012 at 15:23

@phkahler: you'll need the division anyway if you're computing sqrt(x) (you can avoid it if you're computing 1 / sqrt(x)). – Rudie 12/4, 2012 at 15:55

@ja72 a conditional branch is several cycles on the processor I'm using - almost as bad as evaluating a function. Still better than sqrt. – Desexualize 12/4, 2012 at 16:38

In my experiments so far, this does not appear to be stable - for constant input of 1 it works. For 0.5 it seems to diverge. – Desexualize 12/4, 2012 at 17:44

I think repeated use will diverge fast. Maybe a lookup table is a better solution. Oh well, I tried! – Fairway 12/4, 2012 at 18:14

R

1