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