[comp.lang.c] sqroot of a 2**n number ?

oz@yunexus.UUCP (Ozan Yigit) (05/13/89)

given any number that is an exact power of two, finding the square
root: this can be done by counting to the bit position on some
architectures. Does anyone know of a really portable version of
this, assuming the code below will not work on all architectures ?
(no sqrt call)

	if (n > 0)	/* n > 0 && n == 2**X */
		for (i = 0; !(n & 1); i++)
			n >>= 1;

oz
-- 
use the source, luke !!     	        Usenet:    oz@nexus.yorku.ca
uh... we forgot to tell you...          ......!uunet!utai!yunexus!oz
it is unintelligible, but hey, you      Bitnet: oz@[yulibra|yuyetti]
got it, for free (!).                   Phonet: +1 416 736-5257x3976

hascall@atanasoff.cs.iastate.edu (John Hascall) (05/13/89)

In article <1906@yunexus.UUCP> oz@yunexus.UUCP (Ozan Yigit) writes:
>given any number that is an exact power of two, finding the square
>root: this can be done by counting to the bit position on some
>architectures. Does anyone know of a really portable version of
>this, assuming the code below will not work on all architectures ?
>(no sqrt call)
 
>	if (n > 0)	/* n > 0 && n == 2**X */
>		for (i = 0; !(n & 1); i++)
>			n >>= 1;
 
  By an amazing coincidence, I happen to have just written something
  of a similar nature--everyone can feel free to use it (or take potshots
  at it).

  John Hascall

-------------------- cut here ---------------------------------------
/*
** HIBIT        - find the highest set bit in an int, this can be used
**                      as a quick base-2 logarithm.
**
**      John Hascall / 12-May-1989
**
** Description:  This algorithm uses a divide and conquer strategy and is
**      O(log n), n the number of bits.  The routine works by dividing   
**      the int, v,  in half and testing that half, this continues until 
**      the last half (2 bits) is tested.  Note that if no bits are set  
**      zero is returned.
**
**      Hint: m goes: 0xffffffff, 0x0000ffff, 0x00ff00ff, 0x0f0f0f0f,
**                    0x33333333, 0x55555555
*/

#define HALF    16              /* one half of the number of bits in an int */

unsigned int hibit( v )
unsigned int v;
{
        int             i;                      /* size of chunks to test */
        unsigned int    m;                      /* mask used to divide v */ 
        unsigned int    s;                      /* hibit accumulated sum */ 

        s = 0;
        m = ~0;                                 /* all ones (0xffffffff) */
        for ( i = HALF; i > 0; i >>= 1 ) {
                m ^= (m << i);  
                if (~m & v) {   
                        s |= i; 
                        v >>= i;
                }   
        }           
        return( s );
}
-------------------------- end -------------------------------------

oz@yunexus.UUCP (Ozan Yigit) (05/14/89)

In article <4890@umd5.umd.edu> setzer@wam.umd.edu writes:

>int sqrt_of_2_to_the_n(n)
>int n;
>{
>    return (1 << (n / 2));  /* integer division */
>}

Say what ?? (1 << (16 / 2) == 256 !!! I always thought sqrt(16) == 4.

oz
-- 
use the source, luke !!     	        Usenet:    oz@nexus.yorku.ca
uh... we forgot to tell you...          ......!uunet!utai!yunexus!oz
it is unintelligible, but hey, you      Bitnet: oz@[yulibra|yuyetti]
got it, for free (!).                   Phonet: +1 416 736-5257x3976

envbvs@epb2.lbl.gov (Brian V. Smith) (05/14/89)

	In article <1928@yunexus.UUCP> oz@yunexus.UUCP (Ozan Yigit) writes:
	>In article <4890@umd5.umd.edu> setzer@wam.umd.edu writes:
	>
	>>int sqrt_of_2_to_the_n(n)
	>>int n;
	>>{
	>>    return (1 << (n / 2));  /* integer division */
	>>}

	>Say what ?? (1 << (16 / 2) == 256 !!! I always thought sqrt(16) == 4.

Try reading the name of the function.  It returns the square root 
of 2**n, not n!

>use the source, luke !!     	        Usenet:    oz@nexus.yorku.ca

Use the brain, luke !!

_____________________________________
Brian V. Smith    (bvsmith@lbl.gov)
Lawrence Berkeley Laboratory
We don't need no signatures!

ai04@umd5.umd.edu (CMSC 620) (05/15/89)

In article <1928@yunexus.UUCP> oz@yunexus.UUCP (Ozan Yigit) writes:
>In article <4890@umd5.umd.edu> setzer@wam.umd.edu writes:
>
>>int sqrt_of_2_to_the_n(n)
>>int n;
>>{
>>    return (1 << (n / 2));  /* integer division */
>>}
>
>Say what ?? (1 << (16 / 2) == 256 !!! I always thought sqrt(16) == 4.

Apparently my cancel didn't get through.  After posting my message, I realized
that your argument was 2**n instead of n.  The argument of my programs is the
exponent.  BTW, if you send the value 'i' from your patch of code to my
program, everything works fine.  My apologies for the mixup.  I'm going to
find my article and try to kill it again.

--
William Setzer
setzer@wam.umd.edu

trt@rti.UUCP (Thomas Truscott) (05/15/89)

> In article <1906@yunexus.UUCP> oz@yunexus.UUCP (Ozan Yigit) writes:
> >given any number that is an exact power of two, finding the square
> >root: this can be done by counting to the bit position on some
> >architectures. Does anyone know of a really portable version of this ...

I believe the following is portable code:
	/* "n" is type "int", exact positive power of two */
	int i;
	for (i = 0; (1<<i) != n; i++)
		;
	/* "i" is now lg2(n), to get square root we do: */
	sqrt_n = 1 << (i/2);
	if ((i%2) != 0)		/* odd power-of-two */
		sqrt_n *= sqrt(2.0);	/* oops! need sqrt after all */

Portability and performance sometimes clash,
and I think this problem is an example.
Typical C machines only support 31 different exact powers of two
(2^i, 0 <= i <= 30) which suggests a table look up:

unsigned char lg2_hashtable[] = {
	0,0,1,26,2,23,27,0,3,16,24,30,28,11,0,13,4,7,17,
	0,25,22,0,15,29,10,12,6,0,21,14,9,5,20,8,19,18,0,
};
#define lg2(twotoN) lg2_hashtable[(twotoN)%(37)]

I generated this with a simple program.
The problem is: this code is not portable!
We could build a slightly bigger table, to handle "64-bit" ints,
but it would be slightly slower and *still* be unportable
(to a machine with 144-bit integers, for example).
We could generate the table at runtime,
but (a) it would have to be allocated dynamically,
(b) it would be tricky to generate it portably,
and (c) there would be a bigger-yet performance hit.
	Tom Truscott