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