home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #3 / NN_1993_3.iso / spool / comp / lang / c / 20342 < prev    next >
Encoding:
Text File  |  1993-01-28  |  3.0 KB  |  139 lines
  1. Newsgroups: comp.lang.c
  2. Path: sparky!uunet!ferkel.ucsb.edu!taco!gatech!rpi!ghost.dsi.unimi.it!univ-lyon1.fr!scsing.switch.ch!dxcern!marmac2.in2p3.fr!user
  3. From: Meessen@marvax.IN2P3.fr (Meessen Christophe)
  4. Subject: Fixed point sqrt
  5. Message-ID: <Meessen-280193094030@marmac2.in2p3.fr>
  6. Followup-To: comp.lang.c
  7. Sender: news@dxcern.cern.ch (USENET News System)
  8. Organization: Centre de Physique des Particules de Marseille
  9. Date: Thu, 28 Jan 1993 08:35:23 GMT
  10. Lines: 127
  11.  
  12. In article <VJ.93Jan19115116@dlsi.dl.ac.uk>, vj@gserv1.dl.ac.uk (Dave
  13. Hines) wrote:
  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?
  17. > 
  18. > Hmmm - I have a long integer square root routine. 
  19.  
  20. Dave Hines proposed a sqrt function receiving an unsigned long (X) and
  21. returning
  22. an unsigned long. The remainder, an unsigned long, is accessible as
  23. 'result2'.
  24. The algorithm uses 4 unsigned long variables. 
  25. Thus it's a long -> long function.
  26.  
  27.  
  28. I suggest two algorithms using the same variables. The first one will
  29. compute a sqrt of a long, returning a fixed point value. The second will
  30. compute the sqrt of a fixed point value, returning a fixed point value.
  31. Thus it's a long -> fixed or fixed -> fixed function.
  32.  
  33. We will suppose a fixed is equivalent to a long, and longs are 32 bits
  34. wide.
  35.  
  36. /*
  37.  * fixed sqrtL2F( long X );
  38.  *
  39.  * Long to fixed point square roots. 
  40.  * RETURNS the fixed point square root of X (long).
  41.  * REQUIRES X is positive.
  42.  *
  43.  * Christophe MEESSEN, 1993.
  44.  */
  45.  
  46. fixed sqrtL2F( long X ) {
  47.  
  48.   unsigned long t, q, b, r;
  49.  
  50.   r = X;
  51.   b = 0x40000000;
  52.   q = 0;
  53.  
  54.   while( b > 0 ) {
  55.     t = q + b;
  56.     if( r >= t ) {
  57.       r = r - t;
  58.       q = t + b;
  59.     }
  60.     r = r * 2;     /* shift left  1 bit */
  61.     b = b / 2;     /* shift right 1 bit */
  62.   }
  63.   if( r >= q )
  64.     q = q + 1;
  65.  
  66.   return( q );
  67. }
  68.  
  69.  
  70. /*
  71.  * fixed sqrtF2F( fixed X );
  72.  *
  73.  * Fixed to fixed point square roots. 
  74.  * RETURNS the fixed point square root of X (fixed).
  75.  * REQUIRES X is positive.
  76.  *
  77.  * Christophe MEESSEN, 1993.
  78.  */
  79.  
  80. fixed sqrtF2F ( fixed X ) {
  81.  
  82.   unsigned long t, q, b, r;
  83.  
  84.   r = X;
  85.   b = 0x40000000;
  86.   q = 0;
  87.  
  88.   while( b >= 256 ) {
  89.     t = q + b;
  90.     if( r >= t ) {
  91.       r = r - t;
  92.       q = t + b;
  93.     }
  94.     r = r * 2;     /* shift left  1 bit */
  95.     b = b / 2;     /* shift right 1 bit */
  96.   }
  97.   q = q / 256;     /* shift right 8 bits */
  98.  
  99.   return( q );
  100. }
  101.  
  102. for memory I just rewrote the long to long sqrt function.
  103.  
  104. /*
  105.  * long sqrtL2L( long X );
  106.  *
  107.  * Long to long point square roots. 
  108.  * RETURNS the integer square root of X (long).
  109.  * REMAINDER is in the local variable r of type long on return.  
  110.  * REQUIRES X is positive.
  111.  *
  112.  * Christophe MEESSEN, 1993.
  113.  */
  114.  
  115. long sqrtL2L( long X ) {
  116.  
  117.   unsigned long t, q, b, r;
  118.  
  119.   r = X;
  120.   b = 0x40000000;
  121.   q = 0;
  122.  
  123.   while( b >= 256 ) {
  124.     t = q + b;
  125.     q = q / 2;     /* shift right 1 bit */
  126.     if( r >= t ) {
  127.       r = r - t;
  128.       q = q + b;
  129.     }
  130.     b = b / 4;     /* shift right 2 bits */
  131.   }
  132.  
  133.   return( q );
  134. }
  135.  
  136. This functions have been tested. But I'll be pleased to receive any
  137. comments.
  138. Please respond through mail: Meessen@marvax.in2p3.fr
  139.