Channels ▼

Al Williams

Dr. Dobb's Bloggers

Square Root Part 2 (or 1.414, If You Prefer)

May 13, 2011

Last time I talked about an efficient (and old) technique for generating square root estimates. As you might expect, there are many different ways to do this calculation, some of which are more amenable to computers than others. I was browsing http://is.gd/UZgRhd to see if there were any other methods that made sense on restricted CPUs.

One algorithm that caught my eye was the "digit-by-digit" binary algorithm. If you aren't old enough to have done square roots by hand (or, like me, you just forgot), you might start by reading the base 10 version first — if you understand that one, the base 2 version is the same idea but applied to base 10.

This is one of those algorithms that is easier to code than it is to explain. So with that in mind, here's the code slightly adapted from the code in the Wikipedia article:

/* Another square root -- Williams http://www.ddj.com/embedded */
#include <stdio.h>

#ifdef TEST
#include <stdlib.h> /* random numbers */
#endif

/* Work out the square root bit by bit */
#define BITS 32   // set for the number of bits in an int


int bit_isqrt(int num, float *v) 
   {
     int num0=num;
     int result=0; 
     // The second-to-top bit is set: 1L<<14 or 1L<<30 or whatever for your sizeof(int)
     int tbit=1<<(BITS-2); 
     
     // tbit starts at the highest power of four <= the argument.
     while (tbit>num0)
       tbit>>=2;
     
     while (tbit!=0) 
       {
	 if (num0>=result+tbit) 
	   {
	     num0-=result+tbit;
	     result=(result>>1)+tbit;
	   }
	 else
	   result>>=1;
	 tbit>>=2;
       }
     if (v) *v=(float)result;
     return result;
   }


#ifdef TEST

/* Simple routine to drive the calculation and print stats */
void debugsqrt(int x)
{
  float v;
  int iv=bit_isqrt(x,&v);
  printf("Square root of %d = %d  %f",x,iv,v);
#ifdef TEST
  v*=v;
  printf(" Check=%f error=%f (%f%%)\n",v,v-x,((v-x)/x)*100);
#else
  putchar('\n');
#endif
}



/* Simple main */

// generate random numbers from 0 to MAXTESTINT-1
#define MAXTESTINT (10000)
// Random number seed change for different sequence of randoms */
#define SEED 1

void main(int argc, char *argv[])
{
  int i;
  srand(SEED);
 
  debugsqrt(9);  /* known test case ! */
  debugsqrt(512);  /* known test case ! */
  for (i=0;i<100;i++) 
    {
      debugsqrt(rand()%MAXTESTINT);
    }
}


#endif

If you compare the test harness output with last time's program, you'll see that this algorithm is a little less accurate because it rounds down. So the square root of 512 shows up as 22 instead of 23 (it's really about 22.6, so 23 is closer). You could fix that by always calculating the error of the answer and the answer plus one. But that's two extra multiplies. Or you can do fixed-point math to get some data to the left of the decimal point.

You can extend the algorithm to do fixed point math pretty easily. Or if you don't mind a little floating point, you can first multiply the input by a perfect square scaling factor. Then the answer can be divided (that's the floating-point part) by the square root of the scaling factor.

For example, if you wanted the square root of 2 with a few digits of precision, you could multiply by 256 (shift left 8 times) and actually find the square root of 512. Then divide the answer by 16. That gives you 1.375, which is not too far off from 1.414 (about 3%). I think I'll stick with the old EINAC algorithm for now. But if I had to do a fixed-point square root, this algorithm would definitely be one to try.

Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.
 


Video