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