4.3BSD - /usr/src/usr.lib/libom/sqrt.c

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

Defined functions

Defined variables