1 /* 2 sqrt returns the square root of its floating 3 point argument. Newton's method. 4 5 calls frexp 6 */ 7 8 #include <math.h> 9 #include <errno.h> 10 #define _RESEARCH_SOURCE 11 #include <float.h> 12 13 double sqrt(double arg)14sqrt(double arg) 15 { 16 double x, temp; 17 int exp, i; 18 19 if(isInf(arg, 1)) 20 return arg; 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 } 57