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

# Fast integer distance formula

Started by ●October 4, 2006

Reply by ●October 4, 20062006-10-04

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, MichiganI 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

Reply by ●October 4, 20062006-10-04

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 recollection is correct and help you out. 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

Reply by ●October 4, 20062006-10-04

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 your fixed point code. 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

Reply by ●October 4, 20062006-10-04

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

Reply by ●October 5, 20062006-10-05

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 > your fixed point code. > > 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

Reply by ●October 5, 20062006-10-05

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. CurrentlyLook 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.

Reply by ●October 5, 20062006-10-05

"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 > your fixed point code.[snip] Isn't computing the exact square root of an integer faster than the algorithm you show?

Reply by ●October 5, 20062006-10-05

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 >>your fixed point code. > > [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

Reply by ●October 6, 20062006-10-06

"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 > >>your fixed point code. > > > > [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.It would be interesting to know the results of the comparison.