home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsi / squareroot.c < prev    next >
C/C++ Source or Header  |  1992-09-15  |  2KB  |  90 lines
  1. /*
  2.  * A High Speed, Low Precision Square Root
  3.  * by Paul Lalonde and Robert Dawson
  4.  * from "Graphics Gems", Academic Press, 1990
  5.  */
  6.  
  7.  
  8. /*
  9.  * SPARC implementation of a fast square root by table 
  10.  * lookup.
  11.  * SPARC floating point format is as follows:
  12.  *
  13.  * BIT 31     30     23     22     0
  14.  *     sign    exponent    mantissa
  15.  */
  16.  
  17. #include <math.h>
  18.  
  19. static short sqrttab[0x100];    /* declare table of square roots */
  20.  
  21. void
  22. build_table()
  23. {
  24.     unsigned short i;
  25.     float f;
  26.     unsigned int *fi = (unsigned int *)&f;
  27.                 /* to access the bits of a float in  */
  28.                 /* C quickly we must misuse pointers */
  29.  
  30.     for(i=0; i<= 0x7f; i++) {
  31.         *fi = 0;
  32.  
  33.         /*
  34.          * Build a float with the bit pattern i as mantissa
  35.          * and an exponent of 0, stored as 127
  36.            */
  37.  
  38.         *fi = (i << 16) | (127 << 23);
  39.         f = sqrt(f);
  40.  
  41.         /*
  42.          * Take the square root then strip the first 7 bits of
  43.          * the mantissa into the table
  44.          */
  45.  
  46.         sqrttab[i] = (*fi & 0x7fffff) >> 16;
  47.  
  48.         /*
  49.          * Repeat the process, this time with an exponent of
  50.          * 1, stored as 128
  51.          */
  52.  
  53.         *fi = 0;
  54.         *fi = (i << 16) | (128 << 23);
  55.         f = sqrt(f);
  56.         sqrttab[i+0x80] = (*fi & 0x7fffff) >> 16;
  57.     }
  58. }
  59.  
  60. /*
  61.  * fsqrt - fast square root by table lookup
  62.  */
  63. float
  64. fsqrt(float n)
  65. {
  66.     unsigned int *num = (unsigned int *)&n;
  67.                 /* to access the bits of a float in C
  68.                  * we must misuse pointers */
  69.                             
  70.     short e;        /* the exponent */
  71.     if (n == 0) return (0);    /* check for square root of 0 */
  72.     e = (*num >> 23) - 127;    /* get the exponent - on a SPARC the */
  73.                 /* exponent is stored with 127 added */
  74.     *num &= 0x7fffff;    /* leave only the mantissa */
  75.     if (e & 0x01) *num |= 0x800000;
  76.                 /* the exponent is odd so we have to */
  77.                 /* look it up in the second half of  */
  78.                 /* the lookup table, so we set the   */
  79.                 /* high bit                   */
  80.     e >>= 1;        /* divide the exponent by two */
  81.                 /* note that in C the shift */
  82.                 /* operators are sign preserving */
  83.                 /* for signed operands */
  84.     /* Do the table lookup, based on the quaternary mantissa,
  85.      * then reconstruct the result back into a float
  86.      */
  87.     *num = ((sqrttab[*num >> 16]) << 16) | ((e + 127) << 23);
  88.     return(n);
  89. }
  90.