UNIX 7th Ed. - /usr/src/libm/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 <errno.h>
   9: 
  10: int errno;
  11: double frexp();
  12: 
  13: double
  14: sqrt(arg)
  15: double arg;
  16: {
  17:     double x, temp;
  18:     int exp;
  19:     int i;
  20: 
  21:     if(arg <= 0.) {
  22:         if(arg < 0.)
  23:             errno = EDOM;
  24:         return(0.);
  25:     }
  26:     x = frexp(arg,&exp);
  27:     while(x < 0.5) {
  28:         x *= 2;
  29:         exp--;
  30:     }
  31:     /*
  32: 	 * NOTE
  33: 	 * this wont work on 1's comp
  34: 	 */
  35:     if(exp & 1) {
  36:         x *= 2;
  37:         exp--;
  38:     }
  39:     temp = 0.5*(1.0+x);
  40: 
  41:     while(exp > 60) {
  42:         temp *= (1L<<30);
  43:         exp -= 60;
  44:     }
  45:     while(exp < -60) {
  46:         temp /= (1L<<30);
  47:         exp += 60;
  48:     }
  49:     if(exp >= 0)
  50:         temp *= 1L << (exp/2);
  51:     else
  52:         temp /= 1L << (-exp/2);
  53:     for(i=0; i<=4; i++)
  54:         temp = 0.5*(temp + arg/temp);
  55:     return(temp);
  56: }

Defined functions

Defined variables