home *** CD-ROM | disk | FTP | other *** search
/ Tricks of the Mac Game Programming Gurus / TricksOfTheMacGameProgrammingGurus.iso / More Source / C⁄C++ / Peter's Final Project / src / fixed.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-05-10  |  4.4 KB  |  132 lines  |  [TEXT/KAHL]
  1. /*
  2.  *  Peter's Final Project -- A texture mapping demonstration
  3.  *  © 1995, Peter Mattis
  4.  *
  5.  *  E-mail:
  6.  *  petm@soda.csua.berkeley.edu
  7.  *
  8.  *  Snail-mail:
  9.  *   Peter Mattis
  10.  *   557 Fort Laramie Dr.
  11.  *   Sunnyvale, CA 94087
  12.  *
  13.  *  Avaible from:
  14.  *  http://www.csua.berkeley.edu/~petm/final.html
  15.  *
  16.  *  This program is free software; you can redistribute it and/or modify
  17.  *  it under the terms of the GNU General Public License as published by
  18.  *  the Free Software Foundation; either version 2 of the License, or
  19.  *  (at your option) any later version.
  20.  *
  21.  *  This program is distributed in the hope that it will be useful,
  22.  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  23.  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  24.  *  GNU General Public License for more details.
  25.  *
  26.  *  You should have received a copy of the GNU General Public License
  27.  *  along with this program; if not, write to the Free Software
  28.  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  29.  */
  30.  
  31. #include "fixed.h"
  32.  
  33. /*
  34.  * This is Ron Hunsinger's integer square root algorithm adapted for fixed
  35.  * point numbers. This adaption is simple. Since a fixed point number can
  36.  * be represented as p*q, where p is some constant value (65536) in this
  37.  * case. Then the sqrt(q) represented as a fixed point number is p*sqrt(q).
  38.  * The algorithm below will spit out sqrt(p*q), so all that is needed is to
  39.  * multiply by sqrt(p). Since sqrt(p) is constant (256) and is a multiple of
  40.  * 2, this operation can be done with a left shift of 8 bits.
  41.  * 
  42.  * -Peter Mattis
  43.  */
  44.  
  45. FIX_NUM fix_sqrt(FIX_NUM n)
  46. {
  47. #define lsqrt_max4pow (1UL << 30)
  48.     /* lsqrt_max4pow is the (machine-specific) largest power of 4 that can
  49.      * be represented in an unsigned long.
  50.      *
  51.      * Compute the integer square root of the integer argument n
  52.      * Method is to divide n by x computing the quotient x and remainder r
  53.      * Notice that the divisor x is changing as the quotient x changes
  54.      * 
  55.      * Instead of shifting the dividend/remainder left, we shift the
  56.      * quotient/divisor right. The binary point starts at the extreme
  57.      * left, and shifts two bits at a time to the extreme right.
  58.      * 
  59.      * The residue contains n-x^2. (Within these comments, the ^ operator
  60.      * signifies exponentiation rather than exclusive or. Also, the /
  61.      * operator returns fractions, rather than truncating, so 1/4 means
  62.      * one fourth, not zero.)
  63.      * 
  64.      * Since (x + 1/2)^2 == x^2 + x + 1/4,
  65.      * n - (x + 1/2)^2 == (n - x^2) - (x + 1/4)
  66.      * Thus, we can increase x by 1/2 if we decrease (n-x^2) by (x+1/4)
  67.      */
  68.      
  69.     unsigned long residue;    /* n - x^2  */
  70.     unsigned long root;       /* x + 1/4  */
  71.     unsigned long half;       /* 1/2      */
  72.  
  73.     residue = n;              /* n - (x = 0)^2, with suitable alignment */
  74.  
  75.     /*
  76.      * if the correct answer fits in two bits, pull it out of a magic hat
  77.      */
  78.     if (residue <= 12)
  79.         return (0x03FFEA94 >> (residue *= 2)) & 3;
  80.  
  81.     root = lsqrt_max4pow; /* x + 1/4, shifted all the way left */
  82.     /* half = root + root; 1/2, shifted likewise */
  83.  
  84.     /* 
  85.      * Unwind iterations corresponding to leading zero bits 
  86.      */
  87.     while (root > residue) root >>= 2;
  88.  
  89.     /*
  90.      * Unwind the iteration corresponding to the first one bit
  91.      * Operations have been rearranged and combined for efficiency
  92.      * Initialization of half is folded into this iteration
  93.      */
  94.     residue -= root;    /* Decrease (n-x^2) by (0+1/4)             */
  95.     half = root >> 2;   /* 1/4, with binary point shifted right 2  */
  96.     root += half;       /* x=1. (root is now (x=1)+1/4.)           */
  97.     half += half;       /* 1/2, properly aligned                   */
  98.  
  99.     /*
  100.      * Normal loop (there is at least one iteration remaining)
  101.      */
  102.     do {
  103.         if (root <= residue)  /* Whenever we can,                          */
  104.         {
  105.             residue -= root;  /* decrease (n-x^2) by (x+1/4)               */
  106.             root += half;     /* increase x by 1/2                         */
  107.         }
  108.         half >>= 2;          /* Shift binary point 2 places right          */
  109.         root -= half;        /* x{ +1/2 } +1/4 - 1/8 == x { +1/2 } 1/8     */
  110.         root >>= 1;          /* 2x{ +1 } +1/4, shifted right 2 places      */
  111.     } while (half);          /* When 1/2 == 0, bin. point is at far right  */
  112.  
  113.     /* 
  114.      * round up if (x+1/2)^2 < n
  115.      */
  116.     if (root < residue)
  117.         ++root;
  118.  
  119.     /*
  120.      * Multiply the result by 256. (A left shift by 8 will do
  121.      *  the job quite quickly).
  122.      *
  123.      * -Peter Mattis
  124.      */
  125.     root <<= 8;
  126.  
  127.     /* 
  128.      * Guaranteed to be correctly rounded (or truncated)
  129.      */
  130.         return root;
  131. }
  132.