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?

phkahler
  • 5,687
  • 1
  • 23
  • 31
  • One shot of Halley's method (http://www.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. – Alexandre C. Apr 12 '12 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(!?!)? – BlueRaja - Danny Pflughoeft Apr 12 '12 at 20:50
  • Look at this `FastSqrt` function used in gaming http://www.gamedev.net/topic/278840-fast-sqrt/ – John Alexiou Apr 13 '12 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. – phkahler Apr 16 '12 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. – phkahler Apr 26 '12 at 19:19

4 Answers4

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.

Colin D
  • 5,641
  • 1
  • 23
  • 35
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.

deed02392
  • 4,799
  • 2
  • 31
  • 48
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().

John Alexiou
  • 28,472
  • 11
  • 77
  • 133
  • 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? – phkahler Apr 12 '12 at 15:11
  • @phkahler: Shortcut the zero to return `0` duh! if `y_c` is zero you have to evaluate the `sqrt()`. – John Alexiou Apr 12 '12 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)). – Alexandre C. Apr 12 '12 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. – phkahler Apr 12 '12 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. – phkahler Apr 12 '12 at 17:44
  • I think repeated use will diverge fast. Maybe a lookup table is a better solution. Oh well, I tried! – John Alexiou Apr 12 '12 at 18:14
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

Alexandre C.
  • 55,948
  • 11
  • 128
  • 197
  • 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. – phkahler Apr 27 '12 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. – R.. GitHub STOP HELPING ICE May 21 '12 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. – Alexandre C. May 21 '12 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. – R.. GitHub STOP HELPING ICE May 21 '12 at 20:56