Path: utzoo!attcan!utgpu!jarvis.csri.toronto.edu!mailrus!csd4.milw.wisc.edu!cs.utexas.edu!uunet!telxon!scottf
From: scottf@telxon.UUCP (Scott Fluhrer)
Newsgroups: comp.lang.c
Subject: Re: Integer square root routine needed.
Message-ID: <77@telxon.UUCP>
Date: 1 Aug 89 20:51:06 GMT
References: <7415@ecsvax.UUCP> <5392@ficc.uu.net>
Reply-To: scottf@telxon.UUCP (Scott Fluhrer)
Organization: Telxon Corporation, Akron OH
Lines: 56

In article <5392@ficc.uu.net> cliff@ficc.uu.net (cliff click) writes:
>In article <7415@ecsvax.UUCP>, utoddl@ecsvax.UUCP (Todd M. Lewis) writes:
>> I need C source to a reasonably fast integer square root routine.
>Here's mine:
>
>
>long sqrt(a)                    /* Integer square root */
>long a;
>{
>register unsigned long b,c;
>
>  if( a <= 1 ) return a;        /* Can't handle zero, one */
>  c = b = a;                    /* Start at start */
>  while( b ) c>>=1,b>>=2;       /* Get sqrt accurate to 1 bit in c */
>  for( b = 0; b < 4; b++ )      /* 4 iterations of newtons */
>    c = ((a+c*c)>>1)/c;         /* Newtons method */
>/*  c = ((a>>1)+((c*c)>>1))/c;  /* Slightly slower, prevents overflow, loses 1 bit */
>  return c;                     /* Return it */
>}

Here's a 'better' version:

/* 'unsigned' == 'we don't have to worry about negative values' :-) */
long sqrt(a)                    /* Integer square root */
unsigned long a;
{
	unsigned long b,c;

	if ( a <= 1 ) return a; /* Can't handle 0.  Our initial guess */
				/* can't handle 1 either */
  	c = b = a;              /* Start at start */
	while (b) c>>=1,b>>=2;  /* Get sqrt accurate to 1 bit in c */
	b = (c + a/c)>>1;	/* Do 1 iteration of newtons */
	do {			/* Now keep on doing iterations until */
		c = b;		/* the value ceases to decrease */
		b = (c + a/c)>>1;
	} while (b < c);
	return c;		/* Return the minimum value */
}

By 'better', I mean:

* Always returns the 'right' answer, with 'right' being 'the squareroot rounded
  down to the next lowest integer' (it is an interesting algebraic exercise
  to verify this).  The previous version got it wrong occassionally (like at
  values 80 and 99) and could overflow

Now, you don't say if exactness (except for overflow, the other routine is only
off by 1) is required, but since it is not difficult, and doesn't cost much
time (my version may even go faster, since it doesn't do any multiplies), you
might as well get it right...

BTW: followup's probably shouldn't go to comp.lang.c...
------
Poncho, the Faithful Sidekick
aka scottf@telxon