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 <math.h> 9*3e12c5d1SDavid du Colombier #include <errno.h> 10*3e12c5d1SDavid du Colombier #define _RESEARCH_SOURCE 11*3e12c5d1SDavid du Colombier #include <float.h> 12*3e12c5d1SDavid du Colombier 13*3e12c5d1SDavid du Colombier double sqrt(double arg)14*3e12c5d1SDavid du Colombiersqrt(double arg) 15*3e12c5d1SDavid du Colombier { 16*3e12c5d1SDavid du Colombier double x, temp; 17*3e12c5d1SDavid du Colombier int exp, i; 18*3e12c5d1SDavid du Colombier 19*3e12c5d1SDavid du Colombier if(isInf(arg, 1)) 20*3e12c5d1SDavid du Colombier return arg; 21*3e12c5d1SDavid du Colombier if(arg <= 0) { 22*3e12c5d1SDavid du Colombier if(arg < 0) 23*3e12c5d1SDavid du Colombier errno = EDOM; 24*3e12c5d1SDavid du Colombier return 0; 25*3e12c5d1SDavid du Colombier } 26*3e12c5d1SDavid du Colombier x = frexp(arg, &exp); 27*3e12c5d1SDavid du Colombier while(x < 0.5) { 28*3e12c5d1SDavid du Colombier x *= 2; 29*3e12c5d1SDavid du Colombier exp--; 30*3e12c5d1SDavid du Colombier } 31*3e12c5d1SDavid du Colombier /* 32*3e12c5d1SDavid du Colombier * NOTE 33*3e12c5d1SDavid du Colombier * this wont work on 1's comp 34*3e12c5d1SDavid du Colombier */ 35*3e12c5d1SDavid du Colombier if(exp & 1) { 36*3e12c5d1SDavid du Colombier x *= 2; 37*3e12c5d1SDavid du Colombier exp--; 38*3e12c5d1SDavid du Colombier } 39*3e12c5d1SDavid du Colombier temp = 0.5 * (1.0+x); 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 while(exp < -60) { 46*3e12c5d1SDavid du Colombier temp /= (1L<<30); 47*3e12c5d1SDavid du Colombier exp += 60; 48*3e12c5d1SDavid du Colombier } 49*3e12c5d1SDavid du Colombier if(exp >= 0) 50*3e12c5d1SDavid du Colombier temp *= 1L << (exp/2); 51*3e12c5d1SDavid du Colombier else 52*3e12c5d1SDavid du Colombier temp /= 1L << (-exp/2); 53*3e12c5d1SDavid du Colombier for(i=0; i<=4; i++) 54*3e12c5d1SDavid du Colombier temp = 0.5*(temp + arg/temp); 55*3e12c5d1SDavid du Colombier return temp; 56*3e12c5d1SDavid du Colombier } 57