On Sat, 30 Mar 2002, John Dammeyer wrote:

> This works.  Bytecraft C has long equal to 16 bit integer so you can see
> functions have a high 16 bit and a low 16 bit.  EQ32 sets up the
> variables and eq32() does the subtration to test for equality.
>
>
>
> long
> sqrt32(void) {
>   WORD r @ local6;
>   WORD i;
>         // Handle special cases.
>     EQ32(ArgsqrtHi,ArgsqrtLo,0,0)
>     if ( eq32() ) return(0);
>     i = 0;
>     r = 1;
>     do {
>             SUB32(ArgsqrtHi,ArgsqrtLo,0,r);
>         r += 2;
>         i++;
>     }
>     while (ArgsqrtHi >= 0);
>
>     return(i);
> }


I've always liked this algorithm for its simplicity. The sum of the first
N odd integers is equal to N^2... Unfortunately, this can be rather slow
if N is a large number. This though inspired me to dust off and polish
some old square root algorithms. Each of theses makes only one iteration
for each needed bit (e.g. 12 iterations for 12 bits in the result):


unsigned int isqrt1(unsigned int N)
{

  unsigned int s1 = 1<<(6+RESULT_SIZE);
  unsigned s2 = 0;

  do {

    if((s1 + s2) <= N) {

      N = N - s2 -s1;
      s2 = s2 | (s1<<1);
    }

    s1 = s1 >> 1;
    N  = N  << 1;

  } while(s1 > (1<<6));

  return(s2>>(4 + RESULT_SIZE/2));


}

This is the one I developed for the first 16-bit square root routine for
the PIC about 3 or 4 years ago. BTW, I don't claim to be the first to
discover this. Especially since shortly after writing this I came across a
book written in the '60's by Ivan (?) Flores. I don't even recall the
title of the book. But in there I came across the binary square root
algorithm. It was not too unlike the routine above, at least in concept.
This inspired the second square root routine I wrote for the PIC. Here it
is more or less in C:

unsigned int isqrt2(unsigned int N)
{

  unsigned int mask = 3<<(6+RESULT_SIZE);
  unsigned int root = 1<<(6+RESULT_SIZE);

  do {

    if(root <= N) {

      N = N - root;
      root |= mask;
    }

    mask >>= 1;
    root ^= mask;
    N  <<= 1;

  } while(mask > (1<<6));

  return(root>>(4 + RESULT_SIZE/2));


}

For reasons that are not completely obvious, these routines work :). The
hand written assembly code that implements these are fairly well
commented.

In each routine, RESULT_SIZE is a compile time defined constant (that must
be even...) that defines the size of the result in bits. So for a 12-bit
routine:

#define RESULT_SIZE 12


Here's a simple loop to run through all possible test cases.


int
main(int argc, char **argv)
{
  int i;

  printf(" An integer square root routine\n");
  for(i=0; i<((1<<RESULT_SIZE) - 1); i+=1)
    if(isqrt2(i) != isqrt1(i))
       printf(" %d %d %d\n", i,isqrt1(i),isqrt2(i));

  return 0;
}


FWIW, sdcc compiles isqrt1 to 0x47 instructions, and isqrt2 to 0x33
instructions. The hand written version is about 0x1c words. Hmm, I
need to work on sdcc some more...

Scott

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
email listserv@mitvma.mit.edu with SET PICList DIGEST in the body