Fixed point sqrt
Meessen Christophe
Hines) wrote:
>
> > Well? Does anyone know of a quick and dirty way to get a reasonably accurate
> > square root of a fixed point integer number?
>
> Hmmm - I have a long integer square root routine.
Dave Hines proposed a sqrt function receiving an unsigned long (X) and
returning
an unsigned long. The remainder, an unsigned long, is accessible as
'result2'.
The algorithm uses 4 unsigned long variables.
Thus it's a long -> long function.
I suggest two algorithms using the same variables. The first one will
compute a sqrt of a long, returning a fixed point value. The second will
compute the sqrt of a fixed point value, returning a fixed point value.
Thus it's a long -> fixed or fixed -> fixed function.
We will suppose a fixed is equivalent to a long, and longs are 32 bits
wide.
/*
* fixed sqrtL2F( long X );
*
* Long to fixed point square roots.
* RETURNS the fixed point square root of X (long).
* REQUIRES X is positive.
*
* Christophe MEESSEN, 1993.
*/
fixed sqrtL2F( long X ) {
unsigned long t, q, b, r;
r = X;
b = 0x40000000;
q = 0;
while( b > 0 ) {
t = q + b;
if( r >= t ) {
r = r - t;
q = t + b;
}
r = r * 2; /* shift left 1 bit */
b = b / 2; /* shift right 1 bit */
}
if( r >= q )
q = q + 1;
return( q );
}
/*
* fixed sqrtF2F( fixed X );
*
* Fixed to fixed point square roots.
* RETURNS the fixed point square root of X (fixed).
* REQUIRES X is positive.
*
* Christophe MEESSEN, 1993.
*/
fixed sqrtF2F ( fixed X ) {
unsigned long t, q, b, r;
r = X;
b = 0x40000000;
q = 0;
while( b >= 256 ) {
t = q + b;
if( r >= t ) {
r = r - t;
q = t + b;
}
r = r * 2; /* shift left 1 bit */
b = b / 2; /* shift right 1 bit */
}
q = q / 256; /* shift right 8 bits */
return( q );
}
for memory I just rewrote the long to long sqrt function.
/*
* long sqrtL2L( long X );
*
* Long to long point square roots.
* RETURNS the integer square root of X (long).
* REMAINDER is in the local variable r of type long on return.
* REQUIRES X is positive.
*
* Christophe MEESSEN, 1993.
*/
long sqrtL2L( long X ) {
unsigned long t, q, b, r;
r = X;
b = 0x40000000;
q = 0;
while( b >= 256 ) {
t = q + b;
q = q / 2; /* shift right 1 bit */
if( r >= t ) {
r = r - t;
q = q + b;
}
b = b / 4; /* shift right 2 bits */
}
return( q );
}
This functions have been tested. But I'll be pleased to receive any
comments.
Please respond through mail: Mee...@marvax.in2p3.fr