lib/libm/sqrt.c

   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   domain     # trials      peak         rms
  31  *    DEC       0, 10       60000       2.1e-17     7.9e-18
  32  *    IEEE      0,1.7e308   30000       1.7e-16     6.3e-17
  33  *
  34  *
  35  * ERROR MESSAGES:
  36  *
  37  *   message         condition      value returned
  38  * sqrt domain        x < 0            0.0
  39  *
  40  */
  41 \f
  42 /*
  43 Cephes Math Library Release 2.8:  June, 2000
  44 Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
  45 */
  46
  47
  48 #include "mconf.h"
  49 #ifdef ANSIPROT
  50 extern double frexp ( double, int * );
  51 extern double ldexp ( double, int );
  52 #else
  53 double frexp(), ldexp();
  54 #endif
  55 extern double SQRT2;  /*  SQRT2 = 1.41421356237309504880 */
  56
  57 double sqrt(x)
  58 double x;
  59 {
  60 int e;
  61 #ifndef UNK
  62 short *q;
  63 #endif
  64 double z, w;
  65
  66 if( x <= 0.0 )
  67         {
  68         if( x < 0.0 )
  69                 mtherr( "sqrt", DOMAIN );
  70         return( 0.0 );
  71         }
  72 w = x;
  73 /* separate exponent and significand */
  74 #ifdef UNK
  75 z = frexp( x, &e );
  76 #endif
  77 #ifdef DEC
  78 q = (short *)&x;
  79 e = ((*q >> 7) & 0377) - 0200;
  80 *q &= 0177;
  81 *q |= 040000;
  82 z = x;
  83 #endif
  84
  85 /* Note, frexp and ldexp are used in order to
  86  * handle denormal numbers properly.
  87  */
  88 #ifdef IBMPC
  89 z = frexp( x, &e );
  90 q = (short *)&x;
  91 q += 3;
  92 /*
  93 e = ((*q >> 4) & 0x0fff) - 0x3fe;
  94 *q &= 0x000f;
  95 *q |= 0x3fe0;
  96 z = x;
  97 */
  98 #endif
  99 #ifdef MIEEE
 100 z = frexp( x, &e );
 101 q = (short *)&x;
 102 /*
 103 e = ((*q >> 4) & 0x0fff) - 0x3fe;
 104 *q &= 0x000f;
 105 *q |= 0x3fe0;
 106 z = x;
 107 */
 108 #endif
 109
 110 /* approximate square root of number between 0.5 and 1
 111  * relative error of approximation = 7.47e-3
 112  */
 113 x = 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z;
 114
 115 /* adjust for odd powers of 2 */
 116 if( (e & 1) != 0 )
 117         x *= SQRT2;
 118
 119 /* re-insert exponent */
 120 #ifdef UNK
 121 x = ldexp( x, (e >> 1) );
 122 #endif
 123 #ifdef DEC
 124 *q += ((e >> 1) & 0377) << 7;
 125 *q &= 077777;
 126 #endif
 127 #ifdef IBMPC
 128 x = ldexp( x, (e >> 1) );
 129 /*
 130 *q += ((e >>1) & 0x7ff) << 4;
 131 *q &= 077777;
 132 */
 133 #endif
 134 #ifdef MIEEE
 135 x = ldexp( x, (e >> 1) );
 136 /*
 137 *q += ((e >>1) & 0x7ff) << 4;
 138 *q &= 077777;
 139 */
 140 #endif
 141
 142 /* Newton iterations: */
 143 #ifdef UNK
 144 x = 0.5*(x + w/x);
 145 x = 0.5*(x + w/x);
 146 x = 0.5*(x + w/x);
 147 #endif
 148
 149 /* Note, assume the square root cannot be denormal,
 150  * so it is safe to use integer exponent operations here.
 151  */
 152 #ifdef DEC
 153 x += w/x;
 154 *q -= 0200;
 155 x += w/x;
 156 *q -= 0200;
 157 x += w/x;
 158 *q -= 0200;
 159 #endif
 160 #ifdef IBMPC
 161 x += w/x;
 162 *q -= 0x10;
 163 x += w/x;
 164 *q -= 0x10;
 165 x += w/x;
 166 *q -= 0x10;
 167 #endif
 168 #ifdef MIEEE
 169 x += w/x;
 170 *q -= 0x10;
 171 x += w/x;
 172 *q -= 0x10;
 173 x += w/x;
 174 *q -= 0x10;
 175 #endif
 176
 177 return(x);
 178 }