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