sys/src/cmd/map/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 0.;
  20                 return 0;
  21         }
  22         x = frexp(arg, &exp);
  23         while(x < 0.5) {
  24                 x *= 2;
  25                 exp--;
  26         }
  27         /*
  28          * NOTE
  29          * this wont work on 1's comp
  30          */
  31         if(exp & 1) {
  32                 x *= 2;
  33                 exp--;
  34         }
  35         temp = 0.5 * (1.0+x);
  36
  37         while(exp > 60) {
  38                 temp *= (1L<<30);
  39                 exp -= 60;
  40         }
  41         while(exp < -60) {
  42                 temp /= (1L<<30);
  43                 exp += 60;
  44         }
  45         if(exp >= 0)
  46                 temp *= 1L << (exp/2);
  47         else
  48                 temp /= 1L << (-exp/2);
  49         for(i=0; i<=4; i++)
  50                 temp = 0.5*(temp + arg/temp);
  51         return temp;
  52 }