home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / qfloat / qsqrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.1 KB  |  66 lines
  1. /*                            qsqrt.c        */
  2. /* square root check routine */
  3. /* full precision for 9 word mantissa */
  4. #include "qhead.h"
  5. extern short qsqrt2[];
  6.  
  7. /* constants for first approximation polynomial */
  8. static short qsq2[NQ] =
  9. {-1,037776,0,0150517,0142057,0163633,0124000,0,0,0,0,0,};
  10.  
  11. static short qsq1[NQ] =
  12. {0,040000,0,0161743,0141256,046024,063400,0,0,0,0,0,};
  13.  
  14. static short qsq0[NQ] =
  15. {0,037777,0,0120213,0156175,0152777,0161400,0,0,0,0,0,};
  16.  
  17. int qsqrt( x, y )
  18. short *x, *y;
  19. {
  20. short a[NQ], b[NQ];
  21. int i, m;
  22.  
  23. if( x[0] < 0 )
  24.     {
  25.     qclear(y);
  26.     printf( "qsqrt domain error\n" );
  27.     return;
  28.     }
  29. if( x[1] == 0 )
  30.     {
  31.     qclear(y);
  32.     return;
  33.     }
  34.  
  35. qmov( x, a );
  36. m = a[1];    /* save the exponent        */
  37. a[1] = 040000;    /* change range to 0.5, 1 */
  38.  
  39. /* b = ( a * qsq2 + qsq1) * a + qsq0        */
  40. qmul( a, qsq2, b );    /* b = a * qsq2        */
  41. qadd( qsq1, b, b );    /* b += qsq1        */
  42. qmul( a, b, b );    /* b *= a;        */
  43. qadd( qsq0, b, b );    /* b += qsq0;        */
  44.  
  45. /* divide exponent by 2 */
  46. m -= 040000;
  47. b[1] = (m / 2) + 040000;
  48.  
  49.  
  50. /* multiply by sqrt(2) if exponent odd */
  51. if( (m & 1) != 0 )
  52.     qmul( b, qsqrt2, b );
  53.  
  54.  
  55. /* Newton iterations */
  56.  
  57. for( i=0; i<8; i++ )
  58.     {
  59.     qdiv( b, x, a );
  60.     qadd( a, b, b );
  61.     b[1] -= 1;
  62.     }
  63.  
  64. qmov( b, y );
  65. }
  66.