home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / snip9707.zip / ISQRT.C < prev    next >
C/C++ Source or Header  |  1997-07-05  |  3KB  |  90 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. #include <string.h>
  4. #include "snipmath.h"
  5.  
  6. #define BITSPERLONG 32
  7.  
  8. #define TOP2BITS(x) ((x & (3L << (BITSPERLONG-2))) >> (BITSPERLONG-2))
  9.  
  10.  
  11. /* usqrt:
  12.     ENTRY x: unsigned long
  13.     EXIT  returns floor(sqrt(x) * pow(2, BITSPERLONG/2))
  14.  
  15.     Since the square root never uses more than half the bits
  16.     of the input, we use the other half of the bits to contain
  17.     extra bits of precision after the binary point.
  18.  
  19.     EXAMPLE
  20.         suppose BITSPERLONG = 32
  21.         then    usqrt(144) = 786432 = 12 * 65536
  22.                 usqrt(32) = 370727 = 5.66 * 65536
  23.  
  24.     NOTES
  25.         (1) change BITSPERLONG to BITSPERLONG/2 if you do not want
  26.             the answer scaled.  Indeed, if you want n bits of
  27.             precision after the binary point, use BITSPERLONG/2+n.
  28.             The code assumes that BITSPERLONG is even.
  29.         (2) This is really better off being written in assembly.
  30.             The line marked below is really a "arithmetic shift left"
  31.             on the double-long value with r in the upper half
  32.             and x in the lower half.  This operation is typically
  33.             expressible in only one or two assembly instructions.
  34.         (3) Unrolling this loop is probably not a bad idea.
  35.  
  36.     ALGORITHM
  37.         The calculations are the base-two analogue of the square
  38.         root algorithm we all learned in grammar school.  Since we're
  39.         in base 2, there is only one nontrivial trial multiplier.
  40.  
  41.         Notice that absolutely no multiplications or divisions are performed.
  42.         This means it'll be fast on a wide range of processors.
  43. */
  44.  
  45. void usqrt(unsigned long x, struct int_sqrt *q)
  46. {
  47.       unsigned long a = 0L;                   /* accumulator      */
  48.       unsigned long r = 0L;                   /* remainder        */
  49.       unsigned long e = 0L;                   /* trial product    */
  50.  
  51.       int i;
  52.  
  53.       for (i = 0; i < BITSPERLONG; i++)   /* NOTE 1 */
  54.       {
  55.             r = (r << 2) + TOP2BITS(x); x <<= 2; /* NOTE 2 */
  56.             a <<= 1;
  57.             e = (a << 1) + 1;
  58.             if (r >= e)
  59.             {
  60.                   r -= e;
  61.                   a++;
  62.             }
  63.       }
  64.       memcpy(q, &a, sizeof(long));
  65. }
  66.  
  67. #ifdef TEST
  68.  
  69. #include <stdio.h>
  70. #include <stdlib.h>
  71.  
  72. main(void)
  73. {
  74.       int i;
  75.       unsigned long l = 0x3fed0169L;
  76.       struct int_sqrt q;
  77.  
  78.       for (i = 0; i < 101; ++i)
  79.       {
  80.             usqrt(i, &q);
  81.             printf("sqrt(%3d) = %2d, remainder = %2d\n",
  82.                   i, q.sqrt, q.frac);
  83.       }
  84.       usqrt(l, &q);
  85.       printf("\nsqrt(%lX) = %X, remainder = %X\n", l, q.sqrt, q.frac);
  86.       return 0;
  87. }
  88.  
  89. #endif /* TEST */
  90.