home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / os / os2 / programm / 7424 < prev    next >
Encoding:
Internet Message Format  |  1993-01-08  |  2.2 KB
  1. Xref: sparky comp.os.os2.programmer:7424 comp.graphics:13562 comp.os.msdos.programmer:11913 comp.sys.ibm.pc.programmer:755
  2. Newsgroups: comp.os.os2.programmer,comp.graphics,comp.os.msdos.programmer,comp.sys.ibm.pc.programmer
  3. Path: sparky!uunet!usc!cs.utexas.edu!qt.cs.utexas.edu!yale.edu!ira.uka.de!Sirius.dfn.de!olymp!rs1.thch.uni-bonn.de!greve
  4. From: greve@rs1.thch.uni-bonn.de (Thomas Greve)
  5. Subject: Re: Fixed point sqrt anyone?
  6. Message-ID: <1993Jan8.134438.5422@olymp.informatik.uni-bonn.de>
  7. Sender: usenet@olymp.informatik.uni-bonn.de
  8. Organization: University of Bonn, phys. Chemistry
  9. References: <1idsq6INNced@flop.ENGR.ORST.EDU>
  10. Date: Fri, 8 Jan 1993 13:44:38 GMT
  11. Lines: 37
  12.  
  13. In article <1idsq6INNced@flop.ENGR.ORST.EDU> murrayk@prism.CS.ORST.EDU (the Dodger) writes:
  14. >
  15. >Well? Does anyone know of a quick and dirty way to get a reasonably accurate
  16. >square root of a fixed point integer number? My numbers are stored in a 16.16
  17. >bit format. I've heard of a few ways, but I'm open for ideas. I just want to
  18. >keep on avoiding floating point for ethical reasons. 8^)
  19. I have read a method for integers in a computer magazine some time
  20. ago, and it should work for fixed point numbers, also.
  21.  
  22. Here we go: You can interpret the 16.16 number as (32bit int)/(2^16),
  23. true? The result can be up to 8.16 accurate, but this depends on the
  24. upper 16 bits of your fixed point number. I'll show you a fast way for
  25. 8.8 accuracy and leave the remaining 8 bit as an exercise to the
  26. reader ;-)...
  27.  
  28. The way to go is succesive approximation:
  29. root is the (yet to be calculated) root  8, radicant is the radicant
  30. << 16, which is obtained by simply interpreting the fix as long.
  31.  
  32.     for (j = 128, root = 0; j || root*root == radicant; j >>= 1)
  33.       if ((root + j)*(root + j) >= radicant)
  34.         root += j;
  35.     root <<= 8;    /* now you can cast it back to the fix format */
  36.  
  37. If you need performance, you can do some bit optimization by
  38. calculating (root + j)^2 iteratively from the last value:
  39.  
  40.     for (j = 128, jj = 16384, root = 0, rr2 = rr = 0;
  41.              j || rr = radicant; j >>= 1, jj >>= 2, rr2 >>= 1)
  42.       if ((rr1 = rr + rr2 + jj) >= radicant)
  43.         root += j, rr = rr1, /* = root^2 */ rr2 += jj*2 /* = 2*root*j */; 
  44.  
  45. -- 
  46.                 - Thomas
  47.  
  48.    greve@rs1.thch.uni-bonn.de
  49.    unt145@dbnrhrz1
  50.