[INFO] Yet Another sqrt algorithm.

[dohzono]

I've read gmp-2.0.2 in FrreeBSD-4.0, and found that generic-sqrt is
using Newton's method (PROJECTS says it will be a Newton-Raphson
method. Have you done it? Anyway, just read this as an information).

The cost of the algorithm described below is:

  bits_of_result * [1 shift + 1 sub (or add)] 

and nearly equal to 1 division's cost (means faster than Newton's, but
may be slower than NR). I've heard about 10 years ago that this method
had been used in some Pascal libraries, and wrote this sample for my
own multi-preceision library.

SAMPLE CODE DESCRIPTION:

Capitalized variables are multi-preceision variables. 1 bit and/or to
MP vars are easy operations, as you know.

Usually this kind of algorithm goes like:

  if (XX <= X)
    {
      X -= XX;
      ...;

but compare and sub are almost equal costs (especially for MP vars),
so it always tries to sub, and if the result became negative, runs
add-loop (This doesn't mean going backward. Please read the code
carefully). This sub-add method can also be used in MP division-loop, 
when you use a divide algorithm like this.

Any questions and/or suggestions are welcome, but please e-mail me
directly because I'm not a member of the ML.

Thank you for reading my poor English.

/* assume that 0 <= X < 0x40000000 */
int
isqrt (int X)
{
  unsigned y1 = 0x40000000, y2 = 0x80000000;
  unsigned yr = 0x8000;
  int XX = 0, RET = 0;

  for (;;)
    {
      XX |= y1;
      if (0 <= (X -= XX))
        {
          XX |= y2;
          RET |= yr;
        }
      else
        {
          for (;;)
            {
              int y3 = y1;
              X <<= 1;
              y1 >>= 1;
              yr >>= 1;
              XX |= y1;
              if (!yr)
                {
#ifdef NEAREST
                  if (0 <= (X += XX))
                    RET++;
#endif
                  return RET;
                }
              y2 >>= 1;
              if (0 <= (X += XX))
                {
                  XX &= ~y3;
                  XX |= y2;
                  RET |= yr;
                  break;
                }
              XX &= ~y3;
            }
        }
      XX &= ~y1;
      X <<= 1;
      y1 >>= 1;
      yr >>= 1;
      if (!yr)
        {
#ifdef NEAREST
          XX |= y1;
          if (0 <= (X -= XX))
            RET++;
#endif
          return RET;
        }
      y2 >>= 1;
    }
}