home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / qfloat / qsqrta.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.8 KB  |  105 lines
  1. /*    qsqrta.c        */
  2. /* Square root check routine, done by long division. */
  3. /* Copyright (C) 1984-1988 by Stephen L. Moshier. */
  4.  
  5.  
  6. #include "qhead.h"
  7. #include "mconf.h"
  8.  
  9. extern unsigned short qzero[], qone[];
  10. static unsigned short rndbit[NQ+1] = {0};
  11. static short rndset = 0;
  12.  
  13. int qsqrt( x, y )
  14. short *x, *y;
  15. {
  16. unsigned short temp[NQ+1], num[NQ+1], sq[NQ+1], xx[NQ+1];
  17. int i, j, k, n;
  18. long m, exp;
  19.  
  20. if( rndset == 0 )
  21.     {
  22.     qclear( rndbit );
  23.     rndbit[NQ] = 0;
  24.     rndbit[NQ-1] = 1;
  25.     rndset = 1;
  26.     }
  27. /* Check for arg <= 0 */
  28. /*i = qcmp( x, qzero );*/
  29. if( x[1] == 0 )
  30.     goto iszero;
  31. if( x[0] != 0 )
  32.     {
  33.     mtherr( "qsqrt", DOMAIN );
  34. iszero:
  35.     qclear(y);
  36.     return;
  37.     }
  38.  
  39. /* Bring in the arg and renormalized if it is denormal. */
  40. qmovz( x, xx );
  41. m = xx[1]; /* local long word exponent */
  42. /*
  43. if( m == 0 )
  44.     m -= normlz( xx );
  45. */
  46. /* Divide exponent by 2 */
  47. m -= 0x4000;
  48. exp = (unsigned short )( (m / 2) + 0x4000 );
  49.  
  50. /* Adjust if exponent odd */
  51. if( (m & 1) != 0 )
  52.     {
  53.     if( m > 0 )
  54.         exp += 1;
  55.     shdn1( xx );
  56.     }
  57.  
  58. qclear( sq );
  59. sq[NQ] = 0;
  60. qclear( num );
  61. num[NQ] = 0;
  62. n = 8;
  63.  
  64. for( j=0; j<2*NQ-6; j++ )
  65.     {
  66. /* bring in next word of arg */
  67.     if( j < NQ-3 )
  68.         num[NQ] = xx[j+3];
  69. /* Do additional bit on last outer loop, for roundoff. */
  70.     if( j == 2*NQ-7 )
  71.         n = 9;
  72.     for( i=0; i<n; i++ )
  73.         {
  74. /* Next 2 bits of arg */
  75.         shup1( num );
  76.         shup1( num );
  77. /* Shift up answer */
  78.         shup1( sq );
  79. /* Make trial divisor */
  80.         for( k=0; k<=NQ; k++ )
  81.             temp[k] = sq[k];
  82.         shup1( temp );
  83.         addm( rndbit, temp );
  84. /* Subtract and insert answer bit if it goes in */
  85.         if( cmpm( temp, num ) <= 0 )
  86.             {
  87.             subm( temp, num );
  88.             sq[NQ-1] |= 1;
  89.             }
  90.         }
  91.     }
  92.  
  93. exp -= 1; /* Adjust for extra, roundoff loop done. */
  94. sq[1] = exp;
  95.  
  96. /* Sticky bit = 1 if the remainder is nonzero. */
  97. k = 0;
  98. for( i=3; i<=NQ; i++ )
  99.     k |= num[i];
  100.  
  101. /* Renormalize and round off. */
  102. mdnorm( sq, k );
  103. qmov( sq, y );
  104. }
  105.