sys/src/libc/port/sqrt.c

   1 /*
   2         sqrt returns the square root of its floating
   3         point argument. Newton's method.
   4
   5         calls frexp
   6 */
   7
   8 #include <u.h>
   9 #include <libc.h>
  10
  11 double
  12 sqrt(double arg)
  13 {
  14         double x, temp;
  15         int exp, i;
  16
  17         if(arg <= 0) {
  18                 if(arg < 0)
  19                         return NaN();
  20                 return 0;
  21         }
  22         if(isInf(arg, 1))
  23                 return arg;
  24         x = frexp(arg, &exp);
  25         while(x < 0.5) {
  26                 x *= 2;
  27                 exp--;
  28         }
  29         /*
  30          * NOTE
  31          * this wont work on 1's comp
  32          */
  33         if(exp & 1) {
  34                 x *= 2;
  35                 exp--;
  36         }
  37         temp = 0.5 * (1.0+x);
  38
  39         while(exp > 60) {
  40                 temp *= (1L<<30);
  41                 exp -= 60;
  42         }
  43         while(exp < -60) {
  44                 temp /= (1L<<30);
  45                 exp += 60;
  46         }
  47         if(exp >= 0)
  48                 temp *= 1L << (exp/2);
  49         else
  50                 temp /= 1L << (-exp/2);
  51         for(i=0; i<=4; i++)
  52                 temp = 0.5*(temp + arg/temp);
  53         return temp;
  54 }