Integer Square Root (Was Re: # to the nth power)

[Chris Torek]

First, notes on `pow':

   Various people claim that pow(x,y) is `not as good' as FORTRAN's `**'
   operator because it fails when x is negative but y is integral.  This
   is false (see ANSI C standard X3.159-1989, p. 116).  Others claim that
   a function is `not as good' as an operator because it cannot be inlined;
   this is also false.  The Standard gives license for compilers to
   recognise Standard functions as `intrinsics' once the standard headers
   have been `#include'd (this includes eliminating conversion to and from
   double precision whenever the result is the same).  Still others claim
   that it is not as good because the notation sucks.  This is something of
   an opinion and I will not argue either side (though I agree with both).

Now on to integer square root.

The last time this came up there was a bit of a speed war.  The following
is the result.  Note that it depends on 32-bit integers (but the changes
to make it work on some other size are obvious).

/*
 * Integer square root routine, good for up to 32-bit values.
 * Note that the largest square root (that of 0xffffffff) is
 * 0xffff, so the result fits in a regular unsigned and need
 * not be `long'.
 *
 * Original code from Tomas Rokicki (using a well known algorithm).
 * This version by Chris Torek, University of Maryland.
 *
 * This code is in the public domain.
 */
unsigned
root(v)
	register unsigned long v;
{
	register long t = 1L << 30, r = 0, s;

#define STEP(k) \
	s = t + r; \
	r >>= 1; \
	if (s <= v) { \
		v -= s; \
		r |= t; \
	}
	STEP(15); t >>= 2;
	STEP(14); t >>= 2;
	STEP(13); t >>= 2;
	STEP(12); t >>= 2;
	STEP(11); t >>= 2;
	STEP(10); t >>= 2;
	STEP(9); t >>= 2;
	STEP(8); t >>= 2;
	STEP(7); t >>= 2;
	STEP(6); t >>= 2;
	STEP(5); t >>= 2;
	STEP(4); t >>= 2;
	STEP(3); t >>= 2;
	STEP(2); t >>= 2;
	STEP(1); t >>= 2;
	STEP(0);
	return r;
}

#ifdef MAIN
#include <stdio.h>
#include <math.h>

main()
{
	unsigned long l;
	char buf[100];

	for (;;) {
		(void) printf("gimme a number> ");
		if (fgets(buf, sizeof buf, stdin) == NULL)
			break;
		/* should use strtoul here but some do not have it */
		if (sscanf(buf, " 0x%lx", &l) != 1 &&
		    sscanf(buf, " 0%lo", &l) != 1 &&
		    sscanf(buf, "%lu", &l) != 1 &&
		    sscanf(buf, "%lx", &l) != 1)
			(void) printf("that was not a number\n");
		else
			(void) printf("root(%lu) => %u; sqrt(%lu) => %.17g\n",
			    l, root(l), l, sqrt((double)l));
	}
	exit(0);
}
#endif
-- 
In-Real-Life: Chris Torek, Univ of MD Comp Sci Dept (+1 301 405 2750)
Domain:	chris at cs.umd.edu	Path:	uunet!mimsy!chris