home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / os / os2 / programm / 7351 < prev    next >
Encoding:
Internet Message Format  |  1993-01-06  |  2.8 KB
  1. Xref: sparky comp.os.os2.programmer:7351 comp.graphics:13442 comp.os.msdos.programmer:11825 comp.sys.ibm.pc.programmer:745
  2. Newsgroups: comp.os.os2.programmer,comp.graphics,comp.os.msdos.programmer,comp.sys.ibm.pc.programmer
  3. Path: sparky!uunet!munnari.oz.au!manuel.anu.edu.au!csc.canberra.edu.au!echo!eyal
  4. From: eyal@echo.canberra.edu.au (Eyal Lebedinsky)
  5. Subject: Re: Fixed point sqrt anyone?
  6. Message-ID: <eyal.726317041@ise>
  7. Sender: news@csc.canberra.edu.au
  8. Organization: Info Sci & Eng, University of Canberra, AUSTRALIA
  9. References: <1idsq6INNced@flop.ENGR.ORST.EDU> <1993Jan6.093937.7848@canon.co.uk>
  10. Date:  6 Jan 93 10:44:01 GMT
  11. Lines: 108
  12.  
  13. In <1993Jan6.093937.7848@canon.co.uk> ads@canon.co.uk (Adam Billyard) writes:
  14.  
  15. >murrayk@prism.CS.ORST.EDU (the Dodger) writes:
  16.  
  17.  
  18. >> Well? Does anyone know of a quick and dirty way to get a reasonably accurate
  19. >> square root of a fixed point integer number? My numbers are stored in a 16.16
  20.  
  21. >It depends on what "reasonably accurate" means.  I used a Newton-Raphson
  22. >approximation in a 3D pool game I did a while back, and I got enough accuracy 
  23. >and quickly enough for my needs..
  24.  
  25. >    Adam
  26. >    
  27. >NB If you're looking to do vector normalization then you probably need quite a 
  28. >lot of accuracy.
  29. Here is how I do it. This code is from a flight simulator. I use 16 bit
  30. ints but need to take sqrt of 32bit numbers (the result of squaring the
  31. ints).
  32.  
  33. The program uses a table for the starting value. In practice I write the
  34. tables to a file which I later include into the source.
  35.  
  36. #define FASTCALL _fastcall        /* microsoft style */
  37.  
  38. typedef unsigned char    Uchar;
  39. typedef unsigned int    Uint;
  40. typedef unsigned short    Ushort;
  41. typedef unsigned long    Ulong;
  42.  
  43. static Uint    far sqrtab0[256] = {0};
  44. static Uint    far sqrtab1[256] = {0};
  45. static Uint    far sqrtab2[256] = {0};
  46. static Uint    far sqrtab3[256] = {0};
  47.  
  48. static int far
  49. lsqrt (Ulong x)                /* used for initialization only */
  50. {
  51.     long    e;
  52.     Ulong    r, t;
  53.  
  54.     if (x & 0xffff0000L)
  55.         r = 662 + x / 17916;
  56.     else if (x & 0x0000ff00L)
  57.         r = 3 + x / 70;
  58.     else
  59.         r = 2 + x / 11;
  60.  
  61.     do {
  62.         t = x / r;
  63.         e = (long)(r - t) / 2;
  64.         r = (r + t) / 2;
  65.     } while (e);
  66.     return ((int)r);
  67. }
  68.  
  69. static void far
  70. set_lsqrt (void)            /* initialise tables */
  71. {
  72.     int    i;
  73.  
  74.     for (i = 0; i < 256; ++i) {
  75.         sqrtab3[i] = lsqrt (i*256L*256L*256L);
  76.         sqrtab2[i] = sqrtab3[i] >> 4;
  77.         sqrtab1[i] = sqrtab2[i] >> 4;
  78.         sqrtab0[i] = sqrtab1[i] >> 4;
  79.     }
  80. }
  81.  
  82. static Uint far FASTCALL
  83. my_sqrt (Ulong x)
  84. {
  85.     register Uint    r, t;
  86.     register int    e, ind;
  87.  
  88.     ind = 0;
  89.     if (x & 0xff000000L) {
  90.         if (x & 0xc0000000L) {
  91.             ind = 1;
  92.             x /= 4;
  93.         }
  94.         r = sqrtab3[x>>24];
  95.     } else if (x & 0x7fff0000L)
  96.         r = sqrtab2[x>>16];
  97.     else if (x & 0x007fff00L)
  98.         r = sqrtab1[x>>8];
  99.     else
  100.         return (sqrtab0[x]);
  101.  
  102.     do {
  103.         t = (Uint)(x / r);
  104.         e = (int)(r - t) / 2;
  105.         r -= e;
  106.     } while (e);
  107.  
  108.     if (ind)
  109.         return (t+r);
  110.     else
  111.         return ((t+r)/2);
  112. }
  113.  
  114. main ()
  115. {
  116.     set_lsqrt ();
  117.  
  118.     return (0);
  119. }
  120.  
  121.