home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ieee / esqrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.7 KB  |  92 lines
  1. /*    esqrt.c        */
  2. /* Square root check routine, done by long division. */
  3. /* by Stephen L. Moshier. */
  4.  
  5.  
  6. #include "mconf.h"
  7. #include "ehead.h"
  8.  
  9. void emovi(), emovo(), ecleaz(), emdnorm(), eaddm(), esubm();
  10. int ecmpm();
  11. static unsigned short rndbit[NI] = {0,0,0,0,0,0,0,1,0};
  12.  
  13. void esqrt( x, y )
  14. short *x, *y;
  15. {
  16. unsigned short temp[NI], num[NI], sq[NI], xx[NI];
  17. int i, j, k, n;
  18. long m, exp;
  19.  
  20. /* Check for arg <= 0 */
  21. i = ecmp( x, ezero );
  22. if( i <= 0 )
  23.     {
  24.     eclear(y);
  25.     if( i < 0 )
  26.         mtherr( "esqrt", DOMAIN );
  27.     return;
  28.     }
  29.  
  30. /* Bring in the arg and renormalize if it is denormal. */
  31. emovi( x, xx );
  32. m = (long )xx[1]; /* local long word exponent */
  33. if( m == 0 )
  34.     m -= enormlz( xx );
  35.  
  36. /* Divide exponent by 2 */
  37. m -= 0x3ffe;
  38. exp = (unsigned short )( (m / 2) + 0x3ffe );
  39.  
  40. /* Adjust if exponent odd */
  41. if( (m & 1) != 0 )
  42.     {
  43.     if( m > 0 )
  44.         exp += 1;
  45.     eshdn1( xx );
  46.     }
  47.  
  48. ecleaz( sq );
  49. ecleaz( num );
  50. n = 8; /* get 8 bits of result per inner loop */
  51.  
  52. for( j=0; j<(2*NE-2); j++ )
  53.     {
  54. /* bring in next word of arg */
  55.     if( j < NE )
  56.         num[NI-1] = xx[j+3];
  57. /* Do additional bit on last outer loop, for roundoff. */
  58.     if( j == (2*NE-3) )
  59.         n = 9;
  60.     for( i=0; i<n; i++ )
  61.         {
  62. /* Next 2 bits of arg */
  63.         eshup1( num );
  64.         eshup1( num );
  65. /* Shift up answer */
  66.         eshup1( sq );
  67. /* Make trial divisor */
  68.         for( k=0; k<NI; k++ )
  69.             temp[k] = sq[k];
  70.         eshup1( temp );
  71.         eaddm( rndbit, temp );
  72. /* Subtract and insert answer bit if it goes in */
  73.         if( ecmpm( temp, num ) <= 0 )
  74.             {
  75.             esubm( temp, num );
  76.             sq[NI-2] |= 1;
  77.             }
  78.         }
  79.     }
  80.  
  81. exp -= 1; /* Adjust for extra, roundoff loop done. */
  82.  
  83. /* Sticky bit = 1 if the remainder is nonzero. */
  84. k = 0;
  85. for( i=3; i<NI; i++ )
  86.     k |= (int )num[i];
  87.  
  88. /* Renormalize and round off. */
  89. emdnorm( sq, k, 0, exp, 64 );
  90. emovo( sq, y );
  91. }
  92.