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