home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsiii / sqrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-17  |  2.7 KB  |  98 lines
  1. /* fsqrt.c
  2.  *
  3.  * A fast square root program adapted from the code of
  4.  * Paul Lalonde and Robert Dawson in Graphics Gems I.
  5.  * The format of IEEE double precision floating point numbers is:
  6.  *
  7.  * SEEEEEEEEEEEMMMM MMMMMMMMMMMMMMMM MMMMMMMMMMMMMMMM MMMMMMMMMMMMMMMM
  8.  *
  9.  * S = Sign bit for whole number
  10.  * E = Exponent bit (exponent in excess 1023 form)
  11.  * M = Mantissa bit
  12.  */
  13.  
  14. #include <stdio.h>
  15. #include <math.h>
  16.  
  17. /* MOST_SIG_OFFSET gives the (int *) offset from the address of the double
  18.  * to the part of the number containing the sign and exponent.
  19.  * You will need to find the relevant offset for your architecture.
  20.  */
  21.  
  22. #define MOST_SIG_OFFSET 1
  23.  
  24. /* SQRT_TAB_SIZE - the size of the lookup table - must be a power of four.
  25.  */
  26.  
  27. #define SQRT_TAB_SIZE 16384
  28.  
  29. /* MANT_SHIFTS is the number of shifts to move mantissa into position.
  30.  * If you quadruple the table size subtract two from this constant,
  31.  * if you quarter the table size then add two.
  32.  * Valid values are: (16384, 7) (4096, 9) (1024, 11) (256, 13)
  33.  */
  34.  
  35. #define MANT_SHIFTS   7
  36.  
  37. #define EXP_BIAS   1023       /* Exponents are always positive     */
  38. #define EXP_SHIFTS 20         /* Shifs exponent to least sig. bits */
  39. #define EXP_LSB    0x00100000 /* 1 << EXP_SHIFTS                   */
  40. #define MANT_MASK  0x000FFFFF /* Mask to extract mantissa          */
  41.  
  42. int        sqrt_tab[SQRT_TAB_SIZE];
  43.  
  44. void
  45. init_sqrt_tab()
  46. {
  47.         int           i;
  48.         double        f;
  49.         unsigned int  *fi = (unsigned int *) &f + MOST_SIG_OFFSET;
  50.         
  51.         for (i = 0; i < SQRT_TAB_SIZE/2; i++)
  52.         {
  53.                 f = 0; /* Clears least sig part */
  54.                 *fi = (i << MANT_SHIFTS) | (EXP_BIAS << EXP_SHIFTS);
  55.                 f = sqrt(f);
  56.                 sqrt_tab[i] = *fi & MANT_MASK;
  57.  
  58.                 f = 0; /* Clears least sig part */
  59.                 *fi = (i << MANT_SHIFTS) | ((EXP_BIAS + 1) << EXP_SHIFTS);
  60.                 f = sqrt(f);
  61.                 sqrt_tab[i + SQRT_TAB_SIZE/2] = *fi & MANT_MASK;
  62.         }
  63. }
  64.  
  65. double
  66. fsqrt(f)
  67. double f;
  68. {
  69.         unsigned int e;
  70.         unsigned int   *fi = (unsigned int *) &f + MOST_SIG_OFFSET;
  71.  
  72.         if (f == 0.0) return(0.0);
  73.         e = (*fi >> EXP_SHIFTS) - EXP_BIAS;
  74.         *fi &= MANT_MASK;
  75.         if (e & 1)
  76.                 *fi |= EXP_LSB;
  77.         e >>= 1;
  78.         *fi = (sqrt_tab[*fi >> MANT_SHIFTS]) |
  79.               ((e + EXP_BIAS) << EXP_SHIFTS);
  80.         return(f);
  81. }
  82.  
  83. void
  84. dump_sqrt_tab()
  85. {
  86.         int        i, nl = 0;
  87.  
  88.         printf("unsigned int sqrt_tab[] = {\n");
  89.         for (i = 0; i < SQRT_TAB_SIZE-1; i++)
  90.         {
  91.                 printf("0x%x,", sqrt_tab[i]);
  92.                 nl++;
  93.                 if (nl > 8) { nl = 0; putchar('\n'); }
  94.         }
  95.         printf("0x%x\n", sqrt_tab[SQRT_TAB_SIZE-1]);
  96.         printf("};\n");
  97. }
  98.