home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / SQRT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1987-06-20  |  2.6 KB  |  157 lines
  1. /*                            sqrt.c
  2.  *
  3.  *    Square root
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, sqrt();
  10.  *
  11.  * y = sqrt( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the square root of x.
  18.  *
  19.  * Range reduction involves isolating the power of two of the
  20.  * argument and using a polynomial approximation to obtain
  21.  * a rough value for the square root.  Then Heron's iteration
  22.  * is used three times to converge to an accurate value.
  23.  *
  24.  *
  25.  *
  26.  * ACCURACY:
  27.  *
  28.  *
  29.  *                      Relative error:
  30.  * arithmetic   range      # trials      peak         rms
  31.  *    DEC       0, 30        2000       2.1e-17     5.2e-18
  32.  *    IEEE      0,1.7e308   15000       1.6e-16     5.9e-17
  33.  *
  34.  *
  35.  * ERROR MESSAGES:
  36.  *
  37.  *   message         condition      value returned
  38.  * sqrt domain        x < 0            0.0
  39.  *
  40.  */
  41.  
  42. /*
  43. Cephes Math Library Release 2.0:  April, 1987
  44. Copyright 1984, 1987 by Stephen L. Moshier
  45. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  46. */
  47.  
  48. #include "mconf.h"
  49.  
  50. extern double SQRT2;  /*  SQRT2 = 1.41421356237309504880 */
  51.  
  52. double sqrt(x)
  53. double x;
  54. {
  55. short e, *q;
  56. double z, w;
  57. #ifdef UNK
  58. double frexp(), ldexp();
  59. #endif
  60.  
  61. if( x <= 0.0 )
  62.     {
  63.     if( x < 0.0 )
  64.         mtherr( "sqrt", DOMAIN );
  65.     return( 0.0 );
  66.     }
  67. w = x;
  68. /* separate exponent and significand */
  69. #ifdef UNK
  70. z = frexp( x, &e );
  71. #endif
  72. #ifdef DEC
  73. q = (short *)&x;
  74. e = ((*q >> 7) & 0377) - 0200;
  75. *q &= 0177;
  76. *q |= 040000;
  77. z = x;
  78. #endif
  79. #ifdef IBMPC
  80. q = (short *)&x;
  81. q += 3;
  82. e = ((*q >> 4) & 0x0fff) - 0x3fe;
  83. *q &= 0x000f;
  84. *q |= 0x3fe0;
  85. z = x;
  86. #endif
  87. #ifdef MIEEE
  88. q = (short *)&x;
  89. e = ((*q >> 4) & 0x0fff) - 0x3fe;
  90. *q &= 0x000f;
  91. *q |= 0x3fe0;
  92. z = x;
  93. #endif
  94.  
  95. /* approximate square root of number between 0.5 and 1
  96.  * relative error of approximation = 7.47e-3
  97.  */
  98. x = 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z;
  99.  
  100. /* adjust for odd powers of 2 */
  101. if( (e & 1) != 0 )
  102.     x *= SQRT2;
  103.  
  104. /* re-insert exponent */
  105. #ifdef UNK
  106. x = ldexp( x, (e >> 1) );
  107. #endif
  108. #ifdef DEC
  109. *q += ((e >> 1) & 0377) << 7;
  110. *q &= 077777;
  111. #endif
  112. #ifdef IBMPC
  113. *q += ((e >>1) & 0x7ff) << 4;
  114. *q &= 077777;
  115. #endif
  116. #ifdef MIEEE
  117. *q += ((e >>1) & 0x7ff) << 4;
  118. *q &= 077777;
  119. #endif
  120.  
  121. /* Newton iterations: */
  122. #ifdef UNK
  123. x += w/x;
  124. x = ldexp( x, -1 );    /* divide by 2 */
  125. x += w/x;
  126. x = ldexp( x, -1 );
  127. x += w/x;
  128. x = ldexp( x, -1 );
  129. #endif
  130. #ifdef DEC
  131. x += w/x;
  132. *q -= 0200;
  133. x += w/x;
  134. *q -= 0200;
  135. x += w/x;
  136. *q -= 0200;
  137. #endif
  138. #ifdef IBMPC
  139. x += w/x;
  140. *q -= 0x10;
  141. x += w/x;
  142. *q -= 0x10;
  143. x += w/x;
  144. *q -= 0x10;
  145. #endif
  146. #ifdef MIEEE
  147. x += w/x;
  148. *q -= 0x10;
  149. x += w/x;
  150. *q -= 0x10;
  151. x += w/x;
  152. *q -= 0x10;
  153. #endif
  154.  
  155. return(x);
  156. }
  157.