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