Path: chuka.playstation.co.uk!news
From: gil@snsys.com (Gil Jaysmith)
Newsgroups: scee.yaroze.programming.codewarrior
Subject: Re: "No FP... no comment"
Date: Fri, 15 Aug 1997 17:07:46 GMT
Organization: SN Systems
Lines: 126
Message-ID: <5t221t$d3t2@chuka.playstation.co.uk>
References: <01bca987$927743a0$e6567ec2@default> <01bca99a$10782be0$b8567ec2@default>
Reply-To: gil@snsys.com
NNTP-Posting-Host: gil.snsys.com
X-Newsreader: Forte Free Agent 1.0.82

"Jon Prestidge (alias Moose)" <zmoose@premier.co.uk> wrote:

>I've had a bit of an experiment since I submitted to above item and
>produced the listing below.   This seems to give OK-ish results, but I
>would have thought it was not nearly the most efficient way of doing it....
>also you can't just throw any-old figure at it. 

Integral square roots for 16-bit numbers will fit into a 256-byte
lookup table (max 8 lookups and comparisons with binary chop) ...
similarly square roots for 32-bit numbers will fit into a 64K lookup
table (so don't take square roots of numbers that high :) although I
seem to remember there's a shortcut which shrinks that table
significantly. Don't remember exactly what it is, though.

There's the 'fencing' approach which is traditionally used to show
recursion at work but which being tail-recursive is much faster when
uncoiled into iteration:
The square root of N is somewhere between 1 (the initial lower bound)
and N+1 (the initial upper bound). So repeatedly:
- calculate the average of the lower and upper bounds and square it
- if this is greater than N then set the upper bound to the average,
otherwise set the lower bound to the average
- keep going till the lower bound is one less than the upper bound, at
which point the lower bound is the square root of N

NB I wouldn't expect this to be overly efficient either.

Can you tell it's Friday? It's like stretching just before you go home
for the weekend...

Here is "Fermat's method":

/* Integer square root */
/* This is a stupid way to do square root, but it is portable */
long isqrt (x) long x; {
  long s; long remainder;
  if (x<=0) return 0;
  remainder = x; s = 1;
  while ( remainder >= s ) {
          remainder -= s; 
          s += 2;
  }
  return (s-1)>>1;

}
Here is "Halleck's method":

/* Integer square root */
long isqrt (x) long x;{
  long   squaredbit, remainder, root;
   if (x<1) return 0;
  
   /* We really want to load the binary constant 01 00 00 ... 00, 
    *   (There must be an even number of zeros to the right
    *    of the single one bit, and the one bit as far to the
    *    left as can be done within that constraint)
    * but we don't portably know the word size.  But it must
    * be at least as long as the argument we were handed.
    * So we will compute that size.  If you know your word size,
    * Then replace this loop with a load of squaredbit with
    * the constant.
    */
   remainder = x>>2;  squaredbit = 1;
   while (remainder >   0) {
          remainder >>= 2;
          squaredbit <<= 2;
   }
   /* Form bits of the answer. */
   remainder = x;  root = 0;
   while (squaredbit > 0) {
     if (remainder >= (squaredbit | root)) {
         remainder -= (squaredbit | root);
         root >>= 1; root |= squaredbit;
     } else {
         root >>= 1;
     }
     squaredbit >>= 2; 
   }
   return root;
}


And here is "Halleck and Legalize's improved version":

/* Integer square root */
long isqrt (x) long x;{
  long   squaredbit, remainder, root;
   if (x<1) return 0;
  
   /* Load the binary constant 01 00 00 ... 00, where the number
    * of zero bits to the right of the single one bit
    * is even, and the one bit is as far left as is consistant
    * with that condition.)
    */
   squaredbit  = (long) ((((unsigned long) ~0L) >> 1) & 
                        ~(((unsigned long) ~0L) >> 2));
   /* This portable load replaces the loop that used to be 
    * here, and was donated by  legalize@xmission.com 
    */
   /* Form bits of the answer. */
   remainder = x;  root = 0;
   while (squaredbit > 0) {
     if (remainder >= (squaredbit | root)) {
         remainder -= (squaredbit | root);
         root >>= 1; root |= squaredbit;
     } else {
         root >>= 1;
     }
     squaredbit >>= 2; 
   }
  return root;
}}


All the above was merrily plagiarised from
http://www.cc.utah.edu/~nahaj/factoring/ which probably only goes to
show how wonderful the Web is. I also have absolutely no idea about
whether any of the above plagiarised code works, which lack of
knowledge (and interest!) is totally unrelated to the time in this
locality.


Gil Jaysmith
SN Systems Software Ltd, makers of Psy-Q... http://www.snsys.com
Disclaimer: What I say when I post here represents me, not my employers.