[comp.lang.c] Integer Square Root

weaver (11/05/90)

Here is an integer square root function I wrote. It is derived from 
Newton's method, mentioned earlier by Henry Spencer. After trying 
to learn algorithms used in hardware for square root, it seems to
me that none would translate to more efficient C than Newton's
method, unless divide or multiply are really slow.

Caveats:

I am not a mathematician.
This has not been exhaustively tested.

Cut here---------------------------------------- 

#include <stdio.h>
main(argc, argv)
int             argc;
char          **argv;
{
    int             n;
    argv++;
    argc--;
    while (argc-- > 0) {
	n = atoi(*argv++);
	printf("Ans = %d\n", int_sqrt(n));
    }
}
int 
int_sqrt(n)
int             n;
{
    /* Integer version of Newton's method applied to square root. 	*/
    /* Similar to floating point, but need beware of overflows.		*/
    /* Michael Gordon Weaver 1990/11/04	 weaver@weitek.COM		*/
    int             i, j, k, d;
    int             k = 0;


    /* range check */
    if (n <= 0)
	return -1;

    /* initial estimate (j) */
    i = j = n;
    while (i > 0) {
	i >>= 2;
	j >>= 1;
    }

    /* iterate */
    do {
	d = ((j * j) - n) / (2 * j);
	printf("Iter %d, guess = %d, est err = %d, remainder = %d\n",
	       ++k, j, d, (j * j) - n);
	j = j - d;
    } while (d != 0);

    /* we want answer < sqrt(n) when inexact */
    if ((j * j) == n)
	return j;
    else
	return j - 1;
}

chris@mimsy.umd.edu (Chris Torek) (11/05/90)

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@cs.umd.edu	Path:	uunet!mimsy!chris