home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_09_09 / 9n09010a < prev    next >
Text File  |  1991-08-06  |  1KB  |  52 lines
  1. /* sqrt function */
  2. #include <limits.h>
  3. #include "xmath.h"
  4.  
  5. double (sqrt)(double x)
  6.     {   /* compute sqrt(x) */
  7.     short xexp;
  8.  
  9.     switch (_Dunscale(&xexp, &x))
  10.         {   /* test for special codes */
  11.     case NAN:
  12.         errno = EDOM;
  13.         return (x);
  14.     case INF:
  15.         if (DSIGN(x))
  16.             {   /* -INF */
  17.             errno = EDOM;
  18.             return (_Nan._D);
  19.             }
  20.         else
  21.             {   /* +INF */
  22.             errno = ERANGE;
  23.             return (_Inf._D);
  24.             }
  25.     case 0:
  26.         return (0.0);
  27.     default:    /* finite */
  28.         if (x < 0.0)
  29.             {   /* sqrt undefined for reals < 0 */
  30.             errno = EDOM;
  31.             return (_Nan._D);
  32.             }
  33.          {  /* 0 < x, compute sqrt(x) */
  34.         double y;
  35.         static const double sqrt2 =
  36.             {1.41421356237309505};
  37.  
  38.         y = (-0.1984742 * x + 0.8804894) * x
  39.             + 0.3176687;
  40.         y = 0.5 * (y + x / y);
  41.         y += x / y;
  42.         x = 0.25 * y + x / y;
  43.         if ((unsigned int)xexp & 1)
  44.             x *= sqrt2, --xexp;
  45.         _Dscale(&x, xexp / 2);
  46.         return (x);
  47.          }
  48.         }
  49.     }
  50. /* End of File */
  51.  
  52.