lib/math/sqrt.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7 /* $Header$ */
   8
   9 #include        <math.h>
  10 #include        <float.h>
  11 #include        <errno.h>
  12
  13 #define NITER   5
  14
  15 double
  16 sqrt(double x)
  17 {
  18         int exponent;
  19         double val;
  20
  21         if (__IsNan(x)) {
  22                 errno = EDOM;
  23                 return x;
  24         }
  25         if (x <= 0) {
  26                 if (x < 0) errno = EDOM;
  27                 return 0;
  28         }
  29
  30         if (x > DBL_MAX) return x;      /* for infinity */
  31
  32         val = frexp(x, &exponent);
  33         if (exponent & 1) {
  34                 exponent--;
  35                 val *= 2;
  36         }
  37         val = ldexp(val + 1.0, exponent/2 - 1);
  38         /* was: val = (val + 1.0)/2.0; val = ldexp(val, exponent/2); */
  39         for (exponent = NITER - 1; exponent >= 0; exponent--) {
  40                 val = (val + x / val) / 2.0;
  41         }
  42         return val;
  43 }