Computationally efficient sqrt(x^2 + y^2 + z^2)

On 4/9/2014 2:23 PM, robert bristow-johnson wrote:
> On 4/9/14 5:13 PM, Tim Wescott wrote:
>> On Wed, 09 Apr 2014 16:48:58 -0400, robert bristow-johnson wrote:
>>
>>> On 4/9/14 4:29 PM, Rob Gaddi wrote:
>>>>
>>>> Is http://en.wikipedia.org/wiki/Fast_inverse_square_root at all handy?
>>>>
>>>>
>>>        sqrt(x) =  x * (1/sqrt(x))
>>
>> I was about to tell Rob "no, how could the _inverse_ possibly be useful".
>>
>> Thanks for heading me off!
>>
>
>
> when the SHArC first came out, they had an instruction that's a
> primitive iteration for 1/sqrt(x) and i wondered *why* an iterative for
> that instead of for sqrt(x).
>
> the fact that a simply multiplication turns the one into the other was
> pointed out early, but it wasn't until years later that the
> implementation of Newton's Method (to solve for the inverse of
> f(x) = 1/(x^2)) came out to be so clean.  then it all made sense to me.
>

The Power PC has a frsqrte (Floating Reciprocal Square Root Estimate) 
instruction that looks up the initial value of Newton's method.  The 
initial estimate is correct to one part in 32.  Each iteration 
thereafter is a single instruction and a single clock cycle (I think).

There's also an instruction for doing reciprocals by Newton's method.

Rob.

Rob Doyle <radioengr@gmail.com> wrote:

(snip)
> The Power PC has a frsqrte (Floating Reciprocal Square Root Estimate) 
> instruction that looks up the initial value of Newton's method.  The 
> initial estimate is correct to one part in 32.  Each iteration 
> thereafter is a single instruction and a single clock cycle (I think).

The software square root algorithms I know of, from the days
before SQRT in hardware, do an approximation using fixed point
instructions, divide the exponent by two and splice it in,
with a special case for odd exponent, then two cycles (single
precision) or four cycles (double precision) of Newton-Raphson.

The initial approximation is known to be close enough for that
number of cycles. 

-- glen

There are a bunch of algorithms for  integer approximations of square
roots. They target for FPGA implementation and weak uC hardware.

See 
http://en.wikipedia.org/wiki/Methods_of_computing_square_roots    
for overview.

This 
http://www.cc.utah.edu/~nahaj/factoring/isqrt.legalize.c.html
is tricky and fast.

C source for two different square root algos and cubic roots as well: 
http://www.mikrocontroller.net/topic/246322#new

Cheers
Detlef	 

_____________________________		
Posted through www.DSPRelated.com

On Wed, 09 Apr 2014 15:15:00 -0500, Tim Wescott
<tim@seemywebsite.please> wrote:

>Subject line says it:  I need a quick, approximate (+/- 1%) 3D vector 
>length given its Cartesean components.
>

Here is a very fast 2-D method that you might be able to extend to
3-D.  It involves calculation by cases.  Take absolute values and sort
so that y < x.  You just need to cover 45 degrees of the first
quadrant.  That you can do by precomputing the secants of
2,4,6,8,...,44 degrees (22 values).  For each of these angles, use the
pre-computed secant as an approximation to the secant for angles
within 1 degree of the central angle.  So if (x,y) is closest to 36
degrees, use x * secant(36) as the approximation to the length of
(x,y).  The worst errors will be at the boundaries where the real
angle is 35 or 37 degrees.  A little checking shows that
secant(35)/secant(36) = 0.9876 and secant(37)/secant(36) = 1.023.
This is not quite within your tolerance of 1%, but that can be fixed
by using finer subdivisions of angles.  As for finding which angle is
closest for a given x and y, that you can do by a binary search using
precomuted tangents.  If y > tangent(angle) * x, the angle needs to be
higher.  Otherwise lower.  With the 22 angles listed above, the
closest angle can be found in 5 comparisons.  You could double the
resolution at the cost of twice as much precomuted lookup and just one
more comparision in the binary search.  Then you would get your 1%
accuracy.

This method so far does not involve any divides, so is great for
extremely limited micros.  If you do have a divide instruction, then
the method can be vastly improved by pre-computing not only the secant
but also the slope of the secant with respect to the tangent.  We know
the actual tangent is (y/x).  And we can precompute the nominal
tangent for the central angle.  So by subracting the actual tangent
from the nominal tangent and multipling by the slope of the secant
with respect to the tangent, we can add this first order correction
onto the precomputed secant.  So instead of approximating the secant
function by a step function we would be approximating it with a
piecewise linear function.  If this gives you more accuracy than you
need, you can reduce the number of angle subdivisions and save on some
code space.

To extend to 3-D, maybe there is a faster way, but one way is to use
this method twice.  The first time you approximate the length of
(x,y).  Then using that length as "v", apply the 2-D method again to
find the length of (v,z).

I used a techniqe like this to calculate the arctangent of (y/x) in a
very resource-limited micro with tight timing requirements.

Robert Scott
Hopkins, MN

On 09.04.2014 22:15, Tim Wescott wrote:
> Subject line says it:  I need a quick, approximate (+/- 1%) 3D vector 
> length given its Cartesean components.
> 
> Processor is a nuthin-special ARM m4, without the coprocessor.  The 
> vector is in floating point at the moment, but it could be in integers if 
> that would speed things up.
> 
> I haven't tested it yet, but my mind says "sqrt == expensive", so I want 
> to avoid that.
> 

1st: get the leading zero count of your magnitude and the first few bits
of the mantissa. Use this to lookup a good approximation of sqrt(x) from
a table.

2nd: Iterate once or twice using newton/raphson until you've reached
your desired precision. You can save some time here by fusing multiple
iteration into a single high order iteration.

Details on high order iterations for square roots are available here:

  http://numbers.computation.free.fr/Constants/Algorithms/newton.html

As usual you can trade space for speed. If you increase your lookup
table you may be able to save a iteration or two.

/Nils

Hi,

there is a Verilog implementation in iverilog's example directory (haven't
tried it myself yet).
According to the documentation, it'll give the result in 16 clock cycles
using only binary operations.
Of course, Verilog isn't too useful for an ARM processor but the algorithm
might work just the same.

original file: sqrt-virtex.v Copyright (c) 2002 Stephen Williams
(steve@icarus.com), GPL

module sqrt32(clk, rdy, reset, x, .y(acc));
   input  clk;
   output rdy;
   input  reset;

   input [31:0] x;
   output [15:0] acc;


   // acc holds the accumulated result, and acc2 is the accumulated
   // square of the accumulated result.
   reg [15:0] acc;
   reg [31:0] acc2;

   // Keep track of which bit I'm working on.
   reg [4:0]  bitl;
   wire [15:0] bit = 1 << bitl;
   wire [31:0] bit2 = 1 << (bitl << 1);

   // The output is ready when the bitl counter underflows.
   wire rdy = bitl[4];

   // guess holds the potential next values for acc, and guess2 holds
   // the square of that guess. The guess2 calculation is a little bit
   // subtle. The idea is that:
   //
   //      guess2 = (acc + bit) * (acc + bit)
   //             = (acc * acc) + 2*acc*bit + bit*bit
   //             = acc2 + 2*acc*bit + bit2
   //             = acc2 + 2 * (acc<<bitl) + bit
   //
   // This works out using shifts because bit and bit2 are known to
   // have only a single bit in them.
   wire [15:0] guess  = acc | bit;
   wire [31:0] guess2 = acc2 + bit2 + ((acc << bitl) << 1);

   (* ivl_synthesis_on *)
   always @(posedge clk or posedge reset)
      if (reset) begin
	 acc = 0;
	 acc2 = 0;
	 bitl = 15;
      end else begin
	 if (guess2 <= x) begin
	    acc  <= guess;
	    acc2 <= guess2;
	 end
	 bitl <= bitl - 5'd1;
      end

endmodule // sqrt32
	 

_____________________________		
Posted through www.DSPRelated.com