Fast integer distance formula

Started by October 4, 2006
```I have an approximation for the Pythagorean distance formula (magnitude of
vector [x,y]) using integer arithmetic that I would like to improve.  Currently
I am using this:

if(x<0) x=-x;
if(y<0) y=-y;
if(x < y)
{
int t = x;
x = y;
y = t;     	// ensures that x >= y
}
z = (y < ((13107 * x)>>15)) ?        // * (.4)
(x + ((y * 6310)>>15)) :       // * (.192582403)
(((x * 27926)>>15)             // * (.852245894)
+ ((y * 18414)>>15));      // * (.561967668)
//..(linear approximation to within 2% of the Pythagorean distance formula)..

This is for an ARM processor in an application where I cannot afford the time
for any floating point operations.  The integer values of x and y come from an
array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
above formula.  The final result (z) will be shifted right one bit before being
stored back into an array of 16-bit integers, since the distance formula can
potentially extend the range of x and y by one bit.  In case you are wondering,
this is part of calculating the power spectrum at the end of an FFT.

The improvement I am looking for is in accuracy.  I would like to try for a
4-fold improvement in accuracy (.5%) without substantially increasing the
running time of what I have now.  Does anyone know of a better approximation
that is almost as fast?

Robert Scott
Ypsilanti, Michigan
```
```Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve.  Currently
> I am using this:
>
> if(x<0) x=-x;
> if(y<0) y=-y;
> if(x < y)
> {
>   int t = x;
>   x = y;
>   y = t;     	// ensures that x >= y
> }
> z = (y < ((13107 * x)>>15)) ?        // * (.4)
>           (x + ((y * 6310)>>15)) :       // * (.192582403)
>           (((x * 27926)>>15)             // * (.852245894)
>              + ((y * 18414)>>15));      // * (.561967668)
> //..(linear approximation to within 2% of the Pythagorean distance formula)..
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations.  The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula.  The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit.  In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy.  I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now.  Does anyone know of a better approximation
> that is almost as fast?
>
>
>
> Robert Scott
> Ypsilanti, Michigan

I still like the tried and true isqrt((x^2+y^2)<<16).  You can see we
scale up the sum of x^+y^ by shifing left 16 times so we have a fixed
point number with the interger in the upper 16 bits.   When we take the
square root of at 16.16 we get an 8.8 result.  This meets your .5%
criteria.  There are plenty of fast and simple square root routines.
We use a table to find the must significant 8 bits and then do one
Newton iteration to get the lower 8 bits but then I have a divide.
If you don't there are other fast and simple routines on the internet.
I think our isqrt routine takes about 8 microseconds on a 40 Mhz 196.
The squaring and summing would take two extra microseconds.

Peter Nachtwey

```
```Robert,

I do not have the book handy, but my recollection is that methods very
similar to the method you are using, with varying levels of accuracy,
are detailed in 'Digital Signal Processing in Communications Systems'
by Marvin E. Frerking.  You might look there and see if you can find
something useful.

Maybe someone here who has access to the book could see if my

Dirk Bell
DSP Consultant

Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve.  Currently
> I am using this:
>
> if(x<0) x=-x;
> if(y<0) y=-y;
> if(x < y)
> {
>   int t = x;
>   x = y;
>   y = t;     	// ensures that x >= y
> }
> z = (y < ((13107 * x)>>15)) ?        // * (.4)
>           (x + ((y * 6310)>>15)) :       // * (.192582403)
>           (((x * 27926)>>15)             // * (.852245894)
>              + ((y * 18414)>>15));      // * (.561967668)
> //..(linear approximation to within 2% of the Pythagorean distance formula)..
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations.  The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula.  The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit.  In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy.  I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now.  Does anyone know of a better approximation
> that is almost as fast?
>
>
>
> Robert Scott
> Ypsilanti, Michigan

```
```Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve.  Currently
> I am using this:
<snip>
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations.  The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula.  The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit.  In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy.  I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now.  Does anyone know of a better approximation
> that is almost as fast?

I like this algorithm. It is snipped from my OpenLPC fixed point codec.
It uses only simple operations. It is likely more accurate than you
need, but I think you can simply truncate the algorithm sooner for less
precision. At the end PRECISION is the fractional bits you are using in

static fixed32 fixsqrt32(fixed32 x)
{

unsigned long r = 0, s, v = (unsigned long)x;

#define STEP(k) s = r + (1 << k * 2); r >>= 1; \
if (s <= v) { v -= s; r |= (1 << k * 2); }

STEP(15);
STEP(14);
STEP(13);
STEP(12);
STEP(11);
STEP(10);
STEP(9);
STEP(8);
STEP(7);
STEP(6);
STEP(5);
STEP(4);
STEP(3);
STEP(2);
STEP(1);
STEP(0);

return (fixed32)(r << (PRECISION / 2));
}

--
Phil Frisbie, Jr.
Hawk Software
http://www.hawksoft.com
```
```The comp.dsp "DSP Tricks" at
http://www.dspguru.com/comp.dsp/tricks/tricks.htm contains two entries
on magnitude estimation, which is the same computation.  In particular,
this one http://www.dspguru.com/comp.dsp/tricks/alg/mag_est.htm should
tell you all that you want to know (altho IIRC, the C program contains
a typo -- the coefficients Alpha and Beta are printed backwards in the
printed table (which one multiples the max and which the min?)) .

cheers,
jerry

Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve.  Currently
> I am using this:
>
> if(x<0) x=-x;
> if(y<0) y=-y;
> if(x < y)
> {
>   int t = x;
>   x = y;
>   y = t;     	// ensures that x >= y
> }
> z = (y < ((13107 * x)>>15)) ?        // * (.4)
>           (x + ((y * 6310)>>15)) :       // * (.192582403)
>           (((x * 27926)>>15)             // * (.852245894)
>              + ((y * 18414)>>15));      // * (.561967668)
> //..(linear approximation to within 2% of the Pythagorean distance formula)..
>
> This is for an ARM processor in an application where I cannot afford the time
> for any floating point operations.  The integer values of x and y come from an
> array of signed 16-bit numbers, but x and y themselves are 32-bit numbers in the
> above formula.  The final result (z) will be shifted right one bit before being
> stored back into an array of 16-bit integers, since the distance formula can
> potentially extend the range of x and y by one bit.  In case you are wondering,
> this is part of calculating the power spectrum at the end of an FFT.
>
> The improvement I am looking for is in accuracy.  I would like to try for a
> 4-fold improvement in accuracy (.5%) without substantially increasing the
> running time of what I have now.  Does anyone know of a better approximation
> that is almost as fast?
>
>
>
> Robert Scott
> Ypsilanti, Michigan

```
```Phil Frisbie, Jr. wrote:
> Robert Scott wrote:
>> I have an approximation for the Pythagorean distance formula
>> (magnitude of
>> vector [x,y]) using integer arithmetic that I would like to improve.
>> Currently
>> I am using this:
> <snip>
>>
>> This is for an ARM processor in an application where I cannot afford
>> the time
>> for any floating point operations.  The integer values of x and y come
>> from an
>> array of signed 16-bit numbers, but x and y themselves are 32-bit
>> numbers in the
>> above formula.  The final result (z) will be shifted right one bit
>> before being
>> stored back into an array of 16-bit integers, since the distance
>> formula can
>> potentially extend the range of x and y by one bit.  In case you are
>> wondering,
>> this is part of calculating the power spectrum at the end of an FFT.
>>
>> The improvement I am looking for is in accuracy.  I would like to try
>> for a
>> 4-fold improvement in accuracy (.5%) without substantially increasing the
>> running time of what I have now.  Does anyone know of a better
>> approximation
>> that is almost as fast?
>
> I like this algorithm. It is snipped from my OpenLPC fixed point codec.
> It uses only simple operations. It is likely more accurate than you
> need, but I think you can simply truncate the algorithm sooner for less
> precision. At the end PRECISION is the fractional bits you are using in
>
> static fixed32 fixsqrt32(fixed32 x)
> {
>
>     unsigned long r = 0, s, v = (unsigned long)x;
>
> #define STEP(k) s = r + (1 << k * 2); r >>= 1; \
>     if (s <= v) { v -= s; r |= (1 << k * 2); }
>
>     STEP(15);
>     STEP(14);
>     STEP(13);
>     STEP(12);
>     STEP(11);
>     STEP(10);
>     STEP(9);
>     STEP(8);
>     STEP(7);
>     STEP(6);
>     STEP(5);
>     STEP(4);
>     STEP(3);
>     STEP(2);
>     STEP(1);
>     STEP(0);
>
>     return (fixed32)(r << (PRECISION / 2));
> }
>

That only seems to calculate a part of the answer, leaving lots of zeros
at the end. Why not use something that calculates as many bits as it
can, like:

int32_t isqrt32(int32_t h)
{
int32_t x;
int32_t y;
int i;

/* The answer is calculated as a 32 bit value, where the last
16 bits are fractional. */
/* Calling with negative numbers is not a good idea :-) */
x =
y = 0;
for (i = 0;  i < 32;  i++)
{
x = (x << 1) | 1;
if (y < x)
x -= 2;
else
y -= x;
x++;
y <<= 1;
if ((h & 0x80000000))
y |= 1;
h <<= 1;
y <<= 1;
if ((h & 0x80000000))
y |= 1;
h <<= 1;
}
return  x;
}

Of course, both these routines seem like overkill for what the original
poster needs. :-)

Steve
```
```Robert Scott wrote:
> I have an approximation for the Pythagorean distance formula (magnitude of
> vector [x,y]) using integer arithmetic that I would like to improve.  Currently

Look for a paper entitled

Filip, A. E. , "A baker's dozen magnitude approximations and their
detection statistics"  IEEE Trans on Aerospace and Elcectronic Systems,
Jan 1976

It is a classic work with many tricks similar to the one you are using.

With a little trickery, several of them could be implemented with little
or no multiplication.

--
Mark Borgerding
3dB Labs, Inc
Innovate.  Develop.  Deliver.
```
```"Phil Frisbie, Jr." <phil@hawksoft.com> writes:
[snip]
> I like this algorithm. It is snipped from my OpenLPC fixed point codec.
> It uses only simple operations. It is likely more accurate than you
> need, but I think you can simply truncate the algorithm sooner for less
> precision. At the end PRECISION is the fractional bits you are using in
[snip]

Isn't computing the exact square root of an integer
faster than the algorithm you show?
```
```Everett M. Greene wrote:

> "Phil Frisbie, Jr." <phil@hawksoft.com> writes:
> [snip]
>
>>I like this algorithm. It is snipped from my OpenLPC fixed point codec.
>>It uses only simple operations. It is likely more accurate than you
>>need, but I think you can simply truncate the algorithm sooner for less
>>precision. At the end PRECISION is the fractional bits you are using in
>
> [snip]
>
> Isn't computing the exact square root of an integer
> faster than the algorithm you show?

That depends on your target CPU/DSP. I am now targeting ARM most of the
time so I will make a note to profile that code and see.

--
Phil Frisbie, Jr.
Hawk Software
http://www.hawksoft.com
```
```"Phil Frisbie, Jr." <phil@hawksoft.com> writes:
> Everett M. Greene wrote:
>
> > "Phil Frisbie, Jr." <phil@hawksoft.com> writes:
> > [snip]
> >
> >>I like this algorithm. It is snipped from my OpenLPC fixed point codec.
> >>It uses only simple operations. It is likely more accurate than you
> >>need, but I think you can simply truncate the algorithm sooner for less
> >>precision. At the end PRECISION is the fractional bits you are using in