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 } 59